线性代数方程组数值解法及MATLAB实现综述

线性代数方程组数值解法及MATLAB实现综述

廖淑芳 20122090 数计学院 12计算机科学与技术1班（职教本科） 一、分析课题

随着科学技术的发展，提出了大量复杂的数值计算问题，在建立电子计算机成为数值计算的主要工具以后，它以数字计算机求解数学问题的理论和方法为研究对象。其数值计算中线性代数方程的求解问题就广泛应用于各种工程技术方面。因此在各种数据处理中，线性代数方程组的求解是最常见的问题之一。关于线性代数方程组的数值解法一般分为两大类：直接法和迭代法。

直接法就是经过有限步算术运算，可求的线性方程组精确解的方法（若计算过程没有舍入误差），但实际犹如舍入误差的存在和影响，这种方法也只能求得近似解，这类方法是解低阶稠密矩阵方程组级某些大型稀疏矩阵方程组的有效方法。直接法包括高斯消元法，矩阵三角分解法、追赶法、平方根法。

迭代法就是利用某种极限过程去逐步逼近线性方程组精确解的方法。迭代法具有需要计算机的存储单元少，程序设计简单，原始系数矩阵在计算过程始终不变等优点，但存在收敛性级收敛速度问题。迭代法是解大型稀疏矩阵方程组（尤其是微分方程离散后得到的大型方程组）的重要方法。迭代法包括Jacobi法SOR法、SSOR法等多种方法。

二、研究课题-线性代数方程组数值解法 一、 直接法 1、 Gauss消元法

通过一系列的加减消元运算，也就是代数中的加减消去法，以使A对角线以下的元素化为零，将方程组化为上三角矩阵；然后，再逐一回代求解出x向量。

. 可修编.

. -

1.1消元过程

1. 高斯消元法(加减消元)：首先将A化为上三角阵，再回代求解。

a(1)12a(1)13a(1)1nb(1)1a11a12a1nba(1)111a(2)a22a(2)23a(2)2nb(2)221a22a2nb020a(3)33a(3)3nb(3)30 an1an2annbn000a(n)nnb(n)n步骤如下：

第一步：第1行ai1a第i行,i2,,n

11a11a12a1nb1a1nb1a21a22aba11a12(2)2n2a(2)b(2)2n20a22 an1an2annb(2)n0an2a(2)nnb(2)n第二步：2行a(2)第i2a(2)第i行,i3,,n 22a11a12a1ba11a12a13a1nb1n10a(2)22a(2)b(2)0a(2)22a(2)23a(2)2nb(2)22n200a(3)33a(3)3nb(3)3 0a(2)n2a(2)nnb(2)n00a(3)a(3)n3nnb(3)n类似的做下去，我们有：

第k步：(第k行ak)ika(k)第i行,ik1,,n。

kkn－1步以后，我们可以得到变换后的矩阵为：

a11a12a13a1nb10a(2)22a(2)23a(2)2nb(2)200a(3)33a(3)3nb(3)3 000a(n)nnb(n)n . 可修编.

. -

(k)注意到，计算过程中akk处在被除的位置，因此整个计算过程要保证它不为0。

所以，Gauss消元法的可行条件为：a(k)kk0。 就是要求A的所有顺序主子式均不为0，即

a11a1idet0,i1,,n

ai1aii因此，有些有解的问题，不能用Gauss消元求解。

另外，如果某个a(k)kk很小的话，会引入大的误差。

例 用Gauss消去法解方程组：

12334x115（1）18311x11112x15 331116x42（1）对增广矩阵进行初等变换

123415(32E1E2)E2712334153(112E1E3)E3035151831115(1E1E4)E422243219

111160531112443470474074123341233(56E)E372E33051503715(76E2E4)E4220011292(711E3E4)E42236110011300735736000 . 451515292 61191033 可修编.

. -

得等价方程组

12x13x23x34x4153x27x35x415222 1129xx11343691x4033回代得x40，x33，x22，x11。

第一步：将(1)/3使x1的系数化为1,再将(2)、(3)式中x1的系数都化为零，即由(2)-2×(1)得

(1)

x121x2x32......(1)(1)33由(3)-4×(1)(1)得

24x2x30......(2)(1)331410x2x36......(3)(1)33x22x30......(2)(2)第二步：将(2)除以2/3，使x2系数化为1，得

(1)

再将(3)(1)式中x2系数化为零，由(3)(1)-(-14/3)*(2)(2)，得

18x36......(3)(2)3

