第6章 (新)解线性方程组的迭代法(演示)

第6章 解线性方程组的迭代法

直接方法比较适用于中小型方程组。对高阶方程组，即使系数矩阵是稀疏的，但在运算中很难保持稀疏性，因而有存储量大，程序复杂等不足。迭代法则能保持矩阵的稀疏性，具有计算简单，编制程序容易的优点，并在许多情况下收敛较快。故能有效地解一些高阶方程组。

1 迭代法概述 迭代法的基本思想是构造一串收敛到解的序列，即建立一种从已有近似解计算新的近似解的规则。由不同的计算规则得到不同的迭代法。迭代法的一般格式

(k1)(k)(k1)(km)xFk(x,x,,x),k0,1,

式中x(k1)与x(k),x(k1),,x(km)有关，称为多步迭

(k1)(k)xx代法。若只与有关，即 x(k1)Fk(x(k)),k0,1,

称为单步迭代法。再设Fk是线性的，即

(k1)(k)xBkxfk,k0,1,

式中BkRnn，称为单步线性迭代法。Bk称为迭

1

代矩阵。若Bk和fk与k无关，即

x(k1)Bx(k)f,k0,1,

称为单步定常线性迭代法。本章主要讨论具有这种形式的各种迭代方法。

1.1 向量序列与矩阵序列的极限

由于R中的向量可与R的点建立——对应关系，由点列的收敛概念及向量范数的等价性，可得到向量序列的收敛概念。

定义1 设{x(k)}为Rn中的向量序列，xRn，如果

limx(k)x0

nkn其中为向量范数，则称序列{x(k)}收敛于x，

x(k)x。 记为limk定理1 Rn中的向量序列{x(k)}收敛于Rn中的向量x当且仅当

limxi(k)xi (i1,2,,n) k(k)(k)T其中x(k)(x1(k),x2,,xn),x(x1,x2,,xn)T。 证明 由定义1，{x(k)}收敛于x，即

limx(k)x0 k而对任意1in，有

2

0xi(k)ximaxx(jk)xjx(k)x

1jn由极限存在准则得

limxk(k)ixi0 (i1,2,,n)

即limxi(k)xi (i1,2,,n)。 k定义2 设{A}为Rnn中的矩阵序列，ARn，如果

limA(k)A0

(k)k其中为矩阵范数，则称序列{A(k)}收敛于A，记为limA(k)A。

k定理 2 R中的矩阵序列{A}收敛于R中的矩阵A的充要条件为

(k)limaijaij (i,j1,2,,n) k证明留给读者。

定理1和定理2表明，向量和矩阵序列的收敛可以归结为对应分量或对应元素序列的收敛。

A(k)0的充分必要条件是 定理3 limklimA(k)x0，xRn

knn(k)nn其中两个极限的右端分别指零矩阵和零向量。 证明 对任一种算子范数，有A(k)xA(k)x，

3

从而可证必要性。若x依次取n个单位向量

其它ej (j1,2,n,，)其中ej的第j个分量为1，

分量为零。得

A(k)A(k)I(A(k)e1,A(k)e2,,A(k)en)

所以

(k)(k)(k)(k)limA(limAe1,limAe2,,limAen)0 kkkk充分性得证。

定理4 设BRnn，则下面三个命题等价 (1) limBk0；(2) (B)1；(3) 至少存在一k种从属的矩阵范数，使B1。

证明 (1)(2)：用反证法，假设B有一个特征值，满足1，则有特征向量x0，满足

kkBxx。由此可得Bxx。所以当k时向量序列Bkx不收敛于零向量。由定理3可知(1)不成立，从而知1，即(2)成立。

(2)(3)：根据第五章定理18，对任意实数0，存在一种从属的矩阵范数，使B(B)。由(2)有(B)1，适当选择，便可使B1，即(3)成立。

(3)(1)：由(3)给出的矩阵范数B1，由

4

BkB，可得limBkkk 0，从而有limBk0。

k定理 5 设BRnn，为任一种矩阵范数，则limBk1kk(B)。

1kk 证明 因为(B)(B)由第5章定理18，，

对一切k，有

k(B)(B)B1k1kk

另一方面，对任意的0，记

B(B)B

1显然有(B)1。由定理4可知limB0，所

kk以存在NN()，使kN时

BkBk(B)1kkk1

即当kN时有

(B)B(B)

由任意性即得定理结论。

1.2 迭代公式的构造 将非奇异线性方程组

Ax=b

变形为等价方程组

5

x=Bx+f (1.1)

由此构造迭代公式

x(k1)Bx(k)f (k0,1,2,) (1.2)

给定初始向量xR后，按此迭代公式产生向量序列x(k)，当k充分大时，以x(k)作为方程组

(0)n的近似解，这就是求解线性方程组的单步定常线性迭代法。B称为迭代矩阵。

定义3 如果对任意的初始向量x(0)及f，迭代法(1.2)得出的向量序列x(k)都有

limx(k)x* k成立，其中x*为一确定的向量，它不依赖于x(0)的选取，则称迭代法(1.2)是收敛的，否则称迭代法(1.2)是发散的。

(k)x显然，若按式(1.2)产生的向量序列收敛于向量x*，则有

x*limx(k1)lim[Bx(k)f]Bx*f

kk即x*是方程组(1.1)的解。

