數值分析Python實現系列—— 二、逐次超松弛叠代法(SOR)
二、超松弛叠代法(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\)
\[ \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.在根據叠代公式改變各個向量分量時,要註意遍歷範圍.
結果:
數值分析Python實現系列—— 二、逐次超松弛叠代法(SOR)