LU分解法、列主元高斯法、Jacobi迭代法、Gauss-Seidel法的原理及Matlab程序

一、实验目的及题目1.1实验目的：（1）学会用高斯列主元消去法，LU分解法，Jacobi迭代法和Gauss-Seidel迭代法解线性方程组。（2）学会用Matlab编写各种方法求解线性方程组的程序。1.2实验题目：1.用列主元消去法解方程组：1241234123412343421233234xxxxxxxxxxxxxxx2.用LU分解法解方程组,Axb其中4824012242412120620266216A，4422b3.分别用Jacobi迭代法和Gauss-Seidel迭代法求解方程组：1232341231234102118311210631125xxxxxxxxxxxxx二、实验原理、程序框图、程序代码等2.1实验原理2.1.1高斯列主元消去法的原理Gauss消去法的基本思想是一次用前面的方程消去后面的未知数，从而将方程组化为等价形式：1111221122222nnnnnnnnbxbxbxgbxbxgbxg这个过程就是消元，然后再回代就好了。具体过程如下：对于1,2,,1kn，若()0,kkka依次计算南昌航空大学数学与信息科学学院实验报告第1页()()(1)()()(1)()()/,,1,,kkikikkkkkkijijikkjkkkiiikkmaaaamabbmbijkn然后将其回代得到：()()()()()1/()/,1,2,,1nnnnnnnkkkkkkjjkkjkxbaxbaxaknn以上是高斯消去。但是高斯消去法在消元的过程中有可能会出现()0kkka的情况，这时消元就无法进行了，即使主元数()0,kkka但是很小时，其做除数，也会导致其他元素数量级的严重增长和舍入误差的扩散。因此，为了减少误差，每次消元选取系数矩阵的某列中绝对值最大的元素作为主元素。然后换行使之变到主元位置上，再进行销元计算。即高斯列主元消去法。2.1.2直接三角分解法（LU分解）的原理先将矩阵A直接分解为ALU则求解方程组的问题就等价于求解两个三角形方程组。直接利用矩阵乘法，得到矩阵的三角分解计算公式为：1111111111,1,2,,/,2,,,,,1,,,2,3,()/,1,2,,iiiikkjkjkmmjmkikikimmkkkmuainlauinualujkknknlaluuikknkn且由上面的式子得到矩阵A的LU分解后，求解Ux=y的计算公式为11111,2,3,/()/,1,2,,1iiiijjjnnnnniiijjiijiybyblyinxyuxyuxuinn以上为LU分解法。南昌航空大学数学与信息科学学院实验报告第2页2.1.3Jacobi迭代法和Gauss-Seidel迭代法的原理（1）Jcaobi迭代设线性方程组bAx(1)的系数矩阵A可逆且主对角元素nna,...,a,a2211均不为零,令nna,...,a,adiagD2211并将A分解成DDAA(2)从而(1)可写成bxADDx令11fxBx其中bDf,ADIB1111.(3)以1B为迭代矩阵的迭代法(公式)111fxBxkk(4)称为雅可比(Jacobi)迭代法,其分量形式为,...,,k,n,...,ixabaxnijj)k(jjiiii)k(i21021111(5)其中Tnx,...x,xx002010为初始向量.（2）Gauss-Seidel迭代由雅可比迭代公式可知,在迭代的每一步计算过程中是用kx的全部分量来计算1kx的所有分量,显然在计算第i个分量1kix时,已经计算出的最新分量1111kikx,...,x没有被利用。把矩阵A分解成ULDA(6)其中nna,...,a,adiagD2211,U,L分别为A的主对角元除外的下三角和上三角部分,南昌航空大学数学与信息科学学院实验报告第3页于是,方程组(1)便可以写成bUxxLD即22fxBx其中bLDf,ULDB1212(7)以2B为迭代矩阵构成的迭代法(公式)221fxBxkk(8)称为高斯—塞德尔迭代法,用分量表示的形式为,...,,k,n,,ixaxabaxijnij)k(jij)k(jijiii)k(i210211111112.2程序代码2.2.1高斯列主元的代码functionGauss(A,b)%A为系数矩阵，b为右端项矩阵[m,n]=size(A);n=length(b);fork=1:n-1[pt,p]=max(abs(A(k:n,k)));%找出列中绝对值最大的数p=p+k-1;ifpkt=A(k,:);A(k,:)=A(p,:);A(p,:)=t;%交换行使之变到主元位置上t=b(k);b(k)=b(p);b(p)=t;endm=A(k+1:n,k)/A(k,k);%开始消元A(k+1:n,k+1:n)=A(k+1:n,k+1:n)-m*A(k,k+1:n);b(k+1:n)=b(k+1:n)-m*b(k);A(k+1:n,k)=zeros(n-k,1);ifflag~=0南昌航空大学数学与信息科学学院实验报告第4页Ab=[A,b];endendx=zeros(n,1);%开始回代x(n)=b(n)/A(n,n);fork=n-1:-1:1x(k)=(b(k)-A(k,k+1:n)*x(k+1:n))/A(k,k);endfork=1:nfprintf('x[%d]=%f\n',k,x(k));end2.2.2LU分解法的程序functionLU(A,b)%A为系数矩阵，b为右端项矩阵[m,n]=size(A);%初始化矩阵A，b，L和Un=length(b);L=eye(n,n);U=zeros(n,n);U(1,1:n)=A(1,1:n);%开始进行LU分解L(2:n,1)=A(2:n,1)/U(1,1);fork=2:nU(k,k:n)=A(k,k:n)-L(k,1:k-1)*U(1:k-1,k:n);L(k+1:n,k)=(A(k+1:n,k)-L(k+1:n,1:k-1)*U(1:k-1,k))/U(k,k);endL%输出L矩阵U%输出U矩阵y=zeros(n,1);%开始解方程组Ux=yy(1)=b(1);fork=2:ny(k)=b(k)-L(k,1:k-1)*y(1:k-1);endx=zeros(n,1);南昌航空大学数学与信息科学学院实验报告第5页x(n)=y(n)/U(n,n);fork=n-1:-1:1x(k)=(y(k)-U(k,k+1:n)*x(k+1:n))/U(k,k);endfork=1:nfprintf('x[%d]=%f\n',k,x(k));end2.2.3Jacobi迭代法的程序functionJacobi(A,b,eps)%A为系数矩阵，b为后端项矩阵，epe为精度[m,n]=size(A);D=diag(diag(A));%求矩阵DL=D-tril(A);%求矩阵LU=D-triu(A);%求矩阵Utemp=1;x=zeros(m,1);k=0;whileabs(max(x)-temp)epstemp=max(abs(x));k=k+1;%记录循环次数x=inv(D)*(L+U)*x+inv(D)*b;%雅克比迭代公式endfork=1:nfprintf('x[%d]=%f\n',k,x(k));end2.2.4Gauss-Seidel迭代程序functionGauss_Seidel(A,b,eps)%A为系数矩阵，b为后端项矩阵，epe为精度[m,n]=size(A);D=diag(diag(A));%求矩阵DL=D-tril(A);%求矩阵LU=D-triu(A);%求矩阵Utemp=1;南昌航空大学数学与信息科学学院实验报告第6页x=zeros(m,1);k=0;whileabs(max(x)-temp)epstemp=max(abs(x));k=k+1;%记录循环次数x=inv(D-L)*U*x+inv(D-L)*b;%Gauss_Seidel的迭代公式endfork=1:nfprintf('x[%d]=%f\n',k,x(k));end三、实验过程中需要记录的实验数据表格3.1第一题（高斯列主元消去）的数据A=[1103;21-11;3-1-13;-123-1];b=[4;1;-3;4];Gauss(A,b)x[1]=-1.333333x[2]=2.333333x[3]=-0.333333x[4]=1.0000003.2第二题（LU分解法）数据A=[48-240-12;-24241212;06202;-66216];b=[4;4;-2;-2];LU(A,b)L=1.0000000-0.50001.00000000.50001.00000-0.12500.2500-0.07141.0000U=48.0000-24.00000-12.0000012.000012.00006.00000014.0000-1.000000012.9286南昌航空大学数学与信息科学学院实验报告第7页x[1]=0.521179x[2]=1.005525x[3]=-0.375691x[4]=-0.2596693.3第三题Jacobi迭代法的数据A=[10-120;08-13;2-1100;-13-111];b=[-11;-11;6;25];Jacobi(A,b,0.00005)x[1]=-1.467396x[2]=-2.358678x[3]=0.657604x[4]=2.8423973.4第三题用Gauss_Seidel迭代的数据A=[10-120;08-13;2-1100;-13-111];b=[-11;-11;6;25];Gauss_Seidel(A,b,0.00005)x[1]=-1.467357x[2]=-2.358740x[3]=0.657597x[4]=2.842405四、实验中存在的问题及解决方案4.1存在的问题（1）第一题中在matlab中输入Gauss(A,b)（数据省略）得到m=4n=4???UndefinedfunctionorvariableAb.Errorin==Gaussat8[ap,p]=max(abs(Ab(k:n,k)));没有得到想要的结果。（2）第二题中在matlab中输入y=LU(A,b)得到y=4.00006.0000-5.0000-3.3571不是方程组的解。（3）第三题中在用高斯赛德尔方法时在matlab中输入Gauss-Seidel(A,b,eps)结果程序报错???Errorusing==GaussToomanyoutputarguments.得不到想要的结果。南昌航空大学数学与信息科学学院实验报告第8页4.2解决方案（1）针对第一题中由于程序的第二行漏了一个分号导致输出了m和n的值，第8行中将Ab改为A问题就解决了。（2）由于程序后面出现了矩阵y故输出的事矩阵y的值，但是我们要的事x的值，故只需要将y改成x，或者直接把y去掉就解决了问题。（3）在function文件中命名不能出现“-”应该将其改为下划线“_”，所以将M文件名“Gauss-Seidel(A,b,eps)”改成“Gauss_Seidel(A,b,eps)”就解决问题了。五、心得体会本次试验涉及到了用高斯列主元消去法，LU分解法，Jacobi迭代法以及Ga