数值分析实验报告(Matlab实现).doc

1、学 生 实 验 报 告实验课程名称 数值分析 开课实验室 数学与统计学院实验室 学 院 2010 年级 数学与应用数学专业班 01班 学 生 姓 名 学 号 开 课 时 间 2012 至 2013 学年第 一 学期总 成 绩教师签名课程名称数值分析实验项目名 称Gauss消元法实验项目类型验证演示综合设计其他指导教师何光辉成 绩一、 实验目的：（1）高斯列主元消去法求解线性方程组的过程（2）熟悉用迭代法求解线性方程组的过程（3）设计出相应的算法，编制相应的函数子程序二、 实验内容分别用高斯列主元消元法和直接消元法求解线性方程组：三、 实验原理对于线性方程组 （1）常记为矩阵形式 （2）根据高等

2、代数的知识，若，上式的解存在且唯一。（1） Gauss直接消元法考虑上述线性方程组的增广矩阵，对增广矩阵进行行变换，将（2）式化为等价的三角形方阵，然后回代解之，这就是Gauss消元法。具体如下：a) 消元令；对k=1到n-1，若，进行b) 回代，若（2） Gauss列主元消元法设列主元消元法已完成的第k-1（）次消元，的到方程组在进行第k次消元前，先进行2个步骤：a) 在至这一列内选出最大值，即，若，此时方程组无确定解，应给出退出信息。b) 若，则交换第行和行，然后用Gauss消元法进行消元。四、 MATLAB软件实现（1） 写出Gauss消元法和列主元消元法实现的MATLAB函数 根据以上

