丁麗娟《數(shù)值計算方法》五章課后實驗題答案(源程序很詳細,且運行無誤)

丁麗娟《數(shù)值計算方法》五章課后實驗題答案（源程序都是自丁麗娟《數(shù)值計算方法》五章課后實驗題答案（源程序都是自己寫的，很詳細，且保證運行無誤）己寫的，很詳細，且保證運行無誤）我做的五章數(shù)值實驗作業(yè)題目如下：第二章：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)結(jié)果：L=1.0000 0 0 0 01.0000 1.0000 0 0 01.0000 0.5000 1.0000 0 01.0000 0.7500 0.7500 1.0000 01.0000 0.2500 0.7500 -1.0000 1.0000U=1.0000 1.0000 1.0000 1.0000 1.00000 4.0000 14.0000 34.0000 69.00000 0 -2.0000 -8.0000 -20.50000 0 0 -0.5000 -2.37500 0 0 0 -0.2500(2)inv(A)結(jié)果為：ans=5-1010-51-1030-3519-410-3546-276-519-2717-41-46-41（3）檢驗結(jié)果：E=diag([11111])A\Eans=5-1010-51-1030-3519-410-3546-276-519-2717-42:1-46-41程序：functiond=myfun(a,b,c,d,n)fori=2:nl(i)=a(i)/b(i-1);a(i)=l(i);b(i)=u(i);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）Abmatlabfunction[Ab]=myfun(n)fori=1:nX(i)=1+0.1*i;endfori=1:nforj=1:nA(i,j)=X(i)^(j-1);endfori=1:nn=5A1，b1A12-條件數(shù)程序運行結(jié)果如下：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.1051 7.4416 9.0431 10.9456 13.1875cond2=cond(A1,2)cond2=5.3615e+005n=10A2，b2A22-條件數(shù)程序運行結(jié)果如下：n=10;1.00001.10001.21001.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+011n=20A3，b3A32-條件數(shù)程序運行結(jié)果如下：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.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0001 0.0002 0.0004 0.0010Columns11through200.0025 0.0059 0.0132 0.0287 0.0606 0.1246 0.2494 0.4874 0.9316 1.7434cond2=cond(A3,2)cond2=3.2395e+022n的增大，矩陣的病態(tài)變得嚴(yán)重。（2）n=5時：x1=A1\b1'x1=

1.00001.00001.00001.00001.0000n=10x2=A2\b2'x2=1.00001.00001.00001.00010.99991.00001.00001.00001.00001.0000n=20x3=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由運行結(jié)果可見：x1與精確解吻合，x2與精確解稍有差異，x3與精確解差別很大?？梢婋S著n的增大，矩陣病態(tài)越來越嚴(yán)重。（3）n=10A2(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.9997n=10時，系數(shù)矩陣是病態(tài)的。44：（1）A=[10787;7565;86109;75910];b=[32233331]';ans=1ans=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.00771.0157dx=x1-xdx=-0.99230.0157dA=A1-AdA=0 0.2000 0.1000 -0.10000.08000.07000.020000.2000-0.1100-0.04000.0100-0.02000.0400-0.0300-0.0200||x||2

與||2

之間的關(guān)系：x||2 A||2根據(jù)式（2-39）知：當(dāng)dA充分小，使得||A-1||*||δA||<1時,則有：cond(||2||2

||A||2x||2 1cond(||2A||2norm(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||x||2

與||2

之間的關(guān)系：x||2 A||2根據(jù)式（2-39）知：當(dāng)dA充分小，使得||A-1||*||δA||<1時,則有：cond(||2||2

||A||2x||2 1cond(||2A||2norm(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由計算結(jié)果可知dA充分小，使得||A-1||*||δA||=0.0121<1時,有：

=-

cond(||2A2x||2 1cond(||2||A||2第三章11：jacobi迭代法：jacobim文件如下：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.0297 2.9451 -1.9920 0.9620即用jacobi迭代法求得解為：[1.0318 -2.0297 2.94510.9620]';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-1

-1.9920f=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;forsum=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.0046 2.9921 -1.9993 0.9950即用Gauss-Seidel迭代法求得解為：[1.0055 -2.0046-1.9993 0.9950]';用共軛梯度法解：mfunctionx=gonger(A,b,x0,e)r0=b-(A*x0')';d0=r0;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;end

2.9921x=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.0000 3.0000 -2.0000 1.0000即用共軛梯度法求得解為：[1.0000 -2.0000 3.0000 -2.00001.0000]';22：借用上題編寫的共軛梯度法m文件：gonger.m。首先在matlab命令窗口中構(gòu)造矩陣A,向量b以及初值向量x如下：>>n=1e5;>>fori=1:nx(i)=0;A(i,i)=3;if(i~=n)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然后調(diào)用上題編寫的共軛梯度法解題程序：>>x1=gonger(A,b,x,1e-6)第四章11：（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=208x=286286385x=750216944235x=174361x=33674305563581917971x=525026265250262682941265x=1.0e+009*1.59790.23740.9124x=1.0e+010*2.50612.50613.9968x=1.0e+011*7.69554.3965z=7.6955e+011x=1.0e+011*7.69554.3965A的特征值和特征向量，得不到正確結(jié)果。使計算過程出現(xiàn)了溢出。（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]。22：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);3A的按模最大值特征值的近似值作為下面反冪法程序的輸入。第六章22：>>x=[2356791011121416171920];>>y=[106.42108.26109.58109.5109.86110109.93110.59110.6110.72110.9110.76111.1111.3];>>>>X=1./x;>>size(X)ans=1 14>>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.0008 0.00901y

=

1+0.009x>>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.3240 2.6470 2.9707 3.2885 3.6008 3.90904.2147 4.5191 4.8232 5.1275];>>a=lsqcurvefit(@fun,[1112],xi,yi)a=2.6507 1.8686 1.5236 2.0558于是得：a=2.6507，b=1.8686，c=1.5236，d=2.0558作圖顯示：>>x1=0.1:0.02:1;>>size(x1)ans=1 46>>fori=1:46y1(i)=2.6507*x1(i)+1.8686*x1(i)^2*exp(-1.5236*x1(i))+2.0558;end>>plot(xi,yi,'o',x1,y1,'r',xi,yi,'b')可見，擬合函數(shù)f（x）擬合的非常好。第八章1：x<-2x>2f(x)>0x的區(qū)間為:[-2,2]。方程作圖：>>x=-2: