数值分析线性方程组的实验报告包含代码解析

1、线性代数方程组的直接解法实验目的：线性方程组求解的直接法编程实现，实验内容：线性方程组求解的高斯消去法算法实现线性方程组求解的主元素消去法算法实现线性方程组求解的LU分解得方法算法实现线性方程组求解追赶法算法实现x1x2x3=6实验比较：高斯消去、主元素消去、LU分解都用实例12x1-3x2+3x3=15这个进行比较。18x1+3x2-.x3=15知识理论解线性方程组的方法大致分为直接法和迭代法。直接法是指假设计算过程中不产生舍入误差，经过有限次的运算可求的方程组精确解的方法。ai1X1ai2x2万程（2-1）222.-ainXn=6a2nxn=b2记为：AX=b;一、图斯（Gauss）消元法

2、（1）.Gauss消元法是最基本的一种方法。先逐次消去变量，将方程组化成同解的上三角方程组。消元过程：先逐次消去变量，将方程组化成同解的上三角方程组基本思想一l回代过程：按方程的相反顺序求解三角形方程组，得到原方程组的解。Gauss消去法的求解思路为:若先让第一个方程组保持不变，利用它消去其余方程组中的，使之变成一个关于变元的n-1阶方程组。按照（1）中的思路继续运算得到更为低阶的方程组。经过n-1步的消元后，得到一个三角方程。利用求解公式回代得到线性方程组的解。a（1）消元过程：第一次消元，设a：?#0,记l口=（、（i=1,2,3，n）.将（2-1）中第i个方程减去第a¥一个方程

3、乘以li1（i=1,2,3，n）,完成第一次消元,得等价方程组=b(1)(2)a；2X2(2-2)(2)an2X2简记为：Ax=b苴中.N2)aIa(1)八丁,aijaij_lHalj,(2)(1)(1)bi-bi-1i1b1经过k1次消元后：(1)anX1(1)a12x2(2)a22X2(1) (1)(1).(1)-akXka，k.1Xk1.anXn=b1(2) (2)(2).(2)-a2kXk-a2,k1Xk.1'.a2nXn=b2(k)akkXkak,k1xk：(1.Jk).(k)aknXn=bk(k)ak-1,nXn(k)k1(k).(k)ak：；1,kXkak1,k1Xk!；

4、1(k).(k)ankXkan,k：1Xk：；1八(k)annXn(k)n的1*A(k1).(k1)间记为：Ax=b苴中.o(k1)的的aij-aj-likHkjbi(k1)=b(k)TkKk)(i,j=k1,，n)按上述方法完成n-1次消元后得到。同解的三角方程组:A(n)=,"aa11a12为a22.aa1na(2).a2n.-(n)a；n(1)b1b22)(n)bn简记为：A(n)x=b(n)回代过程：按逆序逐步回代得到方程的解。x_b(n).a(n)xnn-annJ)xl)akk(k=n-1,n2,1)(k)xk=(bkn-、akik)i±1(3)算法:1输入A=(

5、aj),b=(4力2,，bn)T,置k=12若akk=0,转3;否则输出失败信息，停机aik3 .又ti=k+1,.n,置naik,bj-aikbknbakk对j=k,1,.n,置ajaikHkj=a04若k=n-1,转5;否则k+1=k,转25若ann=0,输出失败信息，停机；否则，置殳=6annn6 .Xtk=n1,.1，置(bk£akibj/akk=bki=k17 .输出x=b.停机(5)Matlab程序代码:functionRA,RB,n,X=gaus(A,b)B=A,b;n=length(b);RA=rank(A);RB=rank(B);zhicha=RB-RA;if(zh

6、icha=0)disp('因为RAwRB,所以此方程组无解)；return;endifRA=RBifRA=ndisp('因为RA=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);endb=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('因为RA

7、=RB<n,所以此方程有无穷解)；endend运行检测：>>A=111;12-33;-183-1;b=6;15;-15;RA,RB,n,X=gaus(A,b)因为RA=RB=n,所以此方程有唯一解>>RA=3RB=3n=3X=1.00002.00003.0000二、主元素法1.列主元素法消元(1)基本思想：在每次消元前，在要消去未知数的系数中找到绝对值最大的系数作为主元，通过对换行将其换到对角线上，然后进行消元.(2)消元过程：与GausS艮类似，每次对对角线换成最大的值，后面过程与Gaus谟本相同。如：第k次消元之前，在aikk)(i=k,k+1,.n)中找绝对

8、值最大的元,即:成小r。=k,.n)则先对换增广矩阵的行kHik,再消元。(实质：每次消元之前交换方程的顺序，使对角元素最大)如此经过n-1步，(2-1)的增广矩阵A,b被化为上三角矩阵；回代过程：同GaussB法一样回代求解。(3)算法：det<-1对于k=1,2,3,-n-1;|aik,k|=max|aik|(k<i<n)(i)如果aik=0,则停止计算(det(A)=0)(ii)如果aik=k，贝Uzhuan(i)换彳亍akj<->aik,j(j=k,k+1,n)bk<->bik后边就是Gaus©肖元.Matlab程序functionX

