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

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

       一、实验目的及题目

       1.1 实验目的：

       （1）学会用高斯列主元消去法，LU分解法，Jacobi迭代法和Gauss-Seidel迭代法解线性方程组。

       （2）学会用Matlab编写各种方法求解线性方程组的程序。1.2 实验题目：

       1.用列主元消去法解方程组：

       x1x23x442xxxx1123

       4

       2xxx3x34123x12x23x3x442.用LU分解法解方程组Axb,其中

       1248240424241212，b4

       A0262026621623.分别用Jacobi迭代法和Gauss-Seidel迭代法求解方程组：

       10x1x22x3118xx3x11234

       

       2xx10x6312x13x2x311x425

       二、实验原理、程序框图、程序代码等

       2.1实验原理

       2.1.1高斯列主元消去法的原理

       Gauss消去法的基本思想是一次用前面的方程消去后面的未知数，从而将方程组化为等价形式：

       b11x1b12x2b1nxng1b22x2b2nxng2 bnnxngn这个过程就是消元，然后再回代就好了。具体过程如下：

       (k)对于k1,2,,n1，若akk0,依次计算

       南昌航空大学数学与信息科学学院实验报告

       (k)(k)mikaik/akk(k1)(k)(k)aijaijmikakj

       bi(k1)bi(k)mikbk(k),i,jk1,,n然后将其回代得到：

       (n)(n)xnbn/annn (k)(k)(k)x(bax)/a,kn1,n2,,1kkjjkkkjk1以上是高斯消去。

       (k)但是高斯消去法在消元的过程中有可能会出现akk0的情况，这时消元就无法进行了，(k)即使主元数akk0,但是很小时，其做除数，也会导致其他元素数量级的严重增长和舍入误差的扩散。因此，为了减少误差，每次消元选取系数矩阵的某列中绝对值最大的元素作为主元素。然后换行使之变到主元位置上，再进行销元计算。即高斯列主元消去法。2.1.2直接三角分解法（LU分解）的原理

       先将矩阵A直接分解为ALU则求解方程组的问题就等价于求解两个三角形方程组。直接利用矩阵乘法，得到矩阵的三角分解计算公式为：

       u1ia1i,i1,2,,nli1ai1/u11,i2,,nk1 ukjakjlkmumj，jk,k1,,nm1,k2,3,nk1l(alu)/u,ik1,k2,,n且knikikimmkkkm1由上面的式子得到矩阵A的LU分解后，求解Ux=y的计算公式为

       y1b1i1yibilijyj,i2,3,nj1

       xnyn/unnnxi(yiuijxj)/uii,in1,n2,,1ji1以上为LU分解法。

       第 页

       南昌航空大学数学与信息科学学院实验报告

       2.1.3Jacobi迭代法和Gauss-Seidel迭代法的原理（1）Jcaobi迭代

       设线性方程组

       Axb(1)的系数矩阵A可逆且主对角元素a11,a22,...,ann均不为零,令

       Ddiaga11,a22,...,ann并将A分解成

       AADD(2)从而(1)可写成

       DxDAxb 令

       xB1xf1

       11BIDA,fDb.(3)1其中1

       以B1为迭代矩阵的迭代法(公式)

       xk1B1xkf1(4)称为雅可比(Jacobi)迭代法,其分量形式为

       x(k1)i1aiibiaijx(jk)j1jini1,2,...n,Tk0,1,2,...(5)

       0000为初始向量.xx,x,...x12n其中（2）Gauss-Seidel迭代

       由雅可比迭代公式可知,在迭代的每一步计算过程中是用x所有分量,显然在计算第i个分量xi用。

       把矩阵A分解成

       k1k的全部分量来计算xk1k1的时,已经计算出的最新分量x1,...,xi1k1没有被利

       ADLU(6)其中Ddiaga11,a22,...,ann,L,U分别为A的主对角元除外的下三角和上三角部分,第 页

       南昌航空大学数学与信息科学学院实验报告

       于是,方程组(1)便可以写成

       DLxUxb 即

       xB2xf2

       其中

       B2DLU,1f2DLb(7)

       1以B2为迭代矩阵构成的迭代法(公式)

       xk1B2xkf2(8)称为高斯—塞德尔迭代法,用分量表示的形式为

        x(k1)ii1n1(k1)biaijxjaijx(jk)j1ji1aii i1,2,n,k0,1,2,...2.2程序代码

       2.2.1高斯列主元的代码

       function Gauss(A,b)

       %A为系数矩阵，b为右端项矩阵 [m,n]=size(A);n=length(b);for k=1:n-1

       [pt,p]=max(abs(A(k:n,k)));

       %找出列中绝对值最大的数

       p=p k-1;

       if p>k

       t=A(k,:);A(k,:)=A(p,:);A(p,:)=t;

       %交换行使之变到主元位置上

       t=b(k);b(k)=b(p);b(p)=t;

       end

       m=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);if flag~=0

       第 页

       南昌航空大学数学与信息科学学院实验报告

       Ab=[A,b];end end

       x=zeros(n,1);

       %开始回代

       x(n)=b(n)/A(n,n);

       for k=n-1:-1:1

       x(k)=(b(k)-A(k,k 1:n)*x(k 1:n))/A(k,k);

       end

       for k=1:n

       fprintf('x[%d]=%fn',k,x(k));

       end 2.2.2 LU分解法的程序

       function LU(A,b)

       %A为系数矩阵，b为右端项矩阵 [m,n]=size(A);

       %初始化矩阵A，b，L和U n=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);for k=2:n

       U(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);

       end L

       %输出L矩阵 U

       %输出U矩阵

       y=zeros(n,1);

       %开始解方程组Ux=y y(1)=b(1);for k=2:n

       y(k)=b(k)-L(k,1:k-1)*y(1:k-1);end x=zeros(n,1);

       第 页

       南昌航空大学数学与信息科学学院实验报告

       x(n)=y(n)/U(n,n);for k=n-1:-1:1

       x(k)=(y(k)-U(k,k 1:n)*x(k 1:n))/U(k,k);end for k=1:n

       fprintf('x[%d]=%fn',k,x(k));end 2.2.3 Jacobi迭代法的程序

       function Jacobi(A,b,eps)

       %A为系数矩阵，b为后端项矩阵，epe为精度 [m,n]=size(A);D=diag(diag(A));

       %求矩阵D L=D-tril(A);

       %求矩阵L U=D-triu(A);

       %求矩阵U temp=1;x=zeros(m,1);k=0;while abs(max(x)-temp)>eps

       temp=max(abs(x));

       k=k 1;

       %记录循环次数

       x=inv(D)*(L U)*x inv(D)*b;

       %雅克比迭代公式 end

       for k=1:n

       fprintf('x[%d]=%fn',k,x(k));end 2.2.4 Gauss-Seidel迭代程序

       function Gauss_Seidel(A,b,eps)

       %A为系数矩阵，b为后端项矩阵，epe为精度 [m,n]=size(A);D=diag(diag(A));

       %求矩阵D L=D-tril(A);

       %求矩阵L U=D-triu(A);

       %求矩阵U temp=1;

       第 页

       南昌航空大学数学与信息科学学院实验报告

       x=zeros(m,1);k=0;while abs(max(x)-temp)>eps

       temp=max(abs(x));

       k=k 1;

       %记录循环次数

       x=inv(D-L)*U*x inv(D-L)*b;

       %Gauss_Seidel的迭代公式 end

       for k=1:n

       fprintf('x[%d]=%fn',k,x(k));end

       三、实验过程中需要记录的实验数据表格

       3.1第一题（高斯列主元消去）的数据

       >> A=[1 1 0 3;2 1-1 1;3-1-1 3;-1 2 3-1];>> b=[4;1;-3;4];>> Gauss(A,b)x[1]=-1.333333 x[2]=2.333333 x[3]=-0.333333 x[4]=1.000000

       3.2 第二题（LU分解法）数据

       >> A=[48-24 0-12;-24 24 12 12;0 6 20 2;-6 6 2 16];>> b=[4;4;-2;-2];>> LU(A,b)L =

       1.0000

       0

       0

       0

       -0.5000

       1.0000

       0

       0

       0

       0.5000

       1.0000

       0

       -0.1250

       0.2500

       -0.0714

       1.0000 U =

       48.0000-24.0000

       0-12.0000

       0

       12.0000

       12.0000

       6.0000

       0

       0

       14.0000

       -1.0000

       0

       0

       0

       12.9286

       第 页

       南昌航空大学数学与信息科学学院实验报告

       x[1]=0.521179 x[2]=1.005525 x[3]=-0.375691 x[4]=-0.259669

       3.3第三题Jacobi迭代法的数据

       >> A=[10-1 2 0;0 8-1 3;2-1 10 0;-1 3-1 11];b=[-11;-11;6;25];Jacobi(A,b,0.00005)x[1]=-1.467396 x[2]=-2.358678 x[3]=0.657604 x[4]=2.842397

       3.4第三题用Gauss_Seidel迭代的数据

       >> A=[10-1 2 0;0 8-1 3;2-1 10 0;-1 3-1 11];>> b=[-11;-11;6;25];>> Gauss_Seidel(A,b,0.00005)x[1]=-1.467357 x[2]=-2.358740 x[3]=0.657597 x[4]=2.842405

       四、实验中存在的问题及解决方案

       4.1存在的问题

       （1）第一题中在matlab中输入>> Gauss(A,b)（数据省略）得到m =4 n =4 ??? Undefined function or variable “Ab”.Error in ==> Gauss at 8[ap,p]=max(abs(Ab(k:n,k)));没有得到想要的结果。

       （2）第二题中在matlab中输入>> y=LU(A,b)得到y =4.0000

       6.0000

       -5.0000

       -3.3571不是方程组的解。

       （3）第三题中在用高斯赛德尔方法时在matlab中输入>> Gauss-Seidel(A,b,eps)结果程序报错??? Error using ==> Gauss

       Too many output arguments.得不到想要的结果。

       第 页

       南昌航空大学数学与信息科学学院实验报告

       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迭代法以及Gauss-seidel迭代法等四种方法。需要对这些方法的原理都要掌握才能写出程序，由于理论知识的欠缺，我花了很大一部分时间在看懂实验的原理上，看懂了实验原理之后就开始根据原理编写程序，程序中还是出现了很多的低级错误导致调试很久才能运行。

       通过这次试验使我深刻的体会到理论知识的重要性，没有理论知识的支撑是写不出程序来的。写程序时还会犯很多低级的错误，以后一定要加强理论知识的学习，减少编程时低级错误的产生。

       第 页

