数值分析Python实现系列—— 二、逐次超松弛迭代法(SOR)

二、超松弛迭代法(SOR)

1.原理:

​ 回顾:

​ 在一般情况下 : 收敛过慢甚至不收敛的$B$与$f$,经过对系数矩阵$A$分裂成$A=M-N$的形式, 使得迭代公式变为: $x^{k+1}=(I-M^{-1})A{x}^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$,正是***高思-赛德尔***迭代法,下面推导迭代公式:
$${\text{x}}_{k+1}=I-M^{-1}A{\text{x}}_k+M^{-1}b$$

$${\text{x}}_{k+1}=I-\omega (D-\omega L{)}^{-1}A{\text{x}}_k+\omega (D-\omega L{)}^{-1}b$$

$${\text{x}}_{k+1}=(D-\omega L{)}^{-1}((1-\omega )D+\omega U){\text{x}}_k+\omega (D-\omega L{)}^{-1}b$$

​ 推导完毕,我们较为常用的是下式:
$$(D-\omega L){\text{x}}_{k+1}=((1-\omega )D+\omega U){\text{x}}_k+\omega b$$
以及:
$$\{\begin{array}{ccc}{\text{x}}^{(0)}& =& (x_1^{(0)},\text{...},x_n^{(0)}{)}^T,\\ x_i^{(k+1)}& =& x_i^{(k+)}+\Delta x_i\\ \Delta x_i& =& \omega \frac{b_i-\sum _{j=1}^{i-1}{a}_{ij}{x}_j^{(k+1)}-\sum _{j=1}^n{a}_{ij}{x}_j^{(k)}}{a_{ii}}\\ i& =& 1,2,...,n,k=0,1,...,\\ \omega 为松弛因子& & \end{array}$$

当$\omega >1$时为超松弛迭代,当$\omega <1$时为低松弛迭代

迭代终止条件:$\underset{1\le i\le n}{\max }|\Delta x_i|=\underset{1\le i\le n}{\max }|x_i^{(k+1)}-x_i^{(k)}|<\epsilon$,下面我们试试用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.在根据迭代公式改变各个向量分量时,要注意遍历范围.

结果: