《丁丽娟《数值计算方法》五章课后实验题答案(源程序很详细,且运行无误).docx》由会员分享，可在线阅读，更多相关《丁丽娟《数值计算方法》五章课后实验题答案(源程序很详细,且运行无误).docx（42页珍藏版）》请在三一办公上搜索。
1、丁丽娟数值计算方法五章课后实验题答案(源程序都是自己写的，很详细，且保证运行无误)我做的五章数值实验作业题目如下：第二章：1、2、3、4题第三章：1、2题第四章：1、2题第六章：2、3题第八章：1、2题第二章(D对A进行列主元素三角分解：function1u=myfun（八）n=size（八）;fork=l:nfori=k:nsum=0;m=k;forj=k(k-l)sum=sum+A(i,j)*A(j,k);ends(i)=A(i,k)-sum;ifabs(s(m)A=11111;12345;1361015;14102035;15153570;L,U=myfun（八）结果：L=1.00000
2、0001.00001.00000001.00000.50001.0000001.00000.75000.75001.000001.00000.25000.7500-1.00001.0000O4.000014.000034.000069.000000-2.0000-8.0000-20.5000000-0.5000-2.37500000-0.2500(2)求矩阵的逆矩阵A1:inv（八）结果为：ans=5-1010-51-1030-3519-410-3546-276-519-2717-41-46-41(3)检验结果：E=diag(lIlll)AEans=5-1010-51-519-2717-41-
3、46-41程序：functiond=myfun(a,b,c,d,n)fori=2:nl(i)=a(i)b(i-l);a(i)=l(i);u(i)=b(i)-c(i-l)*a(i);b(i)=u(i);y(i)=d(i)-a(i)*d(i-l);d(i)=y(i);endx(n)=d(n)b(n);d(n)=x(n);fori=(n-l):-l:lx(i)=(d(i)-c(i)*d(i+l)b(i);d(i)=x(i);end求各段电流量程序：a(i)=-2;endb=25555555;c=-2-2-2-2-2-2-2;V=220;R=27;d=VROOOOOOO;n=8;I=myfun(a,b
4、,c,d,n)运行程序得：I=8.14784.07372.03651.01750.50730.25060.11940.0477(1)求矩阵A和向量b的matlab程序:functionAb=myfun(n)fori=l:nX(i)=l0.1*i;endfori=l:nforj=l:nendendfori=l:nb(i)=sum(A(i,:);end求n=5时Al,bl及Al的2.条件数程序运行结果如下:n=5；Al,bl=myfun(n)Al=1.00001.10001.21001.33101.46411.00001.20001.44001.72802.07361.00001.30001.69
5、002.19702.85611.00001.40001.96002.74403.84161.00001.50002.25003.37505.0625bl=6.10517.44169.043110.945613.1875cond2=cond(A1,2)求n=10时A2,b2及A2的2.条件数程序运行结果如下:n=10;A2,b2=myfun(n)A2=1.0000.ooo1.21001.33IO1.46411.61051.77161.94872.14362.35791.00001.2(XX)1.44(X)1.72802.07362.4B832.98603.58324.29985.15981.00
6、001.30001.69002.19702.85613.71294.82686.27498.157310.60451.00001.40001.96002.74403.84165.37827.529510.541414.757920.66101.00001.50002.25003.37505.06257.593811.390617.085925.628938.44341.00001.60002.56004.09606.553610.485816.777226.843542.949768.71951.00001.70002.894.91308.352114.198624.137641.033969
7、.7576118.58791.00001.80003.24005.832010.497618.895734.012261.2220110.1996198.35931.00001.90003.61006.859013.032124.761047.045989.3872169.8356322.68771.00002.00004.00008.0000I6.OO32.000064.00I28.OO256.00512.0000b2=1.0e+003*0.01590.02600.04260.06980.11330.18160.28660.44510.68011.0230cond2=cond(A2,2)co
8、nd2=8.6823e+011求n=20时A3,b3及A3的2-条件数程序运行结果如下：n=20;A3,b3=myfun(n)A3=0.00000.00000.00000.00000.000.0000O.(XX)O0.00000.000.000.000.00000.000.00000.00000.00000.000.00000.00000.000.00000.000.000.00000.000.000.00000.00000.00000.00000.00000.000.00000.00000.000.00000.00000.000.00000.00000.000.00000.000.00000
9、.00000.00000.000.00000.00000.00000.0000O.(XX)O0.00000.0000O.O(XX)0.00O.(XX)O0.00000.00000.000.00000.000.000.00000.000.000.00000.00000.00000.00000.00000.000.00000.00000.000.00000.00000.000.00000.00000.00000.00000.00000.00000.000.0000O.O(XX)O.(XX)O0.00000.00000.000.00000.000.00000.00000.00000.000.0000
10、0.00000.000.00000.000.000.00000.000.00000.000.00000.00000.000.00000.000.00000.00000.000.00000.00000.000.00000.00000.00000.00000.00000.00000.000.0000O.O(XX)O.(XX)O0.00000.00000.00000.000.00000.00000.0000O.(XXX)0.0000O.(XX)O0.00O.O(XX)0.00000.000.000.00000.000.000.00000.00000.00000.00000.000.00000.000
11、.00000.00000.00000.000.00000.00000.00000.00000.00000.00000.00000.000.0000O.O(XX)O.(XX)O0.00000.00000.0000O.(XX)O0.00000.00000.00O.(XXX)0.0000O.(XX)O0.000.00000.00000.000.000.00O.O(XX)0.000.000.000.00000.000.000.00000.000.000.00000.00000.000.00000.00000.0000Columns11through200.00000.00000.000.000.000
12、.00000.00000.00000.000.00000.0000O.(XX)O0.00000.00000.00O.(XXX)0.00000.00000.000.00000.00000.000.000.00000.000.000.00000.00000.00000.00000.000.00000.000.00000.000.0000O.O(XX)0.00000.00000.00000.00000.00000.000.00000.000.000.000.00000.000.00000.000.0000O.(XX)O0.00000.00000.000.000.000.00000.00000.000
13、00.00000.000.000.000.00000.000.00000.00000.00000.000.00000.000.00000.000.0000O.O(XX)0.00000.00000.00010.00000.00000.000.00000.000.000.000.0l0.0l0.00020.000.0000O.(XX)O0.00000.000.000.00010.00010.00030.050.00000.00000.000.000.000.0l0.00010.00030.00060.130.000.00000.000.00000.0l0.0l0.00030.00070.150.3
14、20.000.00000.000.00010.0l0.030.0006O.I40.00320.750.000.00O.(XX)O0.00010.020.050.00120.00290.700.01670.00000.00000.0l0.0l0.040.090.00230.00580.01460.03640.000.00000.0l0.00020.00060.00170.440.01130.02950.07660.0000O.(XX)I0.020.00040.00110.00300.00800.02150.05810.15700.0000().(XX)10.00()20.0007().(X)18
15、0.(X)5I0.0143().(X)0.1119().3133O.O(XX)O.(X)O10.00040.00100.00300.860.02500.07260.21050.61030.0l0.00020.(XX)50.00160.00480.01430.04300.12910.38741.1623b3=1.0e+009*Columns1through100.00000.00000.00000.00000.00000.00000.00010.00020.00040.0010Columns11through200.00250.00590.01320.02870.06060.12460.2494
16、0.48740.93161.7434cond2=cond(A3,2)cond2=3.2395e+022由上述运行结果可知：它们是病态的，而且随着n的增大，矩阵的病态变得严重。(2)当n=5时：xl=Albxl=1.00001.00001.0000x2=A2b2,x2=1.00001.00001.00001.00010.99991.00001.00001.00001.00001.0000当n=20时:x3=A3b3,x3=1.0e+005*0.0203-0.17560.7034-3.43422.9927-1.87650.7820-0.1396-0.07200.0745-0.03500.0108-
17、0.00230.0003-0.00000.00000.00000.0000由运行结果可见：Xl与精确解吻合，X2与精确解稍有差异，X3与精确解差别很大。可见随着n的增大，矩阵病态越来越严重。(3)当n=10时:A2(2,2)=A(2,2)+le-8;A2(10,10)=A(10,10)le-8;1.01370.91971.20890.68441.30530.80391.08370.97711.00360.9997比较可见，系数矩阵出现微小变动，导致解出现较大变化。说明n=10时，系数矩阵是病态的。(1)A=10787;7565;86109;75910;b=32233331;det（八）ans=
18、1cond（八）ans=2.9841e+003eig（八）ans=0.01020.84313.858130.2887(2)Al=107.28.16.9;7.085.076.025;8.25.899.969.01;6.985.048.979.98;x=l11I1xl=Albxl=0.00772.31171.02111.0157dx=xl-xdx=-0.99231.31170.0157dA=Al-AdA=C)0.20000.1000-0.10000.08000.07000.020000.2000-0.1100-0.04000.0100-0.02000.0400-0.0300-0.0200l3x2与
19、HAIb之间的关系:IlXll2ia2根据式(2-39)知:当dA充分小,使得IIAjII*3AIlVl时,则有:C小JA2Cond（八）lx2-IlAIl2raIbIIAII2norm(dx)norm(x)ans=0.8225(cond(A)*norm(dA)norm（八）)(1-cond(A)*norm(dA)norm(A)ans=-1.0358norm(inv（八）)*norm(dA)ans=28.8964显然，上式不成立。显然，原因是因为dA较大，使norm(inv（八）)*norm(dA)=28.8964l(3) dA=0.5*le-4*rand(4);Al=A+dAAl=10.00
20、007.00008.00007.00007.00005.00006.00005.00008.00006.000010.00009.00007.00005.00009.000010.0000xl=AlbXl=0.99961.00070.99981.0001dx=x-xldx=1.0e-003*0.4360-0.67430.1508-0.0952norm(dx)ans=8.2256e-0043x2与Il3A2之间的关系:IlXl2ia2根据式(2-39)知:当dA充分小，使得IIAJlI*AKl时,则有:C.JlAIl2Cona（八）lx2IlAIl2IIAII2norm(dx)norm(x)an
21、s=4.1128e-004(cond（八）*norm(dA)norm（八）)(1-cond（八）*norm(dA)norm（八）)ans=0.0122norm(inv（八）)*norm(dA)ans=0.0121由计算结果可知dA充分小，使得IIA-III*A=0.012kl时,有:Cz7JIa2condA)IlRb=4.1128e-0040.0122=IlA2碇1aA)UA.IlaIl2第三章(1)用jacobi迭代法：编写jacobi迭代法的m文件如下：functionx1=jacobi(A,b,n,x,e,N)fork=l:Nfori=l:nsum=0;forj=knif(j=i)con
22、tinue;endSUm=SUm+A(iJ)*x(j);endx1(i)=(b(i)-sum)A(i,i);endif(norm(xl-x)e)break;endend保存为jacobi.m文件。然后在matlab命令窗口中编程计算：A=101234;19-12-3;2-173-5;32312-1;4-3-5-115;b=12-2714-1712;x=00000;X1=jacobi(A,b,5,x,1e-6,15)xl=1.0318-2.02972.9451-1.99200.9620-1.9920即用jacobi迭代法求得解为：1.0318-2.02972.94510.9620;(2)用GaU
23、SS-Seidel迭代法解：编写Gauss-Seidel迭代法的m文件如下：functionx1=gausdel(A,b,n,x,e,N)fork=l:Nsum=0;forj=2:nsum=sum+A(l,j)*x(j);endx1(1)=(b(1)-sum)A(1,1);fori=2:n-lf=O;g=O;forj=ki-lf=f+A(i,j)*xl(j);endforj=i+kng=g+A(iJ)*x(j);endxl(i)=(b(i)-f-g)A(i,i);endsum=0;forj=l:n-lsum=sum+A(nj)*x1(j);endx1(n)=(b(n)-sum)A(n,n);i
24、f(norm(x1-x)e)rl=b-(A*x)t;m=-rl*(A*dO)(dO*(A*dO);dl=rl+m*d;n=rl*dl7(dl*(A*d);x2=xl+n*dl;d0=dl;xl=x2;x=x1;end保存为gonger.m文件。然后在matlab命令窗口中编程计算：A=101234;19-12-3;2-173-5;32312-1;4-3-5-115;b=12-2714-1712;x=00000;x3=gonger(A,b,x,1e-6)x3=1.0000-2.00003.0000-2.00001.0000即用共辗梯度法求得解为：1.0000-2.00003.0000-2.000
25、01.0000;借用上题编写的共朝梯度法m文件：gonger.mo首先在matlab命令窗口中构造矩阵A,向量b以及初值向量X如下:n=le5;fori=l:nx(i)=0;A(i,i)=3;if(i=n)A(ij+1)=-1;A(i+l,i)=-l;endif(i=n/2&i=n/2+1)A(i,nl-i)=0.5;endif(i=li=n)b(i)=2.5;elseif(i=n2i=n2+l)b(i)=l;elseb(i)=1.5;endend然后调用上题编写的共规梯度法解题程序:X1=gonger(A,b,x,1e-6)这样，只用这一条指令即可得到结果。第四章(1)直接用基法计算：先编一
26、个M文件如下：functionz,x=myprounchg(A,x,e,N)k=l;z0=0;z=maxof(x);while(kN)k=k+l;z=z;x=A*xz=maxof(x);end保存为myprounchg.m然后在matlab命令窗口中编程计算:A=6-418;20-6-6;22-2211;x=l111;z,x=myprounchg(A,x,1e-6,10)X=20811X=286286385750216944235114466114466174361336743055635819179715250262652502626829412651.0e+009*1.59790.2374
27、0.91241.0e+010*2.50612.50613.9968X=1.0e+011*7.69551.11044.3965z=7.6955e+011x=1.0e+011*7.69551.11044.3965由以上计算结果可知：直接用嘉法计算矩阵A的特征值和特征向量,得不到正确结果。原因：因为实际计算时每次迭代所求得的向量没有进行归一化处理，使计算过程出现了溢出。(2)用归一化的基法计算：先写一个向量中求按模最大值的程序：functionm=maxof(x)m=x(l);fori=kmax(size(x)ifabs(m)abs(x(i)m=x(i);endend保存为maxof.m文件。然后写
28、一个用归一化的基法计算矩阵特征值与特征向量的程序:functionz,y=mypro(A,x,e,N)k=1;z0=0;y=x/maxof(x);z=maxof(x);while(ke)k=k+1;z=z;x=A*y;z=maxof(x);y=xz;end保存为mypro.m文件。最后在matlab命令窗口编程计算矩阵A的特征值与特征向量:A=6-418;20-6-6;22-2211;x=l1I1;z,x=mypro(A,x,1e-6,10)1.00000.14430.5713即求得矩阵A的特征值为：19.2540；特征向量为10.14430.5713。借助与上题所编mypro.m文件、max
29、of.m文件;首先编用反幕法解题的m文件：functionay=fanmi(A,aO,x,c,N)k=1;a1=1;B=A-aO*eye(size（八）);y=xmaxof(x);x=By;u=maxof(x);a=a+lu;while(ke)y=xmaxof(x);x=By;u=maxof(x);a=a+lu;k=k+l;end然后在matlab命令窗口编程解题:A=126-6;6162;-6216;x=l1lt;ax=mypro(A,x,1e-10,3)；ay=fanmi(A,a0,x,1e-10,20)a=21.5440y=1.000.7838-0.8069得到矩阵A的按模最大特征值的更
30、精确的近似值：21.544o其中程序：IaOX=111)1*0(人丛1匕-10,3);用累法迭代3次来得到A的按模最大值特征值的近似值作为下面反幕法程序的输入。第八早x=2356791011121416171920;y=106.42108.26109.58109.5109.86110109.93110.59110.6110.72110.9110.76111.1111.3;Y=ly;X=Lx;size(X)ans=A=14sum(X);sum(X)sum(X.2);b=sum(Y);X*Y;a=Bba=0.00900.0008或直接用下面指令则更为简便：a=polyfit(X,Y,l)a=0.0
31、0080.0090即得到L=0.(XX)81+0009yX作拟合曲线图：xl=2:0,1:20;Xl=lxl;Yl=polyval(a,X1);plot(X,Y,o,X1,Y1,T,X,Y,b)得到下图：先编一个M文件：functiony=fun(a,xi)y=a(l)*xi+a(2)*(xi.2).*exp(-a(3)*xi)+a(4);保存为fun.m然后在命令窗口编程：xi=0.1:0.1:1;yi=2.32402.64702.97073.28853.60083.90904.21474.51914.82325.1275;a=lsqcurvefit(fun,l112,xi,yi)2.650
32、71.86861.52362.0558于是得：a=2.6507,b=1.8686,c=1.5236,d=2.0558作图显示：X1=0.1:0.02:1;size(xl)ans=146fori=l:46yl(i)=2.6507*xl(i)1.8686*xl(i)2*exp(-1.5236*xl(i)+2.0558;endplot(xi,yi,o,xl,ylr,xi,yi,b,)可见，拟合函数f(x)拟合的非常好。第八章由方程可知，在v-2或x2时，显然方程f(x)0o所以作图时取X的区间为:2,2方程作图：x=-2:0.1:2;y=x.2+sin(10*x)-1;plot(x,y,o,x,y,
33、b,)gridon由图可知，有根区间为-1.5.5。确定根的近似值：functionfy=ygqLj(f,x,a,b)j=l；fori=a:0.1:bxl=subs(f,x,i);x2=subs(f,x,i+0.1);ifxl*x2le-10y(i)=i;j=ji;endend保存为:ygqj.m。确定根的近似值：symsxy=x2+sin(10*x)-l;b=ygqj(y,x,-2,2)b=-1.5000-1.4000-1.0000-0.6000-0.50000.60000.90001.2000编写牛顿迭代的m文件：functionx1,err=newton(f,x,xO,e,N)fori=
34、l:NxI=XO-SUbS(f,x,x)/SUbS(diff(f),x,x);err=min(abs(x1-x),abs(x1-x)abs(x1);iferrereturn;endx=xl;i=i+l;end保存为：newton.m求方程的所有解：fori=kmax(size(b)x1,err=newton(y,x,b(i),1e-6,16)endxl=-1.4142err=1.1969e-009xl=-1.4142err=1.0212e-009xl=err=3.7299e-012xl=-0.5513err=1.5263e-007xl=-0.5513err=3.6513e-007xl=0.68
35、44err=5.5799e-012xl=0.9287err=2.6598e-012xl=err=6.1208e-008则方程的所有解为：xl=-1.4142,x2=-1.4142,x3=-0.9519,x4=-0.5513x5=-0.5513,x6=0.6844,x7=0.9287,x8=1.2087o(1)先编写迭代的m文件：functionx1,err=diedai(f,x,e,N)fori=l:Nxl=feval(f,x);err=abs(xl-x);iferrereturn;endx=x1;i=i+l;end然后在matlab命令窗口中编程计算：f=inline(,exp(-x),x);x1,err=diedai(f,0.5,1e-8,30)xl=0.5671err=得到方程的近似解：x=0.5671,精确度为7.7332e-009o(2)用加速收敛的迭代格式计算：fl=inline(,0.625*exp(-x)+0.375*xx);x1,err=diedai(f1,0.5,1e-8,10)xl=0.5671err=4.7402e-009得到方程的近似解为：x=0.5671,精确度为4.7402e-009o
