Richardson同步迭代法¶

约 303 个字 66 行代码预计阅读时间 2 分钟

计算原理¶

对于线性方程组：

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

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

Richardson 迭代法每一步的分量更新公式为：

\[ x_i^{(k+1)} = x_i^{(k)} - \omega \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\) 个分量
\(\omega\)：松弛因子（步长参数）

同步迭代:

“同步”指的是所有 \(x_i^{(k+1)}\) 都是 用上一轮 \(k\) 的值 \(x_j^{(k)}\) 计算，而不是边算边用新值。
无需提前 \(x2 = x1\)，只需 \(x1 = x0\) 即可，内外循环分明。每轮新值 \(x2\) 算完后整体赋给 \(x1\)。

Fortran数值实现¶

subroutine richardson(a,b,x,x0,n)
! richardson迭代法
    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 ! 默认精度

    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))
        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,"Richardson Iterative",//)
    503 format(i5, 4(3x, f12.8))
end subroutine richardson

数值案例¶

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

主程序：

program main
    implicit none

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

    a = reshape([ &
        0.8, -0.01, 0.12, &
        -0.01, 0.84, 0.05, &
        0.12, 0.05, 0.88], [n,n], order=[2,1])

    ! 定义向量 B
    b = [1.0, 1.0, 1.0]

    x0 = 0.0d0

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

end program main

输出：

                  Richardson Iterative


    1     1.20000000     1.20000000     1.20000000
    2     1.08959999     1.13280004     0.88800001
    3     1.12930559     1.15073283     0.92540161
    4     1.12572317     1.14882174     0.91651354
    5     1.12683682     1.14932733     0.91764181
    6     1.12672497     1.14926895     0.91738792
    7     1.12675635     1.14928331     0.91742175
    8     1.12675291     1.14928154     0.91741448
    9     1.12675380     1.14928195     0.91741549
   10     1.12675369     1.14928190     0.91741528

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

参考文献¶

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