3、的算法，写出如下程序：%Gauss消元法%function y=Gauss1(A,b)m,n=size(A);%检查系数正确性if m=n error(矩阵A的行数和列数必须相同); return;endif m=size(b) error( b的大小必须和A的行数或A的列数相同); return;end%再检查方程是否存在唯一解if rank(A)=rank(A,b) error( A矩阵的秩和增广矩阵的秩不相同,方程不存在唯一解); return;end%这里采用增广矩阵行变换的方式求解c=n+1;A(:,c)=b; %消元过程for k=1:n-1A(k+1:n, k:c)=A(k+1:

4、n, k:c)-(A(k+1:n,k)/ A(k,k)*A(k, k:c); end%回代结果x=zeros(length(b),1); x(n)=A(n,c)/A(n,n);for k=n-1:-1:1x(k)=(A(k,c)-A(k,k+1:n)*x(k+1:n)/A(k,k);end%显示计算结果%disp(x=);%disp(x);y=x;% %高斯列主元消元法求解线性方程组Ax=b%A为输入矩阵系数,b为方程组右端系数%方程组的解保存在x变量中function y=Gauss_line(A,b)format long;% 设置为长格式显示,显示15位小数m,n=size(A);%先检

5、查系数正确性if m=n error(矩阵A的行数和列数必须相同); return;endif m=size(b) error( b的大小必须和A的行数或A的列数相同); return;end%再检查方程是否存在唯一解if rank(A)=rank(A,b) error( A矩阵的秩和增广矩阵的秩不相同,方程不存在唯一解); return;endc=n+1;A(:,c)=b; %(增广)for k=1:n-1r,m=max(abs(A(k:n,k); %选主元 m=m+k-1; %修正操作行的值 if(A(m,k)=0) if(m=k) A(k m,:)=A(m k,:); %换行 end A

6、(k+1:n, k:c)=A(k+1:n, k:c)-(A(k+1:n,k)/ A(k,k)*A(k, k:c); %消去 endendx=zeros(length(b),1); %回代求解x(n)=A(n,c)/A(n,n);for k=n-1:-1:1 x(k)=(A(k,c)-A(k,k+1:n)*x(k+1:n)/A(k,k);endy=x;format short;% 设置为默认格式显示,显示5位（2） 建立MATLAB界面 利用MATLAB的GUI建立如下界面求解线性方程组： 详见程序。五、 计算实例、数据、结果、分析下面我们对以上的结果进行测试，求解：输入数据后点击和，得到如下结

7、果：更改以上数据进行测试，求解如下方程组：得到如下结果：六、 实验中遇到的问题及解决办法在本实验中，遇到的问题主要有两个：（1） 如何将上述的Gauss消元法的算法在MATLAB中实现针对此问题我借鉴了网上以及 课本上的算法的MATLAB实现的程序；（2） 如何将建立界面使得可以随意输入想要求解的相关矩阵后就可以直接求解针对此问题，我通过网上的一些关于MATLAB的GUI设计的相关资料，总结经验完成了此项任务。七、 实验结论 通过以上的测试，我们发现以上算法和程序能够求出线性方程组的比较精确解。八、 参考文献1杨大地，王开荣 .2006.数值分析.北京:科学出版社2何光辉.2008. 数值分析

8、实验. 重庆大学数理学院数学实验教学中心3百度文库，百度知道教师签名年 月 日课程名称数值分析实验项目名 称插值方法实验项目类型验证演示综合设计其他指导教师何光辉成 绩一、 实验目的：（1） 学会拉格朗日插值、牛顿插值等基本方法（2） 设计出相应的算法，编制相应的函数子程序（3） 会用这些函数解决实际问题二、 实验内容（1）设计拉格朗日插值算法，编制并调试相应的函数子程序（2）设计牛顿插值算法，编制并调试相应的函数子程序（3）给定函数四个点的数据如下：X1.12.33.95.1Y3.8874.2764.6512.117 试用拉格朗日插值确定函数在x=2.101，4.234处的函数值。（4）已知

9、用牛顿插值公式求的近似值。三、 实验原理（1） 拉格朗日插值n次拉格朗日插值多项式为：Ln(x)=y0l0(x)+y1l1(x)+y2l2(x)+ynln(x)n=1时，称为线性插值，L1(x)=y0(x-x1)/(x0-x1)+ y1(x-x0)/(x1-x0)=y0+(y1-x0)(x-x0)/(x1-x0)n=2时，称为二次插值或抛物线插值，精度相对高些L2(x)=y0(x-x1)(x-x2)/(x0-x1)/(x0-x2)+y1(x-x0)(x-x2)/(x1-x0)/(x1-x2)+y2(x-x0)(x-x1)/(x2-x0)/(x2-x1)对节点xi(i=0,1,n)中任一点xk(

10、0=k=n)作一n 次多项式lk(xk)，使它在该点上取值为1，而在其余点xi(i=0,1,k-1,k+1,n)上为0，则插值多项式为Ln(x)=y0l0(x)+y1l1(x)+y2l2(x)+ynln(x)上式表明：n 个点xi(i=0,1,k-1,k+1,n)都是lk(x)的零点。（2） 牛顿插值插商公式Newton插值多项式为四、 MATLAB软件实现（1） 分别写出lagrange插值法和Newton插值法的求解函数%lagrange插值法求解函数%x，y为初始数据，z为插值点function z=lagrange(x,y,a)format long;% 显示15位n=length(x

11、);% 取长度%初始计算s=0;%进入公式计算for j=0:(n-1) t=1; for i=0:(n-1) if i=j t=t*(a-x(i+1)/(x(j+1)-x(i+1); end end s=s+t*y(j+1);endz=s; %显示输出结果format short;%Newton插值法求解函数%x，y为初始数据，z为插值点function j=Newton(x,y,z)n=max(size(x);l=1;a=y(1);B=a;s=1; %一次因子的乘积,预设为1dx=y; %差商for i=1:n-1 dx0=dx; for j=1:n-i dx(j)=(dx0(j+1)-d

12、x0(j)/(x(i+j)-x(j); end df=dx(1); s=s*(z-x(i); %一次因子乘积 a=a+s*df; %计算各次Newton插值的值 l=l+1; B=a;% 结果保存在变量B中endj=B; （2） 建立界面利用MATLAB中的GUI编程建立如下界面： 详见程序。五、 计算实例、数据、结果、分析下面我们对以上的问题进行测试：输入数据：计算结果如下：当x=2.101时，x=4.234时，同理可以测试（4）中的的值。六、 实验中遇到的问题及解决办法在本实验中，遇到的问题主要有两个：（3） 如何将上述的插值的算法在MATLAB中实现针对此问题我借鉴了网上以及 课本上的算

13、法的MATLAB实现的程序；（4） 如何将建立界面使得可以随意输入想要求解的相关矩阵后就可以直接求解针对此问题，我通过网上的一些关于MATLAB的GUI设计的相关资料，总结经验完成了此项任务。七、 实验结论 通过以上的测试，我们发现以上算法和程序能够求出插值的比较精确解。八、 参考文献1杨大地，王开荣 .2006.数值分析.北京:科学出版社2何光辉.2008. 数值分析实验. 数理学院数学实验教学中心3百度文库，百度知道教师签名年 月 日课程名称数值分析实验项目名 称数值微积分实验项目类型验证演示综合设计其他指导教师何光辉成 绩一、 实验目的：（1）学会复化梯形、复化辛浦生求积公式的应用（2）

14、设计出相应的算法，编制相应的函数子程序（3）会用这些函数解决实际问题二、 实验内容（1）设计复化梯形公式求积算法，编制并调试相应的函数子程序（2）设计复化辛浦生求积算法，编制并调试相应的函数子程序（4）分别用复化梯形公式和复化辛浦生公式计算定积分三、 实验原理（1）复化梯形求积公式开始定义；给出结束输入输出图1 复化梯形求积公式算法的流程图Step1给出被积函数、区间端点和等分数；Step2求出；Step3计算；Step4得（2）复化辛普森求积公式开始定义；给出结束输入输出图2 复化辛普森求积公式算法的流程图Step1 给出被积函数、区间端点和等分数；Step2求出；Step3计算； Step

15、4得四、 MATLAB软件实现（1） 分别写出复化梯形和复化辛浦生求积的求解函数%复化梯形公式求积分值%function T=trap(f,a,b)%f为积分函数%a,b 为积分区间%n是等分区间份数n=200;h=(b-a)/n;% 步长T=0;for k=1:(n-1) x0=a+h*k; T=T+limit(f,x0);endT=h*(limit(f,a)+limit(f,b)/2+h*T;T=double(T);%Simpson公式求积分值%function S=simpson(f,a,b)%f为积分函数%a,b为积分区间%n是等分区间份数n=200;h=(b-a)/(2*n);% 步

16、长s1=0;s2=0;for k=1:n x0=a+h*(2*k-1); s1=s1+limit(f,x0);endfor k=1:(n-1) x0=a+h*2*k; s2=s2+limit(f,x0);endS=h*(limit(f,a)+limit(f,b)+4*s1+2*s2)/3;S=double(S);（2） 建立界面利用MATLAB中的GUI编程建立如下界面： 详见程序。五、 计算实例、数据、结果、分析下面我们对以上的问题进行测试：输入数据：计算结果如下： 六、 实验中遇到的问题及解决办法在本实验中，遇到的问题主要有两个：（5） 如何将上述的积分算法在MATLAB中实现针对此问题我

17、借鉴了网上以及 课本上的算法的MATLAB实现的程序；（6） 如何将建立界面使得可以随意输入想要求解的相关矩阵后就可以直接求解针对此问题，我通过网上的一些关于MATLAB的GUI设计的相关资料，总结经验完成了此项任务。七、 实验结论 通过以上的测试，我们发现以上算法和程序能够求出积分的比较精确解。八、 参考文献1杨大地，王开荣 .2006.数值分析.北京:科学出版社2何光辉.2008. 数值分析实验. 数理学院数学实验教学中心3百度文库，百度知道教师签名年 月 日课程名称数值分析实验项目名 称常微分方程的数值解法实验项目类型验证演示综合设计其他指导教师何光辉成 绩一、 实验目的：（1）学会欧拉

18、方法和四阶龙格-库塔方法的使用（2）设计出相应的算法，编制相应的函数子程序（3）会用这些函数解决实际问题二、 实验内容用欧拉方法和四阶龙格-库塔方法求解微分方程初值 问题：y=sin(y)*x，y(0)=10，求y(1)。三、 实验原理(1) 欧拉法欧拉法是解初值问题的最简单的数值方法。从（9.2）式由于y (x0) = y0已给定，因而可以算出设x1 = h充分小，则近似地有：记 从而我们可以取作为y (x1)的近似值。利用y1及f (x1, y1)又可以算出y(x2)的近似值：一般地，在任意点xn+1 = (n + 1)h处y(x)的近似值由下式给出这就是欧拉法的计算公式，h称为步长。 在

19、实际计算时，可将欧拉法与梯形法则相结合，计算公式为：（2）四阶龙格-库塔方法四阶龙格-库塔法求解公式如下：四、 MATLAB软件实现（1） 分别写出欧拉方法和四阶龙格-库塔方法求解微分方程的求解函数%Euler法，初值y(a)=c%function y=Euler(f,a,b,c)n=1000;h=(b-a)/n;X=a:h:b;Y=zeros(1,n+1);Y(1)=c;for i=2:n+1 x=X(i-1); y=Y(i-1); Y(i)=Y(i-1)+eval(f)*h;end%-%function y=RK(f,a,b,c)n=1000;h=(b-a)/n;X=a:h:b;Y=zer

20、os(1,n+1);Y(1)=c;for i=1:n x=X(i); y=Y(i); K1=h*eval(f); x=x+h/2; y=y+K1/2; K2=h*eval(f); x=x; y=Y(i)+K2/2; K3=h*eval(f); x=X(i)+h; y=Y(i)+K3; K4=h*eval(f); Y(i+1)=Y(i)+(K1+2*K2+2*K3+K4)/6;end（2） 建立界面利用MATLAB中的GUI编程建立如下界面： 详见程序。五、 计算实例、数据、结果、分析下面我们对以上的问题进行测试：输入数据：计算结果如下：六、 实验中遇到的问题及解决办法在本实验中，遇到的问题主要

21、有两个：（7） 如何将上述的求解微分方程的算法在MATLAB中实现针对此问题我借鉴了网上以及 课本上的算法的MATLAB实现的程序；（8） 如何将建立界面使得可以随意输入想要求解的相关矩阵后就可以直接求解针对此问题，我通过网上的一些关于MATLAB的GUI设计的相关资料，总结经验完成了此项任务。七、 实验结论 通过以上的测试，我们发现以上算法和程序能够求出微分方程组的比较精确解。八、 参考文献1杨大地，王开荣 .2006.数值分析.北京:科学出版社2何光辉.2008. 数值分析实验. 数理学院数学实验教学中心3百度文库，百度知道教师签名年 月 日课程名称数值分析实验项目名 称估计水塔的水流量实

22、验项目类型验证演示综合设计其他指导教师何光辉成 绩一、 实验目的：（1）学会对实际问题的分析方法（2）学会利用所学的知识解决实际问题（3）设计出相应的算法，编制相应的应用程序二、 实验内容某居民区，其自来水是有一个圆柱形水塔提供，水塔高12.2m，塔的直径为17.4m，水塔是由水泵根据水塔中的水位自动加水，一般水泵每天工作两次。按照设计，当水塔中的水位降低至最低水位，约8.2m时，水泵自动启动加水。当水位升至最高水位，约10.8m时，水泵停止工作。表略。三、 实验原理计算中将流量定义为单位时间流出的水的高度乘以水塔横截面积。把时间分成5段：第1未供水段、水泵开启第1段、第2未供水段、水泵开启第

23、2段、第3未供水段。先直接对第1、2、3未供水段进行5次曲线拟合。再对得到的曲线分别求导，取得流速（即单位时间内流出的水的高度）。水泵开启第1、2段，分别在两端各取两个点，用时刻流速进行拟合得到这两段的流速。流速乘以水塔横截面积就得到任何时刻的水流量。对其进行分段积分，求和得到一天的总水流量。四、 MATLAB软件实现（1）程序function pushbutton1_Callback(hObject, eventdata, handles)% hObject handle to pushbutton1 (see GCBO)% eventdata reserved - to be define

24、d in a future version of MATLAB% handles structure with handles and user data (see GUIDATA)figure(1);x=0,3316,6635,10619,13937,17921,21240,25223,28543,32284,39435,43318,46636,49953,53936,57254,60574,64554,68535,71854,75021,85968,89953,93270;y=31.75,31.10,30.54,29.94,29.55,28.92,28.50,27.87,27.52,26.

25、97,35.50,34.45,33.50,32.67,31.56,30.81,30.12,29.27,28.42,27.67,26.97,34.75,33.89,33.40;t=x/3600; %时间单位为小时h=y/3.281; %水位高度为米x1=t(1:10);y1=h(1:10);f1=polyfit(x1,y1,5);t1=0:0.01:t(10);h1=polyval(f1,t1);plot(x1,y1,o,t1,h1,k);xlabel(时间（h） );ylabel(水位(m);title(第一阶段供水的时间水位图)function pushbutton2_Callback(hObject, eventdata, handles)% hObject handle to pushbutton2 (see GCBO)% eventdata reserved - to be defined in a future version of MATLAB% handles structure with handles and user data (see GUIDATA)figure(2);x=0,3316,6635,10619,13937,17921,21240,25223,28543,32284,39435,43318,46