数值分析实验报告.pdf

《数值分析实验报告.pdf》由会员分享，可在线阅读，更多相关《数值分析实验报告.pdf（44页珍藏版）》请在淘文阁 - 分享文档赚钱的网站上搜索。

1、.数值分析实验报告数值分析实验报告（第二章第二章）实验题目实验题目：分别用二分法、牛顿迭代法、割线法、史蒂芬森迭代法求方程的根，观察不同初始值下的收敛性，并给出结论。问题分析问题分析：题目有以下几点要求：1.不同的迭代法计算根，并比较收敛性。2.选定不同的初始值，比较收敛性。实验原理实验原理：各个迭代法简述二分法：取有根区间的重点，确定新的有根区间的区间长度仅为区间长度的一版。对压缩了的有根区间重复以上过程，又得到新的有根区间，其区间长度为的一半，如此反复，可得一系列有根区间，区间收敛到一个点即为根。牛顿迭代法：不动点迭代法的一种特例，具有局部二次收敛的特性。迭代格式为割线法：是牛顿法的改进，

2、具有超线性收敛的特性，收敛阶为 1.618.迭代格式为.学习帮手.史蒂芬森迭代法：采用不动点迭代进行预估校正。至少是平方收敛的。迭代格式为这里可采用牛顿迭代法的迭代函数。实验内容实验内容：1.写出该问题的函数代码如下：functionpy=f(x)syms k;y=(k2+1)*(k-1)5;yy=diff(y,k);py(1)=subs(y,k,x);py(2)=subs(yy,k,x);end2.分别写出各个迭代法的迭代函数m(1)=a;代码如下：while abs(a-b)e二分法：t=(a+b)/2;functiony=dichotomie(a,b,e)s1=f(a);i=2;s2=f

