当前位置: 首页 > >

3.6线性方程组的迭代法

发布时间:

1/35

§3.6

迭代法

? 它的基本思想是将线性方程组 AX=b 化为 等价方程组 ? X=BX+f (k) ? 再由此构造向量序列{X }: ? X(k+1)=BX (k)+f (k) ? 若{X }收敛至某个向量X *,则向量X *就是 所求方程组 AX=b 的准确解. ? 线性方程组的迭代法主要有Jocobi迭代法、 Seidel迭代法和超松弛(Sor)迭代法.
计算方法三⑥

2/35

一、Jacobi 迭代法
Jacobi迭代和Seidel迭代 是求解线性方程组的两种基 本的迭代法,但由于收敛速 度较慢,已经越来越不适应 当前信息时代人们对计算速 度和精度的要求,所以在实 际应用中使用得并不多。但 是由于它们体现了迭代法的 最基本的思想,所以仍是学 *其它迭代法的基础。
计算方法三⑥

3/35

? x2 ?10 x1 ? 例:设AX=b为 ? ? x1 ? 10 x 2 ?? x ? x2 1 ? 化为等价于方程组: X=BX+f ? x1 ? ? ? x2 ? ?x ? ? 3 0 .1 x 2 0 . 1 x1 0 . 2 x1 ? 0 .2 x 2

? 2 x3 ? 2 x3 ? 5 x3

? 7 .2 ? 8 .3 ? 4 .2 ? 0 . 72 ? 0 . 83 ? 0 . 84

? 0 .2 x3 ? 0 .2 x3

?由此构造向量序列{X (k)}: X(k+1)=BX(k)+f

若X

(k+1)

→X*,则X*即为AX=b的解。

计算方法三⑥

4/35
( ( 取迭代初值 x1( 0 ) ? x 2 0 ) ? x 3 0 ) ? 0 , 依次代入迭代公式:X(k+1)=BX(k)+f

k

x1(k)

x2(k)

x3(k)

0

0.00000

0.00000
0.83000 1.07000 1.15710 1.18534 1.19510 1.19834 1.19981 1.19941 1.19994

0.00000
0.84000 1.15000 1.24820 1.28282 1029414 1.29504 1.29934 1.29978 1.29992

1 0.72000 2 0.97100 3 1.05700 4 1.08535 5 1.09510 6 1.09834 7 1.09944 8 1.09981 9 计算方法三⑥ 1.09994

5/35

X

(8 )

? 1 . 09981 ? ? ? 1 . 19941 ? 1 . 29978 ?

? ? ? ? ?

X

(9)

? 1 . 09984 ? ? ? 1 . 19994 ? 1 . 29992 ? ? 1 . 09984 ? ? ? 1 . 19994 ? 1 . 29992 ?

? ? ? ? ? ? ? ? ? ?

||X(9)-X(8)||≈0.0005<10-3 ∴原方程组的**馕 X
(9)

原方程组的精确解为:x1=1.1,x2=1.2,x3=1.3
计算方法三⑥

上一页

6/35

一般地, 设有方程组

a11x1+a12x2+·· a1nxn=b1 ·+ · a21x1+a22x2+·· a2nxn=b2 ·+ ·
..................

an1x1+an2x2+·· annxn=bn ·+ · 用矩阵表示为 AX =b (A 为系数矩阵,非奇异阵;b为右端 常数组成的向量,X为解向量) n个方程n个未知量
计算方法三⑥

? a 11 x 1 ? a 12 x 2 ? ? ? a 1 n x n ? b1 ? ? a 21 x 1 ? a 22 x 2 ? ? ? a 2 n x n ? b 2 ? ? ???????????? ? a n 1 x 1 ? a n 2 x 2 ? ? ? a nn x n ? b n ?

7/35

(1 )

