跳转至

Gauss-Seidel迭代法

约 799 个字 97 行代码 预计阅读时间 4 分钟

计算原理

对于线性方程组：

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

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

Gauss-Seidel 迭代法是一种常用的线性方程组求解迭代算法。它的核心思想是：在每一轮迭代中，每个分量的计算都即时用上本轮已经更新过的新值，而不是等所有分量都算完才整体更新。这和 Jacobi 方法的“完全用上一步的值”是最大不同。

矩阵分解与推导

对于线性方程组 \(\mathbf{A}\mathbf{x} = \mathbf{b}\)，可以将系数矩阵 \(\mathbf{A}\) 分解为三部分：

  • 对角阵 \(\mathbf{D}\)，对角元为 \(a_{ii}\)
  • 严格下三角部分 \(\mathbf{L}\)
  • 严格上三角部分 \(\mathbf{U}\)

具体分解为：

\[ \mathbf{A} = \mathbf{D}(\mathbf{I} - \mathbf{L}) - \mathbf{D}\mathbf{U} \]

其中

  • \(\mathbf{D} = \text{diag}(a_{11}, a_{22}, \ldots, a_{nn})\)
  • \(\mathbf{L}\) 为单位下三角矩阵，但去掉对角线，\(\mathbf{U}\) 为上三角部分

\(\mathbf{ L }\)为：

\[ \mathbf { L } = { \left( \begin{array} { l l l l l l } { 0 } & { 0 } & { 0 } & { \cdots } & { 0 } & { 0 } \\ { - { \frac { a _ { 2 1 } } { a _ { 2 2 } } } } & { 0 } & { 0 } & { \cdots } & { 0 } & { 0 } \\ { - { \frac { a _ { 3 1 } } { a _ { 3 3 } } } } & { - { \frac { a _ { 3 2 } } { a _ { 3 3 } } } } & { 0 } & { \cdots } & { 0 } & { 0 } \\ { \vdots } & { \vdots } & { \vdots } & { \ddots } & { \vdots } & { \vdots } \\ { - { \frac { a _ { n 1 } } { a _ { n n } } } } & { - { \frac { a _ { n 2 } } { a _ { n n } } } } & { \cdots } & { \cdots } & { - { \frac { a _ { n , n - 1 } } { a _ { n n } } } } & { 0 } \end{array} \right) } \]

\(\mathbf { U }\)为：

\[ \mathbf { U } = \left( \begin{array} { c c c c c c } { 0 } & { - \frac { a _ { 1 2 } } { a _ { 1 1 } } } & { - \frac { a _ { 1 3 } } { a _ { 1 1 } } } & { \cdots } & { - \frac { a _ { 1 , n - 1 } } { a _ { 1 1 } } } & { - \frac { a _ { 1 n } } { a _ { 1 1 } } } \\ { 0 } & { 0 } & { - \frac { a _ { 2 3 } } { a _ { 2 2 } } } & { \cdots } & { - \frac { a _ { 2 , n - 1 } } { a _ { 2 2 } } } & { - \frac { a _ { 2 n } } { a _ { 2 2 } } } \\ { \vdots } & { \vdots } & { \vdots } & { \ddots } & { \vdots } & { \vdots } \\ { 0 } & { 0 } & { 0 } & { 0 } & { 0 } & { - \frac { a _ { n - 1 , n } } { a _ { n - 1 , n - 1 } } } \\ { 0 } & { 0 } & { 0 } & { 0 } & { 0 } & { 0 } \end{array} \right) \]

迭代格式

向量（矩阵）形式

高斯-赛德尔迭代的矩阵格式为：

\[ \mathbf{x}_k = \mathbf{L}\mathbf{x}_k + \mathbf{U}\mathbf{x}_{k-1} + \mathbf{D}^{-1}\mathbf{b} \]
  • 其中 \(\mathbf{x}_k\) 表示第 \(k\) 次迭代得到的近似解。
  • \(\mathbf{L}\mathbf{x}_k\)：用的是当前轮已经更新的分量
  • \(\mathbf{U}\mathbf{x}_{k-1}\)：用的是上一轮的分量（尚未更新）

分量形式

对于第 \(i\) 个分量，第 \(k\) 轮迭代时，更新公式为： $$ x_i^{(k)} = \frac{1}{a_{ii}} \left( b_i - \sum_{j=1}^{i-1} a_{ij} x_j^{(k)} - \sum_{j=i+1}^{n} a_{ij} x_j^{(k-1)} \right) $$

  • 第一项 \(b_i\) 是常数向量的第 \(i\) 个分量
  • \(j < i\)：用的是 本轮刚刚更新的新分量
  • \(j > i\)：用的是 上一轮的旧分量

算法步骤

  1. 初始化
    给定初始猜测 \(x^{(0)}\)（通常可用全零或任意向量），记作 x0。

  2. 主循环迭代

对 \(k = 1, 2, \ldots\) 直到收敛：

  • 对每个分量 \(i=1\) 到 \(n\)：

    • 计算当前 \(x_i^{(k)}\)：
    • 先累加 \(j=1\) 到 \(i-1\) 的 \(a_{ij} x_j^{(k)}\)（已经更新过的新解）
    • 再累加 \(j=i+1\) 到 \(n\) 的 \(a_{ij} x_j^{(k-1)}\)（上一轮的旧解）
    • 代入公式得到新的 \(x_i^{(k)}\)

  • 计算新解与旧解的距离（欧氏范数）判断收敛：\(\Delta x = \sqrt{ \sum_{i=1}^n (x_i^{(k)} - x_i^{(k-1)})^2 }\)如果 \(\Delta x < \varepsilon\)，则认为已收敛，停止迭代。

  • 输出结果：最终 \(x^{(k)}\) 即为近似解。

Fortran数值实现

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

    ! ---------- 增加收敛性判据判断 ----------
    real*8 :: s
    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)') 'Gauss-Seidel iteration may not converge!'
        return
    end if
    ! ----------------------------------------

    write(*,501)
    write(*,502)
    ! 迭代之前两值都设为初值
    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)/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,"Gauss_Seidel Iterative",//)
    502 format("  Iter", 5x, "  x1", 12x, "  x2", 12x, "  x3")
    503 format(i5, 3(3x, f12.8))

end subroutine gauss_seidel

数值案例

\[ \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 gauss_seidel(a, b, x, x0, n)

end program main

输出：

                  Gauss_Seidel Iterative


  Iter       x1              x2              x3
    1     0.90000000     0.79000000     0.91600000
    2     0.97900000     0.98110000     0.99244000
    3     0.99811000     0.99829900     0.99931960
    4     0.99982990     0.99984691     0.99993876
    5     0.99998469     0.99998622     0.99999449
    6     0.99999862     0.99999876     0.99999950
    7     0.99999988     0.99999989     0.99999996
    8     0.99999999     0.99999999     1.00000000

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

参考文献

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