矩陣上機(jī)作業(yè)講解

1、矩陣與數(shù)值分析課程數(shù)值實驗報告1. 方程 在 x=3.0 附近有根，試寫出其三種不同的等價形 式以構(gòu)成兩種不同的迭代格式，再用這兩種迭代求根，并繪制誤差下降曲線，觀察這兩 種迭代是否收斂及收斂的快慢newton 迭代法function x=juz1(x1,tol,max)xo=x1;T=;i=1;while imaxxn=xo-(2*xo3-5*xo2-19*xo+42)/(6*xo2-10*xo-19);T(i)=xn;if abs(xn-xo) juz1(3,10-5,20)迭代法收斂3.5000弦截法function x=juz11(x1,x2,tol,max)xo=x1;xp=x2;i

2、=1;while imaxxn=xp-(2*xp3-5*xp2-19*xp+42)/(2*xp2+2*xp*xo+2*xo2-5*xp-5*xo-19); T(i)=xn;if abs(xn-xp) juz11(2.9,3,10-5,20) 迭代法收斂3.500082. 用復(fù)化梯形公式、復(fù)化辛普森公式、龍貝格公式求下列定積分，要求絕對誤差為，并將計算結(jié)果與精確解進(jìn)行比較：1）2）1）用復(fù)化梯形公式求積分function f=Txing1(a,b,n)% a，b 分別為積分下限和上限， while n104n 為劃分的區(qū)間個數(shù)h1=(b-a)/n; %h1 為積分區(qū)間長度 x1=zeros(1,

3、n+1);for i=1:n+1x1(i)=a+(i-1)*h1;end%計算 n+1 個節(jié)點(diǎn)儲存在 x1 中y1=(2/3)*(x1.3).*exp(x1.2);to=0;for i=1:n%計算每個節(jié)點(diǎn)對應(yīng)的函數(shù)值to=to+h1/2*(y1(i)+y1(i+1);end%用梯形公式求積分h2=(b-a)/(n+1); %h2 為積分區(qū)間長度 x2=zeros(1,n+2);for i=1:n+2x2(i)=a+(i-1)*h2;end%計算 n+2 個節(jié)點(diǎn)儲存在 x2 中y2=(2/3)*(x2.3).*exp(x2.2);tn=0;%計算每個節(jié)點(diǎn)對應(yīng)的函數(shù)值for i=1:n+1tn=

4、tn+h2/2*(y2(i)+y2(i+1);end%用梯形公式求積分if abs(tn-to) Txing1(1,2,5)I =54.59825121-5.0604e-06 復(fù)化辛普森求積分 function f=Simpson1(a,b,n) %a,b 分別為積分下限和上限 %n 為分割區(qū)間個數(shù) while n108 T=; for n=n:n+1 h1=(b-a)/n;%將積分區(qū)間x1=zeros(1,n+1); for i=1:n+1%計算 n+1x1(i)=a+(i-1)*h1;endz1=(2/3)*x1.3.*exp(x1.2); y1=zeros(1,n);for i=1:ny

5、1(i)=a+(2*i-1)*h1/2;endr1=(2/3)*y1.3.*exp(y1.2);t=0; for i=1:n;公式求積分t=t+h1/6*(z1(i)+4*r1(i)+z1(i+1); end T=T,t;endif abs(T(2)-T(1) Simpson1(1,2,3)54.5982N =160-2.9586e-08龍貝格公式求積分function f=romberg1(a,b) t(1,1)=(b-a)/2*(2*a3*exp(a2)/3+2*b3*exp(b2)/3);n=1;while n100for k=1:naa=0;for i=0:(power(2,(k-1)

6、-1);aa=aa+2*(a+(2*i+1)*(b-a)/2k)3*exp(a+(2*i+1)*(b-a)/2k)2)/3; endt(1,k+1)=0.5*t(1,k)+(b-a)/(2k)*aa;endfor m=1:nh=1/(4m)-1);for k=m:nt(m+1,k+1)=(4m)*t(m,k+1)-t(m,k)*h;endendif (abs(t(n,n)-t(n+1,1+n) romberg1(1,2)54.59822）復(fù)化梯形公式求積分function f=Txing2(a,b,n)% a，b 分別為積分下限和上限， while n10000n 為劃分的區(qū)間個數(shù)h1=(b-

7、a)/n; %h 為積分區(qū)間長度 x1=zeros(1,n+1);for i=1:n+1x1(i)=a+(i-1)*h1;end%計算 n+1 個節(jié)點(diǎn)儲存在 x1 中y1=2*x1./(x1.2-3);to=0;for i=1:n%計算每個節(jié)點(diǎn)對應(yīng)的函數(shù)值to=to+h1/2*(y1(i)+y1(i+1);end%用梯形公式求積分h2=(b-a)/(n+1); %h2 為積分區(qū)間長度 x2=zeros(1,n+2);for i=1:n+2x2(i)=a+(i-1)*h2;end%計算 n+2 個節(jié)點(diǎn)儲存在 x2 中y2=2*x2./(x2.2-3);tn=0;%計算每個節(jié)點(diǎn)對應(yīng)的函數(shù)值for

8、i=1:n+1tn=tn+h2/2*(y2(i)+y2(i+1);end%用梯形公式求積分if abs(tn-to) Txing2(2,3,4)I =1.79181025r =-1.0576e-06 復(fù)化辛普森公式求積分 function f=Simpson2(a,b,n) %a,b 分別為積分下限和上限 %n 為分割區(qū)間個數(shù) while n10000 T=; %生成空矩陣以儲存計算結(jié)果 for n=n:n+1 h=(b-a)/(2*n); x=zeros(1,2*n+1); for i=1:2*n+1x(i)=a+(i-1)*h;endz=2*x./(x.2-3);t=0;for i=1:n

9、;t=t+h/3*(z(2*i-1)+4*z(2*i)+z(2*i+1); end T=T,t; end%將積分區(qū)間分成%計算 n+1%計算 2n+12n 份個求積節(jié)點(diǎn)值并儲存在 x 中%將分成 n 份和 n+1 份結(jié)果儲存在個求積節(jié)點(diǎn)處的函數(shù)值%用復(fù)化 simpsonT中公式求積分if abs(T(2)-T(1) Simpson2(2,3,4)%進(jìn)行誤差判斷%計算結(jié)果%實際分割區(qū)間個數(shù)%計算誤差);1.791896r =-4.9476e-09龍貝格公示求積分 function f=romberg2(a,b) t(1,1)=(b-a)/2*(2*a/(a2-3)+2*b/(b2-3);n=1;

10、while n100for k=1:naa=0;for i=0:(power(2,(k-1)-1);aa=aa+2*(a+(2*i+1)*(b-a)/2k)/(a+(2*i+1)*(b-a)/2k)2-3); endt(1,k+1)=0.5*t(1,k)+(b-a)/(2k)*aa;endfor m=1:nh=1/(4m)-1);for k=m:nt(m+1,k+1)=(4m)*t(m,k+1)-t(m,k)*h;endendif (abs(t(n,n)-t(n+1,1+n) romberg2(2,3)1.7918，的情況分別求解3. 使用帶選主元的分解法求解線性方程組,其中當(dāng) 時 對于 精確

11、解為 對得到的結(jié)果與精確解的差異進(jìn)行解釋 function f=LU(n) %生成矩陣 AA=zeros(n);for i=1:nfor j=1:nA(i,j)=i(j-1);endend%生成右邊項 bb=zeros(n,1);b(1)=n;for i=2:nb(i)=(in-1)/(i-1);endx=zeros(n,1);y=zeros(n,1);a=zeros(n,1);% 儲存矩陣 A 的換行元素c=0;% 儲存右邊項的換行元素for i=1:n-1e,m=max(abs(A(i:n,i);% 尋找最大元素及其位置a=A(i,:);A(i,:)=A(m+i-1,:);A(m+i-1,

12、:)=a;c=b(i);b(i)=b(m+i-1);b(m+i-1)=c;if A(i,i)=0disp(A 是奇異矩陣，方程組沒有唯一解 );return;endfor j=i+1:nd=A(j,i)/A(i,i);A(j,i)=d;A(j,i+1:n)=A(j,i+1:n)-d*A(i,i+1:n);endend%求解 Ly=by(1)=b(1);for k=2:ny(k)=b(k)-A(k,1:k-1)*y(1:k-1);end %求解 Ux=y x(n)=y(n)/A(n,n); for k=n-1:-1:1 x(k)=(y(k)-A(k,k+1:n)*x(k+1:n)/A(k,k);

13、 end x LU(3)111 LU(7)x =1111111 LU(11)x =1.00010.99971.00040.99981.00011.00001.00001.00001.00001.00001.00004. 用 4 階 Runge-kutta 法求解微分方程2t 1 1 2t 2tu u 2t 2u,u(0) ,u(t) e 2t te 2t10 1040步(1) 令h 0.1，使用上述程序執(zhí)行 20 步，然后令 h 0.05，使用上述程序執(zhí)行(2) 比較兩個近似解與精確解(3) 當(dāng) h減半時， (1)中的最終全局誤差是否和預(yù)期相符？(4) 在同一坐標(biāo)系上畫出兩個近似解與精確解 (

14、提示輸出矩陣 R包含近似解的 x和 y 坐 標(biāo)，用命令 plot(R(:,1),R(:,2)畫出相應(yīng)圖形 )function f=Rungek(x,t)uo=x;tn=t;%h=0.1,計算 20 步，并將計算結(jié)果儲存在 U1 中U1=zeros(20,1);h1=0.1;i=1;while i=20k1=exp(-2*tn)-2*uo;k2=exp(-2*(tn+0.5*h1)-2*(uo+0.5*h1*k1);k3=exp(-2*(tn+0.5*h1)-2*(uo+0.5*h1*k2); k4=exp(-2*(tn+h1)-2*(uo+h1*k3);un=uo+h1/6*(k1+2*k2+

15、2*k3+k4);U1(i)=un;i=i+1;tn=(i-1)*h1;uo=un;end%h=0.05,計算 40 步，并將計算結(jié)果儲存在 U2 中 U3=zeros(40,1);h2=0.05;i=1;uo=x;tn=t;while i Rungek(0.1,0)U1 =0.16370.20110.21950.22470.22070.21080.19730.18170.16530.14890.13300.11790.10400.09120.07970.06930.06010.05190.04470.0385U2 =0.16370.20110.21950.22470.22070.21080.

16、19730.18170.16530.14890.13300.11790.10400.09120.07970.06930.06010.05190.04470.03850.16370.20110.21950.22470.22070.21080.19730.18170.16530.14890.13300.11790.10400.09120.07970.06930.06010.05190.04470.0385R1 =1.0e-005 *0.21140.32500.37320.37910.35900.32420.28250.23890.19660.15750.12260.09240.06670.0454

17、0.02810.01420.0034-0.0048-0.0108-0.0151R2 =1.0e-006 *0.11910.18300.20980.21270.20100.18110.15740.13260.10870.08660.06700.04990.03560.02360.01400.00630.0003-0.0042-0.0075-0.0097R3 =1.0e-005 *-0.1995-0.3067 -0.3522 -0.3578 -0.3389 -0.3061 -0.2667 -0.2256 -0.1857 -0.1488 -0.1159 -0.0874 -0.0632 -0.0430

18、 -0.0267 -0.0136-0.00340.00440.0101 0.01415. 設(shè)為 階的三對角方陣， 是一個 階的對稱正定矩陣其中 為 階單位矩陣。 設(shè) 為線性方程組的真解， 右邊的向量 由這個真解給出。(1) 用 Cholesky 分解法求解該方程 .(2) 用 Jacobi 迭代法和 Gauss-Seidel迭代法求解該方程組，誤差設(shè)為.其中 取值為 4，5， 6.用 Cholesky 分解法求解方程： function f=Cholesky(n)%生成 Tn 矩陣 T=zeros(n);for i=1:nT(i,i)=2;endfor i=1:n-1T(i,i+1)=-1;e

19、ndfor i=2:nT(i,i-1)=-1;end%生成系數(shù)矩陣 AI=eye(n);B=T+2*I;C=-I;A=zeros(n2);for i=1:nA(n*(i-1)+1:i*n,n*(i-1)+1:i*n)=B;endfor i=1:n-1A(i*n+1:(i+1)*n,(i-1)*n+1:i*n)=C;A(i-1)*n+1:i*n,i*n+1:(i+1)*n)=C;end%生成右邊項 b y=ones(n2,1); b=A*y;分解R=chol(A);%對系數(shù)矩陣 A 進(jìn)行 Choleskyx=R(Rb)%求解方程組 Cholesky(4)x =1.00001.00001.0000

20、1.00001.00001.00001.00001.00001.00001.00001.00001.00001.00001.00001.00001.0000 Cholesky(5)x =1.00001.00001.00001.00001.00001.00001.00001.00001.00001.00001.00001.00001.00001.00001.00001.00001.00001.00001.00001.00001.00001.00001.00001.00001.0000 Cholesky(6)x =1.00001.00001.00001.00001.00001.00001.0000

21、1.00001.00001.00001.00001.00001.00001.00001.00001.00001.00001.00001.00001.00001.00001.00001.00001.00001.00001.00001.00001.00001.00001.00001.00001.00001.00001.00001.00001.0000Jacobi 迭代法求解方程： function f=Jacobi5(x0,n) %生成 Tn 矩陣 T=zeros(n);for i=1:nT(i,i)=2;endfor i=1:n-1T(i,i+1)=-1;endfor i=2:nT(i,i-1)

22、=-1;end%生成系數(shù)矩陣 AI=eye(n);B=T+2*I;C=-I;A=zeros(n2);for i=1:nA(n*(i-1)+1:i*n,n*(i-1)+1:i*n)=B;endfor i=1:n-1A(i*n+1:(i+1)*n,(i-1)*n+1:i*n)=C;A(i-1)*n+1:i*n,i*n+1:(i+1)*n)=C;end%生成右邊項 b y=ones(n2,1); b=A*y;r,s=size(A);xold=x0;C=-A;for i=1:rC(i,i)=0;endfor i=1:rC(i,:)=C(i,:)/A(i,i);end d=zeros(n2,1); fo

23、r i=1:rd(i)=b(i)/A(i,i);endi=1;while i=500 xnew=C*xold+d;if norm(xnew-xold) a=zeros(16,1); Jacobi5(a,4)Jacobi 迭代法收斂1.00001.00001.00001.00001.00001.00001.00001.00001.00001.00001.00001.00001.00001.00001.00001.000087 a=zeros(25,1); Jacobi5(a,5) Jacobi 迭代法收斂1.00001.00001.00001.00001.00001.00001.00001.00

24、001.00001.00001.00001.00001.00001.00001.00001.00001.00001.00001.00001.00001.00001.00001.00001.00001.0000126 a=zeros(36,1); Jacobi5(a,6)Jacobi 迭代法收斂x =1.00001.00001.00001.00001.00001.00001.00001.00001.00001.00001.00001.00001.00001.00001.00001.00001.00001.00001.00001.00001.00001.00001.00001.00001.0000

25、1.00001.00001.00001.00001.00001.00001.00001.00001.00001.00001.0000172 用 Gauss-Seidel 迭代法求解 function 發(fā) =Gauss(xo,n) %生成 Tn 矩陣 T=zeros(n); for i=1:nT(i,i)=2; end for i=1:n-1T(i,i+1)=-1; end for i=2:nT(i,i-1)=-1; end %生成系數(shù)矩陣 A I=eye(n); B=T+2*I; C=-I; A=zeros(n2); for i=1:nA(n*(i-1)+1:i*n,n*(i-1)+1:i*n

26、)=B; end for i=1:n-1A(i*n+1:(i+1)*n,(i-1)*n+1:i*n)=C; A(i-1)*n+1:i*n,i*n+1:(i+1)*n)=C; end %生成右邊項 by=ones(n2,1); b=A*y;%用 Gauss-Seidel 法求解 C=-A;L=zeros(n2); U=zeros(n2); D=zeros(n2);for i=1:n2D(i,i)=A(i,i);end for i=2:n2L(i,1:i-1)=C(i,1:i-1); end for i=1:n2-1U(i,i+1:n2)=C(i,i+1:n2); endB=inv(D-L)*U;

27、 f=inv(D-L)*b; i=1;while i500xn=B*xo+f;if norm(xn-xo) a=zeros(16,1); Gauss(a,4) 迭代法收斂x =1.00001.00001.00001.00001.00001.00001.00001.00001.00001.00001.00001.00001.00001.00001.00001.000046 a=zeros(25,1); Gauss(a,5) 迭代法收斂1.00001.00001.00001.00001.00001.00001.00001.00001.00001.00001.00001.00001.00001.00001.00001.00001.