二、超松弛迭代法(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)}|

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.在根据迭代公式改变各个向量分量时,要注意遍历范围.

结果:

原文出处：https://www.cnblogs.com/MyBlog-MrY/p/10848322.html