數(shù)值線性代數(shù)第二版徐樹方高立張平文上機習(xí)題試驗報告供參考

1、百度文庫-讓每個人平等地提升自我上機習(xí)題1.先用你所熟悉的的計算機語言將不選主元和列主元Gauss消去法編寫成通用的子程序；然后用你編寫的程序求解84階方程組；最后將你的計算結(jié)果與方程的精確解進行比較，并就此談?wù)勀銓auss消去法的看法。Sol/(1)先用matlab將不選主元和列主元Gauss消去法編寫成通用的子程序，得到L,U,P：不選主元Gauss消去法：L,UGaussLA(A)得至ijL,U滿足ALU列主元Gauss消去法：L,U,PGaussCol(A)得至ijL,U,P滿足PALU(2)用前代法解LyborPb,得y用回代法解Uxy,得x求解程序為xGaussA,b,L,U,P

2、(P可缺省，缺省時默認為單位矩陣)(3)計算腳本為ex1_1代碼娟法(計算三角分解：Gauss消去法)/functionL,U=GaussLA(A)n=length(A);fork=1:n-1A(k+1:n,k)=A(k+1:n,k)/A(k,k);/A(k+1:n,k+1:n)=A(k+1:n,k+1:n)-A(k+1:n,k)*A(k,k+1:n);end11百度文庫-讓每個人平等地提升自我U=triu(A);L=tril(A);/L=L-diag(diag(L)+diag(ones(1,n);end%T法計算列主元三角分解：列主元Gauss消去法)functionL,U,P=GaussC

3、ol(A)n=length(A);fork=1:n-1s,t=max(abs(A(k:n,k);p=t+k-1;temp=A(k,1:n);A(k,1:n)=A(p,1:n);A(p,1:n)=temp;u(k)=p;ifA(k,k)=0A(k+1:n,k)=A(k+1:n,k)/A(k,k);A(k+1:n,k+1:n)=A(k+1:n,k+1:n)-A(k+1:n,k)*A(k,k+1:n);elsebreak;endendL=tril(A);U=triu(A);L=L-diag(diag(L)+diag(ones(1,n);22百度文庫-讓每個人平等地提升自我P=eye(n);fori=

4、1:n-1temp=P(i,:);P(i,:尸P(u(i),:);P(u(i),:)=temp;endend喻斯消去法解線性方程組functionx=Gauss(A,b,L,U,P)ifnargin<5P=eye(length(A);endn=length(A);b=P*b;forj=1:n-1b(j)=b(j)/L(j,j);b(j+1:n)=b(j+1:n)-b(j)*L(j+1:n,j);endb(n)=b(n)/L(n,n);y=b;forj=n:-1:2y(j)=y(j)/U(j,j);33百度文庫-讓每個人平等地提升自我y(1:j-1)=y(1:j-1)-y(j)*U(1:j

5、-1,j);end/y(1)=y(1)/U(1,1);/x=y;/、end/ex1_1clc;clear;淋一題A=6*eye(84)+diag(8*ones(1,83),-1)+diag(ones(1,83),1);b=7;15*ones(82,1);14;環(huán)選主元Gauss消去法L,U=GaussLA(A);x1_1=Gauss(A,b,L,U);%主元Gauss消去法L,U,P=GaussCol(A);x1_2=Gauss(A,b,L,U,P);%軍的比較subplot(1,3,1);plot(1:84,x1_1,'o-');title('Gauss');

6、/subplot(1,3,2);plot(1:84,x1_2,'.-');title('PGauss');subplot(1,3,3);plot(1:84,ones(1,84),'*-');title('精確解');結(jié)果為(其中Gauss表示不選主元的Gauss消去法，/PGaus續(xù)示列主元Gauss消去法，精確解為1,1184)：/44百度文庫-讓每個人平等地提升自我由圖，顯然列主元消去法與精確解更為接近，不選主元的Gauss消去法誤差比列主元消去法大，且不如列主元消去法穩(wěn)定。Gauss消去法重點在于A的分解過程，無論A如何分解

7、，后面兩步的運算過程不變。2.先用你所熟悉的的計算機語言將平方根法和改進的平方根法編寫成通用的子程序；然后用你編寫的程序求解對稱正定方程組Ax=boSol:(1)先用matlab將平方根法和改進的平方根法編寫成通用的子程序，得到L,(D)：平方根法：L=Cholesky(A)改進的平方根法：L,D=LDLt(A)(2)求解得Lyb/求解得LTxyorDLTxy求解程序為x=Gauss(A,b,L,U,P)(ULTorUDLT,P此時缺省，缺省時默認為單位矩陣)55百度文庫-讓每個人平等地提升自我(3)計算腳本為ex1_2代碼嫡法(計算Cholesky分解：平方根法)functionL=Chol

8、esky(A)n=length(A);fork=1:nA(k,k)=sqrt(A(k,k);A(k+1:n,k)=A(k+1:n,k)/A(k,k);forj=k+1:nA(j:n,j)=A(j:n,j)-A(j:n,k)*A(j,k);endendL=tril(A);end%十算LDL'分解：改進的平方根法functionL,D=LDLt(A)n=length(A);forj=1:nfori=1:nv(i,1)=A(j,i)*A(i,i);endA(j,j)=A(j,j)-A(j,1:j-1)*v(1:j-1,1);A(j+1:n,j)=(A(j+1:n,j)-A(j+1:n,1:j

9、-1)*v(1:j-1,1)/A(j,j);66百度文庫-讓每個人平等地提升自我endL=tril(A);D=diag(diag(A);L=L-diag(diag(L)+diag(ones(1,n);end喻斯消去法解線性方程組functionx=Gauss(A,b,L,U,P)ifnargin<5P=eye(length(A);endn=length(A);b=P*b;forj=1:n-1b(j)=b(j)/L(j,j);b(j+1:n)=b(j+1:n)-b(j)*L(j+1:n,j);endb(n)=b(n)/L(n,n);y=b;forj=n:-1:2y(j)=y(j)/U(j,

10、j);y(1:j-1)=y(1:j-1)-y(j)*U(1:j-1,j);end77百度文庫-讓每個人平等地提升自我y(1)=y(1)/U(1,1);x=y;/end/ex1_2二題澄一問A=10*eye(100)+diag(ones(1,99),-1)+diag(ones(1,99),1);b=round(100*rand(100,1);對方根法L=Cholesky(A);x1_2_1_1=Gauss(A,b,L,L');峨進的平方根法L,D=LDLt(A);x1_2_1_2=Gauss(A,b,L,D*L');蹴二問A=hilb(40);b=sum(A);b=b'對

11、方根法L=Cholesky(A);x1_2_2_1=Gauss(A,b,L,L');峨進的平方根法L,D=LDLt(A);x1_2_2_2=Gauss(A,b,L,D*L');88百度文庫-讓每個人平等地提升自我結(jié)果分別為x1_2_1_1=99百度文庫-讓每個人平等地提升自我1010百度文庫-讓每個人平等地提升自我1111百度文庫-讓每個人平等地提升自我1212百度文庫-讓每個人平等地提升自我x1_2_1_2=1313百度文庫-讓每個人平等地提升自我1414百度文庫-讓每個人平等地提升自我1515百度文庫-讓每個人平等地提升自我x1_2_2_1=1616百度文庫-讓每個人平等地

12、提升自我+07*1717百度文庫-讓每個人平等地提升自我x1_2_2_2=1818百度文庫-讓每個人平等地提升自我1919百度文庫-讓每個人平等地提升自我3.用第1題的程序求解第2題的兩個方程組并比較所有的計算結(jié)果，然后評價各個方法的優(yōu)劣。Sol:Gauss表示不選主元的Gauss消去法，PGaus汰示列主元Gauss消去法。計算腳本為：蹴三題澄一問A=10*eye(100)+diag(ones(1,99),-1)+diag(ones(1,99),1);b=round(100*rand(100,1);環(huán)選主元Gauss消去法L,U=GaussLA(A);x1_3_1_1=Gauss(A,b,L

13、,U);/%主元Gauss消去法L,U,P=GaussCol(A);x1_3_1_2=Gauss(A,b,L,U,P);/蹴二問A=hilb(40);2020百度文庫-讓每個人平等地提升自我b=sum(A);b=b'環(huán)選主元Gauss消去法L,U=GaussLA(A);x1_3_2_1=Gauss(A,b,L,U);%主元Gauss消去法L,U,P=GaussCol(A);x1_3_2_2=Gauss(A,b,L,U,P);ex1_2;平方根法1');改進的平方根法1');Gauss1');'PGauss1');'平方根法2');

14、'改進的平方根法2');y1=1:100;y2=1:40;subplot(4,2,1);plot(y1,x1_2_1_1);title(subplot(4,2,2);plot(y1,x1_2_1_2);title(subplot(4,2,3);plot(y1,x1_3_1_1);title(subplot(4,2,4);plot(y1,x1_3_1_2);title(subplot(4,2,5);plot(y2,x1_2_2_1);title(subplot(4,2,6);plot(y2,x1_2_2_2);title(subplot(4,2,7);plot(y2,x1_3_2

15、_1);title('Gauss2');2121百度文庫-讓每個人平等地提升自我1050-501050-5050-502000-2000平方根法1100-5100Gaussl1050-5100500-50-10050005050平方根法220Gauss2改進的平萬根法1102030405060708090PGauss1100102030405060708090改進的平方根法21005101520PGauss225303540510152025303540subplot(4,2,8);plot(y2,x1_3_2_2);title('PGauss2');平方根法和

16、改進的平法根法計算量更小，計算過程穩(wěn)定，但使用范圍窄；不選主元和列主元的Gauss消去法計算量較大，但適用范圍廣。例題考慮對稱正定線性方程組Ax=b,其中向量b是隨機生成的，其元素是服從區(qū)間0,1上均勻分布的隨機數(shù)，矩陣ALLT,這里L是隨機生成的一個下三角矩陣，其元素是服從區(qū)間1,2上均勻分布的隨機數(shù)。/對n=10,20,500分別應(yīng)用Gauss消去法、列主元Gauss消去法和Cholesky分解法求解該方程組，畫出它們所用的CPU時間，其中“Gauss”表示Gauss消去法、"PGausS表示列主元Gauss消去法，"Cholesky”表示Cholesky分解法。/So

17、l:經(jīng)試驗知，對應(yīng)課本上圖所示的Cholesky分解法應(yīng)為改進后的Cholesky分解法即LDLT分解。2222百度文庫-讓每個人平等地提升自我此處所用的CPU時間利用cputime測量。計算腳本為eg1_3_2clc;clear;/fori=1:50;/n=i*10;b=rand(n,1);L=tril(unifrnd(1,2,n,n);A=L*L't1(i)=cputime;L1,U1=GaussLA(A);x1=Gauss(A,b,L1,U1);t1(i)=cputime-t1(i);t2(i)=cputime;L2,U2,P2=GaussCol(A);x2=Gauss(A,b,L2,U2,P2);t2(i)=cputime-t2(i);t3(i)=cputime;L3=LDLt(A);x3=Gau