高斯赛德尔迭代c语言_逐次超松弛SOR迭代概述

本文使用 Zhihu On VSCode 创作并发布

逐次超松弛（SOR）迭代法概述

一、方法背景

逐次超松弛迭代法是高斯-塞德尔迭代法的一种变种，是为了解决线性方程组的一种迭代方法。由David M. Young, Jr. 和Stanley P. Frankel于1950年提出，主要为了实现在计算机上求解线性方程组。逐次超松弛迭代法从高斯-塞德尔迭代出发，加入了

${\omega}$ 作为松弛因子，加快了迭代的收敛速度。

高斯-塞德尔迭代：

$x_{i}^{(k+1)}=x_{i}^{(k)}+\frac{1}{a_{i i}}\left(b_{i}-\sum_{j=1}^{i-1} a_{i j} x_{j}^{(k+1)}-\sum_{j=i}^{n} a_{i j} x_{j}^{(k)}\right)$

逐次超松弛迭代：

$x_{i}^{(k+1)}=x_{i}^{(k)}+\frac{\omega}{a_{i i}}\left(b_{i}-\sum_{j=1}^{i-1} a_{i j} x_{j}^{(k+1)}-\sum_{j=i}^{n} a_{i j} x_{j}^{(k)}\right)$

二、数学推导

1. 基础迭代法

对于

$Ax=b$ ，

假定A有如下分解

$A=M-N$ ，其中M是非奇异方阵。

有

$Mx=Nx+b$ 或

$x=Bx+g$ ，其中

$B=M^{-1} N$ 且

$g=M^{-1} b$

建立迭代公式：

$x^{(k+1)}=B x^{(k)}+g$

若

$x^{(k)}$ 收敛于确定的向量

$x^{*}$ ，则

$x^{*}=B x^{*}+g$ ，即

$A x^{*}=b$ ，

$x^{*}$ 即为方程组

$Ax=b$ 的解。

令

$A=D-L-U$ ，其中

$D=diag([a_{11}, a_{22}, \cdots, a_{n n}])$

$L=-\left[\begin{array}{cccc}0 & 0 & \cdots & 0 \\ a_{21} & 0 & \cdots & 0 \\ \vdots & \ddots & \ddots & \vdots \\ a_{n 1} & a_{n 2} & \cdots & 0\end{array}\right]$

$U=-\left[\begin{array}{cccc}0 & a_{12} & \cdots & a_{1 n} \\ \vdots & \ddots & \ddots & \vdots \\ 0 & \ddots & \ddots & a_{n-1, n} \\ 0 & 0 & \cdots & 0\end{array}\right]$

2. 迭代形式

考虑高斯-塞德尔迭代法

$\begin{aligned} x_{i}^{(k+1)} &=-\frac{1}{a_{i t}} \sum_{j=1}^{i-1} a_{i j} x_{j}^{(k+1)}-\frac{1}{a_{i i}} \sum_{j=i+1}^{n} a_{i j} x_{j}^{(k)}+\frac{1}{a_{i i}} b_{i} \\ &=\frac{1}{a_{i i}} b_{i}-\frac{1}{a_{i i}} \sum_{j=1}^{i-1} a_{i j} x_{j}^{(k+1)}-\frac{1}{a_{i i}} \sum_{j=i+1}^{n} a_{i j} x_{j}^{(k)} \\ &=\frac{1}{a_{i i}} b_{i}-\frac{1}{a_{i i}} \sum_{j=1}^{i-1} a_{i j} x_{j}^{(k+1)}-\frac{1}{a_{i j}} \sum_{j=i}^{n} a_{i j} x_{j}^{(k)}+x_{i}^{(k)} \\ &=x_{i}^{(k)}+\frac{1}{a_{i i}}\left(b_{i}-\sum_{j=1}^{i-1} a_{i j} x_{j}^{(k+1)}-\sum_{j=i}^{n} a_{i j} x_{j}^{(k)}\right) \end{aligned}$

令

$\begin{aligned} r_{i}^{(k)} &=x_{i}^{(k+1)}-x_{i}^{(k)} \\ &=\frac{1}{a_{i t}}\left(b_{i}-\sum_{j=1}^{i-1} a_{i j} x_{j}^{(k+1)}-\sum_{j=i}^{n} a_{i j} x_{j}^{(k)}\right) \end{aligned}$

$r_{i}^{(k)}$ 为第

$k+1$ 次迭代时

$\boldsymbol{x}_{i}$ 的改变量

因此

$x_{i}^{(k+1)}=x_{i}^{(k)}+r_{i}^{(k)}$

在改变量r

$_{i}^{(k)}$ 前加一个因子

$\omega,(0<\omega<2),$ 得

${x}_{i}^{(k+1)}={x}_{i}^{(k)}+\omega r_{i}^{(k)}$

$\quad={x}_{i}^{(k)}+\frac{\omega}{{a}_{i i}}\left({b}_{i}-\sum_{j=1}^{i-1} {a}_{i j} {x}_{j}^{(k+1)}-\sum_{j=i}^{n} {a}_{i j} {x}_{j}^{(k)}\right)$

合并

$x_{i}^{k}$ 得到

$x_{i}^{(k+1)}=(1-\omega) x_{i}^{(k)}+\frac{\omega}{a_{i t}}\left(b_{i}-\sum_{j=1}^{i-1} a_{i j} x_{j}^{(k+1)}-\sum_{j=i+1}^{n} a_{i j} x_{j}^{(k)}\right)$

$i=1,2, \cdots n, k=0,1,2, \cdots$

3. 矩阵形式

由高斯-塞德尔迭代的矩阵形式：

$x^{(k+1)}=(D-L)^{-1} U x^{(k)}+(D-L)^{-1} b$

故

$\begin{aligned} {x}^{(k+1)} &={D}^{-1} {L} {x}^{(k+1)}+{D}^{-1} {U} {x}^{(k)}+{D}^{-1} {b} \\ &={x}^{(k)}-{x}^{(k)}+{D}^{-1} {L} {x}^{(k+1)}+{D}^{-1} {U} {x}^{(k)}+{D}^{-1} {b} \end{aligned}$

记

$r^{(k)}=x^{(k+1)}-x^{(k)}=-{x}^{(k)}+{D}^{-1} {L} {x}^{(k+1)}+{D}^{-1} {U} {x}^{(k)}+{D}^{-1} {b}$

则SOR迭代为

$x^{(k+1)}=x^{(k)}+\omega r^{(k)}=(1-\omega) x^{(k)}+\omega\left(D^{-1} L x^{(k+1)}+D^{-1} U x^{(k)}+D^{-1} b\right)$

两边同时左乘D，可得

$D x^{(k+1)}=(1-\omega) D x^{(k)}+\omega\left(L x^{(k+1)}+U x^{(k)}+b\right)$

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

即

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

令

$B_{\omega}=(D-\omega L)^{-1}((1-\omega) D+\omega U)$

$f_{\omega}=\omega(D-\omega L)^{-1} b$

则

${x}^{(k+1)}={B}_{\omega} {x}^{(k)}+{f}_{{\omega}}$ 为SOR迭代的矩阵形式，

$B_{\omega}$ 为SOR迭代法的迭代矩阵

三、收敛分析与误差分析

1. 收敛条件的判定

收敛的充要条件是

$\rho(B \omega)<1$
收敛的必要条件是
$0< \omega <2$
若系数矩阵A对称正定，则
$0< \omega <2$ 时SOR迭代收敛
若系数矩阵A为严格对角占优矩阵，则
$0< \omega \leq 1$ 时SOR迭代收敛

2. 收敛速度分析

假设

$0< \omega <2$ ，Jacobi迭代矩阵

$C_{\mathrm{Jac}}=I-D^{-1} A$ 只有实特征值，Jacobi迭代收敛

