数值分析试验报告Matlab实现

1、学生实验报告实验课程名称数值分析开课实验室数学与统计学院实验室学 院 2010 年级数学与应用数学专业班01班学生姓名 学 号开课时 间 2012 至 2013 学年第 一 学期总成绩教师签名根据高等代数的知识，若| a|。0,上式的解存在且唯一。(1)(2)课程名称数值分析实验项目名 称Gauss Mte 法实验项目类型验证演示综合设计其他指导教师何光辉成 绩实验目的：(1) 高斯列主元消去法求解线性方程组的过程(2) 熟悉用迭代法求解线性方程组的过程(3) 设计出相应的算法，编制相应的函数子程序实验内容分别用高斯列主元消元法和直接消元法求解线性方程组:一2100- 3 1 X1 I一 10

2、1_ 3-4-1213|X251 1123-4X31-2 1149T3 一乂 一1 17 一实验原理对于线性方程组811X1+812X2+ +anXn=b821X1+822X2+ +a2nXn=b2an1X1 +2为 * +annXn=bn常记为矩阵形式Ax = b15(1) Gauss直接消元法考虑上述线性方程组的增广矩阵 A：b,对增广矩阵进行行变换，将(2)式化为等价的三角形方阵，然后回代解之，这就是Gauss消元法。具体如下:a)消元令a/=aj，1，j =1,2,，n ; W)= h , i =1,2,，n对k=1到n-1，若a；：)黄0 ,进行1 一a(kk)lik /, akka

3、(k书 aik=0,(k+) a 口 ja(k) = aij -4 >akjk)=b-ik * A ,i =k +1,k +2,，ni =k +1,k +2,，ni, j =k +1,k +2,，ni =k +1,k +2,，nb)回代，若a；n #0，b；Xn =Ebnn1(i) v (i)、Xi = i (bi - ' aj Xj)aiij 占 1(2) Gauss列主元消元法设列主元消元法已完成Ax=b的第k-1 (1玄kn-1)次消元，的到方程组Ax = b A(k)x = b(k)在进行第k次消元前，先进行 2个步骤：(k)(k)(k)(k)l(k)A =0a) 在ak

