(完整word版)數(shù)值線性代數(shù)第二版徐樹方高立張平文上機習(xí)題第一章實驗報告

(完整word版)數(shù)值線性代數(shù)第二版徐樹方高立張平文上機習(xí)題第一章實驗報告(完整word版)數(shù)值線性代數(shù)第二版徐樹方高立張平文上機習(xí)題第一章實驗報告上機習(xí)題先用你所熟悉的的計算機語言將不選主元和列主元Gauss序求解84Gauss消去法的看法。Sol：先用matlab將不選主元和列主元Gauss消去法編寫成通用的子程序，得到L,U,P：GaussGaussLA(AL,UALUGaussPGaussColL,UPPALULyy用回代法解Uxy，得x求解程序為xGaussA,b,L,U,P（P可缺省，缺省時默認(rèn)為單位矩陣）ex1_1代碼%算法1.1.3（計算三角分解：Gauss消去法）function[L,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);endU=triu(A);L=tril(A);L=L-diag(diag(L))+diag(ones(1,n));end%算法1.2.2(計算列主元三角分解：列主元Gauss消去法）function[L,U,P]=GaussCol(A)n=length(A);fork=1:n-1p=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));P=eye(n);fori=1:n-1temp=P(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);y(1:j-1)=y(1:j-1)-y(j)*U(1:j-1,j);endy(1)=y(1)/U(1,1);x=y;endex1_1A=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);%列主元Gaussx1_2=Gauss(A,b,L,U,P);%解的比較subplot(1,3,1);plot(1:84,x1_1,'o-');title('Gauss');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消去法，PGauss表示列主元Gauss消去法，精確解為184

：x108Gauss PGauss6 1

精確解1.8141.61 1.421.210 11-21-4 1

0.80.60.40.2-60 50

10 50

00 50 100由圖，顯然列主元消去法與精確解更為接近，不選主元的Gauss消去法誤差比列主元消去法大，且不如列主元消去法穩(wěn)定。Gauss消去法重點在于A的分解過程，無論A如何分解，后面兩步的運算過程不變。