把方程组写成容易迭代的形式: X=BX+f (1)设aii≠0(i=1,2,...,n)。将线性方程组的第i个方程的第i 个变元xi用其它n-1个变元表示,即解出xi(i=1,2,...,n)
1 ? ? a 12 x 2 ? a 13 x 3 ? ? a 1 n x n ) ? x 1 ? a ( b1 11 ? ? x ? 1 (b ? a x ? ? a 23 x 3 ? ? a 2 n x n ) 2 2 21 1 ? a 22 ? ???????????? ? 1 ( b n ? a n 1 x 1 ? a n 2 x 2 ? ? a nn ? 1 x n ? 1 ) ? xn ? ? a nn ?

计算方法三⑥

(2)写成迭代格式

8/35

1 ? ( k ?1) (k ) (k ) (k ) ? ( b1 ? a 12 x 2 ? a 13 x 3 ? ? a 1 n x n ) ? x1 a 11 ? (k ) (k ) ? x ( k ?1) ? 1 ( b ? a x ( k ) ? 2 ? a 23 x 3 ? ? a 2 n x n ) 2 21 1 ? a 22 ? ???????????? ? ( k ?1) 1 (k ) (k ) (k ) ? ( b n ? a n 1 x1 ? a n 2 x 2 ? ? ? a nn ? 1 x n ? 1 ) ? xn ? a nn ?

(3)取定初始向量X(0)=(x1(0),x2(0),...,xn(0))T,代入上 式计算,可逐次算出向量序列X(0), X(1) ,..., X(k),...,这 里X(k )=(x1(k),x2(k),...,xn(k)) T ,向量序列{X(k )}收敛时, 类似于非线性方程的迭代法,对于事先给定的精度 要求ε,当||X(k+1)-X(k)||<ε时,就可以得到方程组的* *釾*≈X(k+1) 计算方法三⑥

9/35

Jacobi迭代公式
xi
( k ? 1)

?

1 a ii

(bi ?

?a
j ?1 j?i

n

ij

x

(k ) j

)

i ? 1,2, ? , n


( k ? 1)

xi

?

bi a ii

?

?

n

a ij a ii

x

(k ) j

i ? 1,2, ? , n

j ?1 j?i

计算方法三⑥

10/35

例题 用Jacobi 迭代法解方程组AX=b,其中:
? 10 ? A ? ??1 ??1 ? ?1 ? 2 ? ? 10 ? 2 ? , ?1 5 ? ? ? 7.2 ? b ? ? 8.3 ? 4.2 ? ? ? ? ? ?

解: Jacobi 迭代公式

x1(k+1)=
x2(k+1)= 0.1x1(k)

0.1x2(k)+0.2x3(k)+0.72
+0.2x3(k)+0.83 +0.84

x3(k+1)= 0.2x1(k)+0.2x2(k)

取X(0)=(0,0,0)T,用五位小数计算,结果如下表
计算方法三⑥

11/35

k
0 1 2 3 4 5 6 7 8 9
计算方法三⑥

x1(k)
0.00000 0.72000 0.97100 1.05700 1.08535 1.09510 1.09834 1.09944 1.09981 1.09994

x2(k)
0.00000 0.83000 1.07000 1.15710 1.18534 1.19510 1.19834 1.19981 1.19941 1.19994

x3(k)
0.00000 0.84000 1.15000 1.24820 1.28282 1029414 1.29504 1.29934 1.29978 1.29992

12/35

X

(8 )

? 1 . 09981 ? ? ? 1 . 19941 ? 1 . 29978 ?

? ? ? ? ?

X

(9)

? 1 . 09984 ? ? ? 1 . 19994 ? 1 . 29992 ? ? 1 . 09984 ? ? ? 1 . 19994 ? 1 . 29992 ?

? ? ? ? ? ? ? ? ? ?

||X(9)-X(8)||≈0.0005<10-3
∴原方程组的**馕 X
(9)

原方程组的精确解为:x1=1.1,x2=1.2,x3=1.3
计算方法三⑥

13/35

Jacobi 迭代的矩阵格式 X
( k ?1)

? BJ X

(k )

? f

推导如下:Ax=b 将A作如下分解 ~ ~ 设A=(aij) A=L+D+U
? 0 ? ? a 21 ? ? ? ?a ? n1 0 0 ? an2 ? ? ? ? 0? ? 0? ?? ? 0? ?
? a 11 ? ? 0 ? ? ? ? 0 ? 0 a 22 ? 0 ? ? ? ? 0 ? ? 0 ? ? ? ? a nn ? ?

~~ L,U的主对角线上元素均为零.
计算方法三⑥

?0 ? ?0 ?? ? ?0 ?

a 12 0 ? 0

? ? ? ?

a1n ? ? a2n ? ? ? ? 0 ? ?

~ ~ ~ ~ AX=b即为:(L+D+U)X=b, 即 DX= (-L-U) X+b ~ ~ -1,得:X= D-1(-L-U)X+D-1b 两边同乘D 故得雅可比迭代公式的矩阵表示式: ~ ~ (k+1)=D-1(- L-U)X(k) + D-1b X ~ ~ -1(-L-U),f =D-1b,则有: 记BJ = D J X(k+1)=BJX(k) +fJ

14/35

BJ称为雅可比方法 的迭代矩阵

Jacobi 迭代的矩阵格式 X
( k ?1)

计算方法三⑥

? BJ X

(k )

? fJ

15/35

例题

用Jacobi 迭代法解方程组AX=b,其中:
? 10 ? A ? ??1 ??1 ? ?1 ? 2 ? ? 10 ? 2 ? , ?1 5 ? ? ? 7.2 ? b ? ? 8.3 ? 4.2 ? ? ? ? ? ?

X

~ ~ A=L+D+U ~ ~ (k+1) -1

=D (- L-U)X(k) + D-1b
? 1 . 09984 ? ? ? 1 . 19994 ? 1 . 29992 ? ? ? ? ? ?

||X(9)-X(8)||=0.0005<10-3 ∴原方程组的**馕

X

(9)

计算方法三⑥

16/35

二、Gauss-Seidel迭代法
为了加快收敛速度,同时为了节省计算机的内存, 我们作如下的改进:每算出一个分量的*似值,立即 用到下一个分量的计算中去,这样所得的迭代法就称 为Gauss-Seidel迭代法,也称为“异步迭代法”,简 称为GS迭代法.
? a 11 x 1 ? a 12 x 2 ? ? ? a 1 n x n ? b1 ? ? a 21 x 1 ? a 22 x 2 ? ? ? a 2 n x n ? b 2 ? ? ???????????? ? a n 1 x 1 ? a n 2 x 2 ? ? ? a nn x n ? b n ?

(1 )

把方程组写成容易迭代的形式: X=BX+f
计算方法三⑥

(1)设aii≠0(i=1,2,...,n)。将线性方程组的第i个方程的第i个变元xi 用其它n-1个变元表示,即解出xi(i=1,2,...,n)
1 ? x1 ? ( b1 ? a12 x 2 ? a 13 x 3 ? ? a1 n x n ) ? a11 ? ? x ? 1 (b ? a x ? a 23 x 3 ? ? a 2 n x n ) 2 2 21 1 ? a 22 ? ???????????? ? 1 xn ? ( b n ? a n 1 x1 ? a n 2 x 2 ? ? a nn ?1 x n ?1 ) ? a nn ?

17/35

(2)Seidel迭代格式
x
( k ?1) 1

?
?

1 a 11 1
a 22 1

( b1

? a 12 x 2
( k ?1)

(k )

? a 13 x 3 ? ? a 1 n x n )
(k ) (k )

x2

( k ?1)

( b 2 ? a 21 x1

? a 23 x 3 ? ? a 2 n x n )
(k ) (k )

x n ?1
x

( k ?1)

……………………
? ( b n ? a n 1 x1
( k ?1)

a nn 计算方法三⑥

( k ?1) n

? an2 x2

( k ?1)

? ? a nn ? 1 x n ? 1

( k ?1)

)

? ( k ?1) (k ) (k ) (k ) x1 ? ( b1 ? a 12 x 2 ? a 13 x 3 ? ? a 1 n x n ) ? a 11 ? Jacobil (k ) (k ) ? x ( k ?1) ? 1 ( b ? a x ( k ) ? 2 ? a 23 x 3 ? ? a 2 n x n ) 2 21 1 迭代格式 ? a 22 ? ???????????? ? ( k ?1) 1 (k ) (k ) (k ) xn ? ( b n ? a n 1 x1 ? a n 2 x 2 ? ? ? a nn ? 1 x n ? 1 ) ? ? a nn ? 1

18/35

Seidel迭代格式
1 ? ( k ?1) (k ) (k ) (k ) x1 ? ( b1 ? a 12 x 2 ? a 13 x 3 ? ? a 1 n x n ) ? a 11 ? (k ) (k ) ? x ( k ?1) ? 1 ( b ? a x ( k ?1) ? ? a 23 x 3 ? ? a 2 n x n ) 2 2 21 1 ? a 22 ? ???????????? ? ( k ?1) 1 ( k ?1) ( k ?1) ( k ?1) ? ( b n ? a n 1 x1 ? an2 x2 ? ? ? a nn ? 1 x n ? 1 ) ? xn ? a nn ?

计算方法三⑥

19/35

1 ? ( k ?1) (k ) (k ) (k ) x1 ? ( b1 ? a 12 x 2 ? a 13 x 3 ? ? a 1 n x n ) ? a 11 ? (k ) (k ) ? x ( k ?1) ? 1 ( b ? a x ( k ?1) ? ? a 23 x 3 ? ? a 2 n x n ) 2 2 21 1 ? a 22 ? ???????????? ? ( k ?1) 1 ( k ?1) ( k ?1) ( k ?1) xn ? ( b n ? a n 1 x1 ? an2 x2 ? ? ? a nn ? 1 x n ? 1 ) ? ? a nn ?

Seidel迭代公式

xi

( k ?1)

?

1 a ii

( bi ?

?a
j ?1

i ?1

ij

xj

( k ?1)

?

j ? i ?1

?a

n

ij

xj

(k )

)

(i=1,2,...,n)
计算方法三⑥

20/35 Seidel迭代法的矩阵表示形式 ~ ~ ~ ~ 利用AX=b 及A=L+D+U,其中D为对角矩阵,L,U分 别为严格下,上三角矩阵.则有,GS迭代法的矩 ~ ~ 阵形式为:

( L ? D ) X ? ?U X ? b

~ ~ ?1 ~ ?1 X ? ?(L ? D ) UX ? (L ? D ) b
X
( k ?1)

记Bs= (-L-D) U,fs=(L+D)-1b,则有:

~ ?1 ~ ? ?(L ? D ) UX ~ ~ -1 ~
( k ?1)

(k )

~ ?1 ? (L ? D ) b

X

? Bs X

(k )

? fs

Bs称为高斯—塞德 尔方法的迭代矩阵 计算方法三⑥

21/35

例题

用Seidel 迭代法解方程组AX=b,其中:
? 10 ? A ? ??1 ??1 ? ?1 ? 2 ? ? 10 ? 2 ? , ?1 5 ? ? ? 7.2 ? b ? ? 8.3 ? 4.2 ? ? ? ? ? ?

解: Seidel 迭 代 公 式 Jacobi 迭 代 公 式

x1(k+1)=

0.1x2(k)+0.2x3(k)+0.72

x2(k+1)= 0.1x1(k+1)
x1(k+1)=

+0.2x3(k)+0.83
+0.84

x3(k+1)= 0.2x1(k+1)+0.2x2(k+1)

0.1x2(k)+0.2x3(k)+0.72

x2(k+1)= 0.1x1(k)

+0.2x3(k)+0.83
+0.84

x3(k+1)= 0.2x1(k)+0.2x2(k)

计算方法三⑥

22/35

k
0 1 2 3

x1(k)
0.00000 0.72000 1.04308 1.09313

x2(k)
0.00000 0.90200 1.16719 1.19572

x3(k)
0.00000 1.6440 1.282.5 1.29778

4 5
6

1.09913 1.09989
1.09999

1.19947 1.19993
1.19999

1.29972 1.29996
1.30000

Seidel迭代公式比Jacobi迭代公式得到的结果好。 一般来说,前者比后者收敛速度更快。
计算方法三⑥

三、 逐次超松驰迭代法(SOR法)

23/35

Jacobi迭代和Seidel迭代对具体的线性方程组来说, *媸到獾乃俣仁枪潭ú槐涞模粲龅绞樟埠苈 的情况怎样对其做加速处理呢?超松驰迭代法就是 迭代的一种加速方法。它可以看成是Seidel迭代法 的加速,而Seidel迭代是它的一个特例。 设有方程组 AX=b 其中A=(aij)n为非奇异阵, X=(x1,x2,...,xn)T,b=(b1,b2,...,b n)T 记X(k)为第k步迭代的*似值,先用Gauss-Seidel迭代公 式计算一步,并将计算结果记为

~ ( k ?1) X

~ ( k ?1) 并记: ? X ? X ? X
计算方法三⑥

(k )

24/35

然后用某个参数ω乘以△X,再将其加到X(k)上做新的
**釾(k+1) , 即令 X
( k ?1) (k )

( k ?1 )

? X

(k )

? ??X

~ ( k ?1) (k ) X ? X ? ?(X ? X ) ~ ( k ?1) ( k ?1) (k ) X ? (1 ? ? ) X ? ?X

△X
SOR迭代 公式

逐次超松弛迭代公式的分量形式如下: ω称为松驰因子.
? ? ( k ?1) (k ) (k ) (k ) (k ) x1 ? (1 ? ? ) x1 ? ( b1 ? a 12 x 2 ? a 13 x 3 ? ? a 1 n x n ) ? a 11 ? (k ) (k ) ? x ( k ? 1 ) ? (1 ? ? ) x ( k ) ? ? ( b ? a x ( k ? 1 ) ? 2 ? a 23 x 3 ? ? a 2 n x n ) 2 2 21 1 ? a 22 ? ?? ? ? ( k ?1) (k ) ( k ?1) ( k ?1) ( k ?1) ? (1 ? ? ) x n ? ( b n ? a n 1 x1 ? an2 x2 ? ? a nn ? 1 x n ? 1 ) ? xn ? a nn ? 计算方法三⑥

25/35

? ? ( k ?1) (k ) (k ) (k ) (k ) x1 ? (1 ? ? ) x1 ? ( b1 ? a 12 x 2 ? a 13 x 3 ? ? a 1 n x n ) ? a 11 ? (k ) (k ) ? x ( k ? 1 ) ? (1 ? ? ) x ( k ) ? ? ( b ? a x ( k ? 1 ) ? 2 ? a 23 x 3 ? ? a 2 n x n ) 2 2 21 1 ? a 22 ? ?? ? ? ( k ?1) (k ) ( k ?1) ( k ?1) ( k ?1) xn ? (1 ? ? ) x n ? ( b n ? a n 1 x1 ? an2 x2 ? ? a nn ? 1 x n ? 1 ) ? ? a nn ?

SOR迭代 公式
xi
( k ?1)

? (1 ? ? ) x i

(k )

?

?
a ii

( bi ?

?

i ?1

a ij x j

( k ?1)

?

j ?1

j ? i ?1

?

n

a ij x j

(k )

)

i=1,2,…,n
计算方法三⑥

xi

( k ?1)

? (1 ? ? ) x

(k ) i

?

?
a ii
i ?1

( bi ?

?a
j ?1 ( k ?1)

i ?1

ij

xj

( k ?1)

?

j ? i ?1 n

?a

n

26/35
ij

xj

(k )

)

xi

( k ?1)

? x

(k ) i

?

?
a ii

( bi ?

?a
j ?1

ij

xj

? a ii x

(k ) i

?

j ? i ?1

?a

ij

xj

(k )

)

迭代公式也可写为:
? ? ( k ?1) (k ) (k ) x1 ? x1 ? ( ? a 11 x1 ? a 11 ? ? x ( k ?1) ? x ( k ) ? ? ( ? a x ( k ?1) ? 2 2 21 1 ? a 22 ? ? ? ( k ?1) (k ) ( k ?1) ? xn ? ( ? a n 1 x1 ? ? xn ? a nn ?
? a 12 x 2
(k ) (k )

i=1,2,…,n
? a 13 x 3 ? ? a 1 n x n
(k ) (k )

? b1 ) ? b2 )

? a 22 x 2 ?? an2 x2

(k )

? a 23 x 3 ? ? a 2 n x n

(k )

( k ?1)

? ? a nn ? 1 x n ? 1 ? a nn x n ? b n )
(k )

( k ?1)

计算方法三⑥

27/35 ?( k ) ? ( k ) ? ( k ?1) (k ) (k ) (k ) (k ) (k ) (k ) (k ) x1 ? ? ? (1 ? ? ) x1 ( ? a 11 x1 1 ? a 12 xa 12 x 2 a 13 xa 13 x 3 ?? 1?xa 1 n x n b1)) (b ? 2 ? ? 3 ? a n n ? ? x1 a 11 a 11 ? ? ( ) k k ( 1) ( ? k ? xx( (kk?? 1 ) ? (1 (? )??) x?k ) ( ? a x bk ? 1? a ax ( k x1()k ) ? a ? (a ) ?( k ) ? ? ak ) ? (b ) ) ? 22 ? x2 x 3 23 x 3 ? a 2 n x n 2 n x n 2 ? 21 (1 2 ? 21 22 2 23 2 1 ? a 22 a 22 ? ?? ?? ? ? ( k ?1) ?( k ) ( k ?1) ( k (? 1? 1 ) ( k ?1) (k( k ) ?1) (k ) k ) ( k ?1) ? (1n ? ? ) x n (? a n 1 x1 b n ??aa1nx1x 2 ??n ?xa nn ? 1 ??? a? a1nn n ?n1 ? b n)) ( a 2 2 x ? ? x n 1 nn ? x x ? xn n 2 a nn ? a nn ?

当ω=1时,上式就是GS迭代公式。

当ω<1时称为低松驰迭代法, ω称为松驰因子.
当ω>1时称为高松驰迭代法。 在SOR迭代中,松弛因子ω的取值对迭代公 式的收敛速度影响极大。实际计算时,要根据方 程的系数矩阵的性质,或实际经验来选取ω。 计算方法三⑥

28/35

例:用SOR迭代法解方程组:? ? 4 ?
? 1 ? 1 ? ? 1 ?

1 ?4 1 1

1 1 ?4 1

1 ? ? x1 ? ? 1 ? ?? ? ? ? 1 ?? x 2 ? ?1 ? ?? x ? ? ?1 ? 1 ?? 3 ? ? ? ? 4 ?? x 4 ? ?1 ? ?? ? ? ?

迭代公式为:
? ? ( k ?1) (k ) (k ) x1 ? x1 ? ( ? a 11 x1 ? a 11 ? ? x ( k ?1) ? x ( k ) ? ? ( ? a x ( k ?1) ? 2 2 21 1 ? a 22 ? ? ? ( k ?1) (k ) ( k ?1) ? xn ? ( ? a n 1 x1 ? ? xn ? a nn ?

(它的精解为x*=(-1,-1,-1,-1)T)
? a 12 x 2 ? a 13 x 3 ? ? a 1 n x n
(k ) (k )

(k )

(k )

? b1 ) ? b2 )

? a 22 x 2 ?? an2 x2

(k )

? a 23 x 3 ? ? a 2 n x n

(k )

( k ?1)

? ? a nn ? 1 x n ? 1 ? a nn x n ? b n )
(k )

( k ?1)

计算方法三⑥

29/35

解:取X(0)=0,迭代公式为
? ? ( k ?1) (k ) (k ) (k ) (k ) (k ) x1 ? x1 ? ( 4 x 1 ? x 2 ? x 3 ? x 4 ? 1) ? 4 ? ? ( k ?1) (k ) ( k ?1) (k ) (k ) (k ) ? x2 ? x2 ? ( ? x1 ? 4 x 2 ? x 3 ? x 4 ? 1) ? 4 ? ? ( ( ( k ?1) ( k ?1) (k ) (k ) ? x 3 k ?1) ? x 3 k ) ? ( ? x1 ? x2 ? 4 x 3 ? x 4 ? 1) 4 ? ? ? ( k ?1) (k ) ( k ?1) ( k ?1) ( k ?1) (k ) x4 ? x4 ? ( ? x1 ? x2 ? x3 ? 4 x 4 ? 1) ? 4 ?

取ω=1.3,迭代11次后结果为 X(11)=(-1.0000,-1.0000,-1.0000,-1.0000)T 满足精度要求:||X(k)-X*||2<10-5
计算方法三⑥

30/35

ω若取其他值,迭代次数如下表:
松驰因子ω 迭代次数(<10-5) 22 17 13 迭代次数(<10-6) 26 20 15

1.0
1.1 1.2 1.3 1.4 1.5 1.6
计算方法三⑥

最少迭代次数

12
15 19 25

14
18 23 31

ω=1. 3 是最佳松驰因子

? ? ( k ?1) (k ) (k ) (k ) (k ) x1 ? (1 ? ? ) x1 ? ( ? a 12 x 2 ? a 13 x 3 ? ? a 1 n x n ? b1 ) ? a 11 ? (k ) (k ) ? x ( k ? 1 ) ? (1 ? ? ) x ( k ) ? ? ( ? a x ( k ? 1 ) ? 2 ? a 23 x 3 ? ? a 2 n x n ? b 2 ) 2 21 1 ? a 22 ? ?? ? ? ( k ?1) (k ) ( k ?1) ( k ?1) ( k ?1) ? (1 ? ? ) x n ? ( ? a n 1 x1 ? an2 x2 ? ? a nn ? 1 x n ? 1 ? bn ) ? xn ? a nn ? ~

SOR迭代公式的矩阵表示

31/35

~ 将公式写成矩阵形式如下: A=L+D+U ~ ( k ?1) ~ ( k ) ( k ?1) (k ) ?1 X ? (1 ? ? ) X ? ? D (? L X ? UX ? b) ~ ( k ?1) ~ (k ) ( k ?1) (k ) DX ? (1 ? ? ) DX ? ? (? L X ? UX ? b)
~ (D ? ?L ) X
计算方法三⑥
( k ?1)

~ ? [( 1 ? ? ) D ? ? U ] X

(k )

? ?b

32/35

~ (D ? ?L ) X
X
( k ?1)

( k ?1)

~ ? [( 1 ? ? ) D ? ? U ] X
(k )

(k )

? ?b

~ ?1 ~ ? ( D ? ? L ) [( 1 ? ? ) D ? ? U ] X
?

~ ?1 ? ? (D ? ?L ) b
?

S?
SOR 迭代的矩阵格式 X
( k ?1 )

f
(k )

?

? S? X

? f?

X

( k ?1)

? BX

(k )

? f

计算方法三⑥

33/35

四. 相关程序设计

Jacobi 迭代的矩阵格式

( k ?1) (k ) 1)给出矩阵A,b,x0 X ? BJ X ? fJ ~ ~ -1(-L-U), f =D-1b 2)计算B= D J ~ ~ (k+1)= D-1(-L-U)X(k) + D-1b X 3) X1=BX 0+f ~ ~ -1(-L-U), 记BJ = D 4)若 ||X1-X0||<=1.0e-6,则转到 fJ=D-1b, 5),否则,转到3) 则有: X(k+1)=BJX(k) +fJ 5)输出X1和n。