2 基本迭代法 2.1 Jacobi（雅可比）迭代法 考虑方程组Axb，即

6

a11x1a12x2a1nxnb1axaxaxb2112222nn2 (2.1)  an1x1an2x2annxnbn其中A(aij)nn非奇异，故不妨设aii0

(i1,,n)。(2.1)等价变形为

aiixibiaijxjj1i1ji1ax (i1,,n)

ijjn有

i1n1xibiaijxjaijxj (i1,,n) (2.2)

aiij1ji1由此构造迭代公式

xi(k1)i1n1(biaijx(jk)aijx(jk))(i1,,n) (2.3) aiij1ji1记

00a12a1na021,UL0an1n0 an1ann10a11Dann

7

则ADLU。由式(2.1)到式(2.2)的过程用矩阵形式表示为

Axb(DLU)xb

Dx(LU)xbxD1(LU)xD1b 因此式(2.3)的矩阵形式为

x(k1)Jx(k)f (2.4)

其中JD1(LU)ID1A,fD1b。

式(2.3)为迭代法的分量形式，它可用于计算迭代近似解；式(2.4)为迭代法的矩阵形式，它主要用于验证迭代法是否收敛及定性分析。 算法6.1

1．输入A(aij),b(b1,,bn)，维数n，精度,

(0)x(0)(x1(0),,xn)，最大容许迭代次数N。 2．置k1

3．对i1,2,,n

xi(biaijx)/aii

(0)jj1jin4．若xx(0)，输出x，停机，否则转5。

5．若kN，置k1k,xixi(0)(i1,2,,n)，

转3；否则，输出失败信息，停机。

2.2 Gauss-Seidel（高斯-赛德尔）迭代法

8

在Jacobi迭代法中，是用x(k)的全部分量来计算x(k1)的全部分量的，然而在计算分量xi(k1)时，

