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

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

第 页