数值计算方法上机实习题

1、数值计算方法上机实习题（1）由递推公式1=-511+（,从Io=O.1824,Io=O.1823出发，计算1刀；（2）仆=0，LEOOOO,用1=一；1一+，计算I。：(1)由I计算I2e递推子程序：Rinctionf=fib(n,i)ifn>=lf=fib(n-l,i)*(-5)+(l/(n);elseifn=Of=l,end计算和显示程序：1=01824;11=01823.fibl=fib(20,I);fib2=fib(20,Il).fynntfCI_O=O1824Hj.l_20%dn：6bl),4)nntfCL0=01823时.IJOWdSfb?).计算结果显示：I_0=01824

2、时.IJ0=7480927e+09I1=0.1823时,I20=-2.055816e+09(2)由计算Ic程序：n=21,11=0,12=10000.fl=11.f2=i2,whilen-=0tfl=flaGl/5)+(l/(5an);f2=f2s(-l/5)*(l/(5an).n=n-l,endpnntfCIJO=O时，】_0=%4).mnntf(L?0=10000时，I_0=%4.8fiT,f2),计算结果显示：IJ0=0时，L0=018232156I20=10C00时，10=018232156（3）分析结果的可靠性及产生此现象的原因（重点分析原因）答：第一个算法可得出eo=Ro-/ol%

3、=4一/；1=5啜。|易知第一个算法每一步计算都把误差放大了5倍，n次计算后更是放大了5n倍，可靠性低。第二个算法可得出可以看出第二个算法每一步计算就把误差缩小5倍，n次后缩小了5n倍，可靠性高。计算根与步数程序：fylot(x)exp(x)+10*x-2,0,l),gndon.symsx.f=exp(x)+10ax-2troot,n=EFF3(f；0,l),jnntfCroot=%68f,n=%dn：rooLn),计算结果显示：root=009057617,n=ll(2)取初值=0,并用迭代看出2-ex102.求方程eX+10x2=0的近似根，要求,日一%|5乂107,并比较计算量。(1)在