第二篇：列主元高斯消去 LU分解 迭代法

       数学与信息科学学院

       实 验 报 告

       课程名称： 《 数值计算方法 》 实验名称： 数值积分 实验类型： 验证性■综合性□设计性□ 实验室名称： 数学实验室 班级学号： 090721 学生姓名：

       任课教师（教师签名）： 成 绩：

       实验日期： 2022-4-27

       一、实验目的及题目

       实验目的：

       南昌航空大学数学与信息科学学院实验报告

       熟悉线性方程组求解原理，运用列主消元高斯消去法，LU分解法及Jacobi迭代法与Gauss-Seidel迭代法丢出线性方程的解 实验题目：

       1、用列主元消去法解方程组：

       x1x23x442xxxx1123

       4

       3xxx3x32341x12x23x3x44

       2、用LU分解法界方程组Axb，其中

       48240124122412124

       A,b

       062022626216

       3、分别用jacobi迭代法和gauss-seidel迭代法求解方程组：

       10x1x22x3118xx3x1123

       4

       2xx10x6231x13x2x311x425

       二、实验原理、程序框图、程序代码等

       实验原理：

       1、列主元高斯消去原理：

       每次消去中，选取每列的绝对值最大的元素作为消去对象并作为主元素。然后换行使之变到主元位子上，在进行消元计算。设A(k)xb(k)，确定第k列主元所在位置ik，在交换ik行和k行后，在进行消元。根据矩阵理论，交换ik和k两个方程的位置，列主元素的消去过程相当于对交换后的新矩阵进行消元，即

       LkIk,ikAkA(k1)

       同时，右端向量b(k)变化为 LkIk,ikbkb(k1)

       2、LU分解原理：

       第 页

       南昌航空大学数学与信息科学学院实验报告

       直接三角分解：先将A分解为上三角矩阵U和下三角矩阵L（A=LU），则原为题等价于求解两个方程组：Lyb,Uxy

       3、迭代法

       （1）jacobi迭代法（又称简单迭代法）考虑n阶线性代数方程组

       a11x1a12x2.....a1nxnb1axax.....axb2112222nn

       2

       an1x1an2x2.....annxnbn其矩阵形式为

       Axb

       设该方程组的系数矩阵A非奇异且aii0(i1,2,.....,n)，可将A分解为：

       ADLU

       其中Ddiag[a11,a22,.....,ann]；

       0a21La31an1 0a31an10an10a120，U=－0a13a230a1na2n

       an1,n0然后化为如下等价形式

       xD1(LU)xD1b

       简记为

       xBJxfJ

       选取初始向量x(0)x1,x2(0)(0)xn(0)T,通过上式可得到线性代数方程组的迭代格式

       x(k1)BJx(k)fJ,(k0,1,2)

       其分量形式为

       xi即(k1)n1k(biaijxj),(i1,2,n)(k0,1,2)aiii1j1第 页

       南昌航空大学数学与信息科学学院实验报告

       (k1)1(k)(k)(k)x(baxaxax)11221331nn1a111(k1)(k)(k)(k)x(baxaxax)222112332nna 221(k1)(k)(k)(k)x(baxaxax)nnn1n2n,n1n112ann

       以上计算过程称迭代法，矩阵BJ陈为jacobi迭代法的迭代矩阵。

       （2）Gauss-Seidel迭代法 将jacobi迭代格式改为如下形式

       xi(k1)i1n1(k)(k)(biaijxjaijxj),(i1,2,n)(k0,1,2)aiii1ji1其矩阵形式为

       x(k1)(DL)1Ux(k)(DL)1b

       于是得Gauss-Seidel迭代公式为

       x(k)BGx(k)fG,(k0,1,2)

       程序代码： 第一题程序代码： function GEpiv(A,b)[m,n]=size(A);nb=n 1;Ab=[A b];for i=1:m-1 [pivot,p]=max(abs(Ab(i:n,j)));ip=p i-1;if ip~=i Ab([i ip],:)=Ab([ip i],:);end

       第 页 称BG为Gauss-Seidel迭代法的迭代矩阵。

       南昌航空大学数学与信息科学学院实验报告

       pivot=Ab(i,i);for k=i 1:m Ab(k,i:nb)=Ab(k,i:nb)-(Ab(k,i)/pivot)*Ab(i,i:nb);end end x=zeros(n,1);x(n)=Ab(n,nb)/Ab(n,n);for i=n-1:-1:1 x(i)(Ab(i,nb)-Ab(i,i 1:n)*x(i 1:n,1))/Ab(i,i);end for k=1:n fprintf('x[%d]=%fn',x(k));end 运行结果为 输入

       A=[1 1 0 3;2 1-1 1;3-1-1 3;-1 2 3-1];b=[4;1;-3;4];得结果为 x =

       -1.3333 2.3333-0.3333 1.0000

       第二题代码：

       function x=lufj(A,b)%lu分解求值函数 [l u p]=lu(A);y=lx(l,b);x=us(u,y);end

       第 页

       南昌航空大学数学与信息科学学院实验报告

       function y=lx(A,b)%求系数矩阵为下三角方程函数 [n,m]=size(A);y(1)=b(1)/A(1,1);for i=2:n s=0;for k=1:i-1 s=s A(i,k)*y(k);end y(i)=(b(i)-s)/A(i,i);end end

       function x=us(A,b)%求系数矩阵为上三角方程函数 [n,m]=size(A);x(n)=b(n)/A(n,n);for i=0:n-2 s=0;for k=0:i s=s A(n-i-1,n-k)*x(n-k);end x(n-i-1)=(b(n-i-1)-s)/A(n-i-1,n-i-1);end end 输入矩阵

       A=[48-24 0 12;-24 4 12 12;0 6 20 2;-6 6 2 16];b=[4;4;-2;-2];结果为 x =

       第 页

       南昌航空大学数学与信息科学学院实验报告

       0.5212 1.0055-0.3757-0.2597

       第三题代码：（1）Jacobi迭代：

       function x = agui_jacobi(a,b)n=length(b);N=100;e=1e-4;x0=zeros(n,1);x=x0;x0=x 2*e;k=0;d=diag(diag(a));l=-triu(a,-1);u=-triu(a,1);while norm(x0-x,inf)>e&k

       （2）Gauss-seidel迭代 function x=gau(A,b,x)[d u l]=fenjie(A);B=inv(d-l)*u;f=inv(d-l)*b;

       第 页

       南昌航空大学数学与信息科学学院实验报告

       for i=1:9 y=x;x=B*y f;end

       function [d u l]=fenjie(A)[n m]=size(A);for i=1:n for j=1:n d(i,j)=0;end end for i=1:n d(i,i)=A(i,i);end for i=1:n for j=1:i l(i,j)=-A(i,j);u(j,i)=-A(j,i);end for k=i:n l(i,k)=0;u(k,i)=0;end end（1）

       输入矩阵：a=[10-1 2 0;0 8-1 3;2-1 10 0;-1 3-1 11];>> b=[-11;-11;6;25];>> x = agui_jacobi(a,b)结果为

       第 页

       南昌航空大学数学与信息科学学院实验报告

       x=-1.4674-2.3587 0.6576 2.8424

       （2）A=[10-1 2 0;0 8-1 3;2-1 10 0;-1 3-1 11];b=[-11;-11;6;25];x=[-1.3;-2.1;0.55;2.66];gau(A,b,x)结果为 x=-1.4674-2.3587 0.6576 2.8424

       四、实验中存在的问题及解决方案

       命令需多次修改，如Gauss迭代中，m文件中因end的不匹配使得程序无法运行处结果，则需细心调试，去除多余的end，才能使程序正常运行。

       五、心得体会

       在线性代数中，从理论上解决了线性方程组的求解问题，但常用的高斯消元法有时会导致结果产生严重的偏差，因此需要研究数值解法。Jacobi迭代与Gauss-Seidel有一个共同特点：新的近似解是已知近似解的线性函数，并且新的近似解只与已知近似解有关。

       编写程序是一个繁杂的过程，需要不断调试，并且善于发现程序中的细小问题，才能够正常运行处程序。

       第 页