y=X1
计算方法三⑥

作业:P74

11、13、14

开始

输入a,b,x0 由矩阵a,生成D,U,L,B,f y=B*x0+f n=1

34/35

Jacobi 迭代的矩阵格式 X
( k ?1)

? BJ X

(k )

? fJ

X(k+1)=

~ ~ -1(-L-U)X(k) + D ~ ~ -1(-L-U), D

D-1b

N

||y-x0||>=1.0e-6

记BJ =

Y
x0=y y=B*x0+f

f=D-1b,则有:
X(k+1)=BJX(k) +f

N
n=n+1 输出y 结束 计算方法三⑥

MATLAB实现Jacobi迭代法

35/35

function y=Jacobi(a,b,x0) D=diag(diag(a)); %生成对角矩阵 U=triu(a,1); %生成严格上三角矩阵 L=tril(a,-1); %生成严格下三角矩阵 B=inv(D)*(-L-U); % B=D\(-L-U); ~ ~ f=inv(D)*b; % f=D\b; (k+1)= D-1(-L-U)X(k) + D-1b X y=B*x0+f; ~ ~ n=1; -1(-L-U), 记BJ = D while norm(y-x0)>=1.0e-6 x0=y; f=D-1b,则有: y=B*x0+f; n=n+1; X(k+1)=BJX(k) +f end y n 计算方法三⑥

