Jacobi迭代法与Gauss-Seidel迭代法回顾
上一篇博客介绍了J迭代法和G-S迭代法的原理,现在复习一下,J迭代法的迭代公式为
x
i
(
k
+
1
)
=
1
a
i
i
(
b
i
−
∑
j
=
1
i
−
1
a
i
j
x
j
(
k
)
−
∑
j
=
i
+
1
n
a
i
j
x
j
(
k
)
)
i
=
1
,
2
,
⋯
,
n
,
k
=
0
,
1
,
⋯
x_i^{(k+1)}=\frac1{a_{ii}}(b_i-\sum_{j=1}^{i-1}a_{ij}x_j^{(k)}-\sum_{j=i+1}^na_{ij}x_j^{(k)})\\i=1,2,\cdots,n,k=0,1,\cdots
xi(k+1)=aii1(bi−j=1∑i−1aijxj(k)−j=i+1∑naijxj(k))i=1,2,⋯,n,k=0,1,⋯矩阵形式为
x
(
k
+
1
)
=
D
−
1
(
L
+
U
)
x
(
k
)
+
D
−
1
b
x^{(k+1)}=D^{-1}(L+U)x^{(k)}+D^{-1}b
x(k+1)=D−1(L+U)x(k)+D−1b(这里将系数矩阵A分裂为A = D - L - U)
G-S迭代法的迭代公式为
x
i
(
k
+
1
)
=
1
a
i
i
(
b
i
−
∑
j
=
1
i
−
1
a
i
j
x
j
(
k
+
1
)
−
∑
j
=
i
+
1
n
a
i
j
x
j
(
k
)
)
i
=
1
,
2
,
⋯
,
n
,
k
=
0
,
1
,
⋯
x_i^{(k+1)}=\frac1{a_{ii}}(b_i-\sum_{j=1}^{i-1}a_{ij}x_j^{(k+1)}-\sum_{j=i+1}^na_{ij}x_j^{(k)})\\i=1,2,\cdots,n,k=0,1,\cdots
xi(k+1)=aii1(bi−j=1∑i−1aijxj(k+1)−j=i+1∑naijxj(k))i=1,2,⋯,n,k=0,1,⋯矩阵形式为
x
(
k
+
1
)
=
(
D
−
L
)
−
1
U
x
(
k
)
+
(
D
−
L
)
−
1
b
x^{(k+1)}= (D - L)^{-1}Ux^{(k)}+(D-L)^{-1}b
x(k+1)=(D−L)−1Ux(k)+(D−L)−1b(仍将系数矩阵A分裂为A = D - L - U)
SOR法(逐次超松弛迭代法)
超松弛迭代法可视为J迭代法或G-S迭代法的改进版,我们可以把J迭代法的迭代公式做如下变换
x
i
(
k
+
1
)
=
1
a
i
i
(
b
i
−
∑
j
=
1
i
−
1
a
i
j
x
j
(
k
)
−
∑
j
=
i
+
1
n
a
i
j
x
j
(
k
)
)
=
x
i
(
k
)
+
1
a
i
i
(
b
i
−
∑
j
=
1
i
−
1
a
i
j
x
j
(
k
)
−
∑
j
=
i
n
a
i
j
x
j
(
k
)
)
x_i^{(k+1)}=\frac1{a_{ii}}(b_i-\sum_{j=1}^{i-1}a_{ij}x_j^{(k)}-\sum_{j=i+1}^na_{ij}x_j^{(k)})\\ =x_i^{(k)} + \frac1{a_{ii}}(b_i-\sum_{j=1}^{i-1}a_{ij}x_j^{(k)}-\sum_{j=i}^na_{ij}x_j^{(k)})
xi(k+1)=aii1(bi−j=1∑i−1aijxj(k)−j=i+1∑naijxj(k))=xi(k)+aii1(bi−j=1∑i−1aijxj(k)−j=i∑naijxj(k))然后,在
1
a
i
i
\frac1{a_{ii}}
aii1 处加一个参数
ω
\omega
ω
x
i
(
k
+
1
)
=
x
i
(
k
)
+
ω
a
i
i
(
b
i
−
∑
j
=
1
i
−
1
a
i
j
x
j
(
k
)
−
∑
j
=
i
n
a
i
j
x
j
(
k
)
)
x_i^{(k+1)}=x_i^{(k)} + \frac\omega{a_{ii}}(b_i-\sum_{j=1}^{i-1}a_{ij}x_j^{(k)}-\sum_{j=i}^na_{ij}x_j^{(k)})
xi(k+1)=xi(k)+aiiω(bi−j=1∑i−1aijxj(k)−j=i∑naijxj(k))这样就得到的SOR法的迭代公式了,其矩阵格式为
x
(
k
+
1
)
=
x
(
k
)
+
ω
D
−
1
(
b
+
(
L
+
U
−
D
)
x
(
k
)
)
=
D
−
1
(
ω
(
L
+
U
)
+
(
1
−
ω
)
D
)
x
(
k
)
+
ω
D
−
1
b
x^{(k+1)}=x^{(k)}+\omega D^{-1}(b+(L+U-D)x^{(k)})\\ =D^{-1}(\omega (L+U)+(1-\omega )D)x^{(k)}+\omega D^{-1}b
x(k+1)=x(k)+ωD−1(b+(L+U−D)x(k))=D−1(ω(L+U)+(1−ω)D)x(k)+ωD−1b(依然将系数矩阵A分裂为A = D - L - U)
同理,如果把G-S迭代公式做相同的变换,有
x
i
(
k
+
1
)
=
1
a
i
i
(
b
i
−
∑
j
=
1
i
−
1
a
i
j
x
j
(
k
+
1
)
−
∑
j
=
i
+
1
n
a
i
j
x
j
(
k
)
)
=
x
i
(
k
)
+
1
a
i
i
(
b
i
−
∑
j
=
1
i
−
1
a
i
j
x
j
(
k
+
1
)
−
∑
j
=
i
n
a
i
j
x
j
(
k
)
)
x_i^{(k+1)}=\frac1{a_{ii}}(b_i-\sum_{j=1}^{i-1}a_{ij}x_j^{(k+1)}-\sum_{j=i+1}^na_{ij}x_j^{(k)})\\ =x_i^{(k)} + \frac1{a_{ii}}(b_i-\sum_{j=1}^{i-1}a_{ij}x_j^{(k+1)}-\sum_{j=i}^na_{ij}x_j^{(k)})
xi(k+1)=aii1(bi−j=1∑i−1aijxj(k+1)−j=i+1∑naijxj(k))=xi(k)+aii1(bi−j=1∑i−1aijxj(k+1)−j=i∑naijxj(k))类似地,在
1
a
i
i
\frac1{a_{ii}}
aii1 处加一个参数
ω
\omega
ω
x
i
(
k
+
1
)
=
x
i
(
k
)
+
ω
a
i
i
(
b
i
−
∑
j
=
1
i
−
1
a
i
j
x
j
(
k
+
1
)
−
∑
j
=
i
n
a
i
j
x
j
(
k
)
)
x_i^{(k+1)}=x_i^{(k)} + \frac\omega{a_{ii}}(b_i-\sum_{j=1}^{i-1}a_{ij}x_j^{(k+1)}-\sum_{j=i}^na_{ij}x_j^{(k)})
xi(k+1)=xi(k)+aiiω(bi−j=1∑i−1aijxj(k+1)−j=i∑naijxj(k))这样就得到的SOR法的迭代公式了,其矩阵格式为
x
(
k
+
1
)
=
x
(
k
)
+
ω
D
−
1
(
b
+
L
x
(
k
+
1
)
+
(
U
−
D
)
x
(
k
)
)
D
x
(
k
+
1
)
=
D
x
(
k
)
+
ω
(
b
+
L
x
(
k
+
1
)
+
(
U
−
D
)
x
(
k
)
)
(
D
−
ω
L
)
x
(
k
+
1
)
=
D
x
(
k
)
+
ω
(
U
−
D
)
x
(
k
)
+
ω
b
x
(
k
+
1
)
=
(
D
−
ω
L
)
−
1
(
ω
U
+
(
1
−
ω
)
D
)
x
(
k
)
+
ω
(
D
−
ω
L
)
−
1
b
x^{(k+1)}=x^{(k)}+\omega D^{-1}(b+Lx^{(k+1)}+(U-D)x^{(k)})\\ \;\\ Dx^{(k+1)}=Dx^{(k)}+\omega(b+Lx^{(k+1)}+(U-D)x^{(k)})\\ \;\\ (D-\omega L)x^{(k+1)}=Dx^{(k)}+\omega (U-D)x^{(k)}+ \omega b\\ \;\\ x^{(k+1)}=(D-\omega L)^{-1}(\omega U+(1-\omega)D)x^{(k)} + \omega(D-\omega L)^{-1} b
x(k+1)=x(k)+ωD−1(b+Lx(k+1)+(U−D)x(k))Dx(k+1)=Dx(k)+ω(b+Lx(k+1)+(U−D)x(k))(D−ωL)x(k+1)=Dx(k)+ω(U−D)x(k)+ωbx(k+1)=(D−ωL)−1(ωU+(1−ω)D)x(k)+ω(D−ωL)−1b
至此SOR方法的公式推导完毕。
值得注意的是同J迭代法和G-S迭代法一样,SOR方法有一些用来判断是否收敛的条件,这里就不介绍了,有兴趣的同学请查阅相关书籍。
基于Python实现SOR法
下面基于公式
x
(
k
+
1
)
=
(
D
−
ω
L
)
−
1
(
ω
U
+
(
1
−
ω
)
D
)
x
(
k
)
+
ω
(
D
−
ω
L
)
−
1
b
x^{(k+1)}=(D-\omega L)^{-1}(\omega U+(1-\omega)D)x^{(k)} + \omega(D-\omega L)^{-1} b
x(k+1)=(D−ωL)−1(ωU+(1−ω)D)x(k)+ω(D−ωL)−1b实现算法
算法代码
import numpy as np
def SORfunction(A,b,w,k):
n = A.shape[1]#系数矩阵的列数
D = np.eye(n)
D[np.arange(n),np.arange(n)] = A[np.arange(n),np.arange(n)]
LU = D - A
L = np.tril(LU)
U = np.triu(LU)
D_wL = D - w*L
X = np.zeros(n)
print("第",0,"次迭代的结果:",X)
for i in range(k):
D_wL_inv = np.linalg.inv(D_wL)
X = np.dot(np.dot(D_wL_inv,w * U + (1-w) * D),X) + w * np.dot(D_wL_inv,b)
print("第",i+1,"次迭代的结果:",X)
return X
实验
A = np.array([
[4,-2,-4],
[-2,17,10],
[-4,10,9]
])
b = np.array([10,3,-7])
w = 1.46
k=20
SORfunction(A,b,w,k)
实验结果