Jacobi超松弛迭代法¶

约 348 个字 74 行代码预计阅读时间 2 分钟

计算原理¶

对于线性方程组：

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

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

Jacobi 超松弛迭代法对第 \(i\) 个分量的同步迭代公式为：

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

其中：

\(x_i^{(k)}\)：第 \(k\) 次迭代的第 \(i\) 个分量
\(a_{ij}\)：系数矩阵 \(\mathbf{A}\) 的元素
\(b_i\)：右端向量的第 \(i\) 个分量
\(a_{ii}\)：第 \(i\) 行的对角元
\(\omega\)：松弛因子（超松弛 \(\omega>1\)，标准Jacobi为 \(\omega=1\)）
\(k\)：当前迭代步数

有时也写作：

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

所有 \(x_i^{(k+1)}\) 都是 用上一轮 \(x_j^{(k)}\) 计算的（同步更新）
当 \(\omega=1\) 时，退化为普通 Jacobi 迭代
\(\omega>1\) 时为超松弛，加速收敛（但\(\omega\)过大可能发散）

Fortran数值实现¶

subroutine jacobi_sor(a,b,x,x0,n)
! Jacobi超松弛迭代法
    integer, intent(in) :: n
    real*8, intent(in) :: a(n,n), b(n), x0(n)
    real*8, intent(out) :: x(n)

    real*8 :: x1(n), x2(n)
    integer :: i, j, k, iter_max = 1000 ! 默认迭代次数
    real*8 :: eps = 1.0d-7 ! 默认精度
    real*8 :: omiga = 1.04d0 ! 松弛因子

    write(*,501)
    x1 = x0
    do k =1, iter_max
        do i = 1, n
            s = 0
            do j = 1, n
                s = s + a(i,j)*x1(j)
            end do
            ! 迭代更新
            x2(i) = x1(i)-omiga*(s-b(i))/a(i,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,"Jacobi_SOR Iterative",//)
    503 format(i5, 4(3x, f12.8))
end subroutine jacobi_sor

数值案例¶

\[ \begin{array} { r } { \left( \begin{array} { l l l l } { 7 } & { 2 } & { 1 } & { - 2 } \\ { 9 } & { 1 5 } & { 3 } & { - 2 } \\ { - 2 } & { - 2 } & { 1 1 } & { 5 } \\ { 1 } & { 3 } & { 2 } & { 1 3 } \end{array} \right) \mathbf { x } = \left( \begin{array} { l } { 4 } \\ { 7 } \\ { - 1 } \\ { 0 } \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([ &
            7, 2, 1, -2, &
            9, 15, 3, -2, &
            -2, -2, 11, 5, &
            1, 3, 2, 13], [4,4], order=[2,1])

    b = [4.0, 7.0, -1.0, 0.0]

    x0 = 0.0d0

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

end program main

输出：

                  Jacobi_SOR Iterative


    1     0.59428571     0.48533333    -0.09454545     0.00000000
    2     0.44034771     0.11475117     0.11338251    -0.14889558
    3     0.48148565     0.16173590     0.07627040    -0.07495348
    4     0.49336413     0.15215906     0.05946363    -0.08654059
    5     0.49478866     0.14701902     0.06604866    -0.08203986
    6     0.49661801     0.14559013     0.06295508    -0.08215385
    7     0.49739517     0.14513343     0.06320843    -0.08145773
    8     0.49766899     0.14471058     0.06292982    -0.08147868
    9     0.49781886     0.14461168     0.06292268    -0.08135368
   10     0.49788045     0.14454094     0.06287352    -0.08134579
   11     0.49790866     0.14451665     0.06287002    -0.08132619
   12     0.49792109     0.14450347     0.06286164    -0.08132285
   13     0.49792675     0.14449845     0.06286025    -0.08131947
   14     0.49792923     0.14449587     0.06285883    -0.08131863
   15     0.49793035     0.14449484     0.06285847    -0.08131801
   16     0.49793085     0.14449434     0.06285821    -0.08131782
   17     0.49793107     0.14449413     0.06285813    -0.08131771
   18     0.49793117     0.14449403     0.06285809    -0.08131767

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

参考文献¶

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