精通MATLAB科学计算（第3版）(王正林)09-3r

第章线性方程组求解

在自然科学和工程技术中，很多问题可以归结为求解线性方程组。采用MATLAB,

不仅可以利用其提供的相关函数直接解决一些简单的线性方程组，而且可以通过简洁的编

程来解决一些复杂的线性方程组。

通过本章的介绍，读者既能应用MATLAB中相应的函数求解线性方程组，又能通过

编程，灵活使用迭代法和其他的特殊解法来求解线性方程组。

9.1求逆法

MATLAB中求解线性方程最直接的方法是矩阵求逆法，它适用于系数矩阵的数据无

规律且系数矩阵的阶数比较小的情况。如果系数矩阵的阶数太大的话，系数矩阵求逆就需

要花费很长时间。

在MATLAB中用求逆法求解，线性方程可以直接用左除法“\”求解，也可以用求逆函

数inv()求解。

【例9-1】左除法和求逆法求解线性方程组实例。用左除法和求逆法求解如下线性

方程组，

X\+2%2+3x3=1

<-x\+3工2+7工3=4

9xi+3x3=7

解：用左除法求解，在MATLAB命令窗口中输入求解程序：

»A=[l23;-137;903]

»b=[l47]z;

»x=A\b

输出计算结果为：

»X=

0.3333

-1.6667

1.3333

用逆矩阵法求解，在MATLAB命令窗口中输入求解程序：

第9章线性方程组求解

»A=[l23;-137;903];

»b=[l47]z;

>>x=inv(A)*b

输出计算结果为：

>>X=

0.3333

-1.6667

1.3333

由计算结果可知,方程组的解为[的42户3]=[03333,-1.6667,1.3333]。

9.2分解法

矩阵分解法是指根据一定的原理用某种算法将系数矩阵分解成若干个矩阵的代数运

算，常用的分解是乘积分解，而分解后的新矩阵一般是某种特殊矩阵。

常用的分解有LU分解、QR分解、Cholesky分解、Schur分解、Hessenberg分解和奇

异分解等。

9.2.1LU分解法

矩阵的LU分解也叫三角分解，还有的书称之为Doolittle分解，它是将矩阵分解为一

个单位下三角矩阵与上三角矩阵的乘积。只要矩阵非奇异，这种分解总是可以进行的。

MATLAB提供了lu函数来求矩阵A的LU分解，常用格式为：

也,S=lu(Z)：产生一个上三角矩阵U和一个变换形式的下三角矩阵Z(进行了行交

换)，使得A=LU；

