矩陣與數值分析實驗報告任課教師：張宏偉教學班號：02院系：電信（計算機應用技術）學生姓名：學生學號：方程在x=3.0附近有根，試寫出其三種不同的等價形式以構成兩種不同的迭代格式，再用這兩種迭代求根，并繪制誤差下降曲線，觀察這兩種迭代是否收斂及收斂的快慢。解答：代碼如下：clear;symsxty;x=3;%初始迭代點t=0;%中間變量y=0;%繪制下降曲線變化量k=0;%迭代計數變量epx=1;%變量計算差E=1e-20;%精度f1=(2*x^3-5*x^2+42)/19;%迭代1f2=(2*x^3-19*x+42)^(1/2);%迭代2f3=(5*x^2+19*x-42)^(1/3);%迭代3while(epx>E&k<1000)%循環迭代k=k+1;y(k)=x;t=f1;%標記1epx=abs(t-x);x=t;f1=(2*x^3-5*x^2+42)/19;%標記2end;plot(y);title('迭代1誤差下降曲線');迭代公式1收斂結果:x=2迭代公式1誤差變化曲線迭代公式2誤差變化曲線迭代公式3收斂結果：x=6.8750迭代公式3誤差變化曲線結果分析：迭代公式1：f1=(2*x^3-5*x^2+42)/19和迭代公式3：f3=(5*x^2+19*x-42)^(1/3)計算結果收斂，其中迭代公式1收斂速度快于迭代公式3。迭代公式2：f2=(2*x^3-19*x+42)^(1/2)計算結果不收斂。用復化梯形公式、復化辛普森公式、龍貝格公式求下列定積分，要求絕對誤差為，并將計算結果與精確解進行比較：（1）（2）解答：（1）復化梯形公式代碼：clear;symsxsum1sum2a1b1;fun=(2/3)*(x^3)*exp(x^2);sum1=0;%積分結果變量1sum2=0;%積分結果變量2n=1;%迭代計數變量epx=1;%階段誤差E=0.5e-8;%精度a=1;%積分下限b=2;%積分上限while(epx>E)h=(b-a)/n;fori=0:(n-1)%for循環求解一次N點積分結果a1=subs(fun,x,(a+i*h));b1=subs(fun,x,(a+(i+1)*h));t=((a1+b1)*h)/2;sum2=sum2+t;end;epx=abs(sum1-sum2);%計算階段誤差sum1=sum2;%使用sum1暫存上次的計算結果sum2=0;n=n*2;end;disp('積分結果為：'),vpa(sum1)積分結果為：ans=分析：由結果可知計算結果含有8位有效數字，已滿足精度要求。使用數據點數為N=8192。復化辛普森公式代碼：clear;symsxsum1sum2a1b1;fun=(2/3)*(x^3)*exp(x^2);sum1=0;%積分結果變量1sum2=0;%積分結果變量2n=1;%迭代計數變量epx=1;%變量計算差E=0.5e-8;%精度a=1;%積分下限b=2;%積分上限while(epx>E)h=(b-a)/n;fori=0:(n-1)a1=subs(fun,x,(a+i*h));b1=subs(fun,x,(a+(i+1)*h));ab=subs(fun,x,((a+i*h)+(a+(i+1)*h))/2);t=((a1+4*ab+b1)*h)/6;sum2=sum2+t;end;epx=abs(sum1-sum2);sum1=sum2;sum2=0;n=n*2;end;disp('積分結果為：'),vpa(sum1)積分結果為：ans=分析：由結果可知計算結果含有11位有效數字，已滿足精度要求。使用數據點數為N=1024,對比復化梯形公式可以復化辛普森公式計算精度更高。龍貝格公式代碼：clear;symsxa1b1y;fun=(2/3)*(x^3)*exp(x^2);sum1=0;%積分結果變量1sum2=0;%積分結果變量2epx=1;%變量計算差E=0.5e-8;%精度a=1;%積分下限b=2;%積分上限k=0;%迭代計數變量1m=0;%迭代計數變量2while(epx>E)k=k+1;m=m+1;h=(b-a)/(2^k);fori=0:(2^k-1)a1=subs(fun,x,(a+i*h));b1=subs(fun,x,(a+(i+1)*h));t=((a1+b1)*h)/2;sum2=sum2+t;end;sum1=sum2;y(m)=sum2;sum2=0;fori=0:k-2;m=m+1;y(m)=((4^(k-1))*y(m-1)-y(m-k))/(4^(k-1)-1);end;if(m>1)epx=abs(sym(y(m)-y(m-1)));end;end;disp('積分結果為：'),vpa(y(m))積分結果為：ans=分析：由結果可知計算結果含有9位有效數字，已滿足精度要求。（2）積分結果為：ans=真值約為：1.791759469228分析：由結果可知計算結果含有8位有效數字，已滿足精度要求。使用數據點數為N=65536。積分結果為：ans=真值約為：1.791759469228分析：由結果可知計算結果含有10位有效數字，已滿足精度要求。使用數據點數為N=512，可知復化辛普森公式精度明顯高于復化梯形公式。積分結果為：ans=真值約為：1.791759469228分析：由結果可知計算結果含有8位有效數字，已滿足精度要求。3.使用帶選主元的分解法求解線性方程組,其中，，當時,．對于的情況分別求解。精確解為．對得到的結果與精確解的差異進行解釋。解答：程序代碼：functionx=lusolve(A,b)[n,n]=size(A);x=zeros(n,1);y=zeros(n,1);temprow=zeros(n,1);tempconstant=0;Pvector=zeros(n,1);forcol=1:n-1[max_element,index]=max(abs(A(col:n,col)));temprow=A(col,:);A(col,:)=A(index+col-1,:);A(index+col-1,:)=temprow;tempconstant=b(col);b(col)=b(index+col-1);b(index+col-1)=tempconstant;ifA(col,col)==0disp('Aissingular.nouniquesolution');returnendforrow=col+1:nmult=A(row,col)/A(col,col);A(row,col)=mult;A(row,col+1:n)=A(row,col+1:n)-mult*A(col,col+1:n);endendy(1)=b(1);fork=2:ny(k)=b(k)-A(k,1:k-1)*y(1:k-1);endx(n)=y(n)/A(n,n);fork=n-1:-1:1x(k)=(y(k)-A(k,k+1:n)*x(k+1:n))/A(k,k);end實驗結果：n=3時，解得x=[1,1,1]T;n=7時，解得x=[1,1,1,1,1,1,1]T;n=11時，解得x=[1.0001,0.9998,1.0002,0.9999,1,1,1,1,1,1,1]T;實驗結果分析：高斯列主元消去法，避免了小主元做除數，在Gauss消去法中增加選主元的過程，在第K步消元時，首先在第K列主對角元以下元素中挑選絕對值最大的數，并通過初等行變換，使得該數位于主對角上，然后繼續消元。當矩陣維數較小時，精度較高，但隨著矩陣維數增大，精度下降。4.用4階Runge-kutta法求解微分方程(1)令，使用上述程序執行20步，然后令，使用上述程序執行40步(2)比較兩個近似解與精確解(3)當減半時，(1)中的最終全局誤差是否和預期相符？(4)在同一坐標系上畫出兩個近似解與精確解．（提示輸出矩陣包含近似解的和坐標，用命令plot(R(:,1),R(:,2))畫出相應圖形．）解答：程序代碼：functionR=rk4(f,a,b,ya,n)symstuh=(b-a)/n;T=zeros(1,n+1);Y=zeros(1,n+1);T=a:h:b;Y(1)=ya;fork=1:nK1=h*subs(f,{t,u},{T(k),Y(k)});K2=h*subs(f,{t,u},{T(k)+h/2,Y(k)+K1/2});K3=h*subs(f,{t,u},{T(k)+h/2,Y(k)+K2/2});K4=h*subs(f,{t,u},{T(k)+h,Y(k)+K3});Y(k+1)=Y(k)+(K1+2*K2+2*K3+K4)/6;endR=[T'Y'];h=0.1時，得到的解：t近似解精確解00.200000000000000.300000000000000.400000000000000.224660691276870.224664482058610.500000000000000.220724075006090.220727664702870.600000000000000.700000000000000.800000000000000.900000000000001.000000000000001.200000000000001.300000000000001.400000000000001.500000000000000.079659028732490.079659309388581.600000000000000.069295604433720.069295746763221.700000000000000.060071851516820.060071885928591.800000000000001.900000000000000.044741651887010.044741543712332.000000000000000.038462992616640.03846284166634h=0.05時得到的解：t近似解精確解00.050000000000000.200000000000000.250000000000000.300000000000000.350000000000000.223463173139530.223463386706130.400000000000000.224664269349820.224664482058610.450000000000000.223613104618570.223613312857330.500000000000000.220727463675860.220727664702870.550000000000000.600000000000000.650000000000000.204398675243600.204398844775510.700000000000000.750000000000000.800000000000000.850000000000000.900000000000000.950000000000001.000000000000001.050000000000001.200000000000001.250000000000001.300000000000001.350000000000000.097447964174010.097447993472641.400000000000001.450000000000000.085285972550030.085285991087431.500000000000000.079659295419300.079659309388581.550000000000000.074331174050690.074331183949371.600000000000000.069295740471530.069295746763221.650000000000000.064545539837660.064545542952171.700000000000000.060071885594700.060071885928591.750000000000000.055865161413770.055865159331291.800000000000001.850000000000000.048210882562780.048210876617161.900000000000000.044741551162660.044741543712331.950000000000002.000000000000000.038462851405420.03846284166634當h減半時，（1）中的全局誤差縮小，和最終的預期相符。近似解與精確解圖形：設為階的三對角方陣，是一個階的對稱正定矩陣其中為階單位矩陣。設為線性方程組的真解，右邊的向量由這個真解給出。用Cholesky分解法求解該方程。用Jacobi迭代法和Gauss-Seidel迭代法求解該方程組，誤差設為.其中取值為4，5，6.解答：（1）Cholesky分解程序代碼（Matlab）只列出n等于4時的代碼，n等于5和6時代碼類似。clearclc%定義T(n)T4=[2-100;-12-10;0-12-1;00-12];%定義單位矩陣I4=eye(4);%定義0矩陣Z4=zeros(4,4);%定義A(n)%定義A(n)中的對角元素d4=T4+2*I4;A4=[d4-I4Z4Z4;-I4d4-I4Z4;Z4-I4d4-I4;Z4Z4-I4d4];%定義xx4=ones(16,1);%求右邊向量bb4=A4*x4;%(1)用Cholesky分解法求解方程組l4=zeros(16,16);%n為4時求解Ln=16;forj=1:nsum=0;fork=1:j-1sum=sum+l4(j,k)*l4(j,k);endl4(j,j)=sqrt(A4(j,j)-sum);fori=j+1:nsum=0;fork=1:j-1sum=sum+l4(i,k)*l4(j,k);endl4(i,j)=(A4(i,j)-sum)/l4(j,j);endend%n等于4時求解方程的解LL'x=b%Ly=bL'x=y%計算yy4=zeros(16,1);y4(1)=b4(1)/l4(1,1);fori=1:16sum=0;fork=1:i-1sum=sum+l4(i,k)*y4(k);endy4(i)=(b4(i)-sum)/l4(i,i);end%計算x_4x_4=zeros(16,1);x_4(16)=y4(16)/l4(16,16);fori=15:-1:1sum=0;fork=i+1:16sum=sum+l4(k,i)*x_4(k);endx_4(i)=(y4(i)-sum)/l4(i,i);end運行結果：（n等于4,5,6時）x的轉置n=4時：[1111111111111111]n=5時：[1111111111111111111111111]n=6時：[111111111111111111111111111111111111]（2）Jacobi迭代法代碼（只列出n等于4時的代碼，n等于5和6時類似）%取初始向量為0jx40=zeros(16,1);jx41=zeros(16,1);dx4=zeros(16,1);while1fori=1:16sum=0;forj=1:16sum=sum+A4(i,j)*jx40(j);enddx4(i)=(b4(i)-sum)/A4(i,i);jx41(i)=jx40(i)+dx4(i);endifnorm(jx41-jx40)>1e-8%迭代jx40=jx41;elsebreak;endenddisp(jx41');運行結果：（n等于4,5,6時）x的轉置當n=4時：[0.99999999485577 0.999999991676460.999999991676460.999999994855770.999999991676460.999999986532230.999999986532230.999999991676460.999999991676460.999999986532230.999999986532230.999999991676460.99999999485577 0.999999991676460.999999991676460.99999999485577]當n=5時：[0.999999994767790.999999991030500.999999989535580.999999991030500.999999994767790.999999991030500.999999984303370.999999982061000.999999984303370.999999991030500.999999989535580.999999982061000.999999979071160.999999982061000.999999989535580.999999991030500.999999984303370.999999982061000.999999984303370.999999991030500.999999994767790.999999991030500.999999989535580.999999991030500.99999999476779]當n=6時的結果：[0.999999995214840.999999991377430.999999989247830.999999989247830.999999991377430.999999995214840.999999991377430.999999984462670.999999980625260.999999980625260.999999984462670.999999991377430.999999989247830.999999980625260.999999975840100.999999975840100.999999980625260.999999989247830.999999989247830.999999980625260.999999975840100.999999975840100.999999980625260.999999989247830.999999991377430.999999984462670.999999980625260.999999980625260.999999984462670.999999991377430.999999995214840.999999991377430.999999989247830.999999989247830.999999991377430.99999999521484]（3）Gauss-Seidel迭代法代碼（只列出n等于4時的代碼，n等于5和6時類似）%取初始向量為0gsx40=zeros(16,1);gsx41=zeros(16,1);gs_dx4=zeros(16,1);while1fori=1:16sum=0;forj=1:i-1sum=sum+A4(i,j)*gsx40(j);endforj=i:16sum=sum+A4(i,j)*gsx40(j);endgs_dx4(i)=(b4(i)-sum)/A4(i,i);gsx41(i)=gsx40(i)+gs_dx4(i);endifnorm(gsx41-gsx40)>1e-8%迭代gsx40=gsx41;elsebreak;endenddisp(gsx41');運行結果：（n等于4,5,6時）x的轉置當n=4時的結果：[0.999999994855770.999999991676460.999999991676460.999999994855770.999999991676460.999999986532230.999999986532230.999999991676460.999999991676460.999999986532230.999999986532230.999999991676460.999999994855770.999999991676460.999999991676460.99999999485577]當n=5時的結果：[0.999999994767790.999999991030500.999999989535580.999999991030500.999999994767790.999999991030500.999999984303370.999999982061000.999999984303370.999999991030500.999999989535580.999999982061000.999999979071160.999999982061000.999999989535580.999999991030500.999999984303370.999999982061000.999999984303370.999999991030500.999999994767790.999999991030500.999999989535580.999999991030500.99999999476779]當n=6時的結果：[0.999999995214840.999999991377430.999999989247830.999999989247830.999999991377430.9999999952