MATLAB实现Jacobi迭代法

36/35

function y=Jacobi(a,b,x0) D=diag(diag(a)); %生成对角矩阵 U=-triu(a,1); %生成严格上三角矩阵 L=-tril(a,-1); %生成严格下三角矩阵 B=inv(D)*(L+U); % B=D\(L+U); ~ ~ f=inv(D)*b; % f=D\b; (k+1)= D-1(-L-U)X(k) + D-1b X y=B*x0+f; ~ ~ n=1; -1(-L-U), 记BJ = D while norm(y-x0)>=1.0e-6 x0=y; f=D-1b,则有: y=B*x0+f; n=n+1; X(k+1)=BJX(k) +f end y n 计算方法三⑥

MATLAB实现seidel迭代法

原公式

37/35

function y=Seidel(a,b,x0) ( k ? 1 ) ~ ?1 ~ ~ ?1 (k ) X ? ( D ? L ) ( ?U ) X ? (D ? L ) b D=diag(diag(a)); U=-triu(a,1); ( k ?1) ?1 (k ) ?1 L=tril(a,-1); X ? ( D ? L ) UX ? (D ? L) b B=(D+L)\U; 记Bs= (D+L)-1U, f=(D+L)\b; y=B*x0+f; n=1; fs=(D+L)-1b,则有: while norm(y-x0)>=1.0e-6 x0=y; ( k ?1) (k ) y=B*x0+f; X ? Bs X ? fs n=n+1; ~ ~ end ( D ? L ) X ? ?U X ? b y ~ ~ (k ) ( k ?1) n (D ? L ) X ? ?U X ?b 计算方法三⑥

