丁丽娟《数值计算方法》五章课后实验题答案(源程序很详细-且运行无误)

丁丽娟《数值计算方法》五章课后实验题答案（源程序都是自己写的，很详细，且保证运行无误）我做的五章数值实验作业题目如下：：1、2、3、4题：1、2题：1、2题：2、3题第八章：1、2题第二章1：(1)对A进行列主元素三角分解：function[lu]=myfun(A)n=size(A);fork=1:nfori=k:nsum=0;m=k;forj=1:(k-1)sum=sum+A(i,j)*A(j,k);ends(i)=A(i,k)-sum;ifabs(s(m))<abs(s(i))m=i;endendforj=1:nc=A(m,j);A(m,j)=A(k,j);A(k,j)=c;endforj=k:nsum=0;forr=1:(k-1)sum=sum+A(k,r)*A(r,j);endu(k,j)=A(k,j)-sum;A(k,j)=u(k,j);endfori=1:nl(i,i)=1;endfori=(k+1):nsum=0;forr=1:(k-1)sum=sum+A(i,r)*u(r,k);endl(i,k)=(A(i,k)-sum)/u(k,k);A(i,k)=l(i,k);endend求A的列主元素三角分解：>>A=[11111;12345;1361015;14102035;15153570];>>[L,U]=myfun(A)结果：L=1.000000001.00001.00000001.00000.50001.0000001.00000.75000.75001.000001.00000.25000.7500-1.00001.0000U=1.00001.00001.00001.00001.000004.000014.000034.000069.000000-2.0000-8.0000-20.5000000-0.5000-2.37500000-0.2500(2)求矩阵的逆矩阵A-1:inv(A)结果为：ans=5-1010-51-1030-3519-410-3546-276-519-2717-41-46-41（3）检验结果：E=diag([11111])A\Eans=5-1010-51-1030-3519-410-3546-276-519-2717-41-46-412:程序：functiond=myfun(a,b,c,d,n)fori=2:nl(i)=a(i)/b(i-1);a(i)=l(i);u(i)=b(i)-c(i-1)*a(i);b(i)=u(i);y(i)=d(i)-a(i)*d(i-1);d(i)=y(i);endx(n)=d(n)/b(n);d(n)=x(n);fori=(n-1):-1:1x(i)=(d(i)-c(i)*d(i+1))/b(i);d(i)=x(i);end求各段电流量程序：fori=2:8a(i)=-2;endb=[25555555];c=[-2-2-2-2-2-2-2];V=220;R=27;d=[V/R0000000];n=8;I=myfun(a,b,c,d,n)运行程序得：I=8.14784.07372.03651.01750.50730.25060.11940.04773：（1）求矩阵A和向量b的matlab程序：function[Ab]=myfun(n)fori=1:nX(i)=1+0.1*i;endfori=1:nforj=1:nA(i,j)=X(i)^(j-1);endendfori=1:nb(i)=sum(A(i,:));end求n=5时A1，b1及A1的2-条件数程序运行结果如下：n=5；[A1,b1]=myfun(n)A1=1.00001.10001.21001.33101.46411.00001.20001.44001.72802.07361.00001.30001.69002.19702.85611.00001.40001.96002.74403.84161.00001.50002.25003.37505.0625b1=6.10517.44169.043110.945613.1875cond2=cond(A1,2)cond2=5.3615e+005求n=10时A2，b2及A2的2-条件数程序运行结果如下：n=10;[A2,b2]=myfun(n)A2=1.00001.10001.21001.33101.46411.61051.77161.94872.14362.35791.00001.20001.44001.72802.07362.48832.98603.58324.29985.15981.00001.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.89004.91308.352114.198624.137641.033969.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.000016.000032.000064.0000128.0000256.0000512.0000b2=1.0e+003*0.01590.02600.04260.06980.11330.18160.28660.44510.68011.0230cond2=cond(A2,2)cond2=8.6823e+011求n=20时A3，b3及A3的2-条件数程序运行结果如下：n=20;[A3,b3]=myfun(n)A3=1.0e+009*Columns1through100.00000.00000.00000.00000.00000.00000.00000.00000.00000.00000.00000.00000.00000.00000.00000.00000.00000.00000.00000.00000.00000.00000.00000.00000.00000.00000.00000.00000.00000.00000.00000.00000.00000.00000.00000.00000.00000.00000.00000.00000.00000.00000.00000.00000.00000.00000.00000.00000.00000.00000.00000.00000.00000.00000.00000.00000.00000.00000.00000.00000.00000.00000.00000.00000.00000.00000.00000.00000.00000.00000.00000.00000.00000.00000.00000.00000.00000.00000.00000.00000.00000.00000.00000.00000.00000.00000.00000.00000.00000.00000.00000.00000.00000.00000.00000.00000.00000.00000.00000.00000.00000.00000.00000.00000.00000.00000.00000.00000.00000.00000.00000.00000.00000.00000.00000.00000.00000.00000.00000.00000.00000.00000.00000.00000.00000.00000.00000.00000.00000.00000.00000.00000.00000.00000.00000.00000.00000.00000.00000.00000.00000.00000.00000.00000.00000.00000.00000.00000.00000.00000.00000.00000.00000.00000.00000.00000.00000.00000.00000.00000.00000.00000.00000.00000.00000.00000.00000.00000.00000.00000.00000.00000.00000.00000.00000.00000.00000.00000.00000.00000.00000.00000.00000.00000.00000.00000.00000.00000.00000.00000.00000.00000.00000.00000.00000.00000.00000.00000.00000.0000Columns11through200.00000.00000.00000.00000.00000.00000.00000.00000.00000.00000.00000.00000.00000.00000.00000.00000.00000.00000.00000.00000.00000.00000.00000.00000.00000.00000.00000.00000.00000.00000.00000.00000.00000.00000.00000.00000.00000.00000.00000.00000.00000.00000.00000.00000.00000.00000.00000.00000.00000.00000.00000.00000.00000.00000.00000.00000.00000.00000.00000.00000.00000.00000.00000.00000.00000.00000.00000.00000.00000.00000.00000.00000.00000.00000.00000.00000.00000.00000.00000.00010.00000.00000.00000.00000.00000.00000.00000.00010.00010.00020.00000.00000.00000.00000.00000.00000.00010.00010.00030.00050.00000.00000.00000.00000.00000.00010.00010.00030.00060.00130.00000.00000.00000.00000.00010.00010.00030.00070.00150.00320.00000.00000.00000.00010.00010.00030.00060.00140.00320.00750.00000.00000.00000.00010.00020.00050.00120.00290.00700.01670.00000.00000.00010.00010.00040.00090.00230.00580.01460.03640.00000.00000.00010.00020.00060.00170.00440.01130.02950.07660.00000.00010.00020.00040.00110.00300.00800.02150.05810.15700.00000.00010.00020.00070.00180.00510.01430.04000.11190.31330.00000.00010.00040.00100.00300.00860.02500.07260.21050.61030.00010.00020.00050.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.24940.48740.93161.7434cond2=cond(A3,2)cond2=3.2395e+022由上述运行结果可知：它们是病态的，而且随着n的增大，矩阵的病态变得严重。（2）当n=5时：x1=A1\b1'x1=1.00001.00001.00001.00001.0000当n=10时：x2=A2\b2'x2=1.00001.00001.00001.00010.99991.00001.00001.00001.00001.0000当n=20时：x3=A3\b3'x3=1.0e+005*0.0203-0.17560.7034-1.72282.8742-3.43422.9927-1.87650.7820-0.1396-0.07200.0745-0.03500.0108-0.00230.0003-0.00000.00000.00000.0000由运行结果可见：x1与精确解吻合，x2与精确解稍有差异，x3与精确解差别很大。可见随着n的增大，矩阵病态越来越严重。（3）当n=10时：A2(2,2)=A(2,2)+1e-8;A2(10,10)=A(10,10)+1e-8;x=A2\b2'x=1.01370.91971.20890.68441.30530.80391.08370.97711.00360.9997比较可见，系数矩阵出现微小变动，导致解出现较大变化。说明n=10时，系数矩阵是病态的。4：（1）A=[10787;7565;86109;75910];b=[32233331]';det(A)ans=1cond(A)ans=2.9841e+003eig(A)ans=0.01020.84313.858130.2887(2)A1=[107.28.16.9;7.085.076.025;8.25.899.969.01;6.985.048.979.98];x=[1111]'x1=A1\bx1=0.00772.31171.02111.0157dx=x1-xdx=-0.99231.31170.02110.0157dA=A1-AdA=00.20000.1000-0.10000.08000.07000.020000.2000-0.1100-0.04000.0100-0.02000.0400-0.0300-0.0200根据式（2-39）知：当dA充分小，使得||A-1||*||δA||<1时,则有：norm(dx)/norm(x)ans=0.8225(cond(A)*norm(dA)/norm(A))/(1-cond(A)*norm(dA)/norm(A))ans=-1.0358norm(inv(A))*norm(dA)ans=28.8964显然，上式不成立。显然，原因是因为dA较大，使norm(inv(A))*norm(dA)=28.8964>1（3）dA=0.5*1e-4*rand(4);A1=A+dAA1=10.00007.00008.00007.00007.00005.00006.00005.00008.00006.000010.00009.00007.00005.00009.000010.0000x1=A1\bx1=0.99961.00070.99981.0001dx=x-x1dx=1.0e-003*0.4360-0.67430.1508-0.0952norm(dx)ans=8.2256e-004根据式（2-39）知：当dA充分小，使得||A-1||*||δA||<1时,则有：norm(dx)/norm(x)ans=4.1128e-004(cond(A)*norm(dA)/norm(A))/(1-cond(A)*norm(dA)/norm(A))ans=0.0122norm(inv(A))*norm(dA)ans=0.0121由计算结果可知dA充分小，使得||A-1||*||δA||=0.0121<1时,有：第三章1：用jacobi迭代法：编写jacobi迭代法的m文件如下：functionx1=jacobi(A,b,n,x,e,N)fork=1:Nfori=1:nsum=0;forj=1:nif(j==i)continue;endsum=sum+A(i,j)*x(j);endx1(i)=(b(i)-sum)/A(i,i);endif(norm(x1-x)<e)break;endx=x1;end保存为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)x1=1.0318-2.02972.9451-1.99200.9620即用jacobi迭代法求得解为：[1.0318-2.02972.9451-1.99200.9620]';（2）用Gauss-Seidel迭代法解：编写Gauss-Seidel迭代法的m文件如下：functionx1=gausdel(A,b,n,x,e,N)fork=1:Nsum=0;forj=2:nsum=sum+A(1,j)*x(j);endx1(1)=(b(1)-sum)/A(1,1);fori=2:n-1f=0;g=0;forj=1:i-1f=f+A(i,j)*x1(j);endforj=i+1:ng=g+A(i,j)*x(j);endx1(i)=(b(i)-f-g)/A(i,i);endsum=0;forj=1:n-1sum=sum+A(n,j)*x1(j);endx1(n)=(b(n)-sum)/A(n,n);if(norm(x1-x)<e)break;endx=x1;end保存为gausdel.m文件。然后在matlab命令窗口中编程计算：>>A=[101234;19-12-3;2-173-5;32312-1;4-3-5-115];>>b=[12-2714-1712];>>x=[00000];>>x2=gausdel(A,b,5,x,1e-6,15)x2=1.0055-2.00462.9921-1.99930.9950即用Gauss-Seidel迭代法求得解为：[1.0055-2.00462.9921-1.99930.9950]';（3）用共轭梯度法解：编写共轭梯度法的m文件如下：functionx=gonger(A,b,x0,e)r0=b-(A*x0')';d0=r0;z0=r0*d0'/(d0*(A*d0'));x1=x0+z0*d0;r1=b-(A*x1')';while(norm(r1)>e)r1=b-(A*x1')';m=-r1*(A*d0')/(d0*(A*d0'));d1=r1+m*d0;n=r1*d1'/(d1*(A*d1'));x2=x1+n*d1;d0=d1;x1=x2;endx=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.00001.0000]';2：借用上题编写的共轭梯度法m文件：gonger.m。首先在matlab命令窗口中构造矩阵A,向量b以及初值向量x如下：>>n=1e5;>>fori=1:nx(i)=0;A(i,i)=3;if(i~=n)A(i,i+1)=-1;A(i+1,i)=-1;endif(i~=n/2&&i~=n/2+1)A(i,n+1-i)=0.5;endif(i==1||i==n)b(i)=2.5;elseif(i==n/2||i==n/2+1)b(i)=1;elseb(i)=1.5;endend然后调用上题编写的共轭梯度法解题程序：>>x1=gonger(A,b,x,1e-6)这样，只用这一条指令即可得到结果。第四章1：直接用幂法计算：先编一个M文件如下：function[z,x]=myprounchg(A,x,e,N)k=1;z0=0;z=maxof(x);while(k<N)k=k+1;z0=z;x=A*xz=maxof(x);end保存为myprounchg.m然后在matlab命令窗口中编程计算：>>A=[6-418;20-6-6;22-2211];>>x=[111]';>>[z,x]=myprounchg(A,x,1e-6,10)x=20811x=286286385x=750216944235x=114466114466174361x=33674305563581917971x=525026265250262682941265x=1.0e+009*1.59790.23740.9124x=1.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(1);fori=1:max(size(x))ifabs(m)<abs(x(i))m=x(i);endend保存为maxof.m文件。然后写一个用归一化的幂法计算矩阵特征值与特征向量的程序：function[z,y]=mypro(A,x,e,N)k=1;z0=0;y=x/maxof(x);z=maxof(x);while(k<N&&abs(z-z0)>e)k=k+1;z0=z;x=A*y;z=maxof(x);y=x/z;end保存为mypro.m文件。最后在matlab命令窗口编程计算矩阵A的特征值与特征向量：>>A=[6-418;20-6-6;22-2211];>>x=[111]';>>[z,x]=mypro(A,x,1e-6,10)z=19.2540x=1.00000.14430.5713即求得矩阵A的特征值为：19.2540；特征向量为[10.14430.5713]。2：借助与上题所编mypro.m文件、maxof.m文件；首先编用反幂法解题的m文件：function[ay]=fanmi(A,a0,x,e,N)k=1;a1=1;B=A-a0*eye(size(A));y=x/maxof(x);x=B\y;u=maxof(x);a=a0+1/u;while(k<N&&abs(a-a1)>e)y=x/maxof(x);x=B\y;u=maxof(x);a=a0+1/u;k=k+1;end然后在matlab命令窗口编程解题：>>A=[126-6;6162;-6216];>>x=[111]';>>[a0x]=mypro(A,x,1e-10,3);>>[ay]=fanmi(A,a0,x,1e-10,20)a=21.5440y=1.00000.7838-0.8069得到矩阵A的按模最大特征值的更精确的近似值：21.544。其中程序：>>[a0x]=mypro(A,x,1e-10,3);用幂法迭代3次来得到A的按模最大值特征值的近似值作为下面反幂法程序的输入。第六章2：>>x=[2356791011121416171920];>>y=[106.42108.26109.58109.5109.86110109.93110.59110.6110.72110.9110.76111.1111.3];>>Y=1./y;>>X=1./x;>>size(X)ans=114>>A=[14sum(X);sum(X)sum(X.^2)];>>b=[sum(Y);X*Y'];>>a=B\ba=0.00900.0008或直接用下面指令则更为简便：>>a=polyfit(X,Y,1)a=0.00080.0090即得到作拟合曲线图：>>x1=2:0.1:20;>>X1=1./x1;>>Y1=polyval(a,X1);>>plot(X,Y,'o',X1,Y1,'r',X,Y,'b')得到下图：3：先编一个M文件：functiony=fun(a,xi)y=a(1)*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,[1112],xi,yi)a=2.65071.86861.52362.0558于是得：a=2.6507，b=1.8686，c=1.5236，d=2.0558作图显示：>>x1=0.1:0.02:1;