4、0,1上用二分法：(I)0,1上的二分法二分法子程序：functionroot,n=EFF3(txkx2)%第二题(1)二分法fl=subs(£symvar(D，xl),%函数在x=xl的值f2=subs(£symvar(f),x2),%x=x2n=0,%步数er=5T0八4,%误层if(fl=O)ioot=xl,return.elsaf2=0)ioot=x2,return,eham*包>0)dispC两端点函数值乘枳大于0!).return.elsewhile(abs(xlx2)>er)%循环x3=(xl+x2)/2tf3=subs(£symvar(

5、f).x3),n=n+l,if(f3=0)root=x3,break.elsaf(frf3>0)xl=x3.elsex2=x3,endendioot=(xl+x2)/2.%while循环少一步需加上end(2)初值*0迭代迭代法子程序：Hinctionroot>n=DDF(g.xO,err,max)(接F页)计算根与步数程序：symsx,f=(2-exp(x)/10,(接下页)%root根，n+1步数，g函数，xO初值，err误片.max最大迭代次数X(l)=xO,forn=2maxX(n)=subs(g,symvar(g).X(n-l)tc=abs(X(n)-X(n-l);roo

6、t=X(n)tif(c<err)break.endifn=inaxdispC超出预设迭代次数)；endendn=nl,%步数减1xO=O;err=5s10A(-4);max=1001root,n=DDF(CxO,err.max)tfynntf(Yodt=%6,8f,n=%dn*,root,n)»计算结果显示：root=009051262,n=4（3）加速迭代的结果,(3)加速迭代加速迭代计算程序：xO=O,en=510八(Y),max=100>symsx,g=x-(f(x)-x)A2/(f(f(x)-2*f(x)-x),root,Q=DDF(g,xO,err,max)tf

7、ynntf(*root=%6,8Cn=%d*,root,n),程序函数设置：functiong=f(x)g=(2-exp(x)/10.计算结果显示：root=0.09048375,n=2C）取初值5=0,并用牛顿迭代法;(4)牛顿迭代法牛顿迭代法子程序：functionroot,n=NDDDFx(g,x0.err.max)%root根，n+1步数,g函数，xO初值，err误度，max最大迭代次数X(l)=xO;forn=2maxX(n)=subs(g.symvar(g).X(n-l)tc=abs(X(n)-X(n-l)»root=X(n);if(c<err)|(root=0)b

8、reak.endifn=maxdispC超出预设迭代次数)；(接下页)牛顿迭代法计算程序：x0=0;en=5s10A(-4);max=100.symsx,g=exp(x)+10ax-2tg=x-(g/diff(g)trootn=NDDDFx(g,xO»err,inax)tfynntfCroot=%6.81n=%dSjooLn)：计算结果显示：root=009052511,n=2endendn=nl,%步数减1(5)分析绝对误差。答：可以看到，在同一精度卜.，二分法运算了11次，题设迭代算式下运算了4次，加速迭代卜运算了2次，牛顿迭代下运算了2次。因不动点迭代法和二分法都是线性收敛的，

9、但二分法压缩系数比题设迭代方法大，收敛速度较慢。加速迭代速度是超线性收敛，牛顿法是二阶，收敛速度快。3.钢水包使用次数多以后，钢包的容积增大，数据如下:X3456789y642829.58959710993999101112131415161049105910601081061091076ax+b试从中找出使用次数和容枳之间的关系，计算均方差。(用y二拟合)c+x拟合曲线程序：x=216.产64282958959.710993999104910.5910601081061091076,£fval=fininsearch(deltaL0,0,0,optimset.x,y)%finmse

10、arch为求解多元函数的最小值函数%f为多元函数初位x0附近的极小值点%fval为极小值k=2100.yl=(f(l)/k+f(2)./(fi(3)+k),figure1=figureCColor,111),gxt=plot(xy,*xfkylk，).legend原数据？拟合曲线','Location','northwest'),%y为数据点连接曲线，yl为拟合曲线titleC函数y=(ax+b)/(x+c)的拟合,FontSizW14FontWeight?Bold),xlabdC次数？FontSize,14,FontWeight,Bold,),yiab

11、dC容积:FontSize；14F(mtWeight?Bol(r),setCgxt/LmeWidth1.15);gridcnt%计算均方差for1=1:15构造函数子程序：Rmctiondelta=deltal(txy)a=f(l),b=f(2),c=f(3),delta=0.fork=l15delta=delta+(x(k)+c)sy(k)-(a*x(k)+b)A2.end计算结果显示：拟介出的方程为：(x+07110)y=11.2657x+-155024均方差为033165089总结：指标选择，闪虺设方程I=E巴%一为非线性的，要转化为线性方程故需提指标为：y2(1)二(f(l)伊好)/(

12、的七)，(接卜页)endJ=o.for1=1:15j=j+(y(i)-y2(i)r3,endjfc=sqrt(j/15).SnntfC拟合出的方程为(x)%44f)y=%44&+%44fn均方差为48fn4Kl),f(2),jfc),(ex+=ax+b=E：3Kc9+Dy：-ax.-b1其驻点方程为：dlO元=Z2(cxf+l)yf-axi-b(-x)=0；i-iisl£=£2(5+1)比一ax,-b=0；次数-10-100-104.设A=-104-10-10-14-10-1(00-14-10-10-10-104-1T4,0、4J-2，Ax=b4J-2Jacobi迭

13、代计算程序：A=4-10-100?14-10-10,0-14-10-1;-10-14-10t0-10-14-1.00-10-14J.b=05-25-261,x0=0000001,q)s=10A(-4);maxl=200.X,k=Jacobifl2(A>b,xO>eps,maxl)t3nmf('近似解x=n%48fn%48fn%48fn%48fn%48fn%48fn迭代次数n=%dn',X(l),X(2),X(3).X(4).X(5).X(6),k),计算结果显示：近似解x=0999981771.999950180999975091999950180999975091

14、99996353迭代次数n=28GAUSS-SEIDEL迭代计算程序：A=(4-10-100,-14-10-10,0-14-10-1;.10-14-10.0-10-14.1.00.10.14J.b=05-25-26x0=(000000;eps=l0A(-4),maxl=200,X,k=GSDD(A,b,xO,q)s,maxl)t用nmf('近似解x=n%4Sfn%4Sfn%48f分析卜.列迭代法的收敛性，并求410Y的近似解及相应的迭代次数。(1) JACOBI迭代；(1)Jacobi迭代Jacobi迭代子程序：Rinctionxk=Jacobifl2(A.b.x0.eps.maxl)

15、%jacobi按矩阵的分量算法%x0初值，eps误差，maxi最大迭代次数n=length(xO),fork=lmaxix(l)=(b(l)-A(l,2n)-xO(2n)/A(l,l),fori=2:nx(i)=(b(i)-A(i,li-l)*xO(ln)*xO(i+l:n)/A(i.i)tendifnonn(x*-xO»2)<epsreturnelsexO=xendend(2)GAUSS-SEIDEL迭代：(2) GAUSISEIDEL迭代GAUSSSEIDEL迭代子程序：Rinctionx,k=GSDD(A,b,xO,q)s,maxl)n=length(xO),fork=l

16、maxix(l)=(b(l)-A(l,2n)axO(2n)')/A(l,l),fori=2:nx(i)=(b(i)-A(iJxl)*x(l»iyA(u+ln)sx0(i+ln)»A(i,i),endifnorm(x-x0»2)<epsn%48fn%48fn%48fn迭代次数n=%dn；X,X(?),X,X(4),X(5),X(6),k),return计算结果显示：else近似解x=x0=x.099996014end199995676end099996351199996662099997253199998401迭代次数n=15(3) SOR迭代(取3=

17、0.l:01:l.9,找到迭代步数最少的“)(3)SOR迭代SOR迭代子程序：functionx,km=SORl(A,b,x0,qjs.maxl)n=length(xO),a=xO;fors=l19工0=贰w(s)=s/10.x()=x06w，(b.A。/n).x0(Ln»A(LD，fori=2:nx(i)=xO(i)+w(s)a(b(i)-A(i.1i-l)sx(li-iy-AO.in)axO(in)/A(1.0.endk(s)=O;whilenorm(x-xOt2)>=epsxO=x,x(l>xOOw(s)*(bA(1/n)sxO(lnn/ACU),fori=2:nx

18、(i)=xO(i)+w(s)*(b(i)A(i,lil)，x(ln)xO(in)")/A(1.1).endk(s)=k(s)+l;if(k(s)>=maxl)dispC超出迭代次数上限！n),return,SOR迭代计算程序：A=4-10-100,-14-10-10,0-14-10-1,-10-14-10.0-10-14-1.00-10-14jtb=05-25-261;xOHOOOOOO.eps=10A(-4),maxl=1000,Xk,m=SORl(Ab,xO,eps,maxl),fynntf('n近似解x=n%48fn%48fn%4,8fn%48fn%48fn%48

19、f：x),fynntf(1n步数最小omega=%21£步数ndngk),计算结果显示：近似解x=100000977199999782099999567200001071099999662199999923步数最小omega=l20,步数n=9.总结：Jacobi迭代次数为28次：GAUSS-SEIDEL迭代次数为15次：SOR迭代次数在w=l2时达到最小次数9次。系数矩阵严格对角占优，则Jacobi迭代法与Gauss-Seidel迭代法均收敛：乂因系数矩阵对称正定，H0<w<2,则SOR迭代法也收敛.三种迭代法的结果分析及比较：(1)JacobiendendData(s

20、)=x.endminl.index=min(k)tx=Data(index)»k=nunl>m=w(index)tfynntf(1n解x=n%48M%48fn%48fn%48fn%48fn%48f,x),4jrmtf(*n步数最小omega=%2.2£步数n=%d。"mk).法,其设计思想是将线性方程组的求解归结为对角方程组求解过程的重复,计算规则简单而易于编写计算程序.但收敛速度慢,谱半径相对偏大(2)Gauss-Seidel充分利用新信息进行计算,其迭代的效果比Jacobi迭代好，谱半径比Jacobi法小。(3)SOR迭代法，计算公式简单、编制程序容易,

21、是求解大型稀疏方程组的一种有效方法,旦若松弛因子选择得当超松弛法收敛速度比GaussSeidel迭代和Jacobi迭代都要快。63r5.用逆案迭代法求人=321最接近于11的特征值和特征向量，准确到10一3。逆嘉迭代法子程序：functionlanda,Vk=NMSF2(A,v,eps,p)n,n=size(A).A=A-p*eye(n).vO=v,tinax,tindex=max(abs(vO),Iamd0=v0(tindex)tuO=vO/laindO.flag=0,k=0,while(flag=O)V=AuO.tinax,tindex=max(abs(V)>lamdl=V(tmde

22、x),uO=V/lamdltif(abs(lamdO)A(-l)-(laindl)A(-l)<=epsflag=l;endlamdO=lamdl,k=k+l;end(接下页)lamda=(lamdl)A(-l)+ptV=uO逆累迭代计算程序：A=(631,321,111,V=l11,eps=0001,p=ll,lamda.Ufk=NMSF2(A.xcps,P)；用nntfiC最接近11的特征值为:%48Al特征向量:n%4Sfn%48fn%4.8fn',lamda.U(l)tU(2)tU(3)t计算结果显示：最接近于11的特征值为：787313712特征向量：1000000000

23、54922698022545618J1V6 .用经典R-K方法求解初值问题(1)y=-2yi+y?+2siiixy；=Yj-2y2+2cosx-2sinxxgOJO，|Yi(0)=2ly2(o)=3xe0J0,fyi(o)=2ly2(o)=3y=-2y1+y2+2sinxy；=998-999y2+999cosx-999sinx和精确解J%(x)=2e，+smx比较，进行误差分析得到结论，图形显示精确解和数值解。y2(x)=2ex+cosx经典四阶R-K方法R-K程序：住(工y)2*y(l)+y(2)+2*sin(x),y(l)-2*y(2)+2*cos(x).2*sin(x),g=(x,y)-

24、2*y(l)+y(2)+2*sin(x);998*y(l)-999*y(2)+999*cos(x)-999*sm(x),(接下页)x=00110,Yl=2*exp(x)+sm(x),Y2=2*exp(x)+cos(x);xl=010/1000010,Y3=2'exp(xl)+sm(xl),Y4=2#exp(-xl)+cos(xl)tN=100,Nl=10000.al,bl=ode23(00110,23),a2,b2=ode23(g,0:l/100010,23);a3,b3=ode45(C00110t23);a4,b4=ode45(&0:l/100010t23)tn=length

25、(bl),figurefigureCColofJl11);dbt=plot(x>bl(:J)、,wbl(:,2),kHYl,gHY2,y-'),legend(yi：y2;精确解Yl?精确解Y2);titleC四阶RK算法下的方程组的数值解&方程的精确解？FontSize；14FontWeight?Bold),yiabd(*YTontSize*,14,TontWeight1,'Bold;xiabdCXVFontSizeUlTontWeightVBoldO；set(dbt»,LineWidth,>15)»gndon,Q=bl(:,l).YlW

26、=bl(:,2).Y2,Ql=b2(.l)-Y3,.Wl=b2(2)-Y4,tQ2=b3C，l)Y,W2=b3(:,2>Y2,Q3=b4(.J)-Y3W3=b4(>2).Y4,>ERl=max(abs(Q)tER2=max(abs(W),ER3=max(abs(Ql),ER4=max(abs(Wl)tER5=max(abs(Q2).ER6=max(abs(W2).fpnntK,(1)中yl在ode23下求出的最大误差为：%4Sfn(1)中y2在ode23下求出的最大误差为：4Sfn；ERl,ER2),grrntfCW(2)中yl在ode23下求出的大误屋为：%4Sfn(2)中

27、y2在ode23下求出的最大误差为：%48fn：ER3,ER4),SrmtRAn(1)中yl在ode45下求出的最大误圣为：%4Sfn(1)中y2在ode45下求出的最大误差为：%48fn：ER5.ER6);fynnt*n(2)中yl在ode45F求出的最大误星为:%4.8fn(2)中y2在。de45下求出的最大误差为:%48MER7.ER8);计算与显示程序：(1)中yl在。de23F求出的最大误差为:000169850(D中y2在ode23下求出的最大误差为：000207422(2)中yl在。de23卜求出的最大误片为：000000583(2)中y2在。de23下求出的最大误差为：000581329(1)中yl在。de45F求出的最大误差为:000052155(1中y2在ode45下求出的最大误差为：000045793(2)中yl在o加45卜求出的最大误是为：000000276(2)中y2在。de45下求出的最大误差为:000275331总结：对比2种方程在。加23,。加45命令下,与精确值的误差可以发现，在此膻中，。加45