38/35

用Matlab求解线性方程组示例
用Jacobi方法求解下列方程组,设x(0)=0,精度为10-6
10x1-x2 =9

-x1+10x2-2x3 =7 -2x1 +10x3=6 结果:y = 0.9959 0.9594 0.7992 n =10

MATLAB实现 >>a=[10 -1 0;-1 10 -2;-2 0 10];

>>b=[9;7;6];
>>Jacobi(a,b,[0;0;0])
计算方法三⑥

39/35

>> Jacobi(a,b,[0;0;0]) y= 0.9959 0.9594 0.7992 n = 10

>> Jacobi(a,b,[1;1;1]) y= 0.9959 0.9594 0.7992 n= 8

计算方法三⑥

40/35

求解方程组
用seidel方法求解下列方程组,设x(0)=0,精度为10-6 10x1-x2 =9

-x1+10x2-2x3 =7
-2x1 +10x3=6 结果:y = 0.9959 0.9594 0.7992

MATLAB实现 >>a=[10 -1 0;-1 10 -2;-2 0 10]; >>b=[9;7;6]; >>Seidel(a,b,[0;0;0])
计算方法三⑥

n= 7

41/35

>>Seidel(a,b,[0;0;0])

>>Seidel(a,b,[1;1;1])
结果:y = 0.9959 0.9594 0.7992

