二、超松弛迭代法(SOR)

1.原理:

? 回顧:

? 在一般情況下 : 收斂過慢 甚至 不收斂 的 \(B\) 與 \(f\) ,經過對系數矩陣 \(A\) 分裂 成 \(A = M - N\) 的形式, 使得迭代公式變為: \(x^{k+1}=(I-M^{-1})Ax^{k}+M^{-1}f\)

? 雅克比迭代法 選取 : 現將 \(A\) 如下分解 \(A = D-L-U\) , \(D\) 為對角陣, \(L\) 為下三角陣, \(U\) 為上三角陣,取 \(M \equiv D\) ,取 \(N \equiv L+U\) ,

? 在這一章中我們選取下三角矩陣 \(M=\frac{1}{\omega}(D-\omega L),\omega>0\) ,其中 \(\omega\) 為松弛因子,我們可以發現當 \(\omega\) 為1時, \(M=D-L\) ,正是 高思-賽德爾 迭代法,下面推導迭代公式:
\[ \textbf{x}_{k+1}=I-M^{-1}A\textbf{x}_{k}+M^{-1}b \]

\[ \textbf{x}_{k+1}=I-\omega(D-\omega L)^{-1}A\textbf{x}_{k}+\omega (D-\omega L)^{-1}b \]

\[ \textbf{x}_{k+1}=(D-\omega L)^{-1}((1-\omega)D+\omega U)\textbf{x}_{k}+\omega (D-\omega L)^{-1}b \]

? 推導完畢,我們較為常用的是下式:
\[ (D-\omega L)\textbf{x}_{k+1}=((1-\omega)D+\omega U)\textbf{x}_{k}+\omega b \]
以及:
\[ \left\{ \begin{matrix} \textbf{x}^{(0)} &=& (x_1^{(0)},\textbf{...},x_n^{(0)})^{T}, \\ x_i^{(k+1)} &=& x_i^{(k+)}+ \Delta x_{i} \\ \Delta x_{i} &=& \omega \frac{b_{i}- \sum\limits_{j=1}^{i-1}a_{ij}x_{j}^{(k+1)}-\sum\limits_{j=1}^{n}a_{ij}x_{j}^{(k)}}{a_{ii}} \\ i &=&1,2,...,n,k=0,1,...,\\ \omega為松弛因子 \end{matrix} \right. \]

當 \(\omega>1\) 時為超松弛迭代,當 \(\omega<1\) 時為低松弛迭代

迭代 終止條件 : \(\mathop{max}\limits_{1\le i\le n}|\Delta x_{i}| = \mathop{max}\limits_{1\le i\le n}|x_{i}^{(k+1)}-x_{i}^{(k)}|<\varepsilon\) ,下面我們試試用Python實現這一功能.

2.實現:

          
            import numpy as np
import matplotlib.pyplot as plt
MAX = 110                                         # 遍歷最大次數
A = np.array([[-4, 1, 1, 1], [1, -4, 1, 1], [1, 1, -4, 1], [1, 1, 1, -4]])
b = np.array([[1], [1], [1], [1]])                # 注意這里取列向量
omega_list = [1 + 0.005 * i for i in range(100)]  # 取到不同的omega值,觀察趨勢
length = len(A)
count = []                                        # 記錄遍歷的次數
for omega in omega_list:                          # 遍歷每一個omega值
    times = 0
    x_0 = np.zeros((length, 1))
    x_hold = x_0 + np.ones((length, 1))
    while (np.linalg.norm(x_hold - x_0, ord=2) >= 10 ** (-5)) and (times <= MAX):
        # 遍歷停止條件以k+1次與k次迭代的向量差的二范數以及遍歷最大次數為標準
        x_hold = x_0.copy()                      # 這里不要用賦值,要用copy
        x_new = x_0.copy()
        for i in range(length):
            # 根據迭代公式迭代
            x_new[i][0] = x_0[i][0] + omega * (b[i][0] - sum([A[i][j] * x_new[j][0] for j in range(i)]) - sum(
                [A[i][j] * x_0[j][0] for j in range(i, length)])) / A[i][i]
            x_0 = x_new.copy()
        times += 1
    count.append(times)
plt.plot(omega_list, count)                      # 觀察omega與迭代次數的關系
plt.show()
          
        

思路:

? 1.遍歷設限:第一種是到達精度,到達精度停止迭代,第二種是到達規定最大次數,這個可以自己設定.

? 2.在根據迭代公式改變各個向量分量時,要注意遍歷范圍.

結果: