chap3解线性方程组的迭代法

第三 解线性方程组的迭代3‐1：取0,10,1u u f22sin(x)sin(y，使用有限元法求解此4X4169

c1 c

2 000000 000000 000000000000 00000000000 000000 c 4 c 5 c

7

一般情形：对区域进行64X64划分；求所有节点上近似数据，最后化归为求解内部节点满足的稀疏方程例3‐2：现实中的问题大多数是连续的，例如工程中求解结§3.1迭代法直接法适用于中小型方程组对高阶方程组,量大,程序复杂,运算量巨大.矩阵稀疏性,计算简单,程序编制容易.迭代法基本思想：构造向量序列xk,使得它的极限x是Axb的解迭代法的一般其中AM为n阶方阵，xg ：x(k1)Mx(k)

k0,1,2,M为迭代矩阵，{x(k)}为迭代序列结论：若迭 x(k

Mx(k

g产生的迭代序列收敛于AxxMx则必有xAxxMx此类方法称为单步定长线性迭代法向量序列与矩阵序列的收敛定义3.1：设{x{k}}为Rn中的向量序列,xRn,如果：limx(k)xk其中为Rn中的向量范数,则称序列{xk)}收敛于记 limx(k)k定理 limx(k)xlimx(k)x(i1,2, k 其 x(k)(x(k),x(k),...,x(k))T

x(x,x,...x 定义3.2：设Ak为n阶方阵序列,A为n阶方阵,如果 A(k)Ak其中为矩阵范数,则称序列Ak)}收敛于A记作limA(k)k定理3.2 limA(k)Alima(k)a(i,j1,2,k k

其 A(k)(ak A(a定理3.1,3.2说明向量和矩阵序列的收敛等价于对应分量 几种基本的迭 雅可比（Jacobi）a11x1a12x2a1nxna21a

a22

a2n

Ax xMx

annxnx1

m1nx1

g1x

x g2 2n22

xmxmx xxn

mnnxn gn

21

若aii (i1,2,,n)方程组可同解变形为x

xa1n

x

x

ann1

xMx

x(k1)Mx(k) k0,1,nJacobi迭代法的计 n

x(k)

nn

x(k)

22

n n即

i

n

ix(k1)i

ba x(k)

x(k)

(i1,2,,ai ai

ji

记 g

ba

(ij i,j1,2,,(i1,2,,n

g1

gM

2n

g 2bb

0 0

g g nx(k1)Mx(k)Jacobi迭代法的矩阵形式 x(k1)Mx(k)选取初始向量x(0),按以上 产生的迭代序列称为Jacobi迭代（简单迭代法）Jacobi MI g

Ddiag(a 算法3.1(Jacobi迭代、输入Aab(bb,...b)Tx(0)

(x(0),x(0),...x(0))T 误差,最大允许迭代次数N2、取k

i

(0

(03、对i1,2,..., xi (bi

aijxj

ji

aijx 4、若

xx(

,输出,停止.否则转5、若kN,取k1kxx(0i12 转3,否则输出失败信息functionM%用途：Jacobi迭代求解方程组n=length(b);N=500;ep=1e-x=ones(n,1k=0;whilek<Nfor%[kx’] ifnorm(x-x0,inf)<ep,break;endifk==N,Warning(‘已达最大迭代次数');enddisp([‘k=’,num2str(k)]) 10x1x22x3x10 2 2x 5 2310x1x22x3解：Jacobi迭代计

x10 2 3x 5 31x(k112x(k12x (k1x

0.1x(k)0.2x(k 0.1x(k)0.2x(k 0.2x(k)0.2x(k

788x(0)0,0,0)T

x(1)(7.2,8.3,1x(2)0.831.687.212x(2)0.721.688.323x(2)1.441.668.43依次下去收敛于真解x(11, kJacobimethodx

10x1x22x3 10 2 5 0.1

MID1A

0gD1b

令x(k1)Mx(k)

x(0)(0得x(1)Mx(0)g 8.4)Tx(9)(10.9994,11.9994,依次下去，x(k1)Mx(k) 收敛于精确x(11,12,13)T －赛德尔(Gauss－Seidel)Jacobi迭代的计 x(k

x(k

x(k

a a

nx(k n

21x(k

a2

x(k

2x(k2

x(k

an

x(k

ann

x(k

a a

i

n

n

x(k1) baijx(k)

aijx(k)

(i1,2,,ii jii

ji Gauss-Seidel迭代的计 x(k

x(k

x(k

a a

nx(k n

21x(k

x(k

2x(k2

x(k

an

x(k

x(k

a a即

x(k

b

i

ax(k1)

ax(k)

(i1,2,,na na

j

ji 0 0

0La

a 00 0 an an

0

a1n0 a 0

2nU

a3n00

Ddiag(a11,a22,,annADL(k 1

i

(k

(k)a a

aijxj

aijx ji

(i1,2,,x(k1)D1(Lx(k1)Ux(k

(DL)x(k1)Ux(k)x(k1)(DL)1Ux(k)(DL)1迭代矩阵为

Ms(D functionM%用途：Gauss-Seidel迭代求解方程组whilefor ifx(1)=(b(1)-elseifx(i)=(b(i)-A(i,1:i-1)*x(1:i-1)-ifnorm(x-x0,inf)<ep,break;end10x1x22x3 10 2 5 10x1x22x3解：G-S迭代计

x10x2 xx5 x(k 0.1x(k

0.2x(k)

2x(k2

0.1x(k

0.2x(k)13x(k 0.2x(k1)0.2x(k1)13

取x(0)(0,0,0)T x(1)(7.2,9.02,1x(2)0.9022.32887.212x(2)1.043082.32888.323x(2)2.086162.334388.43依次下去收敛于真 x(11, er(A,b,20,1e- 0,1e-Gauss_seidelmethodxJacobimethodconvergedx=10x1x22x3 10 2 写成G-S迭代矩阵

2 5 23

01

0.2

0.22MS(D

g(DL)1b S令xk1)MxkS

x(0)(0,依次下去,x(k)收敛于真解 逐次超松弛法（SOR法换个角度看GaussSeidelx(k1) 1

i

x(k1)

(k)

aijx

j1

i

jin

r(kix(k) bax(ki

ax(k)x(k

aii

iij

j

r(k x(k1)x(k

a通过选取合适的来加速收敛1时即为Gauss-SeidelSOR计 (k

(k

i

(k

(k

(i

(bia

aijx

aijx j

i

j (1)x(k)

(b

ax(k

ax(k)a a j

ji

(k1)

b

i ax(k1)

(k)

aij ji x(k

(1-)x(k

(

x(k

x(k

b1

x(k

(1-)x(k

(

x(k

a2

x(k

b2SOR:

x(k

(1-)x(k

(

x(k1)ann

x(k

bn

n (k

(k)

i (k

(k (1)

(bi

aijxj

aijxj ji

(ix(k1)(1)x(k)D1(bLx(k1)Ux(k)Dx(k1)(1)Dx(k)(bLx(k得

Ux(k)xx(k1)(DL)1(1)DUx(k)(D迭代矩阵为

M(DL)1(1)DUfunctionM%用途：SOR迭代求解方程组while for ifx1(1)=(b(1)-elseifx1(i)=(b(i)-A(i,1:i-1)*x(1:i-1)-ifnorm(x0-x,inf)<ep,break;end松弛因子的选取对收敛速度影响极大,最优松弛因子A为对称正定的三对角矩阵时,有：1 11 1(B)22其中(B)为Jacibi迭代矩阵的谱半径 迭代方法评JacobiG-S:减少 量,要求计算顺SORG-SG-SJacobi

x(k1)D1(LU)x(k

x(k1)(DL)1Ux(k)(D x(k1)(DL)1[(1)DU]x(k)(D矩 再探讨：Ax xMx AE

E1存在且方程组Exd容易求解AxAxbExFx xE1(Fx

a1n

0

a1n

0 a

2n

2n

nn

nn AAD(LU

Axb Dx(LU)xb

x(k1)D1(LU)x(k)A(DL)

Axb (DL)xUxG-S:x(k1)(DL)1Ux(k)(DA1DL1DU AxbAxb

A(DL)[(1)DU(DL)x[(1)DU]x

x(k1)(DL)1[(1)DU]x(k)(DA1D1DA JOR:x(k1)(ID1A)x(k)A1I1IA

Richardson

x(k1)(IA)x(k)GeneralRichardson

diag(1,2,,nx(k1)(IA)x(k)A1DL1DU 1DU1DL

x(k1)Sx(kSM

M(DU)1[(1)DN(DL)1[(1)DUg(2)(DU)1(DL)1 矩阵的谱定义3.3设A为方阵，i(i12,n)为A的特征值，称特征值模的最大值为矩阵的谱半径，记为：(A)max结论：、Ak)、(A) A其中为任意由向量范数诱导出的矩阵范数3、0,存在一种矩阵范数,使得A(A) 2当A为对称矩阵时,有(A) 2

(A)(A)4An

limAk(A)k4、设A为n阶方阵limAkA)k证明

Ak

[(A)]k (Ak)

(A)""对任意存在一种范数，使A(A)(A)1，存在范数，使 A而

A

limAk limAkx,xk k迭代法的收敛定理3.6对任意初始向量x0右端项x(k1)Mx(k

(k0,1,2,)产生的向量序列{xk)}收敛的充要条件是：(M)推论：若 则对任意初始向量x( 右端项x(k1)Mx(k)

(k0,12,产生的序列x(k)收敛定理3.6:xk1收敛(M)

x(k1)Mx(k)limAkx,xRnlimAkk k

limAk(A)k"若xk1)x*

x*Mx*x(k

x*Mx(k)

Mk(x(0)x*Mk (M

(M)

（IM）xg存在唯一解xk1x*Mx(k)Mx*Mkx0)x*(M)1Mk x(k1)x*推论：逐次松弛法收敛的必要条件是：0证明：逐次超松弛法的迭代矩阵为M(DL)-1[(1-)DUdet(M)(DL)- (1-)D (1-)n (1-)na11

det(M

[(M

由定理 (M)det(M

(1)n 得0 x12x22x3x 2x2 2 2A

2L

U 2

0 00DL (DL)1 00

2 ID1A1

0IMJ

3MJ0Jacobi迭代法 0

2 0 2M

0

2

0

2I

2

(2)2定义：若n阶方阵A(aij)满足n

j

(i1,2,,且至少有一个i值，使上式不等号严格成立，则称A为弱对角占优阵。若对所有i,上式中不等号均严格成立，则称A为严格对角占优阵。

5 5 A

65 A65

定义3.5---不可约矩阵如果矩阵A不能通过行的互换和相应列的互换成为形 22其中 A22为方阵，则称矩阵A不可约不存在排列矩阵P,使PTAP

A12 22 0

A

B1 2结论若A没有零元素,则A

n3时A只有一个零元素A不可设A(aij)为n阶矩阵(n2若存在非空集合I1,2,使得当iIjI时，有aij0则A是可约阵。几种常用的判别收敛条件：已知Ax1、若A为严格对角占优阵或不可约弱对角占优阵，则Jacobi迭代法与G-S迭代法收敛。2、若A为严格对角占优阵，01则SOR法收敛3、若A为对称正定矩阵,则SOR法收敛的充要条件是02.、Jacobi:A为严格对角占优或不可约弱对角2、G-SA为严格对角占优不可约弱对角占优3、SORA为对称正定矩阵(0严格对角占优(0

判别迭代法收敛的步骤1、先观察系数矩阵A是否对称正定或对角占3、给出迭 ,讨论迭代矩阵的谱半例：JacobiGauss-SeidelAxb,

A b

-1 A 1

A对称正定，G-S迭代法收 12 22 1Jacobi的迭代矩阵B

1(LU) 21 1 0 IB03－5－3 B)1Jacobi迭代发散。

A为不可约弱对角占优,Jacobi和G-S迭代法均收例

Ax 1 A 0.5 b 2 0.5

2

2 3

1A1Jacobi,G-S01的SORJacobi迭 x(k1)0.9x(k

0.05x(k) x(k

0.25x(k

0.5x(k

x(k1)0.1x(k(k(k

G-S迭 1x(k1)0.9x(k1

0.05x(k)x(k

0.25x(k

0.5x(k)(k

1(k

3(k

2323

01x(k1

(1)x(k

(0.9x(k)0.05x(k

1x(k1

(1)x(k

(0.25x(k

0.5x(k)(k

2(k

1(k

3(k1)1

(1)x3

误差估定理3.7设有迭代格式xk1)Mx(k)Mk收敛于x,且有误差估Mk

g若M1,则x(kx(k

1

x(1)x(

(3 M1(M)3.6,limx(k)xxMxx(k

xM(xk1x)Mk(x(0)x两边取范

x(k

M

x(

x* x(0)

x(0)x*Mk又 Mk

1,所以有：

x(0

x

x x(01 即：xk1

1

x(1)x(0

证毕。注 由定理3.7易得 M越小,{x(k)}收敛越若给定精度(误差限),由(322)(1(1Mx(1)x(0)lnM定理3在定理3.7的条件下Mx(k)

1

x(k)x(k证明：x(k)两边取

xM(x(k1)xx(k)

x(k

M同时由M

x(k

x(k1)x(k

x(k

当 1时

x(k

1

x(k)x(k注2.

x(k)x(k

作为停止准 最速下降法与共轭梯求解Axb,

其中A对称正定可转化为极值问 (x)1xTAxbTx1Ax,xb,x 因 (x)Ax求解Axb求x)的极小值 Ax*b min(x)1x*TAx*bT

(x)1xTAxbTx1Ax,xb,x 函数(x)具有以下性质：梯 (x)Axx,yRn,2(x+y)1A(x+y),(x+y)b,(x+y)2=(x)Axb,y2 Ay,2Ax*，则x*)1bx*1Ax*,x* 且xRn,xx*)1Axxbx1bx* =1A(xx*),(xx*)23.：(x*)min(证明：必要性.Ax*b,由性质(3)及对称正(x)(x*)1A(xx*),(xx*)2即x*)min充分性

i

ai2

ain

)

i1,2,,=Ax若(x*)min( 则(x)在x*处从而Ax*b证最速下

(x)1xTAxbT2(x)Axstep1:从初始点x0)出发找负梯度方向r0)r0)x0bAx0搜索方向step2:在r0方向上求x)的极小值点x(1(x(1))min(x(0)r(0)