9、=liezhu(A,b)B=A,b;n=length(b);RA=rank(A);RB=rank(B);ifRA=RBifRA=ndisp('因为系数矩阵与增广矩阵的秩相同且为n,此方程有唯一解)X=zeros(n,1);forp=1:n-1y,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

10、)/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('因为系数矩阵与增广矩阵的秩相同且小于n,方程组有无穷解)endelsedisp(因为系数矩阵与增广矩阵的秩不相同，所以此方程组无解end测试结果：>>clearallclcA=111;12-33;-183-1;b=615-15'liezhu(A,b)>>因为系数矩阵与增广矩阵的秩相同且为n,此方程有唯一解ans=1.00002.00003.0000例2测试结果>>clearallclcA=0

11、.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.高斯消元的矩阵表示第一次消元二LiA,b.LAb11000-l21100L1=-l310100100li1;整可？i=2,3,，nln1001_j经过K次消元后得到：Lk-10010lik11k1k(k)aikakk-1nk消元过程：LnLn/-L1A,bI=A(n),b(n)lL=LflL21-LLn11=一1l21131l32Jn1ln2ln3记作A=LA(n)=LUDoolittle

12、分解,即:A=LU.则A,b1=LA(n),b(n)I-La,Lb1消元过程uA=LU(Doolittle分解)其中L为单位下三角阵，U为上三角阵.2.定理：设心n阶方阵，若颔勺顺序主子式A(i=1.n-1)均不为零，则矩阵A存在唯一的(1)LU解:a1nIa21a22a2n一1l21U11U12U22an2annJn11n2因为A的各元素已知，可求出L和U的各元素。u1nIunnk看A的第卜行，有ULjEjujkD，再看A的第k列，有1kUkk。这样就可以计算计算出U的第k行和L的第k列lik=_lirUrk)(i=k1,k2,.n)的元素。(3)从k=1到k=n交替计算U和L,就能逐步计算

13、出U和L的全部元素。(4)分两步求解LUx=b,先解方程组Ly=b,因为L是下三角矩阵，求解只要逐步向前回代。第二步是解方程组Ux=y,因为U是上三角矩阵，用逐次向后回代的方法求出x。3.Matlab程序functionL,U,Y,X=Doolittle(A,b)n=length(b);L=zeros(n,n);U=zeros(n,n);%先将A进行LU分解fori=1:n咐巴第一列第一行求出forj=i:nifi=1A(i,j)=A(i,j);ifj>iA(j,i)=A(j,i)/A(i,i);endendend%下面at算u的第i行值forj=i:nifi=1m=0;fork=1:i

14、-1m=m+A(i,k)*A(k,j);endA(i,j)=A(i,j)-m;endend%下面at算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,进行求解Ab=L,bY=zeros(n,1);Y(1)=Ab(1,n+1)/Ab(1

15、,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运行检测：>>A=111;12-33;-183-1;b=6;15;-15;L,U,Y,X=Doolittle(A,b)运行结果：>>Ab=12.00001.000

16、0015.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.设系数矩阵为三对角矩阵七Ga?b?0a3A=I，-0HI000C2HI000b3HI000+q，q+q，q+0Iaaan1bnCn0HI0anbn则方程组Ax=f称为三对角方程组。2.设矩阵A非奇异，A有Crout分解A=LU,其中

17、L为下三角矩阵,上001H00"%久0川000%P3用00+F+，+r+F+0001H九01000川备儿100+*0<0U为单位上三角矩阵，记工0愕001篇2愕0001川00fqrfbdhbfq，f00愕0L00aIII1先依次求出L,U中的元素后，令Ux=y,先求解下三角方程组Ly=f得出y,再求解上三角方程组Ux=y。求解三对角方程组的2追赶法将矩阵三角分解的计算与求解两个三角方程组的计算放在一起，使算法吏为紧凑。其计算公式为：禺=0，m=,力=5PlPl对i=2,3,|,n支i=ai,R=bi-ai彳，i=卷ifiaiyi'"Xn=ynxti=n-1,n

18、-2,111,1x=yi-ixi13.算法：1 .输入a=(&，an),b=(bi,，bn),c=(c,，cn-i),d=(di,，dn),n2 .对i=2,3,，nai/bi-i=>a(li)bi-c-iai=>bi(ui)di-aidi-1=>d(yi)3 .dn/bn=>dn(xn)4 .对i=n-1,n-2，1(di-odi+1)/bi=>di(x)5 .输出d=x,停机4.Matlab程序：clc;clear;%如果需要验证其他矩阵的解，只需更改A、B和n的值即可A=2,1,0,0;1,3,1,0;0,1,4,-1;0,0,2,3B=3;5;4;

19、5n=4;辛彳4A中不为0的数存到相应数组中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%a变换后的值存入相应矩阵中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,进行求解disp('增广矩阵)；AB=L,BB=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);e