迭代法求解线性方程组

第六章 线性方程组的迭代解法
/* iteration methods for the solution of linear systems */
Linear systems： ：
a11 x1 + a12 x 2 + L + a1n x n = b1 a x + a x + L + a x = b 21 1 22 2 2n n 2 L a n1 x1 + a n 2 x2 + L + a nn xn = bn
A = D + A − D ⇒ x = ( I − D−1 A) x + D−1b = x + D−1 (b − Ax) ( aii ≠ 0) A = L + D + U ⇒ x = −D−1 ( L + U ) x + D−1b ( aii ≠ 0)
迭代法思想： 第一步 第二步
将Ax = b转化为等价方程组：x = Bx + f 设x *为精确解，x* = Bx * + f
计算公式
1 b ( k +1) (k (k (k x1 =− a12 x 2 ) + a13 x3 ) + ... + a1n x n ) + 1 a11 a11 1 b2 ( k +1) (k ) (k ) (k ) =− + a 21 x1 + a 23 x3 + ... + a 2 n x n + x2 a 22 a 22 ... bn 1 ( k +1) (k ) (k ) (k ) =− a n1 x1 + a n 2 x 2 + ... + a nn −1 xn −1 + xn a nn a nn
x ( 2 ) = (5,−3,−3)T , x ( 2 ) − x (1)
x (3) = (1,1,1)T , x ( 3) − x ( 2 ) x ( 4 ) = (1,1,1)T , x ( 4 ) − x ( 4 )
∞ ∞
=5
∞
=8
=4 =0
所以x* = x ( 4 ) = (1,1,1)T
x(k), k=0,1,2,… hopefully, limx(k)=x*
Ax=b A x*=b
x(k+1)=f(x( k ))
x(0)
怎样设计迭代格 式? Ax=b x=Bx+f
等价变换
迭代矩 阵
x(k+1)=B x(k) +f
B 并不是 唯一的, 因 此迭代格 式也不是 唯一的!
e.g. A = I + A − I ⇒ x = ( I + A) x − b
计算中判断迭代终止条件的方法：
|| x ( k +1) − x ( k ) ||< ε , x* ≈ x ( k +1)
6.2基本迭代法
A = M − N , M非奇异，称为分裂矩阵 则Ax = b ⇔ ( M − N ) x = b ⇔ Mx = Nx + b ⇔ x = M −1 ( Nx + b) 令B = M −1 N , f = M −1b即可得等价方程组。
b a ij x (jk ) + i , ∑ j = i +1 a ii
n
i = 1, 2, ..., n
comparison
Jacobi iteration 计算x(k+1)时需要x(k)的 所有分量，因此需开 两组存储单元分别存 放x(k)和x(k+1)
Gauss-Seidel iteration 计算xi(k+1)时只需要x(k) 的i+1~n个分量，因 此x(k+1)的前i个分量可 存贮在x(k)的前i个分 量所占的存储单元， 无需开两组存储单元
构造迭代格式，任取初始向量x ( 0) , 按逐次代 B 与k无关，称为
一阶定常迭代法 入方法构造{x ( k ) } 收敛？发散？
x ( k +1) = Bx ( k ) + f
若 lim x ( k ) = x*, 则称迭代法收敛，且x * 为方程的解。
k →∞
否则此迭代法发散。
判断收敛的方法：
How to check if a certain iteration system converges or not?
Conditions of convergence
Not flexible to use actually
1. lim x k = x * ⇔ lim B k = 0
k →∞ k →∞
(
)()() Nhomakorabeaor
xi
( k +1 )
1 =− a ii
j =1, j ≠ i
∑
n
a ij x j
(k )
bi , + a ii
i = 1, 2,..., n
上面例对应的分量形式为 x ( 0) = (0,0,0)T (k ) ( ( x1 = (1 − 2 x2k ) + 2 x3k ) ) / 1 (k ) ( x2 = (3 − x1( k ) − x3k ) ) / 1 ( ( x3k ) = (5 − 2 x1( k ) − 2 x2k ) ) / 1
x ( k +1) 0 −2 2 1 (k ) = − 1 0 − 1 x + 3 , k = 0,1,... − 2 − 2 0 5
∞
x ( 0 ) = (0,0,0)T , x (1) = (1,3,5)T , x (1) − x ( 0 )
选不同的M , 则可得不同的迭代法。
L
a11 a A = 21 an1 a12 L a22 L L L an2 L
D
U
0 0 a12 a13 L a1n a1n a11 a21 0 0 a23 L a2n a2n a22 = a a + + 0 M 0 M 31 32 O M M M O O an−1n ann ann a a a L 0 0 n1 n2 n3
xi
( k +1)
n 1 i −1 (k ) ( k ) bi = − ∑ aij x j + ∑ aij x j + , i = 1, 2,..., n aii j =1 j =i +1 aii
Component form
xi
( k +1)
1 i −1 ( k +1) =− + ∑ a ij x j a ii j =1
Jacobi iteration
取M为D
Matrix form
x ( k +1) = Bx ( k ) + f = − D −1 ( L + U ) x ( k ) + D −1b
x ( 0 )，初始向量
Jacobi迭代法简单，迭代一次只需作 一次矩阵与向量的乘
法即可。
x1 + 2 x2 − 2 x3 = 1 例： x1 + x2 + x3 = 3 ，取初始向量x ( 0 ) = (0,0,0)T , 2 x + 2 x + x = 5 2 3 1 写出Jacobi的矩阵迭代形式并求解。
x ( 0 ) = (0,0,0)T x (1) = (1,3,5)T x ( 2 ) = (5,−3,−3)T
Convenient in programming
( k +1)
k +1)
xi
= xi
(k )
1 + aii
n (k ) bi − ∑ aij x j , i = 1, 2,..., n j =1
Convergence of iteration
Error vector of iteration
B →0
k
Convergence of matrix
{
( A( k ) = aijk ) ∈ R n×n | k ∈ Z
( )
2
}
A = ( aij ) ∈ R n×n
example
10 −2 −1 x1 3 −2 10 −1x = 15 2 −1 −2 5 x3 10
( D + L) −1 难求 原式 ⇔ ( D + L) x ( k +1) = −Ux ( k ) + b ⇔ Dx ( k +1) = − Lx ( k +1) − Ux ( k ) + b
Gauss-Seidel iteration
( k +1)
Jacobi分量形式
xi
1 n ( k ) bi = − ∑ aij x j + , i = 1, 2,..., n aii j =1, j ≠i aii
x = (1
*
3)
t
Jacobi iteration
1 .9 9 9 8 1 .9 9 9 9 2 .9 9 9 8 )
t
x( x(
9) 5)
= ( 0 .9 9 9 8 = ( 0 .9 9 9 7
G-S iteration
t
2 .9 9 9 9 )
In this example, G-S iteration converges faster than Jacobi iteration. It’s not always true. Sometimes the later converges faster. For some linear systems, G-S converges while J diverges, and vice versa.