Jacobi迭代法¶

约 514 个字 125 行代码预计阅读时间 3 分钟

计算原理¶

对于线性方程组：

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

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

逐项展开¶

把第 $i$ 个分量单独写出来：

\[ a_{i1}x_1 + a_{i2}x_2 + \cdots + a_{ii}x_i + \cdots + a_{in}x_n = b_i \]

解 $x_i$ 得：

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

也就是：

\[ x_i = -\sum_{j \neq i} \frac{a_{ij}}{a_{ii}}x_j + \frac{b_i}{a_{ii}} \]

设 $b_{ij} = -a_{ij}/a_{ii}$，$g_i = b_i/a_{ii}$，就变成

\[ x_i = \sum_{j \neq i} b_{ij} x_j + g_i \]

用向量、矩阵表达¶

把所有 $x_i$ 放到一起

\[ \mathbf{x} = \mathbf{B}\mathbf{x} + \mathbf{g} \]

Jacobi 迭代公式¶

把“新一轮未知数”的所有分量都记为 $\mathbf{x}^{(k+1)}$，用“上一轮算出来的解”$\mathbf{x}^{(k)}$ 来近似代入右侧，得 $$ \boxed{ \mathbf{x}^{(k+1)} = \mathbf{B} \mathbf{x}^{(k)} + \mathbf{g} } $$

每一步迭代：

用上一次迭代的结果，按上述公式计算出新的每个分量

具体操作流程¶

初始化
给定初始猜测 $x^{(0)}$（可设置为0），令 $x_1 = x^{(0)}$。
迭代主循环（最多迭代 iter_max 次）
对每一轮迭代 $k=1,2,\ldots$：
1. 逐个分量计算新解 对每个分量 $i=1,2,\ldots,n$，根据 Jacobi 公式计算： $$ x_i^{(k+1)} = \frac{1}{a_{ii}} \left( b_i - \sum_{j=1,\, j\neq i}^n a_{ij} x_j^{(k)} \right) $$
2. 收敛性判断 计算新旧解的欧氏距离： $$ \Delta x = \sqrt{ \sum_{i=1}^{n} \left( x_i^{(k+1)} - x_i^{(k)} \right)^2 } $$
  - 若 $\Delta x < \epsilon$，即变化小于预设容差，终止迭代。
3. 更新迭代解 将 $x^{(k+1)}$ 赋值给 $x^{(k)}$，即 $$ x_1 \leftarrow x_2 $$
输出结果
迭代结束后，$x_2$ 即为所求近似解 $x$。

Fortran数值实现¶

subroutine jacobi(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 ! 默认精度

    x1 = x0 
    write(*,501)
    write(*,502)

    do k = 1, iter_max
        do i = 1, n
            s = 0
            do j = 1, n
                if (j == i) cycle
                s = s + a(i,j)*x1(j)
            end do
            x2(i) = (b(i) - s)/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(//,26x,"Jacobi Iterative",//)
502 format("  Iter", 5x, "  x1", 12x, "  x2", 12x, "  x3")
503 format(i5, 3(3x, f12.8))

end subroutine jacobi

数值案例¶

\[ \begin{pmatrix} 10 & -1 & 0 \\ -1 & 10 & -2 \\ 0 & -4 & 10 \end{pmatrix} \begin{pmatrix} x_1 \\ x_2 \\ x_3 \end{pmatrix} = \begin{pmatrix} 9 \\ 7 \\ 6 \end{pmatrix} \]

主程序：

program main
    implicit none

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

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

    b = [9.0d0, 7.0d0, 6.0d0]

    x0 = 0.0d0

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

end program main

输出：

                          Jacobi Iterative


  Iter       x1              x2              x3
    1     0.90000000     0.70000000     0.60000000
    2     0.97000000     0.91000000     0.88000000
    3     0.99100000     0.97300000     0.96400000
    4     0.99730000     0.99190000     0.98920000
    5     0.99919000     0.99757000     0.99676000
    6     0.99975700     0.99927100     0.99902800
    7     0.99992710     0.99978130     0.99970840
    8     0.99997813     0.99993439     0.99991252
    9     0.99999344     0.99998032     0.99997376
   10     0.99999803     0.99999410     0.99999213
   11     0.99999941     0.99999823     0.99999764
   12     0.99999982     0.99999947     0.99999929
   13     0.99999995     0.99999984     0.99999979
   14     0.99999998     0.99999995     0.99999994

收敛性¶

Jacobi 和高斯赛德尔迭代的收敛性都可通过 对角占优 的方法进行判断，即对每一行：

\[\begin{array} { r } { | a _ { i i } | > \sum _ { j \neq i } | a _ { i j } | } \end{array}\]

在子程序最开始 自动检查每一行是否对角占优，如果发现不满足则打印提示信息，并且可以直接return退出子程序（避免盲目迭代）。若采用 弱对角占优（非严格），可以把判断条件改成 abs(a(i, i)) < s（允许等号）

subroutine jacobi(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 :: s, dx2
    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, 2F12.5)') 'Warning: Row ', i, ' is not strictly diagonally dominant |a_ii|, sum|a_ij|=', abs(a(i, i)), s
        end if
    end do

    if (.not. diagonally_dominant) then
        write(*, '(A)') 'The coefficient matrix is not strictly diagonally dominant. Jacobi iteration may not converge!'
        return
    end if
    ! ----------------------------------------

    x1 = x0

    do k = 1, iter_max
        do i = 1, n
            s = 0
            do j = 1, n
                if (j == i) cycle
                s = s + a(i, j) * x1(j)
            end do
            x2(i) = (b(i) - s) / a(i, i)
        end do
        ! 精度判断
        dx2 = norm2(x2 - x1)
        if (dx2 < eps) exit

        x1 = x2
    end do
    x = x2

end subroutine jacobi

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

参考文献¶

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