[£,(/,P]=lu(/1)：产生一个上三角矩阵U和一个下三角矩阵乙及置换矩阵P,使得

PA=LU；

分解以后，方程组=〃的解可写为x=LA(ZA〃)或x=(A(ZAP*〃)。

【例9-21LU分解法求解线性方程组实例。用LU分解法求解以下线性方程组，

1.5xi+3x2-0.8^3+414=4

2xi+9%3+10%4=0

—lx\+4.8x2—0.6x3+=1

14%|+12.3%2—4M+5X4=—2°

解：用LU分解法求解，在MATLAB命令窗口中输入求解程序：

»A=[1.53-0.84;20910;-74.8-0.61;1412.3-45];

»b=[401-2]';

»[LzU]=lu(A);

1◄◄◄187

精通MATLAB科学计算(第2忡----------------------------

»x=U\(L\b)

输出计算结果为：

>>x=

-0.3721

-0.8291

-1.5336

1.4547

由计算结果可知,方程组的解为即X2,X3,X4]=[-0.3721,-0.8291,-1.5336,1.4547]。

9.2.2QR分解法

矩阵的QR分解就是把矩阵分解为一个正交矩阵和一个上三角矩阵的乘积。MATLAB

中对矩阵A进行QR分解的函数调用格式如下：

[Q,K]=qr(/)：产生一个正交矩阵。和一个上三角矩阵A,使得；

[Q,K,E]=qr(/)：产生一个正交矩阵Q、一个上三角矩阵R以及一个置换矩阵E,使得

AE=QR；

实现分解以后，方程组4r=6的解可写为x=K\(Q功)或x=E(K\(QS))。

【例9-31QR分解法求解线性方程组实例。用QR分解法求解以下线性方程组，

x\+0.5x2+0.3333%3+0.25%4=1

0.5%1+0.3333%2+0.25x3+0.2%4=2

0.3333x1+0.25x2+0.2%3+0.1667x4=2

0.25xi+0.2必+0.1667x3+0.1429x4=1。

解：用QR分解法求解，在MATLAB命令窗口中输入求解程序：

»A=[1.00000.50000.33330.2500;

0.50000.33330.25000.2000;

0.33330.25000.20000.1667;

0.25000.20000.16670.1429];

»b=[l221],;

»[Q,R]=qr(A)

»x=R\(Q\b)

输出。、A的计算结果为：

Q=

-0.83810.5226-0.1540-0.0263

-0.4191-0.44170.72780.3157

-0.2794-0.5288-0.1395-0.7892

-0.2095-0.5021-0.65360.5261

-1.1932-0.6705-0.4749-0.3698

0-0.1185-0.1257-0.1175

188►►►►

第9章线性方程组求解

00-0.0062-0.0096

0000.0002

输出解的计算结果为：

x=

1.0e+003*

0.1160

-1.4400

3.6000

-2.3800

由计算结果可知,方程组的解为[的广2,孙孙]=口的,-1440,3600,-2380]。

9.2.3Cholesky分解法

当系数矩阵N正定且对称时，它能分解为以下的形式:其中K为上三角矩阵，

肝为R的转置，是下三角矩阵，这种分解称为Cholesky分解。

MATLAB中与Cholesky分解有关的函数如下：

R=chol(Z):对《进行Cholesky分解，使得A=RXR;

[K,p]=chol(4):对A进行Cholesky分解，使得A=RrRa

如果4对称正定，则p=0；否则，p为一正整数。如果4满秩，/?是为2-1阶的上三

角矩阵，且有相/?=41:仍-1),l:(p-l))；实现分解后，方程组的解写成上大\(如财

的形式。

【例9-4]Cholesky分解法求解线性方程组实例。用Cholesky分解法求解以下线性

方程组，

9xi-36切+3013=1

*-36xi+192^2一180工3=1

30%1-180%2+180%3=10

解：用Cholesky分解法求解，在MATLAB命令窗口中输入求解程序：

»A=[9-3630;-36192-180;30-180180);

»b=ones(3,1);

>>R=chol(A)

>>x=R\(Rz\b)

输出/?的计算结果为：

R=

3.0000-12.000010.0000

06.9282-8.6603

002.2361

输出解的计算结果为：

◄◄◄◄189

精通MATLAB科学计算(第2版i----------------------------------------

X=

1.8333

1.0833

0.7833

由计算结果可知，方程组的解为［莺/2,必］=［1.8333,1.0833,0.7833］。

起合在MATLAB^中，当N正定但不对称时，用K=chol(N)也能得出一个矩阵，

MATLAB并不报错，但是只要验证一下Rt*R就会看出R「*R并不等于A,因

此用公式x=R\(*M)就会得出错误的结果。

190►►►>

第9章线性方程组求解

9.2.4其他分解法

MATLAB中还提供了矩阵的奇异值分解、Hessenberg分解和Schur分解等函数，可用

来求解线性方程组，如表9-1所示。

表9-1其他的矩阵分解函数

函数及常用语法说明

[USM=svd(j)奇异值分解，使得

[P.H]=hess(^)Hessenberg分解.使得A=P*H*/，其中P为酉矩阵

[U,T]=schur(/l)Schui■分解，使得力。，其中U为正交矩阵，T为Schur矩阵

【例9-5】奇异值分解法求解线性方程组实例。用奇异值分解法求解以下线性方程

组，

X)+0.5x2+0.3333%3+0.25x4+0.2x5=1

0.5%1+0.3333x2+0.25%3+0.2x4+0.1667x5=0

<0.3333汨+0.25x2+0.2啊+0.1667x4+0.1429x5=1

0.25xi+0.2x2+0.1667x3+0.1429x4+0.125x5=0

0.2xi+0.1667x2+0.1429x3+0.125x4+0.111IX5=1o

解：用奇异值分解法求解，在MATLAB命令窗口中输入求解程序：

»A=[1.00000.50000.33330.25000.2000;

0.50000.33330.25000.20000.1667;

0.33330.25000.20000.16670.1429;

0.25000.20000.16670.14290.1250;

0.20000.16670.14290.12500.1111]

>>[UzSrV]=svd(A)

z

>>b=[l0101]r

»x=V*inv(S)*UI*b

输出分解的I/、S、H的计算结果为：

U=

-0.76790.6019-0.21420.04720.0062

-0.4458-0.27590.7241-0.4327-0.1167

-0.3216-0.42490.12050.66740.5062

-0.2534-0.4439-0.30960.2330-0.7672

-0.2098-0.4290-0.5652-0.55760.3762

S=

1.56710000

00.2085000

000.011400

0000..00030

00000.,0000

V=

-0.76790.6019-0.21420.04720.0062

◄◄◄◄191

精通MATLAB科学计算（第2蚪

-0.4458-0.27590.7241-0.4327-0.1167

-0.3216-0.42490.12050.66740.5062

-0.2534-0.4439-0.30960.2330-0.7672

-0.2098-0.4290-0.5652-0.55760.3762

输出的计算结果为：

x=

1.0e+004*

0.0865

-1.2467

4.4659

-5.8404

2.5426

由计算结果可知，方程组的解为，

4

[x1,X2,x3,X4,x5]=10*[0.0865,-1.2467,4.4659,-5.8404,2.5426]。

选W•奇异值分解很有用，将系数矩阵进行奇异值分解以后，A的奇异值都在S

的对角线上，这样就可以粗略估计N的条件数，看看N是否是病态矩阵，

以此决定相应的求解方法。

•如果系数矩阵对称，那么奇异值分解后的。和/都是正交矩阵，而S是对

角矩阵，这样求解方程就十分方便。

【例9-6]Hessenberg分解法求解线性方程组实例。用Hessenberg分解法求解以下

