數(shù)值分析線性方程組實(shí)驗(yàn)報(bào)告包含代碼

線性代數(shù)方程組的直接解法實(shí)驗(yàn)?zāi)康模壕€性方程組求解的直接法編程實(shí)現(xiàn)，實(shí)驗(yàn)內(nèi)容：線性方程組求解的高斯消去法算法實(shí)現(xiàn)線性方程組求解的主元素消去法算法實(shí)現(xiàn)線性方程組求解的LU分解得方法算法實(shí)現(xiàn)線性方程組求解追趕法算法實(shí)現(xiàn)實(shí)驗(yàn)比較：高斯消去、主元素消去、LU分解都用實(shí)例這個(gè)進(jìn)行比較。知識(shí)理論解線性方程組的方法大致分為直接法和迭代法。直接法是指假設(shè)計(jì)算過(guò)程中不產(chǎn)生舍入誤差，經(jīng)過(guò)有限次的運(yùn)算可求的方程組精確解的方法。方程（2-1）記為：AX=b；一、高斯（Gauss）消元法(1).Gauss消元法是最基本的一種方法。先逐次消去變量，將方程組化成同解的上三角方程組。 消元過(guò)程：先逐次消去變量，將方程組化成同解的上三角方程組基本思想 回代過(guò)程：按方程的相反順序求解三角形方程組，得到原方程組的解。(2)Gauss消去法的求解思路為：若先讓第一個(gè)方程組保持不變，利用它消去其余方程組中的,使之變成一個(gè)關(guān)于變?cè)膎-1階方程組。按照（1）中的思路繼續(xù)運(yùn)算得到更為低階的方程組。經(jīng)過(guò)n-1步的消元后，得到一個(gè)三角方程。利用求解公式回代得到線性方程組的解。消元過(guò)程：第一次消元，（i=1,2,3……，n）.將（2-1）中第i個(gè)方程減去第一個(gè)方程乘以（i=1,2,3……，n），完成第一次消元，（2-2）其中：；簡(jiǎn)記為：其中：按上述方法完成n-1次消元后得到。同解的三角方程組：簡(jiǎn)記為：回代過(guò)程：按逆序逐步回代得到方程的解。（3）算法：（5）Matlab程序代碼：function[RA,RB,n,X]=gaus(A,b)B=[A,b];n=length(b);RA=rank(A);RB=rank(B);zhicha=RB-RA;if(zhicha~=0)disp('因?yàn)镽A≠RB，所以此方程組無(wú)解');return;endifRA==RBifRA==ndisp('因?yàn)镽A=RB=n,所以此方程有唯一解');X=zeros(n,1);forp=1:n-1fork=p+1:nm=B(k,p)/B(p,p);B(k,p:n+1)=B(k,p:n+1)-m*B(p,p:n+1);endendb=B(1:n,n+1);A=B(1:n,1:n);X(n)=b(n)/A(n,n);forq=n-1:-1:1X(q)=(b(q)-sum(A(q,q+1:n)*X(q+1:n)))/A(q,q);endelsedisp('因?yàn)镽A=RB<n,所以此方程有無(wú)窮解');endend運(yùn)行檢測(cè)：>>A=[111;12-33;-183-1];b=[6;15;-15];[RA,RB,n,X]=gaus(A,b)因?yàn)镽A=RB=n,所以此方程有唯一解>>RA=3RB=3n=3X=1.00002.00003.0000二、主元素法1.列主元素法消元(1)基本思想：在每次消元前，在要消去未知數(shù)的系數(shù)中找到絕對(duì)值最大的系數(shù)作為主元，通過(guò)對(duì)換行將其換到對(duì)角線上，然后進(jìn)行消元.(2)消元過(guò)程：與Gauss很類似，每次對(duì)對(duì)角線換成最大的值，后面過(guò)程與Gauss基本相同。如此經(jīng)過(guò)n-1步,（2-1）的增廣矩陣[A,b]被化為上三角矩陣；回代過(guò)程：同Gauss算法一樣回代求解。(3)算法：det<-1對(duì)于k=1,2,3,…n-1;|aik,k|=max|aik|(k≤i≤n)(i)如果aik=0,則停止計(jì)算（det(A)=0）(ii)如果aik=k,則zhuan（i）換行akj<-->aik,j(j=k,k+1,…n) bk<-->bik后邊就是Gauss消元(4).Matlab程序functionX=liezhu(A,b)B=[A,b];n=length(b);RA=rank(A);RB=rank(B);ifRA==RBifRA==ndisp('因?yàn)橄禂?shù)矩陣與增廣矩陣的秩相同且為n，此方程有唯一解')X=zeros(n,1);forp=1:n-1[y,j]=max(abs(B(p:n,p)));C=B(p,:);B(p,:)=B(j+p-1,:);B(j+p-1,:)=C;endforp=1:n-1fork=p+1:nm=B(k,p)/B(p,p);B(k,p:n+1)=B(k,p:n+1)-m*B(p,p:n+1);endendb=B(1:n,n+1);A=B(1:n,1:n);X(n)=b(n)/A(n,n);forq=n-1:-1:1X(q)=(b(q)-sum(A(q,q+1:n)*X(q+1:n)))/A(q,q);endelsedisp('因?yàn)橄禂?shù)矩陣與增廣矩陣的秩相同且小于n，方程組有無(wú)窮解')endelsedisp('因?yàn)橄禂?shù)矩陣與增廣矩陣的秩不相同，所以此方程組無(wú)解')endend測(cè)試結(jié)果：>>clearallclcA=[111;12-33;-183-1];b=[615-15]';liezhu(A,b)>>因?yàn)橄禂?shù)矩陣與增廣矩陣的秩相同且為n，此方程有唯一解ans=1.00002.00003.0000例2測(cè)試結(jié)果>>clearallclcA=[0.501.13.1;2.04.50.36;5.00.966.5];b=[6.00.0200.96]';liezhu(A,b)>>ans=-2.60001.00002.0000 三、直接三角分解法1.高斯消元的矩陣表示第一次消元經(jīng)過(guò)K次消元后得到：消元過(guò)程：則（1）LU分解：因?yàn)锳的各元素已知，可求出L和U的各元素。看A的第k行，有，再看A的第k列，有。這樣就可以計(jì)算計(jì)算出U的第k行和L的第k列的元素。從k=1到k=n交替計(jì)算U和L，就能逐步計(jì)算出U和L的全部元素。分兩步求解LUx=b，先解方程組Ly=b，因?yàn)長(zhǎng)是下三角矩陣，求解只要逐步向前回代。第二步是解方程組Ux=y，因?yàn)閁是上三角矩陣，用逐次向后回代的方法求出x。3.Matlab程序 function[L,U,Y,X]=Doolittle(A,b)n=length(b);L=zeros(n,n);U=zeros(n,n);%先將A進(jìn)行LU分解fori=1:n%把第一列第一行求出forj=i:nifi==1A(i,j)=A(i,j);ifj>iA(j,i)=A(j,i)/A(i,i);endendend%下面計(jì)算U的第i行值forj=i:nifi~=1m=0;fork=1:i-1m=m+A(i,k)*A(k,j);endA(i,j)=A(i,j)-m;endend%下面計(jì)算L的第i列值forj=i:nifi~=1&&j>il=0;fork=1:i-1l=l+A(j,k)*A(k,i);endA(j,i)=(A(j,i)-l)/A(i,i);endendend%將LU的值分別加入fori=1:nforj=i:nU(i,j)=A(i,j);endforj=1:iifi==jL(i,j)=1;elseL(i,j)=A(i,j);endendend%LY=B，UX=Y，進(jìn)行求解Ab=[L,b]Y=zeros(n,1);Y(1)=Ab(1,n+1)/Ab(1,1);fori=2:nforj=1:i-1Y(i)=Y(i)+Ab(i,j)*Y(j);endY(i)=(Ab(i,n+1)-Y(i))/Ab(i,i);endAb=[U,Y];X=zeros(n,1);X(n)=Ab(n,n+1)/Ab(n,n);fori=n-1:-1:1forj=i+1:nX(i)=X(i)+Ab(i,j)*X(j);endX(i)=(Ab(i,n+1)-X(i))/Ab(i,i);endend運(yùn)行檢測(cè)：>>A=[111;12-33;-183-1];b=[6;15;-15];[L,U,Y,X]=Doolittle(A,b)運(yùn)行結(jié)果：>>Ab=1.0000006.000012.00001.0000015.0000-18.0000-1.40001.0000-15.0000L=1.00000012.00001.00000-18.0000-1.40001.0000U=1.00001.00001.00000-15.0000-9.0000004.4000Y=6.0000-57.000013.2000X=1.00002.00003.0000四、解三角方程組的追趕法1.設(shè)系數(shù)矩陣為三對(duì)角矩陣則方程組Ax=f稱為三對(duì)角方程組。2.設(shè)矩陣A非奇異，A有Crout分解A=LU，其中L為下三角矩陣，U為單位上三角矩陣，記先依次求出L，U中的元素后，令Ux=y，先求解下三角方程組Ly=f得出y，再求解上三角方程組Ux=y。求解三對(duì)角方程組的2追趕法將矩陣三角分解的計(jì)算與求解兩個(gè)三角方程組的計(jì)算放在一起，使算法吏為緊湊。其計(jì)算公式為：3.算法：1.輸入a=(a2,…，an),b=(b1,…，bn),c=(c1,…，cn-1),d=(d1,…，dn),n2.對(duì)i=2,3，…，n ai/bi-1=>ai(li) bi-ci-1ai=>bi(ui) di-aidi-1=>di(yi)3.dn/bn=>dn（xn）4.對(duì)i=n-1,n-2…，1 (di-cidi+1)/bi=>di(xi)5.輸出d=x,停機(jī)4.Matlab程序：clc;clear;%如果需要驗(yàn)證其他矩陣的解，只需更改A、B和n的值即可A=[2,1,0,0;1,3,1,0;0,1,4,-1;0,0,2,3]B=[3;5;4;5]n=4;%獲得A中不為0的數(shù)存到相應(yīng)數(shù)組中fori=1:nforj=1:nifi==jb(i)=A(i,j);endifi==j-1c(j-1)=A(i,j);endifj==i-1a(i)=A(i,j);endendend%將變換后的值存入相應(yīng)矩陣中u(1)=b(1);fori=2:n;l(i)=a(i)/u(i-1);u(i)=b(i)-l(i)*c(i-1);endfori=1:nforj=1:nifi==jL(i,j)=1;U(i,j)=u(i);endifi==j-1U(i,j)=c(i);endifj==i-1L(i,j)=l(i);endendend%LY=B，UX=Y，進(jìn)行求解disp('增廣矩陣');AB=[L,B]B=zeros(n,1);B(1)=AB(1,n+1)/AB(1,1);fori=2:nforj=1:i-1B(i)=B(i)+AB(i,j)*B(j);endB(i)=(AB(i,n+1)-B(i))/AB(i,i);enddisp('解方程LY=B得');Y=BAB=[U,B];B=zeros(n,1);B(n)=AB(n,n+1)/AB(n,n);fori=n-1:-1:1forj=i+1:nB(i)=B(i)+AB(i,j)*B(j);endB(i)=(AB(i,n+1)-B(i))/AB(i,i);enddisp('方程的解為：');X=B運(yùn)行結(jié)果：A=21001310014