$\mu=\rho\left(C_{\mathrm{Jac}}\right)<1$ ，方程组存在唯一解即

$det A \neq 0$

收敛速度可以如下表示

$\rho\left(C_{\omega}\right)=\left\{\begin{array}{ll}\frac{1}{4}(\omega \mu+\sqrt{\omega^{2} \mu^{2}-4(\omega-1)})^{2}, & 0<\omega \leq \omega_{\text {opt }} \\ \omega-1, & \omega_{\text {opt }}<\omega<2\end{array}\right.$

最佳松弛参数

$\omega _ {opt}$ 为

$\omega_{\mathrm{opt}}=1+\left(\frac{\mu}{1+\sqrt{1-\mu^{2}}}\right)^{2}$

3. 误差分析

给定初值

$x^{(0)}$ ，

$x^{*}$ 是真解，若

$\|B\|<1$ ，则有如下估计

$\left\|x^{(k)}-x^{*}\right\| \leqslant \frac{\|B\|^{k}}{1-\|B\|}\left\|x^{(1)}-x^{(0)}\right\|$

$\left\|x^{(k)}-x^{*}\right\| \leqslant \frac{\|B\|}{1-\|B\|}\left\|x^{(k)}-x^{(k-1)}\right\|$

四、理论分析

对于线性方程组，可以通过高斯消去、LU分解等直接法进行求解，但是直接法对于系数矩阵要求高，且程序复杂，对于求解大型线性方程组有较多限制。而像SOR迭代法的基础迭代法，只需按照迭代格式进行迭代，便于编程求解，可以迭代到任意精度，节省计算量和空间。

SOR迭代法由高斯-塞德尔迭代法改进而来，通过加入了松弛因子

$\omega$ ，调整每次迭代的增量，在设定合适的

$\omega$ 时，SOR迭代速度明显快于高斯－塞德尔迭代。

编程时常使用SOR迭代法的矩阵形式，当对角矩阵D中对角线上元素有0时，下三角矩阵

$(D+\omega L)$ 对角线上元素为0，此时无法求逆，进而无法进行迭代。故在编程实现时，进行线性方程组的预处理。由SOR收敛条件可知，当系数矩阵

$A$ 对称正定且

$0<\omega <2$ 时，SOR迭代收敛，对于非奇异方针

$A$ ，

$A^{T}A$ 对阵正定，故若将求解

$Ax=b$ 的线性方程组转化为求解

$A^{T}Ax=A^{T}b$ ，新的系数矩阵为

$A^{T}A$ ，此时使用SOR迭代的矩阵形式可以避免无法求逆的情况出现。

五、MATLAB实现与实验

function [x,iter] = sor(A,b,x0,omega,tol)
if det(A) == 0error('无解或无唯一解');
end%方程组预处理
b=A'*b;
A=A'*A;D = diag(diag(A));
L = D-tril(A);
U = D-triu(A);
x=x0;
for iter=1:500B = (D-L*omega)((1-omega)*D+omega*U);f = (D-L*omega)b*omega;x = B*x+f;eerror = norm( b-A*x ) / norm(b);if ( eerror < tol )break;end
end

不同松弛因子的迭代速度对比

不同迭代方式的迭代速度对比

六、总结

逐次超松弛(SOR)迭代法从高斯-塞德尔迭代法出发，通过加入松弛因子

$\omega$ 来加快收敛速度，对编程实现友好，常用于求解大型线性方程组。SOR迭代提出的目的是为在计算机上进行线性方程组的求解，通过SOR迭代的矩阵形式，可以很简洁地实现程序编写，并能通过简单地线性方程组预处理解决可能出现的无法求逆的情况。从实验中能够看出，SOR迭代的速度比高斯-塞德尔迭代快，但需要基于合适的松弛因子

$\omega$ 。

七、参考资料

[1] https://en.wikipedia.org/wiki/Successive_over-relaxation

[2] https://zhuanlan.zhihu.com/p/31066592

[3] https://www.jianshu.com/p/e14d9e910984