逐次超松弛迭代法SOR

约 372 个字 99 行代码 预计阅读时间 2 分钟

计算原理

对于线性方程组：

\[ \mathbf{A}\mathbf{x} = \mathbf{b} \]

其中 \(\mathbf{A}\) 是 \(n \times n\) 的系数矩阵，\(\mathbf{x}\) 是未知向量，\(\mathbf{b}\) 是已知常数向量。

SOR（Successive Over-Relaxation）超松弛迭代法是在 Gauss-Seidel 方法基础上的推广，通过引入松弛因子 \(\omega\)，加速线性方程组 \(\mathbf{A}\mathbf{x} = \mathbf{b}\) 的迭代收敛过程。

• 当 \(\omega = 1\) 时，SOR 退化为 Gauss-Seidel 方法。
• 当 \(\omega > 1\) 时，称为 逐次超松弛迭代 （Over-Relaxation），可加速收敛。
• 当 \(0 < \omega < 1\) 时，称为 逐次低松弛迭代

分量形式

SOR 方法对第 \(i\) 个分量的迭代公式为：

\[ x_i^{k+1} = \frac{1}{a_{ii}} \left[ (1 - \omega) x_i^k + \omega \left( -\sum_{j=1}^{i-1} a_{ij} x_j^{k+1} - \sum_{j=i+1}^{n} a_{ij} x_j^k + b_i \right) \right] \]

其中：

• \(k\) 表示第 \(k\) 次迭代；
• \(x_j^{k+1}\) 表示本轮已更新的新分量（\(j < i\)）；
• \(x_j^k\) 表示上一轮的旧分量（\(j > i\)）；
• \(b_i\) 是方程组右端项；
• \(\omega\) 是松弛因子。

松弛因子经验值推荐：

问题类型	推荐 \(\omega\)
对称正定线性方程	1.2 ~ 1.8
一般稀疏线性方程	1.0 ~ 1.7
高维稀疏系统	1.1 ~ 1.6

算法步骤

1. 初始化
2. 给定初始猜测 \(x^0\)（如全零向量或任意向量），记作 x0。
3. 主迭代循环
4. 对 \(k = 1, 2, \ldots\)，直到满足收敛判据为止：
   • 对每个分量 \(i = 1, 2, \ldots, n\)：
   • 累加 \(j=1\) 到 \(i-1\) 的 \(a_{ij} x_j^{k+1}\)（已更新的新值）；
   • 累加 \(j=i+1\) 到 \(n\) 的 \(a_{ij} x_j^k\)（上一轮的旧值）；
   • 按上式用 \(\omega\) 加权，求出新的 \(x_i^{k+1}\)；
   • 计算新旧解的范数差 \(\|x^{k+1} - x^k\|\) 判断是否收敛。
5. 输出结果
6. 当收敛判据（如 \(\|x^{k+1} - x^k\| < \varepsilon\)）满足时，停止迭代，输出 \(x^{k+1}\) 作为近似解。

Fortran数值实现

subroutine sor(a,b,x,x0,n)
! sor迭代法
    integer, intent(in) :: n
    real*8, intent(in) :: a(n,n), b(n), x0(n)
    real*8, intent(out) :: x(n)
    real*8 :: omiga = 1.2d0 ! 松弛因子

    real*8 :: x1(n), x2(n)
    integer :: i, j, k, iter_max = 1000 ! 默认迭代次数
    real*8 :: eps = 1.0d-7 ! 默认精度
    logical :: diagonally_dominant

    ! ---------- 增加收敛性判据 ----------
    diagonally_dominant = .true.
    do i = 1, n
        s = 0.0d0
        do j = 1, n
            if (j /= i) s = s + abs(a(i, j))
        end do
        if (abs(a(i, i)) <= s) then
            diagonally_dominant = .false.
            write(*, '(A, I3, A)') 'Warning: Row ', i, &
                ' is not strictly diagonally dominant:'
            write(*, '(A, F12.5, A, F12.5)') '  |a_ii| = ', abs(a(i, i)), &
                ', sum|a_ij| = ', s
        end if
    end do
    if (.not. diagonally_dominant) then
        write(*, '(A)') 'The coefficient matrix is not strictly diagonally dominant.'
        write(*, '(A)') 'SOR iteration may not converge!'
        return
    end if
    ! -------------------------------------

    write(*,501)
    ! 迭代之前两值都设为初值
    x1 = x0
    x2 = x1
    do k = 1, iter_max
        do i = 1, n
            s = 0
            do j = 1, n
                ! 若j < i，则表示这些量已经更新过了，则下一个元素就用最新的量计算
                ! 若j > i，则表示还没有计算到这些量，所以就用上一次迭代的结果
                if (j < i) then
                    s = s + a(i,j)*x2(j)
                else if (j > i) then
                    s = s + a(i,j)*x1(j)
                end if
            end do
            x2(i) = (b(i)-s)*omiga/A(i,i)+(1-omiga)*x1(i)
        end do
        ! 精度判断
        dx2 = norm2(x2 - x1) ! 可使用内置函数计算二范数
        if (dx2 < eps) exit

        x1 = x2
        write(*,503) k, (x2(i), i=1,n)

    end do
    x = x2

    501 format(//,18x,"SOR Iterative",//)
    503 format(i5, 4(3x, f12.8))
end subroutine sor

数值案例

\[ \begin{array} { r } { \left( \begin{array} { l l l l } { 5 } & { - 1 } & { - 1 } & { - 1 } \\ { - 1 } & { 1 0 } & { - 1 } & { - 1 } \\ { - 1 } & { - 1 } & { 5 } & { - 1 } \\ { - 1 } & { - 1 } & { - 1 } & { 1 0 } \end{array} \right) \mathbf { x } = \left( \begin{array} { l } { - 4 } \\ { 1 2 } \\ { 8 } \\ { 3 4 } \end{array} \right) } \end{array} \]

主程序：

program main
    implicit none

    integer, parameter :: n = 4
    real*8 :: a(n,n), b(n), x(n), x0(n)

    a = reshape([5.0d0, -1.0d0, -1.0d0, -1.0d0, &
                -1.0d0, 10.0d0, -1.0d0, -1.0d0, &  
                -1.0d0, -1.0d0, 5.0d0, -1.0d0, &  
                -1.0d0, -1.0d0, -1.0d0, 10.0d0], [n, n], order=[2,1]) 

    b = [-4.0d0, 12.0d0, 8.0d0, 34.0d0]

    x0 = 0.0d0

    call sor(a, b, x, x0, n)

end program main

输出：

                  SOR Iterative


    1    -0.96000000     1.32480000     2.00755200     4.36468224
    2     1.07928822     2.06922269     3.32165596     3.98348358
    3     1.07398929     2.03165092     2.95705852     4.01082713
    4     0.98509092     1.98802700     3.00473511     3.99517694
    5     1.00008719     2.00239451     2.99849105     4.00108134
    6     1.00045462     1.99952434     3.00055626     3.99984796
    7     0.99989193     2.00013067     2.99985768     4.00001604
    8     1.00002267     1.99996143     3.00002850     3.99999830
    9     0.99999264     2.00001005     2.99999454     4.00000001
   10     1.00000257     1.99999764     3.00000115     4.00000016
   11     0.99999923     2.00000054     2.99999975     3.99999991
   12     1.00000020     1.99999988     3.00000005     4.00000003
   13     0.99999995     2.00000003     2.99999999     3.99999999

注意：子程序中默认：implicit real*8(a-z)

参考文献

1. 宋叶志,茅永兴,赵秀杰.Fortran 95/2003科学计算与工程[M].清华大学出版社,2011.