结果:y = 0.9959 0.9594 0.7992
n= 7

n= 6

计算方法三⑥

42/35

作业:P74 作业:P75
实验 令
? 25 ? ? ? 41 A?? 10 ? ? ?6 ? ? 41 68 ? 17 10

11、13、14 15、16、17
10 ? 17 5 ?3 ? 6? ? 10 ? ? 3? ? 2 ? ? ? 32 ? ? ? ? 23 ? b ?? 33 ? ? ? ? 31 ? ? ?

1 、给出Jacobi迭代、Seidel迭代与超松驰迭代的算 法流程和程序,并求方程组AX=b的**狻 2 、给出超松驰迭代的算法流程和程序,并求方程 组AX=b的**饧白罴阉沙垡蜃印
计算方法三⑥

迭代法的特点

43/35

? 方法简单,每次迭代都是简单的重复运算,易 于编制程序;与求解线性方程的精确法相比, Jacoai迭代法对于字长位数较少的计算机更为 适用,它可以用增加迭代次数来弥补字长位数 少的不足. ? 初值可以任取,因而中间结果偶然错误不影响 最后结果的获得。慢。用计算机计算时,。 ? 就其收敛性而言,某些用Jacoai 迭代法不能收 敛,而无法得出结果的线性代数方程组,用 Seidel迭代法却能进行收敛计算,反之亦然.
计算方法三⑥

