1、上機(jī)習(xí)題1.先用你所熟悉的的電腦語言將不選主元和列主元Gauss 消去法編寫成通用的子程序； 然后用你編寫的程序求解 84 階方程組；最后將你的計(jì)算結(jié)果與方程的精確解進(jìn)行比擬，并就此 談?wù)勀銓?duì) Gauss 消去法的看法。Sol：1先用matlab將不選主元和列主元Gauss消去法編寫成通用的子程序，得到L,U,P：不選主元Gauss消去法：L,U GaussLA(A)得到L,U滿足A LU列主元 Gauss消去法：L,U , P GaussCol (A)得到 L,U , P 滿足 PA LU(2)用前代法解 Ly b or Pb ，得 y用回代法解 Ux y ，得 x求解程序?yàn)閤 Gauss代

2、b,L,U,P P可缺省，缺省時(shí)默認(rèn)為單位矩陣 3 計(jì)算腳本為 ex1_1代碼%|法1.1.3 計(jì)算三角分解：Gauss消去法function L,U=GaussLA(A)n=length(A);for k=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);endU=triu(A);L=tril(A);L=L-diag(diag(L)+diag(ones(1,n);end%算法 1.2.2( 計(jì)算列主元三角分解：列主元 Gauss 消去法function L,U,P=Gau

3、ssCol(A)n=length(A);for k=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;if A(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);P=eye(n);for i=1:n-1temp=P

4、(i,:);P(i,:)=P(u(i),:);P(u(i),:)=temp;endend%高斯消去法解線性方程組function x=Gauss(A,b,L,U,P)if nargin<5P=eye(length(A);endn=length(A);b=P*b;for j=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;for j=n:-1:2y(j)=y(j)/U(j,j);y(1:j-1)=y(1:j-1)-y(j)*U(1:j-1,j);endy(1)=y(1)/U(1,

5、1);x=y;end ex1_1 clc;clear;%第一題A=6*eye(84)+diag(8*ones(1,83),-1)+diag(ones(1,83),1);b=7;15*ones(82,1);14;%不選主元 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);%解的比擬'o-' );title('Gauss' );'.-' );title('PGauss' );'*-

6、' );title( ' 精確解 ' );subplot(1,3,1);plot(1:84,x1_1,subplot(1,3,2);plot(1:84,x1_2,subplot(1,3,3);plot(1:84,ones(1,84),結(jié)果為其中Gauss表示不選主元的Gauss消去法，PGauss表示列主元Gauss消去法，精確解為 1, ,1184 ：由圖，顯然列主元消去法與精確解更為接近,不選主元的Gauss消去法誤差比列主元消去法大，且不如列主元消去法穩(wěn)定。Gauss消去法重點(diǎn)在于 A的分解過程，無論A如何分解，后面兩步的運(yùn)算過程不變。2先用你所熟悉的的電腦語言將

7、平方根法和改進(jìn)的平方根法編寫成通用的子程序；然后用你編寫的程序求解對(duì)稱正定方程組Ax=b。Sol:1先用matlab將平方根法和改進(jìn)的平方根法編寫成通用的子程序，得到L, (D)：平方根法：L=Cholesky(A)改進(jìn)的平方根法：L,D=LDLt(A)(2)求解得Ly b求解得 LTx y or DLTxy求解程序?yàn)閤=Gauss(A,b,L,U,P)ULT or U DLT，P此時(shí)缺省，缺省時(shí)默認(rèn)為單位矩陣(3)計(jì)算腳本為ex1_2 代碼%算法 1.3.1 計(jì)算 Cholesky 分解：平方根法 function L=Cholesky(A) n=length(A);for k=1:nA(k

8、,k)=sqrt(A(k,k);A(k+1:n,k)=A(k+1:n,k)/A(k,k);for j=k+1:nA(j:n,j)=A(j:n,j)-A(j:n,k)*A(j,k);endendL=tril(A);end%十算LDL '分解：改進(jìn)的平方根法function L,D=LDLt(A) n=length(A);for j=1:nfor i=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-1)*v(1:j-1,1)/A(j,j);en

9、dL=tril(A); D=diag(diag(A);L=L-diag(diag(L)+diag(ones(1,n);end%高斯消去法解線性方程組function x=Gauss(A,b,L,U,P)if nargin<5P=eye(length(A);endn=length(A);b=P*b;for j=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;for j=n:-1:2y(j)=y(j)/U(j,j);y(1:j-1)=y(1:j-1)-y(j)*U(1:j-1,j)

10、;endy(1)=y(1)/U(1,1);x=y;endex1_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');%改進(jìn)的平方根法L,D=LDLt(A);x1_2_1_2=Gauss(A,b,L,D*L');%第二問A=hilb(40);b=sum(A);b=b'%平方根法L=Cholesky(A);x1_2_2_1=Gauss(A,b,L,L');%改進(jìn)

11、的平方根法L,D=LDLt(A);x1_2_2_2=Gauss(A,b,L,D*L');結(jié)果分別為x1_2_1_1 =7.25868.4143-0.40138.59845.41770.22492.33364.43898.27728.7890-0.16678.87848.38243.29787.64010.30143.34578.24186.23688.39065.8575-0.96567.79817.98425.36016.41436.49666.30240.35637.1342-0.69852.85060.19270.22277.58015.97621.65839.4409-1.06

12、774.23622.70606.70377.25700.72654.47783.49595.56375.86756.76141.51806.05825.90020.93974.02949.00321.93825.61500.91207.26521.43604.37495.81467.47918.39424.57890.81691.25231.66038.14480.89157.94010.70758.98492.44371.57771.77905.63193.90182.35064.72454.16278.64831.35436.80876.55892.60275.41400.25770.00

13、904.65226.46858.6626-0.09485.28564.2385-0.67063.4671 x1_2_1_2 =7.25868.4143-0.40138.59845.41770.22494.43898.27728.7890-0.16678.87848.38243.29787.64010.30143.34578.24186.23688.39065.8575-0.96567.79817.98425.36016.41436.49662.62016.30240.35637.1342-0.69852.85060.22277.58015.97621.65839.4409-1.06774.23

14、622.70606.70377.25700.72654.47783.49595.56375.86756.76141.51806.05825.90020.93970.70314.02949.00321.93825.61500.91201.43604.37495.81467.47918.39424.57890.81691.25231.66038.14480.89157.94010.70758.98492.44371.57771.77905.63193.90182.35067.59254.72454.16278.64831.35436.80872.60275.41400.25770.00904.65

15、226.46858.6626-0.09485.28564.2385-0.67063.4671x1_2_2_1 =1.0e+07 *0.0000-0.00000.0001-0.0004-0.00140.0424-0.29801.1419-2.73354.2539-4.30182.7733-1.19890.5406-0.36880.3285-0.44380.4621-0.25130.05650.0000-0.00510.0071-0.0027-0.00310.0036-0.00190.00090.0002-0.0002-0.00060.00040.0001-0.00020.00010.0000-0

16、.0000-0.0000-0.0000x1_2_2_2 =1.00001.00000.99981.00111.00640.86811.8034-1.56935.5763-2.5315-1.769310.4883-6.28070.5882-4.715722.8299-19.91348.703210.3265-25.214010.0282-1.942514.1891-12.0552-0.5803-12.47918.56529.8724-10.550216.3871-5.813213.421611.1767-64.315446.383712.6957-21.755612.1204-1.9342然后評(píng)

17、價(jià)各個(gè)方法的優(yōu)Gauss消去法3.用第 1 題的程序求解第 2題的兩個(gè)方程組并比擬所有的計(jì)算結(jié)果， 劣。Sol：Gauss表示不選主元的Gauss消去法，PGauss表示列主元計(jì)算腳本為：%第三題%第一問A=10*eye(100)+diag(ones(1,99),-1)+diag(ones(1,99),1);b=round(100*rand(100,1);%不選主元 Gauss 消去法L,U=GaussLA(A); x1_3_1_1=Gauss(A,b,L,U);%列主元 Gauss 消去法 L,U,P=GaussCol(A); x1_3_1_2=Gauss(A,b,L,U,P);%第二問A=

18、hilb(40);b=sum(A);b=b'%不選主元 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;y1=1:100;y2=1:40;subplot(4,2,1);plot(y1,x1_2_1_1);title(' 平方根法 1' );subplot(4,2,2);plot(y1,x1_2_1_2);title(' 改進(jìn)的平方根法 1' );subplot(4,2,3);p

19、lot(y1,x1_3_1_1);title('Gauss1' );subplot(4,2,4);plot(y1,x1_3_1_2);title('PGauss1' );subplot(4,2,5);plot(y2,x1_2_2_1);title(' 平方根法 2' );subplot(4,2,7);plot(y2,x1_3_2_1);title('Gauss2');10平方根法150-5Gauss11050200-200subplot(4,2,8);plot(y2,x1_3_2_2);title(1M燭A7幾 i k.L1i i

20、 v、' 11 . I 1： -Ulfn rJl? ?/w./ 氷 J; J I.1 r ¥血上-J:I' f f ; /.rg-rraEr50-53 f 飛 /'亠A A A .vx/V/VVf y ii a- J v VlJ7 v y*j V1 1irr.r.r10改進(jìn)的平方根法1I.0-52050100105/ I加040PGauss1-55060100改進(jìn)的平方根法2'PGauss2');平方根法和改進(jìn)的平法根法計(jì)算量更小，計(jì)算過程穩(wěn)定，但使用范圍窄;不選主元和列主元的 Gauss消去法計(jì)算量較大，但適用范圍廣。例題考慮對(duì)稱正定線性方程組Ax=b，其中向量b是隨機(jī)生成的，其元素是服從區(qū)間0,1上均勻分布的隨機(jī)數(shù)，矩陣A LLT,這里L(fēng)是隨機(jī)生成的一個(gè)下三角矩陣，其元素是服從區(qū)間1,2上均勻分布的隨機(jī)數(shù)。對(duì)n=10, 20，.,500分別應(yīng)用 Gauss消去法、列主元Gauss消去法和Cholesky分解法求解該 方程組，畫出它們所用的 CPU時(shí)間，其中"Gauss"表示Gauss消去法、