4、k至ank这一列内选出取大值，即礼*| = PQaX aik | ,右aik,k = 0 ,方程组无确定解，应给出退出信息。b) 若a(kkk #0，则交换第ik行和k行，然后用Gauss消元法进行消元。四、MATLA歌件实现(1) 写出Gauss消元法和列主元消元法实现的MA TLAB函数根据以上的算法，写出如下程序：%Gauss元法 %function y=Gauss1(A,b)m,n=size(A);%检查系数正确性if m=nerror( '矩阵A的行数和列数必须相同);return ;endif m=size(b)error( ' b的大小必须和A的行数或A的列数相同

5、');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:n, k:c)-(A(k+1:n,k)/ A(k,k)*A(k, k:c); end%U代结果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)

6、*x(k+1:n)/A(k,k); end% 显示计算结果%disp('x=');%disp(x);y=x;% %U主元消元法求解线，性方程组Ax=b%隅输入矩阵系数,b为方程组右端系数 %方程组的解保存在x变量中function y=Gauss_line(A,b) format long ; %设置为长格式显示，显示15位小数 m,n=size(A);% 先检查系数正确性if m=nerror('矩阵A的行数和列数必须相同);return;end if m=size(b)error(' b的大小必须和A的行数或A的列数相同');return;end%再

7、检查方程是否存在唯一解if rank(A)=rank(A,b)error(' A矩阵的秩和增广矩阵的秩不相同，方程不存在唯一解');return;endc=n+1;A(:,c)=b;%(增广)for k=1:n-1 r,m=max(abs(A(k:n,k);%i 主元m=m+k-1;%修正操作行的值if (A(m,k)=0) if (m=k)A(k m,:)=A(m k,:);濒行endA(k+1:n, k:c)=A(k+1:n, k:c)-(A(k+1:n,k)/ A(k,k)*A(k, k:c);唯肖去end end x=zeros(length(b),1);%回代求解x(

8、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); end y=x; format short ; %设置为默认格式显示，显示5位(2) 建立MATLAB界面利用MATLAB的GUI建立如下界面求解线性方程组：详见程序。:、计算实例、数据、结果、分析下面我们对以上的结果进行测试，求解:一2100-3 | x 1一 10-3-4-1213|X251123-4X31-214149一13 一乂 一7 一G-L求君京WG求程*受输入数据后点击,得到如下结果：Gauss直接法求解结果Gauss列主元法求解结果

9、更改以上数据进行测试，求解如下方程组:一4得到如下结果:求解线性方程组Ax=b清按matlab中输入矩阵方式输入矩阵A4 33 1；3432343;1 234请按matlab中输入矩阵方式输入b11 -1 -1Gauss直接法求解结果-1-1.3323fl-0T62775Gf-0l7Gauss列主元法求解结果'福L术霍元-t. 3323e-016Copyright: yfg Tel、实验中遇到的问题及解决办法在本实验中，遇到的问题主要有两个：(1) 如何将上述的Gauss消元法的算法在 MATLAB中实现针对此问题我借鉴了网上以及课本上的算法的 MATLAB实

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

11、编制相应的函数子程序(3) 会用这些函数解决实际问题实验内容(1) 设计拉格朗日插值算法，编制并调试相应的函数子程序(2) 设计牛顿插值算法，编制并调试相应的函数子程序(3) 给定函数四个点的数据如下：X1.12.33.95.1Y3.8874.2764.6512.117试用拉格朗日插值确定函数在x=2.101, 4.234处的函数值。(4)已知1 =1,w'4 = 2,廿9 =3,用牛顿插值公式求 .5的近似值。三、实验原理(1) 拉格朗日插值n 次拉格朗日插值多项式为：Ln(X)=y0l0(X)+y1l1(X)+y2l2(X)+- +ynln(X)n=1 时，称为线性插值,L1(x)

12、=y0(x-x1)/(x0-x1)+ y1(x-x0)/(x1-x0)=y0+(y1-x0)(x-x0)/(x1-x0)n=2时，称为二次插值或抛物线插值，精度相对高些1,而在其+ynln(x)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(0<=k<=n)作一 n次多项式lk(xk),使它在该点上取值为余点 xi(i=0,1,-k,k+1, -,n止为 0,则插值多项式为Ln(x)=y0

13、l0(x)+y1l1(x)+y2l2(x)+-上式表明：n个点xi(i=0,1, - -k,k+1, -,n廊是lk(x)的零点。(2) 牛顿插值插商公式f(Xj)k和3*,.冬:j(Xj -X1).(Xj-Xj)(Xj -Xj).(Xj -Xk)Newton插值多项式为NJX) = f(Xi) fX,X2(X-Xi)fX,X2，X3(X-Xi)(X -X2)f Xi,X2,X3,.Xn(X X)(X X2).(X XnQ；四、MATLABC件实现(1) 分别写出lagrange插值法和Newton插值法的求解函数%lagran 阿法求解函数 % %x, y为初始数据，z为插值点 functi

14、on z=lagrange(x,y,a) format long ; % 显示 15 位 n=length(x); % 取长度 %初始计算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);endends=s+t*y(j+1);endz=s;%显示输出结果format short ;%NeWt姒求解函数 %x, y为初始数据，z为插值点function j=Newton(x,y,z)n=max(size(x);l=1;a=y(1);B=a;s=1;%一次因子的乘积，预设为1dx=y;%li

15、商for i=1:n-1dx0=dx;for j=1:n-idx(j)=(dx0(j+1)-dx0(j)/(x(i+j)-x(j);enddf=dx(1);s=s*(z-x(i);%一次因子乘积a=a+s*df;%计算各次Newton插值的值l=l+1;B=a; %结果保存在变量 B中endj=B;(2)建立界面利用MATLAB中的GUI编程建立如下界面：请按matlab中的矩阵方式输入Xi请按matlab中的矩阵方式输入Yi输入插值点x|-Lagrange插值结果:执行插应LNewton插值结果曲虾边立网Copyright: yfg Tel:见程序。计算实例、数据、

16、结果、分析 下面我们对以上的问题进行测试： 输入数据：2 何光辉.2008.数值分析实验.数理学院数学实验教学中心3 百度文库，百度知道教师签名年 月曰课程名称数值分析实验项目名 称数值微积分实验项目类型验证演示综合设计其他指导教师何光辉成 绩一、实验目的：(1) 学会复化梯形、复化辛浦生求积公式的应用(2) 设计出相应的算法，编制相应的函数子程序(3) 会用这些函数解决实际问题二、实验内容(1) 设计复化梯形公式求积算法，编制并调试相应的函数子程序(2) 设计夏化辛浦生求积算法，编制并调试相应的函数子程序(4) 分别用复化梯形公式和复化辛浦生公式计算定积分,i sin x , dx0 x三、

17、实验原理(1)复化梯形求积公式:开始定义f (x)；给出a,b/输入 n7, b a -hu ;xu a;T u 0nXi =a 十i h;Tu T 十 f (%);(i =1,n)Tu gf(a)+2T + f(b) 2*/输出Tn/.，.吉束图1复化梯形求积公式算法的流程图Stepl给出被积函数f (x)、区问a,b端点a,b和等分数n ;b aStep2 求出 xk = kh, h=;nStep3 计算 f(a), f(b)* f(xQ ;k=0Step4 得Tn =1hf(a)+iff(xk) + f(b) 2k =0(2)复化辛普森求积公式图2复化辛普森求积公式算法的流程图Stepl

18、给出被积函数f(x)、区问a,b端点a,b和等分数n ；b aStep2 求出 xk = kh, h =;nStep3 计算 f(a), f(b),成 f(xQ,成 f (x)；k=sk=0k2n -1n -1Step4 得&=hf(a)+4£ f(x)+2£ f(xQ + f(b)6ek2 E四、 MATLA歌件实现(1) 分别写出复化梯形和复化辛浦生求积的求解函数%«W 公式求积分值%function T=trap(f,a,b)%f为积分函数 %a,b为积分区间%门是等分区间份数 n=200; h=(b-a)/n;% 步长T=0; for k=1:(n

19、-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=200;h=(b-a)/(2*n);% 步长s1=0; s2=0; for k=1:n x0=a+h*(2*k-1); s1=s1+limit(f,x0); end for k=1:(n-1) x0=a+h*2*k; s2=s2+limit(f,x0); end S=h*(limit(f,a)+li

20、mit(f,b)+4*s1+2*s2)/3; S=double(S);(2) 建立界面利用MATLAB中的GUI编程建立如下界面：在本实验中，遇到的问题主要有两个：(5) 如何将上述的积分算法在MATLAB中实现针对此问题我借鉴了网上以及课本上的算法的 MATLAB实现的程序；(6) 如何将建立界面使得可以随意输入想要求解的相关矩阵后就可以直接求解针对此问题，我通过网上的一些关于MATLAB的GUI设计的相关资料，总结经验完成了此项任务。七、实验结论通过以上的测试，我们发现以上算法和程序能够求出积分的比较精确解。八、参考文献1 杨大地，王开荣.2006.数值分析.北京:科学出版社2 何光辉.2

21、008.数值分析实验.数理学院数学实验教学中心3 百度文库，百度知道教师签名年 月曰课程名称数值分析实验项目名 称常微分方程的数值解法实验项目类型验证演示综合设计其他指导教师何光辉成 绩一、实验目的：(1) 学会欧拉方法和四阶龙格-库塔方法的使用(2) 设计出相应的算法，编制相应的函数子程序(3) 会用这些函数解决实际问题二、实验内容用欧拉方法和四阶龙格-库塔方法求解微分方程初值1可题：y -sin(y)*x , y(0)=10,求y(1)。三、实验原理(1)欧拉法欧拉法是解初值问题的最简单的数值方法。从(9.2)式由于y (X0) = y°已给定，因而可以算出y'(x

22、76;) = f (xo, y°)设x = h充分小，则近似地有：y(x1)y(x0)顷xo) = f(x0,yo) h记 yi =y(xi) i =0,1，,n从而我们可以取y = y(o =hf(x°,y。)作为y (Xi)的近似值。利用yi及f (Xi, yi)又可以算出y(X2)的近似值:V2=V hf (Xi,yi)-般地，在任意点 Xn+i = (n + i)h处y(x)的近似值由下式给出yn d = yn hf ( xn，yn )这就是欧拉法的计算公式，h称为步长。在实际计算时，可将欧拉法与梯形法则相结合，计算公式为：y：0i =yn hf(Xn,yn)k =

23、 0,i,2,.ynk：)= yn +?f(Xn,yn)+ f(Xz,yn%)】(2)四阶龙格-库塔方法四阶龙格-库塔法求解公式如下：h 一 一一 一yn# =yn + 二(ki +2k2 +2k3 +k4) 6ki = f (Xn,yn)(ih'«2 =f Xn+ 二 h, y- kiI22)(ih'k3 = f Xn +h, yn +k2 I k4 =f Xnh, ynhk3四、 MATLA歌件实现(i)分别写出欧拉方法和四阶龙格-库塔方法求解微分方程的求解函数%Euiefc,初值 y(a)=c%function y=Euler(f,a,b,c)n=i000;h=

24、(b-a)/n;X=a:h:b;Y=zeros(i,n+i);Y(i)=c;for i=2:n+iX=X(i-i);y=Y(i-i);Y(i)=Y(i-i)+eval(f)*h;end%?-?a?t - ? - - - %function y=RK(f,a,b,c)i6n=1000;h=(b-a)/n;X=a:h:b;Y=zeros(1,n+1);Y(1)=c;for i=1:nx=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

25、*eval(f);Y(i+1)=Y(i)+(K1+2*K2+2*K3+K4)/6;end(2)建立界面利用MATLAB中的GUI编程建立如下界面:微分方程数值解(V=f(x,v)请输入函数f 3v)请输入初值点X。和所求点X1请输入初值点X。的函数值y(X0)EulerS求解结果四阶R-K公式求解结果详见程序。97卸3Eu®法求解结果Ejler求程四阶R-K公式求解结果六、实验中遇到的问题及解决办法在本实验中，遇到的问题主要有两个：(7) 如何将上述的求解微分方程的算法在MATLAB中实现针对此问题我借鉴了网上以及课本上的算法的 MATLAB实现的程序；(8) 如何将建立界面使得可以

26、随意输入想要求解的相关矩阵后就可以直接求解针对此问题，我通过网上的一些关于MATLAB的GUI设计的相关资料，总结经验完成了此项任务。七、实验结论通过以上的测试，我们发现以上算法和程序能够求出微分方程组的比较精确解。八、参考文献1 杨大地，王开荣.2006.数值分析.北京:科学出版社2 何光辉.2008.数值分析实验.数理学院数学实验教学中心3 百度文库，百度知道教师签名年 月曰课程名称数值分析实验项目名 称估计水塔的水流量实验项目类型验证演示综合设计其他指导教师何光辉成 绩一、实验目的：(1) 学会对实际问题的分析方法(2) 学会利用所学的知识解决实际问题(3) 设计出相应的算法，编制相应的

27、应用程序二、实验内容某居民区，其自来水是有一个圆柱形水塔提供，水塔高12.2m，塔的直径为17.4m,水塔是由水泵根据水塔中的水位自动加水，一般水泵每天工作两次。按照设计，当水塔中的水位降低至最低水 位，约8.2m时，水泵自动启动加水。当水位升至最高水位，约 10.8m时，水泵停止工作。表略。三、实验原理计算中将流量定义为单位时间流出的水的高度乘以水塔横截面积。把时间分成5段：第1未供水段、水泵开启第1段、第2未供水段、水泵开启第2段、第3未 供水段。先直接对第1、2、3未供水段进行5次曲线拟合。再对得到的曲线分别求 导，取得流速(即单位时间内流出的水的高度)。水泵开启第1、2段，分别在两端

28、各取两个点，用时刻流速进行拟合得到这两段的流速。流速乘以水塔横截面积就得 到任何时刻的水流量。对其进行分段积分，求和得到一天的总水流量。四、MATLA歌件实现(1)程序function pushbutton1_Callback(hObject, eventdata, handles)% hObject handle to pushbutton1 (see GCBO)% eventdata reserved - to be defined in a future version of MATLAB% handles structure with handles and user data (se

29、e GUIDATA)figure(1);x=0,3316,6635,10619,13937,17921,21240,25223,28543,32284,39435,43318,46636,4995;y=31.75,31.10,30.54,29.94,29.55,28.92,28.50,27.87,27.52,26.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;%M间单位为小时h=y/3.281;%水位高度为米x1=t(1:10);y1=h(1:10)

30、;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

31、- 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,46636,4995;y=31.75,31.10,30.54,29.94,29.55,28.92,28.50,27.87,27.52,26.97,35.50,34.45,33.50,32.67,31.56,30.81,30.12,29

32、.27,28.42,27.67,26.97,34.75,33.89,33.40;t=x/3600;%M间单位为小时h=y/3.281;%水位高度为米x2=t(11:21);y2=h(11:21);f2=polyfit(x2,y2,5);t2=t(11):0.01:t(21);h2=polyval(f2,t2);plot(x2,y2,'o' ,t2,h2,'r');xlabel( '时间(h)');ylabel( '水位(m)');title('第二阶段供水时段水位图')% - Executes on button

33、press in pushbutton4.function pushbutton4_Callback(hObject, eventdata, handles)% hObject handle to pushbutton4 (see GCBO)% eventdata reserved - to be defined in a future version of MATLAB% handles structure with handles and user data (see GUIDATA)figure(3);x=0,3316,6635,10619,13937,17921,21240,25223

34、,28543,32284,39435,43318,46636,4995;y=31.75,31.10,30.54,29.94,29.55,28.92,28.50,27.87,27.52,26.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;x3=t(22:24);y3=h(22:24);f3=polyfit(x3,y3,5);t3=t(22):0.01:t(24);h3=polyval(f3,t3);plot(x3,y3, 'o' ,t3,h3,'r');xlabel( '时间(h)');ylabel( '水位(m)');title( '第三阶段的时间水位图')% - Executes on button press in pushbutton7.functionpus