x31......(3)(3)第三步：将(3)除以18/3，使x3系数化为1，得

(2)

经消元后，得到如下三角代数方程组：

. 可修编.

. -

1.2回代过程

由(3)(3)得 x3=1,将x3代入(2)(2)得x2=-2,将x2 、x3代入(1)(1)得x2=1,

所以，本题解为[x]=[1,2,-1]T

1.3高斯消元的公式

综合以上讨论，不难看出，高斯消元法解方程组的公式为 第一步，消元

（1） 令

aij(1) = aij , （i,j=1,2,3,…,n） bi(1) =bi , （i=1,2,3,…,n） （2） 对k=1到n-1，若akk(k)≠0，进行

lik = aik / akk , （i=k+1,k+2,…,n） aij(k+1) = aij(k) - lik * akj(k), （i,j= k+1,k+2,…,n） bi(k+1) = bi(k) - lik * bk(k), （i= k+1,k+2,…,n）

第二步，回代

若ann(n) ≠ 0

xn = bn(n) / ann(n)

xi = (bi(i) – sgm(aij(i) * xj)/- aii(i) ,（i = n-1,n-2,…,1）,( j = i+1,i+2,…,n )

2 、LU分解法

求解线性代数方程组除了高斯消元法外，还常用LU分解法（三角形分解法）。LU分解法的优点是当方程组左端系数矩阵不变，仅仅是方程组右端列向量改变，

. 可修编.

(k)

(k)

. -

即外加激励信号变化时，能够方便地求解方程组。矩阵的三角分解法可分为直接三角分解法,列主元三角分解法,平方根法,三对角方程组的追赶法。下面讨论直接三角分解法。

设n阶线性方程组Ax=b 。假设能将方程组左端系数矩阵A，分解成两个三角阵的乘积，即A=LU ，式中，L为主对角线以上的元素均为零的下三角矩阵, 且主对角线元素均为1的上三角矩阵；U为主对角线以下的元素均为零

所以有，LUx=b令 Ux=y， 则 Ly=b

由A=LU,由矩阵的乘法公式： a1j = u1j , j=1,2,…,n

ai1 = li1u11 , i=1,2,…,n

推出 u1j = a1j, j=1,2,…,n

li1 = ai1/u11, i=1,2,…,n

这样就定出了U的第一行元素和L的第一列元素。

设已定出了U的前k-1行和L的前k-1列，现在确定U的第k行和L的第k列。由矩阵乘法：

. 可修编.

. -

akjlkrurjn当r>k时，lkr=0, 且lkk=1，因为

r1akjukjlkrurjn

ukjakjlkrurjnr1r1jk,k1,...,nn所以，

lik(aiklkrurj)/ukkik,k1,...,n同理可推出计算L的第k列的公式：

r1因此得到如下算法——杜利特（Doolittle）算法:

(1) 将矩阵分解为A=LU,对k=1,2,…,n；j=k,k+1,…n； i=k,k+1,…n；

nukiakjlkiuyj r1nljk(aiklkyurj)/ukkr1公式1  l  1

kk

(2) 解Ly=b

公式2:ykbklkryrr1k1k1,2,...,n

(3) 解Ux=y

公式3：xk(ykrk1unkrrx)/ukkkn,n1,...,1

对大规模稀疏问题，如果能够通过调整方程及未知量的顺序使得方程组的系数矩阵成带状结构，则对系数矩阵使用通常的LU分解，可以保障单位下三角矩阵L及上三角矩阵U仍为带状结构.

. 可修编.

. -

3、直接三角分解法

Gauss消去法还有许多变形，有些变形是为了利用特殊技巧减少误差，把Gauss消去法改写为更紧凑的形式，还有一些变形时根据某类矩阵的特性作一些修正和简化，这些方法可统称为直接三角分解法。

矩阵的三角分解

设A的顺序主子式k0(k1,2,,n)，则可建立线性方程组Axb的Gauss

消去法与矩阵分解的关系，即矩阵A的LU分解。这个问题前面已经讲的比较详细了，此处不再赘述。

Doolittle分解法

首先假设A的顺序主子式k都不为零，则A可作Doolittle分解，即ALU，其中L是单位下三角阵，有lii1，ij时lij0；U是上三角阵，ij时uij0。仔细写出为

a11a12a21a22Aan1an2a1n1a2nl211annln1ln2u11u12u221u1nu2nLU（2.11） unn(k)

在前面逐步推导L和U的元素公式都要借助于有关的aij来表示。现在强调

指出，只要从给定的A通过比较（2.11）式的两边就可能逐步地把L和U构造出来，而不必利用Gauss消去法的中间结果，这种方法称为Gauss消去法的紧凑格式。

根据矩阵的乘法规则，比较（2.11）式的两边可得

min(i,j)aijk1likukj，i,j1,2,,n （2.12）

先算出U的第1行，在（2.12）令i1，由l111，l1j0(1jn)可得

. 可修编.

. -

u1ja1j，j1,2,,n

接着在（2.12）中令j1，由ai1li1u11，从而算出L的第1列

li1ai1/u11ai1/a11，i2,,n

若U的第1至k1行和L的第1至k1列已经算出，则由（2.12）可计算出

U的第k行元素

ukjakjlkrurj，jk,k1,r1k1,n （2.13）

接着再计算出L的第k列元素

likaiklirurkr1k1ukk，jk1,,n （2.14）

解方程组Axb就化为求解LUxb，分两步求解。第一步解Lyb，其中L为下三角阵，只要用逐次向前代入的方法；第二步解Uxy，其中U为上三角阵，用逐次向后代入方法即可解除x。

例 用Doolittle方法求解：

62111x16410x21 141x35013x4521解：第1步算出L和U：

1113L11561161062103，U1913712337101139 1019174第2步解Lyb得：

. 可修编.

. -

23191y6,3,,

574第3步解Uxy得：

x(1,1,1,1)T

T

二、 迭代法

迭代法具有的特点是速度快。与非线性方程的迭代方法一样，需要我们构造一个等价的方程，从而构造一个收敛序列，序列的极限值就是方程组的根。 对方程组：Axb做等价变换：xGxg 如：令AMN，则

Axb(MN)xbMxbNxxM1NxM1b

则，我们可以构造序列x(k1)G x(k)g 若x(k)x*x*G x*gAx*b 则，我们可以构造序列x(k1)G x(k)g 若x(k)x*x*G x*gAx*b 同时：x(k1)x*Gx(k)Gx*G(x(k)x*)

Gk1(x(0)x*)

所以，序列收敛Gk0 与初值的选取无关 1. Jacobi迭代

a11x1axn11a1nxnb1

annxnbn . 可修编.

. -

1x1a(a12x2a1nxnb1)111(a21x1a23x3a1nxnb2)x2a22  1xna(an1x1an n1xn1bn)nn(k1)1(k)(k)x(axaxb1)11221nna11(k1)1(a21x1(k)a23x3(k)a1nxn(k)b2)x2a22  (k1)1(k)(k)x(axaxbn)nn11n n1n1ann格式很简单：

xi(k1)1i1(aijxj(k)aiij1ji1axijn(k)jbi)

迭代矩阵

记ADLU

a11D00 ann0a21Lan100ann100a120U00a1n 0an1n0易知，Jacobi迭代有

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

 B0D1(LU)ID1A , f0D1b，B0是简单迭代的迭代矩阵。

. 可修编.

. -

看上述公式和过程很抽象，我们来举个简单例子：

a11x1a12x2a13x3b1a21x1a22x2a23x3b2 （aii0） axaxaxb31132233331x(b1a12x2a13x3)1a111得x2(b2a21x1a23x3)

a221(b3a31x1a32x2)x3a33(0)(0)T,x3)作初始近似，代入， 得上述变换的方程后，我们任取一向量x(0)(x1(0),x2(1)(1)Tx(1)(x1(1),x2,x3),(m)(m)T,x(m)(x1(m),x2,x3)

假设上述方程的准确解是a(a1,a2,a3)T

那么如果limxma,我们就说上述迭代是收敛的。

m2. Gauss－Seidel迭代

在Gauss－Seidel迭代中，使用最新计算出的分量值

(k1)1(k)(k)x(axaxb1)1221nn1a11(k1)1(a21x1(k1)a23x3(k)a1nxn(k)b2)x2a22  (k1)1(an1x1(k1)an n1xn1(k1)bn)xnannxi(k1)1i1(aijxj(k1)aiij1ji1axijn(k)jbi)

易知，Gauss－Seidel迭代有

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

. 可修编.

. -

 G0D1(LU)ID1A , gD1b，G0是简单迭代的迭代矩阵。

二种方法都存在收敛性问题。 收敛条件

迭代格式收敛的充要条件是BG的谱半径<1。

有例子表明：Gauss-Seidel法收敛时，Jacobi法可能不收敛；而Jacobi法收敛时， Gauss-Seidel法也可能不收敛。

1若为严格对角占优的话则是收敛的，如12130.10.52，0.0250.3 10.030.14假如方程组不满足收敛，有时候对其进行变换，可以改变收敛性。 如求下述方程组的解：

1807A109,b8,x(0)(0,0,0)T

91179117A180,b7

1098格式

(k)(k)x1(k1)1/9(x2x37)(k1)1/8(x1(k1)7) x1x(k1)1/9(x(k1)8)11结果

x(1)(0.7778,0.9722,0.9753)Tx(2)(0.9942,0.9993,0.9994)Tx(3)(0.9999,0.9999,0.9999)T

x(4)(1.0000,1.0000,1.0000)T(B0)1，（G）<1也满足收敛。

. 可修编.

. -

122如A111

22122Jacobi，B0的特征方程为110

22 解得1230，所以用Jacobi迭代法收敛。

22Gauss-Seidel，G0的特征方程为10 22 解得10,228,328，(B0)281, 所以Gauss-Seidel迭代法是不收敛的。

3. 松弛迭代 记x(k)x(k1)x(k) 则x(k1)x(k)x(k)

可以看作在前一步上加一个修正量。若在修正量前乘以一个因子，有x(k1)x(k)x(k) x(k1)x(k)(x(k1)x(k))

对Gauss－Seidel迭代格式x(k1)D1(Lx(k1)Ux(k)b)

x(k1)x(k)(D1Lx(k1)D1Ux(k)D1bx(k)) x(k1)(1)x(k)(D1Lx(k1)D1Ux(k)D1b)

写成分量形式，有

.

可修编.

. -

1(k1)(k)(k)(k)x(1)x(axaxb1)111221nna11(k1)(1)x2(k)(a21x1(k)a23x3(k)a1nxn(k)b2)x2a22 x(k1)n(1)x(k)n1a(an1x(k)1an n1x(k)n1bn)nn关于直接法和迭代法的例题： 1用选主元三角分解法求解下列方程组

3x13x22x34x462x14x22x3103x 2x12x24x34174x12x22x36x483324x12420332462420102243x62x1034226x17224317列主

484226826422684226842242010113224317123361123833246213013121331343212120341201100042261100这里L20313811,U2310008,y6 341111320100013再通过Uxy，得x的解。

x1x2422x3 3x342用Jacobi迭代法求解下列方程组，并且使精度为103。

. 8611

3 可修编.

. -

40.240.08x180.0930.15x29以4，3，4分别除3方程两边 0.040.084x32010.060.02x120.0310.05x23 0.010.021x35x10.060.02x12x020.0300.05x23 x30.010.020x35x(m1)00.060.02x(m)11Jacobi迭代式2x(m1)20.0300.05x(m)23 x(m1)30.010.020x(m)35

40.240.080.060.020.0930.150.030.040.0840.050.010.020.060.02020.00180.050.0006 0(0.06)(320.010.0009)03(20.0018)0.010.010910.06,26,0.010.010936

Jacobi迭代法是收敛的。 x(0)(2,3,5)T为初始值，迭代得，

x(1)100.060.02221.92x(1)20.0300.05x(1)0.010.020333.19

3555.04x(2)100.060.021.92x(2)20.0300.053.19231.909x(2)0.010.0203.1945.045…

35.045 . 可修编.

即有可证明此迭代格式是收敛的，极限值即为解。解得由此可知，用取 . -

x(3)00.060.021.90921.9091(3)00.053.19433.194 其实x20.0355.045x(3)0.010.0205.0453即在103，x(3)x(2)，可以停止计算了。所以x(3)即为近似解。

分析：此题题目可以直接说是解方程组，然后求解过程我们可以先验算Jacobi和Gauss－Seidel迭代的收敛情况，再选择算法，一般来说Gauss－Seidel迭代要比Jacobi迭代快。

三、用Matlab 软件求解线性方程组 一．高斯消去法 1.顺序高斯消去法 直接编写命令文件 a=[] d=[]' [n,n]=size(a); c=n+1 a(:,c)=d; for k=1:n-1

a(k+1:n, k:c)=a(k+1:n, k:c)-(a(k+1:n,k)/ a(k,k))*a(k, k:c); 消去 end

x=[0 0 0 0]' 回带 x(n)=a(n,c)/a(n,n); for g=n-1:-1:1

x(g)=(a(g,c)-a(g,g+1:n)*x(g+1:n))/a(g,g)

. 可修编.

. -

end

2.列主高斯消去法

*由于“[r,m]=max(abs(a(k:n,k)))”返回的行是“k:n,k”的第几行，所以要通过修正来把m 改成真正的行的值。该程序只是演示程序，真正机器计算不需要算主元素所在列以下各行应为零的值。 直接编写命令文件 a=[] d=[] ' [n,n]=size(a); c=n+1

a(:,c)=d; (增广) for k=1:n-1

[r,m]=max(abs(a(k:n,k))); 选主

m=m+k-1; (修正操作行的值) if(a(m,k)~=0)

if(m~=k)

a([k m],:)=a([m k],:); 换行 end

a(k+1:n, k:c)=a(k+1:n, k:c)-(a(k+1:n,k)/ a(k,k))*a(k, k:c); end end

x=[0 0 0 0]' 回带

. %消去 可修编.

. -

x(n)=a(n,c)/a(n,n); for g=n-1:-1:1

x(g)=(a(g,c)-a(g,g+1:n)*x(g+1:n))/a(g,g) end 二迭代法 J迭代公式

a=[5 2 1;-1 4 2;2 -3 10] d=[-12;20;3]

x=[0;0;0]; stop=1.0e-4 L=-tril(a,-1) U=-triu(a,1) D=inv(diag(diag(a)))

X=D*(L+U)*x+D*d; n=1;

while norm(X-x,inf)>=stop x=X;

X=D*(L+U)*x+D*d; n=n+1; end X n

初始向量 迭代的精度 J迭代公式 . 可修编.

. -

G-S迭代公式

a=[5 2 1;-1 4 2;2 -3 10]

x=[0;0;0]; 初始向量 d=[-12;20;3] stop=1.0e-4 L=-tril(a,-1) U=-triu(a,1) D=(diag(diag(a)))

X=inv(D-L)*U*x+inv(D-L)*d; n=1;

while norm(X-x,inf)>= stop x=X;

X=inv(D-L)*U*x+inv(D-L)*d; n=n+1; end X n

SOR迭代公式

a=[5 2 1;-1 4 2;2 -3 10]

x=[0;0;0]; d=[-12;20;3] stop=1.0e-4

G-S迭代公式 时迭代中止否则继续 . 可修编.

初始向量 . -

w=1.44; 松弛因子 L=-tril(a,-1) U=-triu(a,1) D=(diag(diag(a)))

X=inv(D-w*L)*((1-w)*D+w*U)*x+w*inv(D-w*L)*d SOR迭代公式 n=1;

while norm(X-x,inf)>=stop 时迭代中止否则继续 x=X;

X=inv(D-w*L)*((1-w)*D+w*U)*x+w*inv(D-w*L)*d; n=n+1;end X n

3. a*x=d a=[5 2 1;-1 4 2;2 -3 10] d=[-12;20;3] (1)考察用J、G-S及SOR解此方程组的收敛性.

(2)分别用雅可比迭代法、高斯-赛德尔迭代法及逐次超松弛迭代法解此方程组。要求时迭代中止，观察收敛速度。

(1) J迭代矩阵的谱半径：max(abs(eig(D*(L+U))))= 0.506 G-S迭代矩阵的谱半径：max(abs(eig(inv(D-L)*U)))= 0.2000

SOR迭代矩阵的谱半径：max(abs(eig(inv(D-w*L)*((1-w)*D+w*U))))=0.9756 所以用J、G-S及SOR均收敛（且有）。 (2)

J： X =-4.0000 3.0000 2.0000 n =18

. 可修编.

. -

G-S： X =-4.0000 3.0000 2.0000 n =8 SOR()： X =-4.0000 3.0000 2.0000 n =482

分析：可见高斯-赛德尔迭代法要比雅可比迭代公式的收敛速度快。同时，如果超松弛迭代法的选取不当不但不会加速收敛还会对收敛速度造成很恶劣的结果。 四、结论

在科学研究和工程技术中有许多问题可归结为求解线性代数方程组的问题。本文主要讨论了直接法解线性方程组及迭代法解线性方程，讲解了各种直接法，如高斯消元法，三角分解法，追赶法等的基本思想和原理，介绍了病态方程和误差分析，了解了它们各自的优缺点及适用围。可以通过计算方法进行一种比较完善的构造用计算机进行计算，使之更普遍化，能够有举一反三的思想，利用直接法解线性方程组解决一些实际中难解的问题。

. 可修编.