(k1)k1)x1(k1),x2,,xi(1都已经算出，如果Jacobi迭代法收敛，试想用多迭代一次的x1(k1),x2(k1),,xi(k11)(k)k)(k1)代替x1(k),x2来计算，可望取得更好的,,xi(x1i结果。这就是Gauss-Seidel迭代法的基本思想。其迭代公式为

xi(k1)i1n1(biaijx(jk1)aijx(jk)) (2.5) aiij1ji1式(2.5)的矩阵形式为

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

因此迭代法的矩阵形式为

x(k1)Gx(k)f (2.6)

其中G(DL)1U,f(DL)1b。

算法6.2

1．输入A(aij),b(b1,,bn)，维数n，

(0)(0)x(0)(x1(0),x2,,xn),，最大容许迭代次数N。 2．置k1 3．计算

x1(b1a1jx(0)j)/a11j2nxi(biaijxjj1

i1ji1(0)axijj)/aii (i2,,n1)9

n

xn(bnanjxj)/ann

j1n14．若xx(0)，输出x，停机，否则转5。 5．若kN，置k1k,xixi(0)(i1,2,,n)，转3；否则，输出失败信息，停机。

2.3 SOR迭代法

解线性方程组的超松弛法，也叫SOR法，是目前求解大型方程组的一种最常用的方法。是Gauss-Seidal迭代法的一种加速方法。 对一个收敛的Gauss-Seidel迭代法，第k1次的迭代结果一般要比第k次的好。第k1次的迭代结果可看作第k次基础上的修正，现在我们引入一个参数，来改变这个修正量。这就是SOR方法的基本思想。

记x(x1,x2,,xn)Tx(k1)x(k)，其中x(k1)由Gauss-Seidel迭代公式算得。于是

xixi(k1)xi(k)i1n 1(k1)(k)(k)(biaijxjaijxj)xi(i1,,n)aiij1ji1可以把x看作Gauss-Seidel迭代的修正项，即

第k次近似解x(k)以此项修正后得到近似解

10

x(k1)x(k)x。松弛法是将x乘上一个参数因

子作为修正项而得到新的近似解，其具体公式为

x(k1)x(k)x

即

xi(k1)xi(k)xii11xi(k)(biaijx(jk1)j1aiiaijx)x ji1n(k)j(k)iji1(k)axijj)n(1)xi(k)aii(biaijx(jk1)j1i1 (i1,,n) (2.7)

式(2.7)也可看成是Gauss-Seidel迭代法第k次和

按式(2.7)计算方程组k1次迭代解的加权平均。

的近似解序列的方法称为松弛法，称为松弛因子。当1为低松弛，1是Gauss-Seidel迭代，当1时称为超松弛法，简称SOR法。

式(2.7)的矩阵形式为

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

即

1(k)1(k1)1x(k1)(1)IDUxDLxDb 注意到ID1L1，有

(2.8)

其中L(DL)1(1)DU,f(DL)1b。

x(k1)Lx(k)f11

算法6.3

1．输入A(aij),b(b1,,bn)，维数n，精度

(0)要求,x(0)(x1(0),,xn最大容许迭代次数N。 )，2．置k1 3．计算

x1(1)x(0)1(b1a1jx)/a11(0)jj2i1(0)axijj)/aiinnxi(1)xi(0)(biaijxjj1n1j1ji1

(0)xn(1)xn(bnanjxj)/ann (i2,,n1)4．若xx(0) ，输出x，停机，否则转5。

(0)i5．若kN，置k1k,xix(i1,2,,)n，转3；否则，输出失败信息，停机。

例1 用Jacobi和Gauss-Seidel迭代法求解线性方程组

10x1x22x372x110x22x383 xx5x42312解 Jacobi迭代法的迭代公式为(算法语言)

12

取x(0)(k1)1(k)(k)x(72x2x23)110(k1)1(k)(k)x(83x2x213) 10(k1)1(k)(k)x3(42x1x2)5(0,0,0)T，迭代9次的近似解

(9)Tx(10.9994,11.9994,12.9992)

容易验证，方程组的精确解为x*(11,12,13)T。 G-S法的迭代公式为

(k1)1(k)(k)x(72x2x23)110(k1)1(k1)(k) x(83x2x213)10(k1)1(k1)(k1)x(42xx)3125迭代6次的近似解

x(6)(10.9999,11.9999,13.0000)T。

例2 取1.4,x(0)(1,1,1)T，用超松弛法求解方

程组

2x1x2 1x12x2x30  x2x1.823

13

解 迭代公式为

(k)x1(k1)0.4x1(k)0.7(1x2)(k1)(k)(k1)(k)x0.4x0.7(xx2213) (k1)(k)(k1)x30.4x30.7(1.8x2) (k0,1,2,)迭代9次后的近似解

x(9)(1.2000,1.3996,1.6001)T

与精确解x*(1.2,1.4,1.6)T比较，误差为

x(9)

13x0.000410

2*3 迭代法的收敛性

3.1 一阶定常迭代法的基本定理

定理 6 对任意的初始向量x(0)和右端项f，由迭代格式

x(k1)Bx(k)f (k0,1,2,) (3.1)

产生的向量序列{x(k)}收敛的充要条件(B)1。 证明 必要性。设存在向量x*，使得(k)*limxx，则 kx*Bx*f

由此及迭代公式(3.1)，有

x(k)x*Bk(x(0)x*)

14

于是

limBk(x(0)x*)lim(x(k)x*)0

kk因为x(0)为任意n维向量，因此上式成立必须

limBk0

k由定理 4，即得(B)1。

充分性。若(B)1，则1不是B的特征值，因而有I-B0，于是对任意n维向量f，方程组(IB)xf有唯一解，记为x*，即

x*Bx*f

Bk0，又因为 并且limkx(k)x*B(x(k1)x*)Bk(x(0)x*)

所以，对任意初始向量x(0)，都有

lim(xk(k)x)limB(xk*k(0)x)0

*即迭代公式(3.1)产生的向量序列{x(k)}收敛。 推论1 对任一种矩阵范数，若B1，则(k){x}收敛。

推论2 SOR法收敛的必要条件是02。 证明 设L有特征值1,2,,n。因为

ndet(L)1,2,,n[(L)]

由定理 6，SOR法收敛必有det(L)1，又因

15

L(DL)1(1)DU 1(1)na11a22ann(1)na11a22ann于是有

det(L)(1)n1

所以02。

定理6表明，迭代法收敛与否只决定于迭代矩阵的谱半径，与初始向量以方程组的右端项无关。

例3 研究用Jacobi和Gauss-Seidel迭代法解方程组Axb的敛散性。

12213434111； (2) A34134 (1) A22134341解 (1) Jacobi法的迭代矩阵为

JD1(LU)

其特征方程为

1ID(LU)0D(LU)0 由已知

D(LU)0

得1330，所以(J)01，因此Jacobi迭代法收敛。

16

3

如果用Gauss-Seidel法，迭代矩阵为

G(DL)1U

特征方程

I(DL)1U0(DL)U0 由已知

(DL)U(44)0

2得

10,32(12),32(12)

所以(G)2(12)1，因此Gauss-Seidel迭代

法发散。

(2) J法的特征方程为

2727D(LU)0

163227273令f()，因为

163249121f(1)0,f(2)0

32323所以f()在[2,1]中有一根，于是 (J)maxi1

1i3因此Jacobi迭代法发散。 G法的特征方程为

8127(DL)U()0

2

17

得特征根10,2,30.63330.146i，所以

(G)maxi2,30.6501

1i3因此G-S迭代法收敛。 注：1）在求G-S法和SOR法的特征方程时，注意到事实

IGI(DL)1U0

(DL)U01ILI(DL)(1)DU0

(DL)(1)DU0 可避免求逆矩阵。

2）J法和G法的收敛性没有必然联系。 3）定理6虽然给出了判别迭代法收敛的充要条件，但实际使用时很不方便。因为求谱半径(B)却是比较困难的事，因此常利用(B)B，作为(B)上界的一种估计式。需要指出的是：矩阵范数可以是任何一种矩阵范数，不限于算子范数，常用的范数有1,2,, F；其次，使用矩阵范数判别只是充分条件，而非必要条件。

3.2特殊方程组迭代法的收敛性

18

定义 4 若A(aij)nn满足

aiiaij (i1,,n)

j1jin且至少有一个i，使上式中不等式严格成立，则称A为弱对角占优。若所有的i不等式均成立，则称A为严格对角占优。

定义5 设A(aij)nn，若存在置换矩阵P使

A11A12TPAP (3.2) 0A22其中A11,A22均为方阵，则称为可约矩阵，否A．．则，称A为不可约。

Definition5 An nn matrix A(ajk) is called

reducible if there exist two nonempty N,M{1,,n} such that

NM,NM{1,,n}

and

ajk0, jN,kM Otherwise the matrix is called irreducible.

sets

A为可约矩阵意即A经过若干行列重排可化(3.2)式或Axb可化为两个低阶线性方程组求解（如果A经过两行交换的同时进行相应两列的交换，称对A进行一次行列重排）。

19

事实上，由Axb可化为

PTAP(PTx)PTb y1Td1T且记yPx,Pb，其中y1,d1为r维向

y2d2量。于是，求解Axb化为求解

A11y1A12y2d1 A22y2d2由上式第2个方程组求出y2，再代人第1个方

程组求出y1。

显然，如果A所有元素都非零，则A为不可约矩阵。

例 设有矩阵

b1c1abc222，ai,bi,ci都不为零 Aabcn1n1n1bn40B20131120401411021401,C

1041150114则A,C为不可约矩阵，B为可约矩阵。因为将B

的第二行和第三行对换，同时也把第二列和第

20

三列对换，就可化为

42002400113111 25定理 7 若A(aij)nn严格对角占优或为不可约弱对角占优矩阵，则aii0,i1,,n，且A非奇异。

证明 从严格对角占优可知aii0,i1,,n。

T若A奇异，则有x(x1,,xn)0满足Ax0。设xkx，则方程组的第k个方程为

akkxkakjxj

j1jkn由此得

akkakjj1jknxjxkakj

j1jkn与A为严格对角占优矛盾。 我们有如下结论： 1．若A为严格对角占优阵或不可约弱对角占优，则Jacobi迭代法和Gauss-Seidel均收敛。 2．若A为严格对角占优阵或不可约弱对角占

21

优，01，则SOR法收敛。

3．设A对称，且对角元aii0 (i1,,n)，则Jacobi迭代法收敛的充要条件是A及2DA均正定，其中Ddiag(a11,a22,,ann)。

4．若A对称正定，则SOR收敛的充要条件是02。

（4．的特殊情形）若A对称正定，则GS4．

迭代法收敛。

例4 考虑系数矩阵

1012210A11102,A2121 115012的方程组。显然，A1为严格对角占优，故J和

G-S法均收敛。A2非严格对角占优，但为对称正定矩阵，故取(0,2)，SOR法收敛。 例5 若

310 A943015,(G)，计算可得(J)所以J和G-S法22均发散。

如果改变方程的次序，有

22

94 A310显然A为严格对角占优，故J和G-S法均收敛。表明，改变方程组中方程的次序，即将系数矩阵作行交换会改变迭代法的收敛性。

定理 8 若A为严格对角占优阵或不可约弱对角占优，则J法和GS法均收敛。

证明 我们只给出严格对角占优情形的证明。只需证明(J)1和(G)1。 J法的迭代阵

0a211JIDAa11an1anna12a110an2anna1na11a2na22 0显然有J1，故(J)1，从而J法收敛。

GS法的特征方程为

I(DL)1U0(DL)U0 令C(DL)U，有C0。

现在证明所有特征值1。采用反证法，若

23

有一1，则由A为严格对角占优阵有

aiiaj1jinijaj1i1ijji1naij (i1,,n)

因此C为严格对角占优阵，故C0，矛盾，因此只能1，即(G)1，从而GS法收敛。 定理9 若A为严格对角占优阵或不可约弱对角占优，01，则SOR法收敛。

证明 我们只就A为严格对角占优阵的情形的给出证明。SOR法的特征方程为

1detI(DL)((1)DU)0det(DL)(1)DU0

即

det(1)DLU0

设是上述方程的任一根，我们只需证明采用反证法，假设1，由已知01，1。

则10，于是

det(DLU)0

11LU，有det(C)0。 令CD11由于

24

1

1(1)(1)所以C也为严格对角占优矩阵，故det(C)0，与已知矛盾。因而只能1，即(L)1，于是SOR法收敛。

定理10 设A对称，且对角元aii0 (i1,,n)，则Jacobi迭代法收敛的充要条件是A及2DA均正定，其中Ddiag(a11,a22,,ann)。

证明参考关治，陆金甫编《数值分析基础》。 定理11 若A为对称正定矩阵，则解Axb的SOR法收敛的充要条件为02。 证明 只需证明充分性。设是L的任一特征值（可能是复数），若能证明1，则定理得证。设y是L的相应于的特征向量（可能是复向量），即

(DL)1((1)DU)yy

也就是

((1)DU)y(DL)y

上式两边与y作复内积，有

(((1)DU)y,y)((DL)y,y)

则

(1)(Dy,y)(Uy,y)

(Dy,y)(Ly,y)

25

利用正定阵的对角元素大于0，有

(Dy,y)aiiyi2p0

i1n记

(Ly,y)i

由于AAT，所以LTU，故

(Uy,y)yT(UT)Ty=(Ly)Ty=(y,Ly)=(Ly,y)i

于是

(1)pi pi222(pp)2 222(p)而

(pp)2(p)2p(2)(2p)

注意到02及A的正定性

2(Ay,y)((DLU)y,y)p20

2可知的分子小于分母，即1，从而(L)1，SOR法收敛。

推论 若A对称正定，则GS迭代法收敛。

3.3 误差估计 定理 12 设有方程组

26

xBx+f

及一阶定常迭代法

x(k1)Bx(k)f

如果有B的某种算子范数Bq1，则

(1) 迭代法收敛，即对任意的x(0)有

limx(k)x*且x*Bx*f

k(2) x*x(k)qkx*x(0) (3) xx*(k)(4) x*x(k)qx(k)x(k1) 1qqkx(1)x(0) 1q证明 由定理 6，结论(1)是显然的事实。 (2) 由关系式

*(k)*(k1)xxB(xx)

及

(k1)(k)(k)(k1)xxB(xx)

有

(a) x(k1)x(k)qx(k)x(k1) (b) xx*(k)qxx*(k1)

反复利用(b)即得(2)。 (3) 利用(b)得

27

x(k1)x(k)x*x(k)x*x(k1)(1q)xx*(k)

于是

xx*(k)1q(k1)(k)(k)(k1) xxxx1q1q(4) 反复利用(a)，则得到(4)。

注：1）利用定理11作误差估计，一般可取1，2或范数。结论(3)是近似解x(k)的误差事后估计式，对于给定的精度（当然应当选得恰当，小于或接近于机器精度可能会造成死循环），只要q不是很接近1，则可用x(k)x(k1)来控制迭代终止。若q1，即使x小，也不能判定x(k)x*很小。

2）结论(4)可用作迭代次数的估计。根据事先给定的精度，可以估算出迭代的次数：

lnk(k)x(k1)很

(1q)x(1)x(0)lnq

迭代法是否收敛虽与初始向量x(0)的选取无关，但由上面的公式看出对迭代次数却有很大的影响，因而应重视初始向量x(0)的选取。

28

3.4 迭代法的收敛速度及最佳松弛因子 对于迭代法，我们不但要知道它是否收敛，而且对它收敛的快慢感兴趣。本节介绍迭代收敛速度的概念。

nn*AR为非奇异矩阵，设x是(1.1)的解，即

**xBxf

以迭代公式(1.2)减去上式，并记误差向量e(k)x(k)x*，则有e(k1)Be(k)(k0,1,,n)。递推得

e(k)Bke(0)

对任何算子范数有e(k)Bke(0)，于是

ee(k)(0)B

k根据算子范数定义，有

k(0)(k)BeekBmaxmax (0)(0)(0)(0)e0e0ee所以B给出了迭代k次后误差向量范数与初始误差向量范数之比的最大值。这样，迭代k次

后，平均每次迭代误差范数的压缩率就可以看成是B

k1kk。如果要求迭代k次后有

29

ee1kk(k)(0)，有

ee(k)(0)B (3.3)

k其中1，可取10。因为(B)1，所以B1，由B1kks两边取对数，得到

1kk1klnB1ln k即

klnlnB1kk1kksln10lnB1kk (3.4)

表明迭代次数k与lnB成反比。因此，可以

1kk用这个量来度量迭代格式的收敛速度。 定义 6 称Rk(B)lnB为迭代法的平均

收敛速度。

值得注意的是平均收敛速度Rk(B)依赖于迭代次数及所取范数，给计算分析带来不便，由于limBk1kk(B)，所以limRk(B)ln(B)。

k定义7 称R(B)ln(B)为迭代法的渐近收

敛速度。

30

显然R(B)与迭代次数及范数选取无关，它反映了迭代次数趋于无穷时迭代法的渐近性质。 为了达到(3.3)的要求，可以用

lnsln10k

R(B)R(B)代替(3.4)作为所需迭代次数的估计。 当SOR方法收敛时，通常希望选择一个最佳的值opt使SOR方法的收敛速度最快。然而遗憾的是，目前尚无确定最佳超松弛因子opt的一般理论。实际计算时，大都由计算经验或通过试算来确定opt的近似值。目前仅对某些特殊类型的矩阵有确定opt的公式。例如针对一类椭圆微分方程数值解得到的线性方程组Ax=b，当矩阵A具有所谓“性质A”及相容次序时，Young于1950年给出了一个最佳松弛因子计算公式

opt211(J)2 (3.5)

其中(J)是Jacobi迭代公式的迭代矩阵J的谱半径。

然而，在实际应用中，一般来说(J)计算也是比较困难。因此人们多使用计算经验或试算

31

的方法确定opt的一个近似值。

所谓试算法就是从同一初始向量出发，取不同的松弛因子迭代相同的次数k（注意：迭代次数k不能太少），然后比较其相应的剩余

(k)(k)(k1)rkbAx（或xx），并选取使其范数最小的松弛因子作为最佳松弛因子opt近似值。实践证明，此方法虽然简单，但往往是行之有效的。

下面定理给出了三对角形正定矩阵A的最佳松弛因子。

定理 13 如果A为三对角形对称正定矩阵，J和G分别是Jacobi和Gauss-Seidel迭代法的迭代矩阵，则

2(1) (G)(J) (2) 最佳松弛因子

opt211(J) 2(3) 松弛迭代矩阵L的谱半径(L)1。

32

例 10 （模型问题）考虑泊松（Possion）方程的边值问题

2u2u(22)f(x,y),(x,y), (3.10) yxu(x,y)0, (x,y), (3.11)其中{(x,y)|0x,y1}，为的边界，用差分法求解边值问题(3.10)式和(3.11)式。

如图6-1所示，用直线xxi,yyi在打上网格，其中

1xiih,yijh,h,i,j1,2,N.

N1分别记网格内点和边界点的集合为

h{(xi,yi)|i,j1,2,,N},h{(xi,0),(xi,1),(0,yi),(1,yi)|i,j0,1,,N1}

在点(xi,yi)上用差商表示二阶偏导数，即

2u12|[u(x,y)2u(x,y)u(x,y)]o(h),i1iiii1i2(xi,yi)2xh2u12|[u(x,y)2u(x,y)u(x,y)]o(h),ii1iiii12(xi,yi)2yh

略去余项o(h2)，用uij表示u(xi,yi)的近似值，由微分方程(3.10)就可以得到差分方程

33

( y 1 yj y1O x1 图6-1 xi 1 x ui1,j2uijui1,j2hh其中fijf(xi,yi)，再整理成

ui,j12uijui,j12)fij,

4uijui1,jui1,jui,j1ui,j1hfij (3.6)

2其中(i,j)对应的点(xi,yi)h。(3.6)式成为泊松方程的五点差分格式。(3.6)式左端若有某项uij对应的点(xi,yi)h，则uij0，为将差分方程写成矩阵形式，我们把网格点逐行按由左至右和由上至下的自然次序排列，记向量

u(u11,u21,,uN1,u12,u22,,uN2,,u1N,u2N,,uNN)Tbh(f11,f21,,fN1,f12,f22,,fN2,,f1N,f2N,,fNN)2T

则(3.6)式可以写成

Au=b (3.7)

34

其中

A11IAIA22IIAN1,N1IN2N2R IANN(3.8)

41141AiiRNN,i1,2,,N 14114(3.9)

通常N是个大数但A的每一I为NN单位矩阵，

行最多只有5个零元素，所以A是一个稀疏矩阵，故线性方程组(3.7)是一个大型系数方程组，它可用SOR迭代法求解，可算出JD1(LU)的特征值为

1ij(cosihcosjh),i,j1,2,,N.

2当ij1时得到J的谱半径

35

122(J)cosh1ho(h4)

2由于A对称正定，故SOR迭代法收敛，且可利用(3.5)式求出最优松弛因子

opt2 1sinh且

cos2h (Lopt)opt12(1sinh)根据(1.7)式及收敛速度定义可得

122R(J)ln(J)ho(h4),

2R(Lopt)ln(opt1)2[lncoshln(1sinh)]32ho(h)可见R(Lopt)比R(J)大一个h的数量级，若取

Th0.05,f(x,y)0.初值取u(0)=(1,1,,1)，计算到

u(k)u(k1)106停止，则雅可比法需要1154次

迭代，而SOR迭代法若取1.73则只需59次迭代。h0.05时，opt1.72945。

3.5 块迭代法 块迭代法是用于型稀疏线性方程组求解。 对方程组Ax=b，设ARnn可写成分块形式

36

A1NA2N

ANN其中Aii为nini的方阵，向量xn1n2nNn。

A11AA21AN1A12A22AN2及b同样分块

xx1x2xN,bb1b2bN

TT其中 xi和bi都是ni维的向量。令

A=DLU

T其中DA11A22ANN，

0AL21AN10AN2(k1)i0,U0N(k)jA120A1N

AN1,N0块雅可比迭代法为

AiixbiAijx, i1,,N (3.10)

j1ji块雅可比迭代法的矩阵形式

x(k1)Jx(k)f (3.10)

其中JID1AD1(LU),fD1b。 块超松弛法（BSOR法）：

Aiixi(k1)(1)Aiixi(k)(biAijx(jk1)j1

37

i1ji1NAijx(jk))

i1,2,,N (3.11)

BSOR法的矩阵形式

x(k1)Lx(k)f (3.12)

其中 L(DL)1((1)DU),f(DL)1b。当1时，(3.11)和(3.12)就是块GS迭代法。 实际计算中，对每个i，(3.10)和(3.11)是一个nini的方程组，一般用直接方法求解。对大型方程组，n是个大数，而ni相对来说较小的数。 我们给出下述结果。

定理 14 设Ax=b，其中ADLU（分块形式）。

（1） 如果A为对称正定矩阵， （2） 02。

则解Ax=b的BSOR迭代法收敛。

例 10的模型问题中(3.8)式和(3.9)式所表示的分块形式和一般形式相比，有niN，例中(3.8)式的分块对应于图6-1的一条条网格线，按分块形式写出的迭代公式也称线迭代法。在BSOR迭代法的收敛性和最优松弛因子的理论分析中，一类特殊的三对角块矩阵有很多好的性质，它就是T-矩阵，其形式为

38

F2 (3.14)

DN1FN1ENDN的块三对角矩阵，其中对角块Di(i1,2,,N)均

D1E2AF1D2EN1为对角矩阵。

记Ddiag(D1,D2,,DN)，块雅可比矩阵JID1A。设块SOR（BSOR）方法的迭代矩阵为L，则有以下结论。

定理 15 设A为非奇异的形如(3.14)式的T-JID1A，矩阵，且D非奇异。则当(J)1时，

对02有(L)1及最优松弛因子

opt211[(J)]2,(L)opt1

且

1222[4(1)], 0<opt,(L)4

1, opt2,(3.15)

其中(J)。

证明见[33]。根据定理有

39

(L)min(L)

opt

02如图6-2所示。由(3.15)式可知，当1时，(G)(L)22(B)，则得高斯-赛德尔迭代法的收敛速度为

R(G)ln(G)2ln(J)2R(J)。

说明此时的高斯-赛德尔迭代法比雅可比迭代法快一倍。由于T-矩阵的特殊形式就是三对角矩阵，因此，当A为对称正定的三对角矩阵时SOR法的最优松弛因子就是(3.5)式给出的。 注意对例10的模型问题得到的(3.8)式的矩阵A是按自然排序得到的，它不是T-矩阵。如果改变网格点的排序，通常成为红-黑排序，则

A可变成T-矩阵。这里不作介绍，参见[2,33]。

(L)

1



O  1 opt 2 图 6-2 40

4共轭梯度法

4.1与方程组等价的变分问题

共轭梯度法简称CG（conjugate gradient）方法，又称共轭斜量法，它是一种变分方法，对应于求一个二次函数的极值。 设A(aij)Rnn对称正定，b(b1,,bn)T，求解的线性方程组为

Axb (4.1)

考虑如下定义的二次函数:RnR:

['kɔndʒugit]

['greidiənt]

1(x)(Ax,x)(b,x)2 (4.2) nnn1aijxixjbjxj2i1j1j1函数有如下性质： n(1) 对一切xR，(x)的梯度

(x)Ax-b (4.3)

事实上，

41

n(x)1nn[(aijxj)xixkakjxj]bkxk2xki1j1j1ikn1n(aikxiakjxjakkxk)bk2i1j1ikn1n(aikxiakjxj)bk2i1j1

akjxjbk (k1,2,,n)j1n(2) 对一切x,yRn及R，

1(xy)(A(xy),xy)(b,xy)2(x)(Axb,y)2 (4.4)

2(Ay,y)(3)设x*A1b是线性方程组(4.1)式的解，则有

111(x)(b,Ab)(Ax*,x*)

22*实际上

1(x)(Ax*,x*)(b,x*)21(b,A1b)(b,A1b) 21(b,A1b)2*

42

且对一切xRn，有

11***(x)(x)(Ax,x)(Ax,x)(Ax,x)22 (4.5) 1(A(xx*),xx*)2*以上性质可根据定义(4.2)式直接运算验证。 定理14 设A对称正定，则x*为线性方程组(4.1)解的充分必要条件是x*满足

(x*)min(x) (4.6)

xRn证明 设x*为(4.1)的解，则x*A1b。由(4.5)式及A的正定性有

1(x)(x)(A(xx*),xx*)0

2**nx所以(x)(x),xR，即使(x)达到最小。

*反之，设x*为(x)的极小点，则

(x)xkxx*1nn[(aijxj)xibkxk]xx*2xki1j1akjx*jbk0 (k1,2,,n)j1n

即Ax*b。

这个定理把求解方程组(4.1)等价于一个多元二次函数（也称泛函）极小化的问题，这种等价性开辟了设计新算法的一个途径。定理14

43

称为求解线性方程组的变分原理或Ritz原理。求解方法是构造一个向量序列{x(k)}使(x(k)) (x*)。

4.2最速下降法

最速下降法的基本思想是从问题(4.6)的近似解x(k)出发，沿使函数(x)下降最快的方向，寻找新的近似解x(k1)，使得(x(k1))在该方向上达

*到极小值，这样一步步逼近问题(4.6)的解x。这就是最速下降法的基本思想。由多元微分学知，函数值下降最快的方向是负梯度方向。 具体做法是：从某个初始点x(0)出发，沿(x)在x(0)处的负梯度方向r(0)(x(0))bAx(0)(1)（称为搜索方向）求得(x)的极小点x，即

(x(1))min(x(0)r(0)) (4.7) R由于

(x(0)r(0))(x(0))(Ax(0)b,r(0))22(Ar,r)(0)(0)

为确定，由

d(x(0)r(0))(Ax(0)b,r(0))(Ar(0),r(0))0

d解得

44

0(r(0),r(0))(Ar(0),r(0))

于是得到

x(1)x(0)0r(0)(0)(0) rbAx(r(0),r(0))(Ar(0),r(0))0然后从x(1)出发，重复上述过程得到点x(2)。如

(k)x此继续下去，得到序列。最速下降法的一般公式

x(k1)x(k)kr(k)(k)(k) (4.8) rbAx(r(k),r(k))(Ar(k),r(k))k实际上二次函数(4.2)的几何意义是一族超椭球面(x)(x(k))（(x(k))(x(k1))），x*为它的中心，若n2就是二维空间的椭圆曲线，我们从(k)x出发，先找一个使函数值(x)减少最快的方向，这就是正交于椭球面的函数(x)的负梯度方向 由于

(r(k1),r(k))(bA(x(k)kr(k)),r(k))(r(k),r(k))k(Ar(k),r(k))0

说明相邻两次的搜索方向是正交的。

45

可看出{(x(k))}是单调下降有下界的序列，因此它存在极限。实际上

(x(k1))(x)k(Ax(k)(k)b,r)(k)(k)k22(Ar,r)(k)(k)k(r,r)(k)2k

(k)(k)1(r(k),r(k))2 0（A正定） (k)(k)2(Ar,r)说明由(4.8)式得到的{(x(k))}是单调下降有下

2(Ar,r)界的序列，它存在极限，满足

limx(k)x*A1b

k而且

1n(0)*xxxx AA1n其中1,n分别为对称正定矩阵A的最大与最小

(k)*k特征值。uA(Au,u)，当1n时收敛是很慢

(k)的，而且当r很小时，由于舍入误差影响，计算将出现不稳定，所以这个算法实际中很少使用，需要寻找对整体而言下降最快的算法。

4.3共轭梯度法（CG方法）

共轭梯度法简称CG法，是一种求解大型稀

46

12

疏对称正定方程组十分有效的方法。

定义8 设A为n阶对称正定阵，x,y是Rn的向量，如果有(x,Ay)(Ax,y)0，则称x与y关于A-共轭正交，简称为A-正交。

(k)共轭梯度法的基本步骤是在点x处选取搜索方向p(k)，使其与前一次的搜索方向p(k1)关于A共轭，即

(p(k),Ap(k1))0 (k1,2,) (4.9)

然后从点x(k)出发，沿方向p(k)求得(x)的极小点x(k1)，即

(x(k1))min(x(k)p(k)) (4.10) 由此得到方程组(4.1)的近似解序列x(k)。 下面推导其计算公式：

第一步与最速下降法相同，即

x(1)x(0)0p(0)(0)(0)(0) prbAx(r(0),p(0))(p(0),Ap(0))0一般地，令

x(k1)x(k)kp(k)

取p(k)r(k)k1p(k1),r(k)bAx(k)。使

(p(k),Ap(k1))0

即

47

(r(k)k1p(k1),Ap(k1))(k)(k1)(k1)(k1)(r,Ap)k1(p,Ap)0(r(k),p(k1))故得k1(k1)(k1)。令

(p,Ap)g()(x(k)p(k))min(x(k)p(k))



dg()0，可解得k(r(k),p(k))(p(k),Ap(k))。利用 d综合在一起，得到共轭梯度法的一般公式

(k)(k1)(k1)(k1)k1(r,Ap)(p,Ap)(k)(k)(k1)prk1p (k1) (4.11)

k(r(k),p(k))(p(k),Ap(k))x(k1)x(k)kp(k)r(k)bAx(k)经过简化可得

k1(r(k),r(k))(r(k1),r(k1))k(r(k),r(k))(p(k),Ap(k)) (4.12)

事实上，因为

r(k1)bAx(k1)bA(x(k)kp(k))r(k)kAp(k)

有

(r(k1),p(k))(r(k),p(k))k(Ap(k),p(k))0,(r,p)(r,r

(k)(k)(k)(k)k1p48

(k1))(r,r)(k)(k)

定理17 对任意i,j且ij，有

(r(i),r(j))0

可用数学归纳法证明。

CG算法

))Rn，计算r(0bAx(1) 任取x(0)0 p(0r(。

(2) 对k0,1,,计算

(r,r)k(k)(k)(p,Ap)x(k1)x(k)kp(k)r(k1)(k)(k)(0，取

r(k)(r,r)kAp,k(r(k),r(k))(k)(k1)(k1)

p(k1)r(k1)kp(k)(k)(k)(k)r0(3) 若，或(p,Ap)0，计算停止，(k)*(k)(k)xx则。由于A正定，故当(p,Ap)0时，

p(k)0，而(r(k),r(k))(r(k),p(k))0，也即r(k)0。

例11 用CG法解线性方程组

3x1x25 x12x2531解显然A是对称正定的，取12(0)T(0)(0)(0)Tx(0,0)prbAx(5,5)， ，则

49

(r,r)20,(0)(0)(Ap,p)7

1010T(1)(0)(0)xx0p(,),7755T(1)(0)(0)rr0Ap(,),77(r(1),r(1)) 0(0)(0),(r,r)3040T(1)(1)(0)pr0p(,)49497(2)类似可计算出1,x(1,2)T为方程组的精

10(0)(0)确解。

评论：1、由于{r(k)}互相正交，故在

(k)(0)(1)nr,r,,r中至少有一个零向量。若r0，则x(k)x*。所以用CG算法求解n维线性方程组，理论上最多n步便可求得精确解，从这个意义上讲CG算法是一种直接法。但在舍入误差

(k)存在的情况下，很难保证{r}的正交性，此外当n很大时，实际计算步数kn，即可达到精度要求而不必计算n步。从这个意义上讲，它是一个迭代法，所以也有收敛性问题。 2、可以证明CG算法有估计式

50

x(k)K1(0)* (4.13) x2xxAAK1*12k其中xA(x,Ax),Kcond(A)2(证明见[34])。 3、共轭梯度法的优点是存储量小，计算简便。它能充分利用矩阵的稀疏性，计算时只需存储A的非零元素。共轭梯度法具有超收敛性。大量的数值经验表明，当A的条件数很小时，用此方法仅需迭代很少几步就能得到高精度的解。

4、由估计式(4.13)看出K1，即A为病态矩阵时，CG法收敛很慢。为改善收敛性，可采用预处理方法降低矩阵的条件数，从而可得到各种预处理共轭梯度法，此处不再介绍，可参见文献[2,8].

51