五、收敛性及误差估计
Jacobi迭代和GS迭代格式可表述为统一形式:
X
( k ?1)

44/35

? BX

(k )

? f

定理1:对任意初始向量X(0)及任意向量 f,由此产生的迭 代向量序列{X(k)}收敛的充要条件是矩阵B的谱半径
? ?B ? ? 1

证明:设X* 为方程的准确解,则有: 又∵
? X
(k )

X

*

? BX

*

? f

X
*

(k )

? BX

( k ?1 )

? f
*

? X

? B X
(K )

?

( k ?1 )

? X

?

?? ? B
k

k

?X

(0)

? X

*

?

两边取极 lim X 限,得: k ? ?
计算方法三⑥

? X * ? lim B ? 0
k??

? ρ?B ? ? 1

注:用上述定理证明迭代公式的收敛性只要计算B的 谱半径ρ(B),看它是否<1? 例:考察下列迭代公式的敛散性
X
( k ? 1)

45/35

? BX

(k )

? f

其中

? 0 ? B ? ? ? 4/11 ? ? 6/12 ?

3/8 0 ? 3/12

? 2/8 ? ? 1/11 ? 0 ? ?

解:计算B的特征值
Det(λE-B)=λ3+0.034090909λ+0.039772727=0
λ1=-0.3082 λ2=-0.1541+0.3245i λ3=- 0.1541-0.3245i 计算方法三⑥ |λ1|=0.3082<1 |λ2|=0.3592<1 |λ3|= 0.3592<1 即ρ(B)<1

所以上述迭代公式对任意的初始值都是收敛的。

46/35

定理2 若||B||<1,则迭代公式X(k+1)=BX(k)+f 收敛. ρ(B) ? ||B||
推论1 (1)
B

若B=(bij)n×n 满足下列条件之一:
?

? max
i

?

n

b ij

? 1;

j ?1

(2)

B

1

?

max ?
j
n

n

b ij ? 1
n

i ?1

0 ? ? B ? ? ? 4 / 11 ? ? 6 / 12 ?
2

3/8 0 ? 3 / 12

? 2 /8? ? 1 / 11 ? 0 ? ?

(3)

B

2

?

? ?
i?1

b ij

? 1
(k )

||B||<1 收敛。

j?1

则迭代公式 X
计算方法三⑥

( k ? 1)

? BX

? f

47/35

推论2:设AX=b,A=L+D+U,其中A为非奇异矩阵, D也为非奇异矩阵, 则: (1) Jacobi迭代法收敛的 ? ? ( B J ) ? 1
BJ = ~ -1(-L-U) D ~

~

~

(2) Seidel迭代法收敛的 ? ? ( B S ) ? 1
~ -1 ~ BS = (D+L) (-U)

(3) SOR迭代法收敛的

? ? (S? ) ? 1 ~ -1 ~ Sw = (D+wL) [(1-w)D-wU]

计算方法三⑥

48/35