3、(b);.学习帮手.s3=f(t);if s1(1)*s3(1)=es=f(x);t=x-s(1)/s(2);en=t-x;x=t;m(i)=t;.i=i+1;endy=x,i+1,m;end牛顿割线法：functiony=Secant(x1,x2,e)i=3;m(1)=x1,m(2)=x2;while abs(x2-x1)=es1=f(x1);s2=f(x2);t=x2-(x2-x1)*s2(1)/(s2(1)-s1(1);x1=x2;x2=t;m(i)=t;i=i+1;endy=x2,i+1,m;end史蒂芬森迭代法：Function p=StephensonIterative(x,e)i

4、=2;学习帮手.m(2)=x;en=2*e;whileabs(en)=ey=fai(x);z=fai(y);t=x-(y-x)2/(z-2*y+x);en=t-x;x=t;m(i)=t;i=i+1;endp=x,i+1,m;end3.因为经常被使用，故可以写一个函数。代码如下：function y=fai(x)s=f(x);y=x-s(1)/s(2);end4.可以绘制不同的图形来比较不同迭代法的收敛性和不同初值下的收敛性。代码如下：clear all;%相同初始值，不同迭代法下的收敛x1=dichotomie(0,3,1e-10);x2=NewtonIterative(0,1e-10);x3

5、=Secant(0,2,1e-10);x4=StephensonIterative(0,1e-10);x1(2),x2(2),x3(2),x4(2).学习帮手.figure,subplot(2,2,1),plot(x1(3:x1(2),title(二分法);subplot(2,2,2),plot(x2(3:x2(2),title(牛顿迭代法);subplot(2,2,3),plot(x3(3:x3(2),title(牛顿割线法);subplot(2,2,4),plot(x4(3:x4(2),title(史蒂芬森迭代法);figure,subplot(2,2,1),plot(x1(4:x1(2)

6、-1)-x1(1)./(x1(3:x1(2)-2)-x1(1),title(二分法);subplot(2,2,2),plot(x2(4:x2(2)-1)-x2(1)./(x2(3:x2(2)-2)-x2(1),title(牛顿迭代法);subplot(2,2,3),plot(x3(4:x3(2)-1)-x3(1)./(x3(3:x3(2)-2)-x3(1),title(牛顿割线法);subplot(2,2,4),plot(x4(4:x4(2)-1)-x4(1)./(x4(3:x4(2)-2)-x4(1),title(史蒂芬森迭代法);%不同初始值，相同迭代法下的收敛性x5=dichotomie

7、(-1,1,1e-10);x6=dichotomie(-2,3,1e-10);x7=dichotomie(0,4,1e-10);x8=dichotomie(-4,4,1e-10);x9=NewtonIterative(-2,1e-10);x10=NewtonIterative(-4,1e-10);x11=NewtonIterative(4,1e-10);x12=NewtonIterative(6,1e-10);figure,subplot(1,2,1),plot(1:x1(2)-2,x1(3:x1(2),1:x5(2)-2,x5(3:x5(2),1:x6(2)-2,x6(3:x6(2),1:x

8、7(2)-2,x7(3:x7(2),1:x8(2)-.学习帮手.2,x8(3:x8(2),title(二分法);subplot(1,2,2),plot(1:x2(2)-2,x2(3:x2(2),1:x9(2)-2,x9(3:x9(2),1:x10(2)-2,x10(3:x10(2),1:x11(2)-2,x11(3:x11(2),1:x12(2)-2,x12(3:x12(2),title(牛顿迭代法);x13=Secant(-1,1,1e-10);x14=Secant(-4,5,1e-10);x15=Secant(0,7,1e-10);x16=Secant(-8,2,1e-10);x17=St

9、ephensonIterative(-1,1e-10);x18=StephensonIterative(-4,1e-10);x19=StephensonIterative(4,1e-10);x20=StephensonIterative(6,1e-10);figure,subplot(1,2,1),plot(1:x3(2)-2,x3(3:x3(2),1:x13(2)-2,x13(3:x13(2),1:x14(2)-2,x14(3:x14(2),1:x15(2)-2,x15(3:x15(2),1:x16(2)-2,x16(3:x16(2),title(牛顿割线法);subplot(1,2,2),

10、plot(1:x4(2)-2,x4(3:x4(2),1:x17(2)-2,x17(3:x17(2),1:x18(2)-2,x18(3:x18(2),1:x19(2)-2,x19(3:x19(2),1:x20(2)-2,x20(3:x20(2),title(史蒂芬森迭代法);实验结果实验结果：1.各个迭代值分布.学习帮手.二 分 法1.510.50.51牛 顿 迭 代 法00102030400050史 蒂 芬 森 迭 代 法100牛 顿 割 线 法21.510.500501001500.51.51002468图 1.1不同迭代法下的得到的迭代值迭代值的情况如下：二分法01.50000000000

11、.75000000001.12500000000.93750000001.03125000000.98437500001.00781250000.9960937500牛顿迭代法00.20000000000.37049180320.50764420760.61461894470.69738690980.76155380910.81154111860.8506763857牛顿割线法02.00000000000.33333333330.38071968010.49828334190.57049963330.63938062440.69427858790.7411692653史蒂芬森迭代法01.355

12、55555560.98161652830.99994600030.9999999995.学习帮手.1.00195312500.99902343750.88144821230.90572974000.78027159970.8132927871当二分法的初始区间选为，误差限为，误差限为，牛顿迭代法初值选为，牛顿割线法初始点为，误差限为，史蒂芬森迭代法初始点选为 ，误差限为，迭代情况如图所示。迭代次数分别为 38 次，100 次，140 次，9 次。故而，史蒂芬森迭代法速度最快，效果最好。2.收敛情况50 x 109二 分 法0.80.70.6牛 顿 迭 代 法-50.50102030400.40

13、50史 蒂 芬 森 迭 代 法0.40.20100-10牛 顿 割 线 法10-1-0.2010203040-0.412345-2图 1.2不同迭代法下迭代值得收敛情况二分法收敛效果较差，牛顿迭代法和牛顿割线法相近，史蒂芬森迭代法收敛次数高于 1，效果最好3.不同初值的收敛情况.学习帮手.二 分 法2651432-110-2-1-2-3-4010203040-40牛 顿 迭 代 法0-350100150图 1.3二分法，牛顿迭代法下不同初值的收敛情况牛 顿 割 线 法862420-2-40.5-6-8012.5史 蒂 芬 森 迭 代 法1.502040608002468图 1.4牛顿割线法，史

14、蒂芬森迭代法下不同初值的收敛情况.学习帮手.1.二分法的五个初始区间分别为；2.牛顿迭代法的五个初始值分别为；3.牛顿割线法的五个初始区间分别为；4.史蒂芬森迭代法的五个初始值分别为；由图可知，它们最终均达到收敛。收敛性分析及结论收敛性分析及结论：1.二分法收敛较慢且不能求解崇根，但算法简单；此处牛顿法具有了平方收敛；从迭代次数上看，牛顿割线法较牛顿法的多，所以收敛性较差，是超线性收敛；史蒂芬森迭代法收敛效果最好。2.因为牛顿迭代法是局部的二次收敛，所以要注重初值的选取，本次实验中选择的初值均得到了收敛，效果比较好。牛顿割线法也应注意初值的选取。（第三章第三章）实验题目实验题目：1.区间作等距

15、划分：以为结点对函数进行插值逼近。（1）分别取用牛顿插值对进行逼近，并在同一坐标系下做出函数的图形，进行比较。写出插值函数对的逼近程度与节点个数的关系，并分析原因；（2）试用三次样条插值对进行逼近，在同一坐标下画出图形，观察样条插值函数对的逼近程度与节点个数的关系；（3）整体插值有何局限性？如何避免？.学习帮手.2.已知一组数据如下，求其拟合曲线.表 2.1数据表02132437485106117148169181019106.42108.2109.5110109.93110.49110.59110.6110.76111111.2（1）求以上数据形如的拟合曲线及其平方误差；（2）求以上数据形如

16、的拟合曲线及其平方误差；（3）通过观察（1）（2）的结果，写出你对数据拟合的认识.问题分析问题分析：题目除上述要求之外还有以下几点：1.明确整体插值和分段插值的不同。牛顿插值多项式属于整体插值，三次样条插值属于低次分段插值。2.将结果在同一坐标下绘制出。但是为了方便分析节点个数对于插值效果的影响，也可以单独绘制。3.第二题中为了确定各个参数的大小，可以进行适当变换，转化为线性，运用最小二乘法，得到拟合。实验原理实验原理：牛顿插值多项式：对于给定的插值节点构造次数不超过 的插值多项式使其满足插值条件这样得到的插值多项式插值多项式。系数称为为差商，可以通过构造差商.学习帮手.表得到。三次样条插值：

17、三次样条插值函数在每个小区间上为三次多项式；在全进上存在二阶连续导数；其次符合插值条件。中存在内置的三次样条插值函数，命令为实验内容实验内容：第一题：1.牛顿插值函数的构造代码如下：functionf=Newton(x0,y0)%牛顿多项式插值函数syms x;SZ=size(x0,2);a(1)=y0(1);y(:,1)=y0;for j=2:SZnx1=1;for i=1:SZ-j+1;nx2=nx1+j-1;y(i,j)=(y(i,j-1)-y(i+1,j-1)/(x0(nx1)-x0(nx2);nx1=nx1+1;endend.学习帮手.f=y(1,1);for j=2:SZff=y(

18、1,j);for i=1:j-1ff=ff*(x-x0(i);endf=f+ff;endend2.牛顿和三次样条插值情况及比较：代码如下：clear all;clc;syms x;fx=1/(5+x2);%牛顿多项式插值x0=-1:0.01:1;y0=subs(fx,x0);n1=1,5,10,20,25;n2=5,55,60,62,67;h1=2./n1;h2=2./n2;.学习帮手.x1=-1:h1(1):1;y1=subs(fx,x1);f1=Newton(x1,y1);y01=subs(f1,x0);x2=-1:h1(2):1;y2=subs(fx,x2);f2=Newton(x2,y

19、2);y02=subs(f2,x0);x3=-1:h1(3):1;y3=subs(fx,x3);f3=Newton(x3,y3);y03=subs(f3,x0);x4=-1:h1(4):1;y4=subs(fx,x4);f4=Newton(x4,y4);y04=subs(f4,x0);x5=-1:h1(5):1;y5=subs(fx,x5);f5=Newton(x5,y5);y05=subs(f5,x0);figure,plot(x0,y0,x0,y01,x0,y02,x0,y03,x0,y04,x0,y05),title(所有结果);x6=-1:h2(1):1;y6=subs(fx,x6);

20、f6=Newton(x6,y6);y06=subs(f6,x0);x7=-1:h2(2):1;y7=subs(fx,x7);f7=Newton(x7,y7);y07=subs(f7,x0);x8=-1:h2(3):1;y8=subs(fx,x8);f8=Newton(x8,y8);y08=subs(f8,x0);x9=-1:h2(4):1;y9=subs(fx,x9);f9=Newton(x9,y9);y09=subs(f9,x0);x10=-1:h2(5):1;y10=subs(fx,x10);f10=Newton(x10,y10);y010=subs(f10,x0);figure,plot

21、(x0,y0,x0,y06,x0,y07,x0,y08,x0,y09,x0,y010),title(龙格现象);%三次样条插值spline自带命令x0=-5:0.01:5;y0=subs(fx,x0);figure,subplot(2,3,1),plot(x0,y0),title(精确图);y11=spline(x1,y1,x0);subplot(2,3,2),plot(x0,y11),title(n=1);y12=spline(x2,y2,x0);subplot(2,3,3),plot(x0,y12),title(n=5);y13=spline(x3,y3,x0);subplot(2,3,4

22、),plot(x0,y13),title(n=10);.学习帮手.y14=spline(x4,y4,x0);subplot(2,3,5),plot(x0,y14),title(n=20);y15=spline(x5,y5,x0);subplot(2,3,6),plot(x0,y15),title(n=25);figure,plot(x0,y0,x0,y11,x0,y12,x0,y13,x0,y14,x0,y15),title(所有结果);第二题：代码如下：clear all;clc;x=2 3 4 7 8 10 11 14 16 18 19;y=106.42 108.2 109.5 110 1

23、09.93 110.49 110.59 110.6 110.76 111 111.2;%Agenda 1A=(x.*x)x ones(11,1);A1=A*A;y1=A*y;c1=A1y1x1=2:0.1:19;y1=polyval(c1,x1);plot(x,y,k*,x1,y1,r),gtext(曲线一),hold ony01=polyval(c1,x);delta1_2=sum(y-y01).2)%Agenda 2xx=-1./x;yy=log(y);.学习帮手.B=ones(11,1)xx;B1=B*B;yy1=B*yy;c2=B1yy1;syms s;a=exp(c2(1);b=c2

24、(2);a bf=a*exp(-b/s);x2=x1;y2=subs(f,x2);plot(x2,y2,b),gtext(曲线二);y02=subs(f,x);delta2_2=sum(y-y02).2)实验结果实验结果：第一题：1.牛顿插值结果.学习帮手.所 有 结 果0.20.1950.190.1850.180.1750.170.165-1-0.8-0.6-0.4-0.200.20.40.60.81图 2.1 不同节点数的牛顿插值多项式综合图龙 格 现 象0.20.190.180.170.160.15-1-0.8-0.6-0.4-0.200.20.40.60.81图 2.2 龙格现象.学习

25、帮手.由图可知，2.三次样条插值结果精 确 图n=1n=50.220.80.1510.60.10.40.0500.20-505-1-5050-505n=10n=20n=250.80.80.80.60.60.60.40.40.40.20.20.20-5050-5050-505图 2.3 不同节点下的三次样条插值结果.学习帮手.所 有 结 果0.70.60.50.40.30.20.10-5-4-3-2-1012345图 2.4 不同节点数的三次样条的综合图由图可知，牛顿插值多项式在较小的时候如差值效果良好，当变大时如就出现了龙格现象，三次样条在各个子区间内为三次多项式，拟合效果好.第二题1.第一问

26、的多项式拟合得到的拟合曲线为平方误差为第二问的拟合曲线为平方误差为.学习帮手.2.拟合曲线如图所示112曲 线 一111110曲 线 二1091081071062468101214161820112111110曲 线 二曲 线 一1091081071062468101214161820图 2.5 拟合情况.学习帮手.从图中可以看出曲线二大体符合黑点的分布情况，拟合效果较好。结论结论：整体插值随着节点个数的增多，多项式的次数也在升高。高次多项式的插值会出现龙格现象，震荡剧烈，逼近效果不理想。但是当节点很多时，低次插值的精度又不够，所以为了避免这一局限性采用分段低次插值。其中三次样条插值有良好的收

27、敛性和光滑性，效果较好。数据拟合时，可以选择趋势相似的函数形式，在求解相关参量时，将函数进行适当变换，从而运用最小二乘法得到拟合结果。（第四章第四章）实验题目实验题目：设计区间分半求积算法、龙贝格求积算法和自适应辛普森求积算法的程序，观察时，积分的结果，并给出相应的评价.问题分析问题分析：1.分半求积算法即为复合梯度求积。2.因为 越大，在内震荡地越严重，自适应辛普森求积算法很难适用，计算复杂度会很高。所以选取辛普森求积算法代替。实验原理实验原理：分半求积算法：将区间等分为 个小区间，分点为利用定积分性质，有每个小区间上均用梯形求积公式，有.学习帮手.令即为复合梯形公式。龙贝格求积算法：将步长

28、为的复合梯形公式进行査尔逊外推，就得到龙贝格求积公式。自适应辛普森求积算法：按照被奇函数在区间上的变化来安排求积结点的辛普森算法。实验内容实验内容：1.写得各个求积算法的函数。代码如下：（1）分半求积算法：functionf=CompositeTrapezoida(n)syms x;y=x*x*cos(n*x);m=1:100;a=-1;b=1;for i=1:100h=(b-a)/m(i);k=0:m(i);xk=a+k*h;yk=subs(y,xk);f(i)=0;for j=1:m(i)f(i)=f(i)+(yk(j)+yk(j+1);endf(i)=h/2*f(i);endend（2）

29、辛普森算法.学习帮手.function f=CompositeSimpson(n)syms x;y=x*x*cos(n*x);m=1:100;a=-1;b=1;for i=1:100h=(b-a)/(2*m(i);k=0:(2*m(i);xk=a+k*h;yk=subs(y,xk);f(i)=0;for j=1:m(i)f(i)=f(i)+(yk(2*j-1)+4*yk(2*j)+yk(2*j+1);endf(i)=h/3*f(i);endm;f;end（4）龙贝格求积算法function f=Romberg(n)epsilon=0.0001;.syms x;y=x*x*cos(n*x);a=

30、-1;b=1;h=b-a;fa=subs(y,a);fb=subs(y,b);T(1)=h/2*(fa+fb);m=1;while 1h=h/2;k=1:(2m-1);xk=a+k*h;yk=subs(y,xk);su=sum(yk);S(1)=h/2*(fa+fb)+h*su;for i=1:mS(i+1)=S(i)+(S(i)-T(i)/(4i-1);endif(abs(S(m+1)-T(m)epsilon)break;end学习帮手.T=S;m=m+1;endf=S(m+1);end（5）自适应辛普森求积算法functionf=AdaptiveSimpson(n)epsilon=0.00

31、1;syms x;y=x*x*cos(n*x);a=-1;b=1;h=b-a;p=a b;p0=p;ep=epsilon;m=0;q=0;f=0;while 1n1=length(ep);n=length(p0);.if(n=1)break;endh=p0(2)-p0(1);k=0:4;yk=subs(y,p0(1)+k*h/4);s0=h/6*(yk(1)+4*yk(3)+yk(5);s1=h/12*(yk(1)+4*yk(2)+yk(3);s2=h/12*(yk(3)+4*yk(4)+yk(5);if(abs(s0-s1-s2)=2ep=ep(2:n1);endq=q+1;elsem=m+

32、1;p0=yk(1),yk(3),p0(2:n);if n1=1ep=ep(1)/2 ep(1)/2;elseep=ep(1)/2 ep(1)/2 ep(2:n1);end学习帮手.if q=0endp=p0;endelseendp=p(1:q)p0;end2.各个求积公式下的积分结果代码如下:clear all;clc;%龙贝格法%复合梯形法y31=Romberg(1);y11=CompositeTrapezoida(1);y32=Romberg(10);y12=CompositeTrapezoida(10);y33=Romberg(100);y13=CompositeTrapezoida(

33、100);y34=Romberg(500);y14=CompositeTrapezoida(500);y31 y32 y33 y341:100;y11;y12;y13;y14%自适应辛普森法%复合辛普森法y41=AdaptiveSimpson(1)y21=CompositeSimpson(1);%以下复杂度太高不出结果y22=CompositeSimpson(10);%y42=AdaptiveSimpson(10);y23=CompositeSimpson(100);%y43=AdaptiveSimpson(100);y24=CompositeSimpson(500);%y44=Adaptiv

34、eSimpson(500);2:2:200;y21;y22;y23;y24%y41 y42 y43 y44实验结果实验结果：表 3.1不同积分法得到的结果.学习帮手.n10.47828319910-0.139940191100500复合梯形法40.478267253复合辛普森法00.478267257龙贝格法457-0.1401909980-0.140190999-0.00601370080.0023823806-0.00985113540.0072247738-0.2492135610.61122123339结果评价结果评价：龙贝格法的精度高，最接近精确解。因为被积分函数为随着的增大，震荡越

35、剧烈。所以积分面积越来越小。（第五章第五章）实验题目实验题目：1.分别取步长用龙格-库塔法求解计算到并与精确解比较，观察龙格-库塔法的收敛性.2.分别取步长用显示欧拉法和隐式欧拉法求解由结果分析算法的稳定性.3.选择某常微分方程初值问题的数值方法计算有效数字.的近似值，并保证有 4 位问题分析问题分析：1.龙格-库塔有很多的格式，第一题以常用的四级四阶龙格库塔方法为例。2.假设，由于.学习帮手.第三题可以转化为求解常微分方程的问题。实验原理实验原理：四级四阶龙格-库塔格式实验内容实验内容：1.求解第一题代码如下：clear all;clc%以四级四阶龙格库塔格式为例t0=1;h=0.5 0.1

36、 0.05 0.01;n=1./h+1;y1(1)=1;y2(1)=1;y3(1)=1;y4(1)=1;t1=1:0.5:2;t2=1:0.1:2;t3=1:0.05:2;t4=1:0.01:2;for j=2:n(1)k1=t1(j-1)*y1(j-1)(1/3);.学习帮手.k2=(t1(j-1)+h(1)/2)*(y1(j-1)+(h(1)/2)*k1)(1/3);k3=(t1(j-1)+h(1)/2)*(y1(j-1)+(h(1)/2)*k2)(1/3);k4=t1(j)*(y1(j-1)+h(1)*k3)(1/3);y1(j)=y1(j-1)+h(1)/6*(k1+2*k2+2*k3

37、+k4);endfor j=2:n(2)k1=t2(j-1)*y2(j-1)(1/3);k2=(t2(j-1)+h(2)/2)*(y2(j-1)+(h(2)/2)*k1)(1/3);k3=(t2(j-1)+h(2)/2)*(y2(j-1)+(h(2)/2)*k2)(1/3);k4=t2(j)*(y2(j-1)+h(2)*k3)(1/3);y2(j)=y2(j-1)+h(2)/6*(k1+2*k2+2*k3+k4);endfor j=2:n(3)k1=t3(j-1)*y3(j-1)(1/3);k2=(t3(j-1)+h(3)/2)*(y3(j-1)+(h(3)/2)*k1)(1/3);k3=(t

38、3(j-1)+h(3)/2)*(y3(j-1)+(h(3)/2)*k2)(1/3);k4=t3(j)*(y3(j-1)+h(3)*k3)(1/3);y3(j)=y3(j-1)+h(3)/6*(k1+2*k2+2*k3+k4);endfor j=2:n(4)k1=t4(j-1)*y4(j-1)(1/3);k2=(t4(j-1)+h(4)/2)*(y4(j-1)+(h(4)/2)*k1)(1/3);.学习帮手.k3=(t4(j-1)+h(4)/2)*(y4(j-1)+(h(4)/2)*k2)(1/3);k4=t4(j)*(y4(j-1)+h(4)*k3)(1/3);y4(j)=y4(j-1)+h(

39、4)/6*(k1+2*k2+2*k3+k4);endx=1:0.01:2;syms tyy=(t2+2)/3)(3/2);yy1=subs(yy,x);figure,subplot(2,3,1),plot(x,yy1),title(精确值);subplot(2,3,2),plot(x,yy1,t1,y1),title(h=0.5);subplot(2,3,3),plot(x,yy1,t2,y2),title(h=0.1;subplot(2,3,4),plot(x,yy1,t3,y3),title(h=0.05);subplot(2,3,5),plot(x,yy1,t4,y4),title(h=

40、0.01);subplot(2,3,6),plot(x,yy1,t1,y1,t2,y2,t3,y3,t4,y4),title(所有结果);2.求解第二题代码如下：clear all;t0=0;a=0.25;%通过更改a，更改取值区间syms tyy=exp(-50*t);.学习帮手.x=0:0.001:a;yy1=subs(yy,x);h=0.1 0.05 0.01 0.005;n=a./h+1;t1=0:0.1:a;t2=0:0.05:a;t3=0:0.01:a;t4=0:0.005:a;%显式欧拉法y1(1)=100;y2(1)=100;y3(1)=100;y4(1)=100;for j=

41、2:n(1)y1(j)=y1(j-1)+h(1)*(-50)*y1(j-1);endfor j=2:n(2)y2(j)=y2(j-1)+h(2)*(-50)*y2(j-1);endfor j=2:n(3)y3(j)=y3(j-1)+h(3)*(-50)*y3(j-1);endfor j=2:n(4)y4(j)=y4(j-1)+h(4)*(-50)*y4(j-1);endfigure,plot(x,yy1,t1,y1,t2,y2,t3,y3,t4,y4),.gtext(h=0.1),gtext(h=0.05),gtext(h=0.01),gtext(h=0.005);%隐式欧拉法y5(1)=10

42、0;y6(1)=100;y7(1)=100;y8(1)=100;for j=2:n(1)y5(j)=y5(j-1)/(50*h(1)+1);endfor j=2:n(2)y6(j)=y6(j-1)/(50*h(2)+1);endfor j=2:n(3)y7(j)=y7(j-1)/(50*h(3)+1);endfor j=2:n(4)y8(j)=y8(j-1)/(50*h(4)+1);end.学习帮手.figure,plot(x,yy1,t1,y5,t2,y6,t3,y7,t4,y8),.gtext(h=0.1),gtext(h=0.05),gtext(h=0.01),gtext(h=0.005

43、);%相同步长下的稳定性比较figuresubplot(2,2,1),plot(x,yy1,k,t1,y1,b,t1,y5,r),title(h=0.1);subplot(2,2,2),plot(x,yy1,k,t2,y2,b,t2,y6,r),title(h=0.05);subplot(2,2,3),plot(x,yy1,k,t3,y3,b,t3,y7,r),title(h=0.01);subplot(2,2,4),plot(x,yy1,k,t4,y4,b,t4,y8,r),title(h=0.005);3.求解第三题代码如下：%避免t=0作分母，所以用隐式欧拉法clear all;t0=0

44、;h=0.0005 0.0001 0.00005;n=1./h+1;y1(1)=5;y2(1)=5;y3(1)=5;y4(1)=5;t1=0:0.0005:1;t2=0:0.0001:1;t3=0:0.00005:1;for j=2:n(1)y1(j)=y1(j-1)+h(1)*sin(t1(j)/t1(j);endfor j=2:n(2)y2(j)=y2(j-1)+h(2)*sin(t2(j)/t2(j);endfor j=2:n(3)y3(j)=y3(j-1)+h(3)*sin(t3(j)/t3(j);end.学习帮手.format longy(1)=y1(n(1)-y1(1);y(2)=

45、y2(n(2)-y2(1);y(3)=y3(n(3)-y3(1);y实验结果实验结果：1.第一题结果精 确 值h=0.5h=0.13332.52.52.52221.51.51.5111.52111.52111.52h=0.05h=0.01所 以 结 果3332.52.52.52221.51.51.5111.52111.52111.52图 4.1 不同步长下的四级龙格-库塔的解2.显示欧拉法求解结果.学习帮手.20001500h=0.11000500h=0.05h=0.010h=0.005-500-100000.050.10.150.20.25图 4.2 不同步长下的显式欧拉法的解隐式欧拉法求解

46、结果100908070h=0.1605040302010000.050.10.150.20.25h=0.005h=0.05h=0.01.学习帮手.图 4.3 不同步长下的隐式欧拉法的解显式与隐式的稳定性比较结果h=0.120001000100050000-50000.10.2h=0.011001000.30.4-100000.10.2h=0.0050.30.4h=0.05-10005050000.10.20.30.4000.10.20.30.4图 4.4 不同步长下显式与隐式的稳定性比较3.第三题结果表 4.1不同步长下的积分结果步长0.00050.00010.00005积分结果0.94604

47、34318390180.9460751436654500.946079107079003题目要求保证有4位有效数字，所以所求结果为：0.9461结果分析结果分析：四级四阶龙格库塔法具有大范围的收敛性。四级四阶龙格库塔法每计算一步需要计算四次的值，计算有一定的复杂性。当步长为 0.1 和 0.05 时显式欧拉法不具有稳定性。隐式欧拉法在各个步长.学习帮手.下的稳定情况相似，均有较好的稳定性。（第六章第六章）实验题目实验题目：使用列主元高斯消去法解希尔伯特矩阵方程组考察给定的及相应的 时结果的变化，分析其中的原因并给出结论，其中问题分析问题分析：1.希尔伯特矩阵是病态的，的微小变化会造成解发生很大

48、变化。2.刻画矩阵病态程度的量为条件数。条件数越大，矩阵性能越差。可以计算条件数来分析解变化很大的现象。实验原理实验原理：高斯列主元素消元法：将矩阵按列取绝对值最大项进行行行调换，将其以下部分消为，不断循环直到化为上三角的过程。实验内容实验内容：1.高斯列主元素消元法函数代码如下：function A b=GaussElimination(A,b)n=size(A,1);for i=1:n-1q=max(abs(A(i:n,i);for j=i:nif(abs(A(j,i)=q).学习帮手.k=j;endendu=A(i,i:n);v=b(i);A(i,i:n)=A(k,i:n);b(i)=b

49、(k);A(k,i:n)=u;b(k)=v;for l=i+1:nA(l,i)=A(l,i)/A(i,i);b(l)=b(l)-A(l,i)*b(i);for s=i+1:nA(l,s)=A(l,s)-A(l,i)*A(i,s);endendendend2.计算不同阶数下的希尔伯特矩阵的线性方程的解，并通过更改，对其做微小调整，观察解的变化代码如下：clear all;clc;n=5 10 20 50 100;%5阶a=1;H=zeros(n(a);.学习帮手.for i=1:n(a)for j=1:n(a)H(i,j)=1/(i+j-1);endendb=ones(1,n(a);b1=b+0

50、.0001;HH,bb=GaussElimination(H,b);HH1,bb1=GaussElimination(H,b1);for i=n(a):-1:1s=0;s1=0;for j=i+1:n(a)s=s+HH(i,j)*x1(1,j);s1=s1+HH1(i,j)*xx1(1,j);endx1(1,i)=(bb(1,i)-s)/HH(i,i);xx1(1,i)=(bb1(1,i)-s1)/HH1(i,i);end%10阶a=2;H=zeros(n(a);for i=1:n(a)for j=1:n(a)H(i,j)=1/(i+j-1);endendb=ones(1,n(a);b1=b+