基于Python实现的SOR法(逐次超松弛迭代法)

---

Jacobi迭代法与Gauss-Seidel迭代法回顾

上一篇博客介绍了J迭代法和G-S迭代法的原理，现在复习一下，J迭代法的迭代公式为
$$\begin{gathered} x_i^{(k+1)}=\frac{1}{a_{ii}}(b_i-\sum _{j=1}^{i-1}{a}_{ij}{x}_j^{(k)}-\sum _{j=i+1}^n{a}_{ij}{x}_j^{(k)}) \\ i=1,2,\cdots ,n,k=0,1,\cdots \end{gathered}$$矩阵形式为
$$x^{(k+1)}=D^{-1}(L+U)x^{(k)}+D^{-1}b$$(这里将系数矩阵A分裂为A = D - L - U)

G-S迭代法的迭代公式为
$$\begin{gathered} x_i^{(k+1)}=\frac{1}{a_{ii}}(b_i-\sum _{j=1}^{i-1}{a}_{ij}{x}_j^{(k+1)}-\sum _{j=i+1}^n{a}_{ij}{x}_j^{(k)}) \\ i=1,2,\cdots ,n,k=0,1,\cdots \end{gathered}$$矩阵形式为
$$x^{(k+1)}=(D-L{)}^{-1}U{x}^{(k)}+(D-L{)}^{-1}b$$(仍将系数矩阵A分裂为A = D - L - U)

---

SOR法(逐次超松弛迭代法)

超松弛迭代法可视为J迭代法或G-S迭代法的改进版，我们可以把J迭代法的迭代公式做如下变换
$$\begin{gathered} x_i^{(k+1)}=\frac{1}{a_{ii}}(b_i-\sum _{j=1}^{i-1}{a}_{ij}{x}_j^{(k)}-\sum _{j=i+1}^n{a}_{ij}{x}_j^{(k)}) \\ =x_i^{(k)}+\frac{1}{a_{ii}}(b_i-\sum _{j=1}^{i-1}{a}_{ij}{x}_j^{(k)}-\sum _{j=i}^n{a}_{ij}{x}_j^{(k)}) \end{gathered}$$然后，在 $\frac{1}{a_{ii}}$ 处加一个参数 $\omega$
$$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}^n{a}_{ij}{x}_j^{(k)})$$这样就得到的SOR法的迭代公式了，其矩阵格式为
$$\begin{gathered} 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 \end{gathered}$$(依然将系数矩阵A分裂为A = D - L - U)

同理，如果把G-S迭代公式做相同的变换，有
$$\begin{gathered} x_i^{(k+1)}=\frac{1}{a_{ii}}(b_i-\sum _{j=1}^{i-1}{a}_{ij}{x}_j^{(k+1)}-\sum _{j=i+1}^n{a}_{ij}{x}_j^{(k)}) \\ =x_i^{(k)}+\frac{1}{a_{ii}}(b_i-\sum _{j=1}^{i-1}{a}_{ij}{x}_j^{(k+1)}-\sum _{j=i}^n{a}_{ij}{x}_j^{(k)}) \end{gathered}$$类似地，在 $\frac{1}{a_{ii}}$ 处加一个参数 $\omega$
$$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}^n{a}_{ij}{x}_j^{(k)})$$这样就得到的SOR法的迭代公式了，其矩阵格式为
$$\begin{gathered} x^{(k+1)}=x^{(k)}+\omega D^{-1}(b+L{x}^{(k+1)}+(U-D)x^{(k)}) \\ D{x}^{(k+1)}=D{x}^{(k)}+\omega (b+L{x}^{(k+1)}+(U-D)x^{(k)}) \\ (D-\omega L)x^{(k+1)}=D{x}^{(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 \end{gathered}$$
至此SOR方法的公式推导完毕。

值得注意的是同J迭代法和G-S迭代法一样，SOR方法有一些用来判断是否收敛的条件，这里就不介绍了，有兴趣的同学请查阅相关书籍。

---

基于Python实现SOR法

下面基于公式
$$x^{(k+1)}=(D-\omega L{)}^{-1}(\omega U+(1-\omega )D)x^{(k)}+\omega (D-\omega L{)}^{-1}b$$实现算法

算法代码

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)

实验结果