数值分析实验报告

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

2、性。迭代 格式为?)?+1 = ?- ?(?), ?= 0,1,2,割线法：是牛顿法的改进，具有超线性收敛的特性，收敛阶为1.618.迭代格式为?)?+1=?- ?-?-1)(?- ?)，?i2,?)(?- ?2? =? =史蒂芬森迭代法：采用不动点迭代进行预估校正。至少是平方收敛的。迭 代格式为?+1 = ?-'"+1- ?- 2?+ ?这里？?可采用牛顿迭代法的迭代函数。实验内容:1. 写出该问题的 ?(?) 函数 代码如下： function py= f(x) syms k;y=(kA2+1)*(k-1)A5;yy=diff(y,k);py(1)=subs(y,k,x)

3、; py(2)=subs(yy,k,x);end2. 分别写出各个迭代法的迭代函数 代码如下：二分法：function y=dichotomie(a,b,e)i=2;m(1)=a;while abs(a-b)>et=(a+b)/2;s1=f(a);s2=f(b);s3=f(t);if s1(1)*s3(1)<=0b=t;elsea=t;endm(i)=t;i=i+1;endy=t,i+1,m;end牛顿迭代法：functiony=NewtonIterative(x,e)i=2;m(1)=x;while abs(en)>=es=f(x);t=x-s(1)/s(2);en=t-x

4、;x=t;m(i)=t;i=i+1;endy=x,i+1,m;end牛顿割线法：function y=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;enden=2*e;y=x2,i+1,m;endz=fai(y);t=x-(y-x)A2/(z-2*y+x);史蒂芬森迭代法：en=t-x;Functionx=t;p=StephensonIterative(x,e)i=2;m(i)=t

5、; i=i+1;m(2)=x;enden=2*e;p=x,i+1,m;while abs(en)>=eendy=fai(x);3. 因为 ?(?)经常被使用，故可以写一个 代码如下：?(?)函数。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=Secant(0,2, 1e-10);x4=Stephens

6、onIterative(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,

7、1),plot(x1(4:x1(2)-1)-x1(1)./(x1(3:x1(2)-2)-x1(1),tit le( ' 二分法 ' );subplot(2,2,2),plot(x2(4:x2(2)-1)-x2(1)./(x2(3:x2(2)-2)-x2(1),tit le( ' 牛顿迭代法 ' ); subplot(2,2,3),plot(x3(4:x3(2)-1)-x3(1)./(x3(3:x3(2)-2)-x3(1),tit le( ' 牛顿割线法 ' ); subplot(2,2,4),plot(x4(4:x4(2)-1)-x4(1)./(

8、x4(3:x4(2)-2)-x4(1),tit le( ' 史蒂芬森迭代法 ' );%不同初始值，相同迭代法下的收敛性x5=dichotomie(-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,subp

9、lot(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:x7(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( ' 牛顿迭

10、代法 ' );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=StephensonIterative(-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

11、,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),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.各个迭代值分布0.5牛顿迭代法0 050

12、100史蒂芬森迭代法1.5 0.5 -I f I 0 EE102468图1.1不同迭代法下的得到的迭代值迭代值的情况如下:二分法牛顿迭代法牛顿割线法史蒂芬森迭代法00001.50000000000.20000000002.00000000001.35555555560.75000000000.37049180320.33333333330.98161652831.12500000000.50764420760.38071968010.99994600030.93750000000.61461894470.49828334190.99999999951.03125000000.697386909

13、80.57049963330.98437500000.76155380910.63938062441.00781250000.81154111860.69427858790.99609375000.85067638570.74116926531.00195312500.88144821230.78027159970.99902343750.90572974000.8132927871当二分法的初始区间选为0,3，误差限为1 X1Q-10，牛顿迭代法初值选为0,误差限为1 X1O-10 ,牛顿割线法初始点为0,2,误差限为1 X10-5 ,史蒂芬森迭代法初始点选为0,误差限为1 X1O-10，迭

14、代情况如下列图。迭代次数分别为38次，100次，140次，9次。故而，史蒂芬森迭代法速度最快，效果最好2收敛情况牛顿迭代法【IJi5r9x 10分法0-5-100.8 0.7，0.6 -0.5 -0.4 c050100010203040史蒂芬森迭代法0.4 0.2-0'-0.2-0.4 *c"12345图1.2不同迭代法下迭代值得收敛情况二分法收敛效果较差，牛顿迭代法和牛顿割线法相近，史蒂芬森迭代法收敛 次数高于1,效果最好3. 不同初值的收敛情况543210-1-2-3-450100150图1.3二分法，牛顿迭代法下不同初值的收敛情况图1.4牛顿割线法，史蒂芬森迭代法下不同

15、初值的收敛情况1. 二分法的五个初始区间分别为0,3,-1,1 ,-2,3 ,0,4,-4,4 ;2. 牛顿迭代法的五个初始值分别为0,-2, -4,4,6 ;3. 牛顿割线法的五个初始区间分别为0,2,-1,1 ,-4,5 ,0,7,-8,2 ;4. 史蒂芬森迭代法的五个初始值分别为0,-1, -4,4,6 ;由图可知，它们最终均到达收敛。收敛性分析及结论：1. 二分法收敛较慢且不能求解崇根，但算法简单；此处牛顿法具有了平方收敛； 从迭代次数上看，牛顿割线法较牛顿法的多，所以收敛性较差，是超线性收 敛；史蒂芬森迭代法收敛效果最好。2. 因为牛顿迭代法是局部的二次收敛，所以要注重初值的选取，本

16、次实验中选择的初值均得到了收敛，效果比拟好。牛顿割线法也应注意初值的选取。第三章实验题目：1. 区间-1,1作等距划分：2 ?= -1 + ? ? = 一 ?= 0 1 ？ ? ?以？?为结点对函数?= 梟进行插值逼近。(1) 分别取？ = 1,5,10,20,25用牛顿插值对？?进行逼近，并在同一坐标系下做 出函数的图形，进行比拟。写出插值函数对？?的逼近程度与节点个数的 关系，并分析原因；(2) 试用三次样条插值对？?进行逼近，在同一坐标下画出图形，观察样条插 值函数对？?的逼近程度与节点个数的关系；(3) 整体插值有何局限性？如何防止？2. 一组数据如下，求其拟合曲线.表2.1数据表?0

17、12345678910?23478101114161819?106.42108.2109.5110109.93110.49110.59110.6110.76111111.2(1)求以上数据形如？?？= ?+ ?1?+ ?的拟合曲线及其平方误差;?(2) 求以上数据形如？?？=?的拟合曲线及其平方误差;(3) 通过观察12的结果，写出你对数据拟合的认识问题分析：题目除上述要求之外还有以下几点：1. 明确整体插值和分段插值的不同。牛顿插值多项式属于整体插值，三次样条 插值属于低次分段插值。2. 将结果在同一坐标下绘制出。但是为了方便分析节点个数对于插值效果的影 响，也可以单独绘制。3. 第二题中为

18、了确定各个参数的大小，可以进行适当变换，转化为线性，运用 最小二乘法，得到拟合。实验原理：牛顿插值多项式：对于给定的插值节点2? ?<?<? < ? < ?构造次数不超过 ?的插值多项式? ?( ?)? = ?0?+ ?1?( ? ?- ?0?)+ ?2(?- ?0?)( ?- ?1?)+?+ ?(? ?)(? ?)？ ?-1) 使其满足插值条件?( ?)? = ?= ?( ?)? ?=? 0,1, , ? 这样得到的插值多项式 ? ?( ? ?) 称为 ?插?值?多?项式。系数 ?为差商，可以 通过构造差商表得到。三次样条插值：三次样条插值函数 S(x) 在每个小区间

19、 ?-1 , ? ? 上为三次多 项式；在全进 ( ?, ?) 上存在二阶连续导数； 其次符合插值条件。 Matlab 中存在内 置的三次样条插值函数，命令为 ? .实验内容：第一题：1. 牛顿插值函数的构造代码如下： function f=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;endendf=y(

20、1,1);for j=2:SZff=y(1,j);for i=1:j-1ff=ff*(x-x0(i);endf=f+ff;endend2. 牛顿和三次样条插值情况及比拟： 代码如下： clear all ;clc;syms x;fx=1/(5+xA2);%牛顿多项式插值 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=su

21、bs(fx,x2); f2=Newton(x2,y2); 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(

22、' 所有结果 ' );x6=-1:h2(1):1; y6=subs(fx,x6); 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=

23、subs(fx,x10); f10=Newton(x10,y10); y010=subs(f10,x0);figure, plot(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( y12=spline(x2,y2,x0

24、);subplot(2,3,3),plot(x0,y12),title( y13=spline(x3,y3,x0);subplot(2,3,4),plot(x0,y13),title( y14=spline(x4,y4,x0);subplot(2,3,5),plot(x0,y14),title( y15=spline(x5,y5,x0);subplot(2,3,6),plot(x0,y15),title( figure,plot(x0,y0,x0,y11,x0,y12,x0,y13,x0,y14,x0,y15),title( 第二题：代码如下：clear all ;clc;x=2 3 4 7

25、8 10 11 14 16 18 19;y=106.42 108.2 109.5 110 109.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 on y01=polyval(c1,x);delta1_2

26、=sum(y-yO1).A2)龙格现象 ' );'n=1' );'n=5' );'n=10' );'n=20' );'n=25' );所有结果 ' );%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(2);a bf=a*exp(-b/s);x2=x1;y2=subs(f,x2);plot(x2,y2, 'b' )

27、,gtext( ' 曲线二 ' );y02=subs(f,x);delta2_2=sum(y-yO2).A2)实验结果：第一题：1. 牛顿插值结果所有结果图2.1不同节点数的牛顿插值多项式综合图龙格现象图2.2龙格现象由图可知，2.三次样条插值结果0.20.150.10.050-505精确图-1n=110-505-505-5-5-5图2.3不同节点下的三次样条插值结果所有结果图2.4不同节点数的三次样条的综合图由图可知，牛顿插值多项式在？较小的时候如5,10,20,25差值效果良好，当 ?变大时如60,62,65,67,100,200就出现了龙格现象，三次样条在各个子区间内 为

28、三次多项式，拟合效果好第二题1.第一问的多项式拟合得到的拟合曲线为?= 106.2927 + 0.6264?- 0.0205?平方误差为备=2.7796第二问的拟合曲线为0.0903? = 111.4940?平方误差为& = 0.47192.拟合曲线如下列图1121111101091081071061121111101091081071062468 10121416 1820*曲线二曲线8 10121416 1820图2.5拟合情况从图中可以看出曲线二大体符合黑点的分布情况，拟合效果较好。结论：整体插值随着节点个数的增多，多项式的次数也在升高。高次多项式的插值 会出现龙格现象，震荡剧烈

29、，逼近效果不理想。但是当节点很多时，低次插值的 精度又不够，所以为了防止这一局限性采用分段低次插值。 其中三次样条插值有 良好的收敛性和光滑性，效果较好。数据拟合时，可以选择趋势相似的函数形式， 在求解相关参量时，将函数进 行适当变换，从而运用最小二乘法得到拟合结果。第四章实验题目：设计区间分半求积算法、龙贝格求积算法和自适应辛普森求积算法的程序，观察？?= 1,10,100,500时，积分?cos(?)?的结果，并给出相应的评价问题分析:1. 分半求积算法即为复合梯度求积。2. 因为？?越大，?= ?在 -1,1 内震荡地越严重，自适应辛普森求 积算法很难适用，计算复杂度会很高。所以选取辛普

30、森求积算法代替。实验原理：分半求积算法：将区间?等分为？?个小区间，分点为? ?= ?+ ? ? ?= 0 1 ？?'利用定积分性质，有99 4? -1 ?+1广? X 广 ?(?)? ?=0 - ?每个小区间上均用梯形求积公式，有?212 ?(?o?-1?-1?/ ? 刀?(?) + ?(?知)+ 刀-? ?=0 ?=0?-1?-1?= 2 X ?(?)+ ?c?+1) = ?(?+ ?(?)+ ? X ?(?=0 ?=0即为复合梯形公式。龙贝格求积算法：将步长为？的复合梯形公式进行査尔逊外推，就得到龙贝格求积公式。自适应辛普森求积算法：按照被奇函数在区间上的变化来安排求积结点的辛普

31、森算法。实验内容：1. 写得各个求积算法的函数代码如下： 1分半求积算法：xk=a+k*h;functionyk=subs(y,xk);f=CompositeTrapezoida(n)f(i)=0;syms x;for j=1:m(i)y=x*x*cos(n*x);f(i)=f(i)+(yk(2*j-1)+4*yk(2*j)+ym=1:100;k(2*j+1);a=-1;endb=1;f(i)=h/3*f(i);for i=1:100endh=(b-a)/m(i);m;f'k=0:m(i);endxk=a+k*h;4龙贝格求积算法yk=subs(y,xk);f(i)=0;functio

32、n f=Romberg(n) epsilon=0.0001;for j=1:m(i)f(i)=f(i)+(yk(j)+yk(j+1);syms x;endy=x*x*cos(n*x);a=-1;f(i)=h/2*f(i);b=1;endh=b-a;end2辛普森算法fa=subs(y,a);function f=CompositeSimpson(n)fb=subs(y,b);T(1)=h/2*(fa+fb);syms x;m=1;y=x*x*cos(n*x);while 1m=1:100;h=h/2;a=b=1;k=1:(2Am-1);xk=a+k*h;for i=1:100h=(b-a)/(

33、2*m(i);yk=subs(y,xk);k=0:(2*m(i);su=sum(yk);S(1)=h/2*(fa+fb)+h*su;for i=1:mS(i+1)=S(i)+(S(i)-T(i)/(4"-1);endif (abs(S(m+1)-T(m)<epsilon)break ;endT=S;m=m+1;endf=S(m+1);end5自适应辛普森求积算法 function f=AdaptiveSimpson(n) epsilon=0.001;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=

34、0;f=0;while 1n1=length(ep);n=length(p0);if (n=1) break ; endh=p0(2)-p0(1);2. 各个求积公式下的积分结果 代码如下 :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)<15*ep(1)f=f+s1+s2;p0=p0(2:n);if n1>=2ep=ep(2:n1);endq=q+1;elsem=m

35、+1; p0=yk(1),yk(3),p0(2:n);if n1=1ep=ep(1)/2 ep(1)/2;elseep=ep(1)/2 ep(1)/2ep(2:n1);endif q=0p=p0;elsep=p(1:q) p0;endendendend%复合梯形法clear all ;clc;y11=CompositeTrapezoida(1);y12=CompositeTrapezoida(10);y13=CompositeTrapezoida(100);y14=CompositeTrapezoida(500);1:100;y11;y12;y13;y14'%复合辛普森法y21=Com

36、positeSimps on( 1);y22=CompositeSimpso n(10);y23=CompositeSimpso n(100);y24=CompositeSimpso n(500);2:2:200;y21;y22;y23;y24'%龙贝格法实验结果：y31=Romberg(1);y32=Romberg(10);y33=Romberg(100);y34=Romberg(500);y31 y32 y33 y34%自适应辛普森法y41=AdaptiveSimps on (1)%以下复杂度太高不岀结果%y42=AdaptiveSimpso n(10);%y43=Adaptive

37、Simpso n(100);%y44=AdaptiveSimpso n(500);%y41 y42 y43 y44表3.1不同积分法得到的结果n110100500复合梯形法0.4782831994-0.1399401910-0.00601370080.0023823806复合辛普森法0.4782672530-0.1401909997-0.00985113540.0072247738龙贝格法0.4782672574-06112212333-0.2492135619结果评价：龙贝格法的精度高，最接近精确解。因为被积分函数为?= ?随着？的增大，震荡越剧烈。所以积分面积越来

38、越小。第五章实验题目:1. 分别取步长？= 0.5,0.1,0.05,0.01,用龙格-库塔法求解y' = ? ?1) = 1 ,3计算到t= 2 ,并与精确解?)?=警2比拟，观察龙格-库塔法的收敛性.2. 分别取步长？= 0.1,0.05,0.01,0.005,用显示欧拉法和隐式欧拉法求解y'=-50y , y(0) = 100 ,由结果分析算法的稳定性.3. 选择某常微分方程初值问题的数值方法计算驾？?的近似值，并保证有4位有效数字.问题分析：1.龙格-库塔有很多的格式，第一题以常用的四级四阶龙格库塔方法为例。» 一 ? ? 一2假设,?= -?, 由 于sin

39、 t第三题可以转化为求解常微分方程？-?= ?-的问题实验原理:四级四阶龙格-库塔格式?+1 = ?+ (? + 2? + 2? + ?)6? = ?(?-?(? 2，?+2?)-?(?+ 2，?+ 2?)?(?+ ?, ?+ ?)?=?=?=实验内容:1.求解第一题代码如下：clear all ;clc %以四级四阶龙格库塔格式为例 to=i； h=0.5 0.1 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.0 5:2;t4=1:0.01:2;for j=2: n(1)k1=t

40、1(j-1)*y1(j-1)A(1/3);k2=(t1(j-1)+h(1)/2)*(y1(j-1)+(h(1)/2)*k1)A(1/3);k3=(t1(j-1)+h(1)/2)*(y1(j-1)+(h(1)/2)*k2)A(1/3);k4=t1(j)*(y1(j-1)+h(1)*k3)A(1/3);y1(j)=y1(j-1)+h(1)/6*(k1+2*k2+2*k3+k4);endfor j=2:n(2)k仁 t2(j-1)*y2(j-1)A(1/3);k2=(t2(j-1)+h(2)/2)*(y2(j-1)+(h(2)/2)*k1)A(1/3);k3=(t2(j-1)+h(2)/2)*(y2

41、(j-1)+(h(2)/2)*k2)A(1/3); k4=t2(j)*(y2(j-1)+h(2)*k3)A(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)A(1/3);k2=(t3(j-1)+h(3)/2)*(y3(j-1)+(h(3)/2)*k1)A(1/3);k3=(t3(j-1)+h(3)/2)*(y3(j-1)+(h(3)/2)*k2)A(1/3); k4=t3(j)*(y3(j-1)+h(3)*k3)A(1/3);y3(j)=y3(j-1)+h(3)/6*(k1+2*k2+2*

42、k3+k4);endfor j=2:n(4)k1=t4(j-1)*y4(j-1)A(1/3);k2=(t4(j-1)+h(4)/2)*(y4(j-1)+(h(4)/2)*k1)A(1/3);k3=(t4(j-1)+h(4)/2)*(y4(j-1)+(h(4)/2)*k2)A(1/3); k4=t4(j)*(y4(j-1)+h(4)*k3)A(1/3);y4(j)=y4(j-1)+h(4)/6*(k1+2*k2+2*k3+k4);end x=1:0.01:2;syms tyy=(tA2+2)/3)A(3/2); yy1=subs(yy,x);figure,');'h=0.5

43、9;);'h=0.1''h=0.05');'h=0.01');所有结果 ');subplot(2,3,1),plot(x,yy1),title(' 精确值subplot(2,3,2),plot(x,yy1,t1,y1),title(subplot(2,3,3),plot(x,yy1,t2,y2),title( subplot(2,3,4),plot(x,yy1,t3,y3),title( subplot(2,3,5),plot(x,yy1,t4,y4),title(subplot(2,3,6),plot(x,yy1,t1,y1,t

44、2,y2,t3,y3,t4,y4),title(2. 求解第二题代码如下：clear all t0=0;a=0.25;%通过更改 a ，更改取值区间syms t yy=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;gtext( 'h=0.1'),gtext('h=0.05'),gte

45、xt('h=0.01'),gtext('h=0.005');forj=2:n(1)y1(j)=y1(j-1)+h(1)*(-50)*y1(j-1);endforj=2:n(2)y2(j)=y2(j-1)+h(2)*(-50)*y2(j-1);endforj=2:n(3)y3(j)=y3(j-1)+h(3)*(-50)*y3(j-1);endforj=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),%隐式欧拉法y5(1)=100;y6(1)=1

46、00;y7(1)=100;y8(1)=100;forj=2:n(1)y5(j)=y5(j-1)/(50*h(1)+1);endforj=2:n(2)y6(j)=y6(j-1)/(50*h(2)+1);endforj=2:n(3)y7(j)=y7(j-1)/(50*h(3)+1);endforj=2:n(4)y8(j)=y8(j-1)/(50*h(4)+1);endfigure,plot(x,yy1,t1,y5,t2,y6,t3,y7,t4,y8),gtext( 'h=0.1' ),gtext('h=0.05'),gtext('h=0.01'),g

47、text('h=0.005'%相同步长下的稳定性比拟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,

48、'r' ),title('h=0.01'););););'h=0.005'；subplot(2,2,4),plot(x,yy1,'k',t4,y4,'b',t4,y8,'r' ),title(3. 求解第三题代码如下：%防止t=0作分母，所以用隐式欧拉法clear all ;t0=0;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;f

49、orj=2:n (1)y1(j)=y1(j-1)+h(1)*si n(t1(j)/t1(j);endforj=2:n (2)y2(j)=y2(j-1)+h (2) *si n(t2(j)/t2(j);endforj=2:n (3)y3(j)=y3(j-1)+h (3) *si n(t3(j)/t3(j);endformat longy(1)=y1( n(1)-y1(1);y(2)=y2( n( 2)-y2(1);y(3)=y3( n(3)-y3(1);y实验结果：1. 第一题结果3-3/3/2.5- -2.52.5/jV/2P/J2 /2/1.5/|1.51.5/1/: 11/精确值h=0.5

50、h=0.11.51.51.5h=0.01所以结果h=0.0511.52图4.1不同步长下的四级龙格-库塔的解2.显示欧拉法求解结果20001500h=0.11000500-500h=0.05h=0.01 h=0.005-1000 L00.050.10.150.20.25图4.2不同步长下的显式欧拉法的解隐式欧拉法求解结果90807060504030h=0.1h=0.05h=0.01I1h=0.00520 -10 -0匚00.050.10.150.20.25图4.3不同步长下的隐式欧拉法的解显式与隐式的稳定性比拟结果1 w1000' 1- 5000+ f1V-500r-1000,L-h=

51、0.1h=0.05200010000-10000.10.20.30.40.10.20.30.4h=0.01100 h=0.005100 0 L00.10.20.30.450 -000.10.20.30.4图4.4不同步长下显式与隐式的稳定性比拟3. 第三题结果表4.1不同步长下的积分结果步长？0.00050.00010.00005积分结果0.9460434318390180.9460751436654500.946079107079003题目要求保证有4位有效数字，所以所求结果为：0.9461 结果分析：四级四阶龙格库塔法具有大范围的收敛性。四级四阶龙格库塔法每计算一步需要计算四次的值，计算有

52、一定的复杂性。当步长？为0.1和0.05时显式欧拉法不具有稳定性。隐式欧拉法在各个步长 下的稳定情况相似，均有较好的稳定性。第六章实验题目：使用列主元高斯消去法解希尔伯特矩阵？?= (? X方程组？?= ?考察给定的？?= 5,10,20,50,100及相应的？寸结果的变化，分析其中的原因并给出结 论，其中?= 1/(?+ ? 1)问题分析：1. 希尔伯特矩阵是病态的，？的微小变化会造成解发生很大变化。2. 刻画矩阵病态程度的量为条件数。条件数越大，矩阵性能越差。可以计算条 件数来分析解变化很大的现象。实验原理：高斯列主元素消元法：将矩阵？?按列取绝对值最大项进行行行调换，将其以 下局部消为0,不断循环直到？?化为上三角的过程。实验内容：1. 高斯列主元素消元法函数代码如下：fun cti onA b=Gauss