线性方程组，

闲+工2+工3+工4=1

X]+2X2+3汹+414=4

X[+3X2+6x3+1014=7

X[+4X2+1013+20^4=-2o

解：用Hessenberg分解法求解，在MATLAB命令窗口中输入求解程序：

»A=[1111;

1234;

13610;

141020;];

»[P,H]=hess(A)

»b=[l47-2]1;

>>x=P*inv(H)*P1*b

输出分解的P、〃的计算结果为：

P=

0.7754-0.6247-0.09250

-0.6092-0.7015-0.36980

192►►►>

第9章线性方程组求解

0.16620.3431-0.92450

0001.,0000

0.2147-0.265400

-0.26541.08441.08020

01.08027.7009-10.8167

00-10.816720.0000

输出的计算结果为：

10.0000

-33.0000

36.0000

-12.0000

由计算结果可知,方程组的解为,[xi,x2,X3]=[10,-33,36,-12]o

起合•Hessenberg矩阵指的是第一子对角线以下的元素为0的矩阵。

•Hessenberg分解法分解后产生正交矩阵和Hessenberg矩阵，正交矩阵的逆

就是其转置；如果系数矩阵对称的话，产生的Hessenberg矩阵就是三对角

矩阵，这种矩阵的求逆速度也很快，因此用这种方法求解线性方程组也很快。

【例9-7】Schur分解法求解线性方程组实例。Schur分解法求解以下线性方程组，

汨+工2+工3+工4=1

X]+2X2+3冷+414=4

X1+3x2+6x3+10x4=7

X\+4^2+10工3+20%4=-2o

解：用Schur分解法求解，在MATLAB命令窗口中输入求解程序:

»A=[1111；

1234;

13610;

141020;]

»b=[l47-2]

»[UzT]=schur(A)

»x=U*inv(T)*Uf*b

输出分解的U、7的计算结果为：

U=

0.3087-0.78730.53040.0602

-0.72310.16320.64030.2012

0.59460.53210.39180.4581

-0.1684-0.2654-0.39390.8638

T=

0.0380000

00.453800

1◄◄◄193

精通MATLAB科学计算（第2版i---------------------------------------------

002.20340

00026.3047

输出的计算结果为：

x=

10.0000

-33.0000

36.0000

-12.0000

由计算结果可知,方程组的解为[XI,X2,X3,X4]=[10,-33,36,-12]。

通多・Schur矩阵是上三角矩阵，且其对角元素为被分解矩阵的特征值。

•Schur分解法在分解后产生正交矩阵和Schur矩阵，如果系数矩阵对称的话，

产生的Schur矩阵就是对角阵，显然对角阵的逆矩阵非常容易得到，因此

用这种方法求解线性方程组也很快。

9.3迭代法

在实际应用中（例如在运筹学、图论等领域中）,往往会出现这样一种情况：系数矩

阵阶数很高（例如ICT）,但系数矩阵含零元素相对较多，如果用前面的几种分解法去求解

的话，非零元素反而增多了，这就有点得不偿失了。

本节介绍另一类常用的求解方法——迭代法。迭代法是将求一组解转换为求一个近似

解序列的过程，并用最终的近似解来逼近真实解。迭代法需要考虑以下3个重要的问题。

（1）迭代的初始值

考虑初始值的选取是否有范围限制，不同的初始值对最终的迭代结果是否有影响。

（2）迭代算法

考虑怎样由当前的迭代结果得出下次迭代的初始值；由于不同的算法会带来不同的误

差，所以应考虑算法使误差放大还是减少。

（3）迭代的收敛性

考虑迭代是否收敛，收敛的过程是快还是慢。

以上问题的具体分析请参考数值分析的教材，下面讲述采用MATLAB编程实现的几

种常见迭代法。

9.3.1逐次逼近法

对于n阶线性方程组，Ax=b的系数矩阵N（假设N为非奇异）进行如下分解

A=Q-C

194►►►►

第9章线性方程组求解

其中。为非奇异。

则方程组可变换为

x=Bx+r

其中B=Q'C,r=Q：'b,那么迭代过程可以写成

xk+l=Bxk+r

这种迭代过程称为逐次逼近法。取不同的。和C,就得到不同的迭代算法。通常情况

下，逐次逼近法收敛的充分条件是迭代矩阵的谱半径小于lo

9.3.2里查森迭代法

里查森迭代法可以说是最简单的迭代法了，它采用的迭代公式如下：

xk+i=(I-A)xk+b

在MATLAB中编程实现的里查森迭代法函数为：richasono

功能：用里查森迭代法求线性方程组的解。

调用格式：[x,n]=richason(^,b,xO,eps,M)

其中，A为线性方程组的系数矩阵；

分为线性方程组中的常数向量；

xO为迭代初始向量；

eps为解的精度控制(此参数可选)；

M为迭代步数控制(此参数可选)；

x为线性方程组的解；

〃为求出所需精度的解实际的迭代步数。

理查森迭代法的MATLAB程序代码如下：

function[x,n]=richason(A,b,xO,eps,M)

%采用里查森迭代法求线性方程组Ax=b的解

W线性方程组的系数矩阵：A

告线性方程组中的常数向量：b

%迭代初始向量：xO

3解的精度控制：eps

为迭代步数控制：M

年线性方程组的解：x

若求出所需精度的解实际的迭代步数：n

if(nargin==3)

eps=1.Oe-6;%eps表示迭代精度

1◄◄◄195

精通MATLAB科学计算(第2忡-----------------------------

M=200;*M表示迭代步数的限制值

elseif(nargin==4)

M=200;

end

I=eye(size(A));

xl=x0;

x=(I-A)*x0+b;

n=l;

告迭代过程

while(norm(x-xl)>eps)

xl=x;

x=(I-A)*xl+b;

n=n+1;会n为最终求出解时的迭代步数

if(n>=M)

disp(•Warning:迭代次数太多，可能不收敛！,);

return;

end

end

【例9・8】理查森迭代法求解线性方程组实例。用理查森迭代法求解以下线性方程

组，其中取初始值为［0,0,0］。

1.0170汨一0.0092M-0.0095x3=1

,-0.0092X］+0.9903x2+0.0136x3=0

-0.0095xi+0.0136x2+0.9898x3=10

解：用理查森迭代法求解，在MATLAB命令窗口中输入求解程序：

»A=[1.0170-0.00920.0095;

-0.00920.99030.0136;

0.00950.01360.9898];

»b=[l01—;

»x0=[000]/_;

>>[x,n]=richason(A,b,xO)

输出的计算结果为：

X=

0.9738

-0.0047

1.0010

输出的迭代次数为：

n=g

身

可见，经过5步迭代，理查森迭代法求出了方程组的解为：

196►►►>

第9章线性方程组求解

[阳,xz,X3]=[0.9738；6.0047,1.0010]l

[巧，巧，演]=[0.9738,-0.0047,1.001]。

对上述迭代计算结果进行验证，在MATLAB命令窗口中输入如下程序：

»A*x

输出结果为：

ans=

1.0000

0.0000

1.0000

经过验证，也确定了解的正确性。一般情况下，理查森迭代法难以收敛，因为要使/-/

的谱半径小于1,通常要求力是严格对角占优的矩阵。

◄◄◄◄197

精通MATLAB科学计算(第2版i----------------------------------------

9.3.3Jacobi迭代法

如果系数矩阵”的主对角元全不为0,在上节4的分解中取

Q=D

C=D-A

其中O是由/的主对角元素组成的对角阵，则有B=I-D'A,r=D'b,迭代公式为：

x^=(I-D'A]xk+D'b

这种迭代方法称为Jacobi迭代法。

在MATLAB中编程实现的Jacobi迭代法函数为：jacobio

功能：用Jacobi迭代法求线性方程组4r的解

调用格式：[x,w]=jacobi(A,h,xO,eps,varargin)

其中，A为线性方程组的系数矩阵；

8为线性方程组中的常数向量；

xO为迭代初始向量；

eps为解的精度控制(此参数可选)；

varargin为迭代步数控制(此参数可选)；

x为线性方程组的解；

«为求出所需精度的解实际的迭代步数。

Jacobi迭代法的MATLAB程序代码如下：

function[x,n]=jacobi(A,b,xO,eps,varargin)

为采用Jacobi迭代法求线性方程组Ax=b的解

电线性方程组的系数矩阵：A

年线性方程组中的常数向量：b

卷迭代初始向量：xO

优解的精度控制：eps

带迭代步数控制：varargin

为线性方程组的解：x

生求出所需精度的解实际的迭代步数：n

ifnargin==3

eps=1.Oe-6;

M=200;

elseifnargin<3

error

return

elseifnargin==5

198►►►>

第9章线性方程组求解

M=varargin{1};

end

D=diag(diag(A));务求A的对角矩阵

L=-tril(Az-1);%求A的下三角阵

U=-triu(A,1);彩求A的上三角矩阵

B=D\(L+U);

f=D\b;

x=B*xO+f;

n=l;%迭代次数

多迭代过程

whilenorm(x-xO)>=eps

xO=x;

x=B*xO+f;

n=n+l;

if(n>=M)

disp(，Warning:迭代次数太多，可能不收敛！,)；

return;

end

end

【例9-9]Jacobi迭代法求解线性方程组实例。用Jacobi迭代法求解以下线性方程

组，其中取初始值为[1,1,1]。

0.9889xj-0.0005x2-0.0002x3=1

<-0.0046x1+0.9946x2+0.0077x3=0

-0.0002x14-0.0092x2+0.9941X3=1。

解：用Jacobi迭代法求解，在MATLAB命令窗口中输入求解程序：

»A=[0.9889-0.0005-0.0002;

-0.00460.99460.0077;

-0.00020.00920.9941];

»b=[l01]匕

>>x0=ones(3Z1);

>>[xzn]=jacobi(A,b,xO)

输出的计算结果为：

X=

1.0114

-0.0031

1.0062

输出的迭代次数为：

n=£

4

Y◄◄◄199

精通MATLAB科学计算（第2版i---------------------------------------

可见，经过4步迭代，Jacobi法求出了方程组的解为：

[XI,X2,X3]=[1,0114,-0.0031,1.0062]。

•Jacobi迭代法对任意初始向量M）收敛的充要条件为B的谱半径小于1,即8

»的所有特征值的绝对值都小于lo

200►►►►

第9章线性方程组求解

9.3.4Gauss-Seidel迭代法

仍然设系数矩阵N的主对角元全不为零，如果对N作以下分解：

A=(D-L)-U

其中。的意义同Jacobi迭代法，，为下三角矩阵，。为上三角矩阵，它的迭代公式为:

xk+i=(D-£)-'Uxk+(D-L)-'b

在MATLAB中编程实现的Gauss-Seidel迭代法函数为：gauseidelo

功能：用Gauss-Seidel迭代法求线性方程组的解。

调用格式：gauseidel(44的,eps,防。

其中，A为线性方程组的系数矩阵；

b为线性方程组中的常数向量；

xO为迭代初始向量；

eps为解的精度控制(此参数可选)；

M为迭代步数控制(此参数可选)；

x为线性方程组的解；

〃为求出所需精度的解实际的迭代步数。

Gauss-Seidel的MATLAB程序代码如下：

function[x,n]=gauseidel(A,b,xO,eps,M)

会采用Gauss-Seidel迭代法求线性方程组Ax=b的解

年线性方程组的系数矩阵：A

%线性方程组中的常数向量：b

他迭代初始向量：xO

软解的精度控制：eps

多迭代步数控制：M

会线性方程组的解：x

务求出所需精度的解实际的迭代步数：n

ifnargin==3

eps=1.Oe-6;

M=200;

elseifnargin==4

M=200;

elseifnargin<3

error

return;

end

D=diag(diag(A));%求A的对角矩阵

◄◄◄◄201

精通MATLAB科学计算（第2回

L=-tril(A,-l);务求A的下三角矩阵

U=-triu(A,1);%求A的上三角矩阵

G=(D-L)\U；

f=(D-L)\b；

x=G*xO+f;

n=l;它迭代次数

国迭代过程

whilenorm(x-xO)>=eps

xO=x;

x=G*xO+f;

n=n+l;

if(n>=M)

disp('Warning:迭代次数太多，可能不收敛！，)；

return;

end

end

【例9-10】GaussSeidel迭代法求解线性方程组实例。用Gauss・Seidel迭代法求解

以下线性方程组，其中取初始值为[1,1,1]。

1.4449xi+0.7948x2+0.880卜3=1

,0.6946XI+1.9568x2+0.1730x3=0

0.6213为+0.5226x2+1.9797x3=1o

解：用Gauss-Seidel迭代法求解，在MATLAB命令窗口中输入求解程序：

»A=[1.44490.79480.8801;

0.69461.95680.1730;

0.62130.52261.9797];

»b=[l01]已；

»x0=zeros(3,1);

>>[x,n]=gauseidel(A,b,xO)

输出的计算结果为：

X=

0.5929

-0.2443

0.3835

输出的迭代次数为：

n=n.

a

可见，经过11步迭代，Gauss-Seidel迭代法求出方程组的解为：

[Xi,x2,x3]=[0.5929,-0.2443,0.3835]。

9.3.5超松弛迭代法

202►►►►

---------------------------------------------------------第9章线性方程组求解

如果对系数矩阵”作以下分解：

A=(D-a)L)-((l-a))D+a)U)

其中。、L、。的意义与Gauss-Seidel迭代法相同。

。是一个事先选好的常数，称为松弛因子，当时叫超松弛法，也叫SOR迭代法；

当。＜1时叫低松弛法。其迭代公式为：

x*+1=(£＞-尸[(1-co)D+a)U]Xk+(o(D-(oL)~'b

关于超松弛迭代法的收敛性有如下结论：

若系数矩阵/对称正定，当0＜。＜2,SOR迭代法收敛。

在MATLAB中编程实现的超松弛迭代法函数为：SORo

功能：用超松弛迭代法求线性方程组的解。

调用格式：[x,n]=SOR(A,b,xO,w,eps,M)

其中，A为线性方程组的系数矩阵；

6为线性方程组中的常数向量；

xO为迭代初始向量；

w为松弛因子；

eps为解的精度控制(此参数可选)；

M为迭代步数控制(此参数可选)；

x为线性方程组的解；

n为求出所需精度的解实际的迭代步数。

超松弛迭代法的MATLAB程序代码如下：

function[x,n]=SOR(A,b,xO,w,eps,M)

带采用超松弛迭代法求线性方程组Ax=b的解

传线性方程组的系数矩阵：A

省线性方程组中的常数向量：b

与迭代初始向量：xO

%松弛因子：w

为解的精度控制：eps

超迭代步数控制：M

%线性方程组的解：x

省求出所需精度的解实际的迭代步数：n

ifnargin==4

eps=1.Oe-6;

M=200;

elseifnargin＜4

error

◄◄◄◄203

精通MATLAB科学计算(第2版i-----------------------------------------

return

elseifnargin==5

M=200;

end

if(w<=0||w>=2)与收敛条件要求

error;

return;

end

D=diag(diag(A));%求A的对角矩阵

L=-tril(Az-1);考求A的下三角矩阵

U=-triu(Azl);%求A的上三角矩阵

B=inv(D-L*w)*((1-w)*D+w*U);

f=w*inv((D-L*w))*b;

x=B*xO+f;

n=l;恭迭代次数

带迭代过程

whilenorm(x-xO)>=eps

x0=x;

x=B*xO+f;

n=n+l;

if(n>=M)

disp(1Warning:迭代次数太多，可能不收敛！，)；

return;

end

end

【例9-11]超松弛迭代法求解线性方程组实例。用超松弛迭代法求解以下线性方程

组，其中取初始值为[0,0,0]。

4xi+3x2=24

<3xi+4x2一%3=30

—X2+4%3=-24

解：用超松弛迭代法求解，在MATLAB命令窗口中输入求解程序：

»A=[430;34-1;0-14];

»b=[2430-24]';

»x0=[000]';

>>[xzn]=SOR(Azb,xOz1.25)

输出的计算结果为：

X=

3.0000

4.0000

-5.0000

输出的迭代次数为：

204►►►>

---------------------------------------------------------第9章线性方程组求解

n=14

可见，经过14步迭代，SOR迭代法求出了方程组的解为[R,X2,X3]=[3,4,-5]。

超松弛迭代法还有一种改进形式，叫做对称逐次超松弛迭代法(SSOR),它采用的是两

步迭代公式：

(D-a)L)x*+i/2=0(fJxk+b)+(1-a))Dxi；

(D-coU)Xk+i=(o(Lxk+b)+(1-<v)Z)x*.+i/2

在MATLAB中编程实现的对称逐次超松弛迭代法函数为：SSORo

功能：对称逐次超松弛迭代法求线性方程组的解。

调用格式：[x,〃]=SSOR(4b,xO,iv,eps,M)o

其中，A为线性方程组的系数矩阵；

b为线性方程组中的常数向量；

xO为迭代初始向量；

w为松弛因子；

eps为解的精度控制(此参数可选)；

M为迭代步数控制(此参数可选)；

x为线性方程组的解；

〃为求出所需精度的解实际的迭代步数。

SSOR算法用MATLAB实现如下：

function[xzn]=SSOR(A,b,xO,w,eps,M)

务采用对称逐次超松弛迭代法求线性方程组Ax=b的解

金线性方程组的系数矩阵：A

告线性方程组中的常数向量：b

超迭代初始向量：xO

%松弛因子：w

三解的精度控制：eps

告迭代步数控制：M

与线性方程组的解：x

%求出所需精度的解实际的迭代步数：n

ifnargin==4

eps=1.Oe-6;

M=200;

elseifnargin<4

error

return

elseifnargin==5

M=200;

◄◄◄◄205

精通MATLAB科学计算（第2蝌

end

if(w<=0|Iw>=2)

error;

return;

end

D=diag(diag(A));务求A的对角矩阵

L=-tril(A,-l);%求A的下三角矩阵

U=-triu(A,1);£求A的上三角矩阵

Bl=inv(D-L*w)*((1-w)*D+w*U);楮第一步的迭代矩阵

B2=inv(D-U*w)*((1-w)*D+w*L);带第二步的迭代矩阵

fl=w*inv((D-L*w))*b;

f2=w*inv((D-U*w))*b;

xl2=Bl*xO+fl;3第一步的迭代得到的中间值

x=B2*xl2+f2;

n=l;称迭代次数

%迭代过程

whilenorm(x-xO)>=eps

xO=x;

xl2=Bl*xO+fl;

x=B2*xl2+f2;

n=n+l;

if(n>=M)

disp('Warning:迭代次数太多，可能不收敛！D;

return;

end

end

【例9-12］对称逐次超松弛迭代法求解线性方程组实例。用SSOR迭代法求解以下

线性方程组，其中取初始值为［0Q0］。

4为+3%2=24

<3xi+4必一M=30

-X2+4%3=-24

解：用SSOR迭代法求解，在MATLAB命令窗口中输入求解程序：

»A=[430;34-1;0-14];

»b=[2430-24]';

»x0=[000]';

»[x,n]=SSOR(A,b,xO,1.25)鼻松弛因子取1.25

输出的计算结果为：

x=

3.0000

4.0000

-5.0000

206►►►>

第9章线性方程组求解

输出的迭代次数为：

n=34

%

可见,经过34步迭代,SSOR迭代法求出方程组的解为[M,X2,X3]=[3,4,-5]。

上面的两个例子采用的是同一个线性方程组，且采用的松弛系数都是1.25,从求解所

需的步骤来看,SOR迭代法（14步迭代）比SSOR迭代法（34步迭代）要快一些。

9.3.6两步迭代法

两步迭代法的迭代公式如下：

(D-L)Xk+MkUxk+b

(D-U)xk+1=Lxk+i.2+b

其中D、4。的定义和前面一样，即D