推论3:如果线性方程组AX=b的系数矩阵A为严 格对角占优阵,则其 Jacobi迭代和Seidel迭代对 任何初始向量X(0)都收敛。
定理3. 设解方程组AX=b的SOR迭代法收敛,

计算方法三⑥

49/35

定理3. 设解方程组AX=b的SOR迭代法收敛, 则0<ω<2。 ( k ?1) (k ) 证: SOR 迭代的矩阵格式 X ? S X ? f
? ?



~ ?1 ~ ? ( D ? ω L ) [(1 ? ω ) D ? ω U ]

设 S ω 的特征值为 λ1 , λ 2 , λ 3 , ? , λ n 则
det ? S ?

?

? λ1 λ 2 ? λ n

? ?ρ?S?

??

n

? det ? S ?

?n
n

1

? ρ?S ω ? ? 1

~ ?1 ~ det ? S ω ? ? det ( D ? ω L ) ? det[(1 ? ω ) D ? ω U ]
?

?

?

?
i ?1

n

1 a ii

? ? ?(1 ? ω ) a ii ? ? (1 ? ω )
i ?1

n

计算方法三⑥

50/35

det S ω ?

? ? ?
i ?1

n

1 a ii

? ? ?(1 ? ω ) a ii ? ? (1 ? ω )
i ?1

n

n

? det ? S ?

?n

1

? 1 - ? ? ρ?S? ? ? 1

即 1-? ? 1

故有:0<ω<2 .

定理3. 设解方程组AX=b的SOR迭代法收敛, 则0<ω<2。 反之不成立

定理4. 设方程组AX=b,如果A为对称正定矩阵,
且0<ω<2,则SOR迭代法收敛。
计算方法三⑥

51/35

如果 矩阵 A=(aij)满足
a ii ?
j ?1, j ? i

?

n

a ij

i ? 1, 2 , ? n

则称方阵A是严格(行)对角占优阵.
?? 4 ? A ?? 1 ? 2 ? 2 ?9 ?6 1 ? ? 7 ? 10 ? ?

例 矩阵

——严格(行)对角占优阵
计算方法三⑥

52/35

U

A=

a11 a21 … an1

a12 a13 … a1n a22 a23 … a2n … … … … =L+D+U an3 an4 … ann
D L

D是满秩的 D+L是满秩的 D+U是满秩的

严格(行)对角占优阵

计算方法三⑥

53/35

举例检验Jacoai迭代的收敛性

三阶线性方程组:

8x1 -3x2 +2x3 =20 4x1 +11x2 -x3 =33 6x1 +3x2 +12x3 =36
?8 ? A ? ?4 ?6 ? ?3 11 3 2 ? ? ? 1? 12 ? ?

解:该方程组的系数矩阵为

∵ 8>|-3|+2 ,11>4+|-1| ,12>6+3 i ∴A为严格对角占优阵

故该方程组采用Jacoai迭代法计算是收敛的。
计算方法三⑥

例:设方程组为

? 5 x 1 ? 2 x 2 ? x 3 ? ? 12 ? ? ? x 1 ? 4 x 2 ? 2 x 3 ? 20 ? 2 x ? 3 x ? 10 x ? 3 1 2 3 ?

54/35

试分别写出其Jacobi迭代格式和Seidel迭代格式以及 (k ) (k ) 相应的迭代矩阵。 ? x ( k ? 1) ? 2 1 ? x ? x ? 解:Jacobi迭代格式为
? ( k ? 1) (k) x2 ? 1 x1 ? 4 ? x ( k ? 1) ? ? 1 x ( k ) ? 1 5 ? 3
1 5 2 5 3

12 5

?
3 10

1 2

x3

(k )

?5 ?
3 10

x2

(k )

故Jacobi迭代矩阵为

BJ =

? 0 ? 1 ? 4 ?? 1 ? 5

? 0
3 10

2 5

? 1? 5 ? ? 1 2? 0 ? ?

BJ的主对角线上元素均为0
计算方法三⑥

Seidel迭代格式为
? 从式中解出xik+1 ? .? ,i=1,2,3得: ? ?

? x1 ? ? x2 ? x3 ? 12 5 ? ( k ? 1) (k ? 1) (k ) 1 1 ? 4 x1 ? 2 x3 ? 5 ? x2 ? x ( k ? 1) ? ? 1 x ( k ? 1) ? 3 x ( k ? 1) 3 ? 10 1 2 5 10 ? 3
2 5 (k ) 1 5 (k )

( k ? 1)

55/35

x1 x2 x3

(k ? 1)

? ? ?

? ?

2 5 1 10

x2

(k) (k)

? ? ?
?

1 5

x3
11 20

(k)

?
(k)

12 5

(k ? 1) (k ? 1)

x2 x2
?

x3
(k)

?
21 10

22 5

1 20

(k)

1 8
1 5

x3

?

故可得Seidel迭代矩阵为

?0 ? Bs ? 0 ? ?0 ?

2 5 1 10

?

1 20

? ? ? 11 20 ? ? 1 ? 8 ?

从例中可以看出Jacobi迭代矩阵Bj的主对角线为零, 而Seidel迭代矩阵Bs的第1列都是零,这对一般情况也 是成立的。
计算方法三⑥




友情链接: