广义 Richardson 迭代法¶

约 364 个字 67 行代码预计阅读时间 2 分钟

计算原理¶

对于线性方程组：

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

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

广义 Richardson 迭代法允许每个分量有不同的松弛因子 \(\omega_i\)，其同步迭代的分量格式为：

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

\[ \begin{array} { r } { \Omega = \left( \begin{array} { l l l l } { \omega _ { 1 } } & { 0 } & { \cdots } & { 0 } \\ { 0 } & { \omega _ { 2 } } & { \cdots } & { 0 } \\ { \vdots } & { \vdots } & { \ddots } & { \vdots } \\ { 0 } & { 0 } & { \cdots } & { \omega _ { n } } \end{array} \right) } \end{array} \]

其中：

\(x_i^{(k)}\)：第 \(k\) 次迭代的第 \(i\) 个分量
\(a_{ij}\)：系数矩阵 \(\mathbf{A}\) 的元素
\(b_i\)：右端向量的第 \(i\) 个分量
\(\omega_i\)：第 \(i\) 个分量对应的松弛因子（步长参数）
\(k\)：当前迭代步数

2. 迭代特性¶

所有 \(x_j^{(k)}\) 都是用 同一轮 的上一迭代结果计算（同步更新）。
每个分量 \(i\) 可以选择 独立的步长 \(\omega_i\)，相当于因子矩阵 \(\mathbf{\Omega}\) 是对角矩阵。

Fortran数值实现¶

subroutine Grichardson(a,b,x,x0,omiga,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 :: x1(n), x2(n), omiga(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(i)*(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,"GRichardson Iterative",//)
    503 format(i5, 4(3x, f12.8))
end subroutine Grichardson

数值案例¶

\[ { \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), omiga(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]
    ! 松弛向量
    omiga = [1.18,1.17,1.23]

    x0 = 0.0d0

    call Grichardson(a, b, x, x0, omiga, n)

end program main

输出：

                  GRichardson Iterative


    1     1.17999995     1.16999996     1.23000002
    2     1.08571799     1.13197503     0.88252501
    3     1.12919196     1.15054519     0.92741151
    4     1.12548970     1.14874738     0.91615404
    5     1.12685522     1.14933171     0.91773867
    6     1.12671420     1.14926503     0.91737061
    7     1.12675763     1.14928377     0.91742585
    8     1.12675246     1.14928137     0.91741374
    9     1.12675386     1.14928197     0.91741565
   10     1.12675367     1.14928189     0.91741525

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

参考文献¶

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