学习笔记 - 高斯-赛德尔迭代法

最新推荐文章于 2024-05-01 07:58:56 发布

K5niper

最新推荐文章于 2024-05-01 07:58:56 发布

阅读量2.1w

点赞数 16

CC 4.0 BY-SA版权

本文链接：https://blog.youkuaiyun.com/zhaoyin214/article/details/88619827

高斯-赛德尔迭代法是数值代数中用于求解线性方程组的迭代方法，尤其适用于大型稀疏方程组。在对角占优或对称正定矩阵情况下保证收敛。该方法通过下三角矩阵的前向替换计算更新值，仅需保存一个向量。收敛性依赖于矩阵特性，满足特定条件时确保收敛。

摘要生成于 C知道，由 DeepSeek-R1 满血版支持，前往体验 >

高斯-赛德尔迭代法

高斯-赛德尔迭代法（Gauss–Seidel method，李伯曼法（Liebmann method），逐次位移法（successive displacement））是数值代数中的一种迭代法，用于求解线性方程组（linear system of equations）。

当矩阵为对角占优（diagonally dominant）或者对称且正定时，高斯-赛德尔迭代法能够保证收敛。
求解大型稀疏方程组时，通常采用迭代法。迭代法能够在内存和运算两方面充分利矩阵的稀疏性。

描述

高斯-赛德尔法用于迭代求解 $n$ 维线性方程组：

$Ax=b\mathbf{Ax} = \mathbf{b}$

$A=[a11a12⋯a1na21a22⋯a2n⋮⋮⋱⋮an1an2⋯ann],x=[x1x2⋮xn],b=[b1b2⋮bn]\mathbf{A} = \begin{bmatrix} a_{11} & a_{12} & \cdots & a_{1n} \\ a_{21} & a_{22} & \cdots & a_{2n} \\ \vdots & \vdots & \ddots & \vdots \\ a_{n1} & a_{n2} & \cdots & a_{nn} \end{bmatrix}, \quad \mathbf {x} = \begin{bmatrix} x_{1} \\ x_{2} \\ \vdots \\ x_{n} \end{bmatrix}, \quad \mathbf {b} = \begin{bmatrix} b_{1} \\ b_{2} \\ \vdots \\ b_{n} \end{bmatrix}$

迭代更新定义为：

$L∗x(k+1)=b−Ux(k)\mathbf{L}_{\ast}\mathbf{x}^{(k+1)} = \mathbf{b} - \mathbf{U} \mathbf{x}^{(k)}$

其中， $A\mathbf{A}$ 分解为下三角矩阵（lower triangular component） $L∗\mathbf{L}_{\ast}$ 和严格上三角矩阵（strictly upper triangular component） $U\mathbf{U}$ ：

$A=L∗+U\mathbf{A} = \mathbf{L}_{\ast} + \mathbf{U}$

$L∗=[a110⋯0a21a22⋯0⋮⋮⋱⋮an1an2⋯ann],U=[0a12⋯a1n00⋯a2n⋮⋮⋱⋮00⋯0]\mathbf{L}_{\ast} = \begin{bmatrix} a_{11} & 0 & \cdots & 0 \\ a_{21} & a_{22} & \cdots & 0 \\ \vdots & \vdots & \ddots & \vdots \\ a_{n1} & a_{n2} & \cdots & a_{nn} \end{bmatrix}, \quad \mathbf{U} = \begin{bmatrix} 0 & a_{12} & \cdots & a_{1n} \\ 0 & 0 & \cdots & a_{2n} \\ \vdots & \vdots & \ddots & \vdots \\ 0 & 0 & \cdots & 0 \end{bmatrix}$

线性方程组可改写为：

$L∗x=b−Ux\mathbf{L}_{\ast} \mathbf{x} = \mathbf{b} - \mathbf{U} \mathbf{x}$

高斯-赛德尔法用第 $k$ 次迭代 $x\mathbf{x}$ 的值（等式右边）计算第 $k + 1$ 次迭代 $x\mathbf{x}$ 的值（等式左边），即

$x(k+1)=L∗−1(b−Ux(k))\mathbf{x}^{(k+1)} = \mathbf{L}_{\ast}^{-1} \left( \mathbf{b} - \mathbf{U} \mathbf{x}^{(k)} \right)$

注意 $L∗\mathbf{L}_{\ast}$ 为三角矩阵， $x(k+1)\mathbf{x}^{(k+1)}$ 的元素通过前向替换（forward substitution）串行计算：

$xi(k+1)=1aii(bi−∑j=1i−1aijxj(k+1)−∑j=i+1naijxj(k)),i=1,2,…,nx_{i}^{(k+1)} = \frac{1}{a_{ii}} \left( b_{i} - \sum_{j=1}^{i-1} a_{ij} x_{j}^{(k+1)} - \sum_{j=i+1}^{n} a_{ij} x_{j}^{(k)} \right), \quad i=1, 2, \dots, n$

当更新增量 $Δx=x(k+1)−x(k)\Delta \mathbf{x} = \mathbf{x}^{(k+1)} - \mathbf{x}^{(k)}$ 小于某个容限（残差，sufficiently small residual）时，迭代停止。

计算 $x_{i}^{(k+1)}$ 使用的是第 $k + 1$ 次迭代中已更新的 $x1(k+1),⋯ ,xi−1(k+1)x_{1}^{(k+1)}, \cdots, x_{i - 1}^{(k+1)}$ 和尚未更新的 $xi+1(k),⋯ ,xn(k)x_{i + 1}^{(k)}, \cdots, x_{n}^{(k)}$ ，因此计算过程只需要保存一个向量。

敛散性

高斯-赛德尔法的敛散性与矩阵 $A\mathbf{A}$ 有关，当 $A\mathbf{A}$ 满足以下两个条件之一时收敛：

（1） $A\mathbf{A}$ 是对称正定的

（2） $A\mathbf{A}$ 严格或不可约对角占优

但即使不满足这些条件，高斯-塞德尔法有时也会收敛。

算法

伪代码：

输入： $A\mathbf{A}$ ， $b\mathbf{b}$

输出： $ϕ\phi$

给定解的猜测值 $ϕ\phi$ 作为初始值：

REPEAT UNTIL 收敛：

FOR i FROM 1 UNTIL n DO

    $\sigma \leftarrow 0$
    
    FOR j FROM 1 UNTIL n DO
        IF j ≠ i THEN
            $\sigma \leftarrow \sigma + a_{ij} \phi_{j}$
        END IF
    END (j-loop)
    
    $\phi_{i} \leftarrow \frac{1}{a_{ii}} ( b_{i} - \sigma )$
    
END (i-loop)

验证是否收敛

END (REPEAT)

实现

$xi(k+1)=1aii(bi−∑j<iaijxj(k+1)−∑j>iaijxj(k)),i,j=1,2,…,nx_{i}^{(k+1)} = \frac {1}{a_{ii}} \left( b_{i} - \sum_{j<i} a_{ij} x_{j}^{(k+1)} - \sum_{j>i} a_{ij} x_{j}^{(k)} \right), \quad i, j = 1, 2, \ldots, n$

import numpy as np

def gauss_seidel(A, b, x0, iters):
    
    m, n = A.shape
    x = x0
    if m != n:
        raise ValueError("A must be a square matrix.")
    
    for k in range(iters):
        
        tol = 0
        for i in range(n):
            
            x_i_old = x[i]
            sum_new = (A[i, : i] * x[: i]).sum()
            sum_old = (A[i, i + 1 :] * x[i + 1 :]).sum()
            x[i] = 1 / A[i, i] * (b[i] - sum_new - sum_old)
            
            tol += abs(x[i] - x_i_old)
            
        if tol / n < 1e-8:
            break
        print("Iteration {0}: x = {1}, tol = {2}".format(k, x, tol))
        
    return x

ITERATION_LIMIT = 1000

# initialize the matrix
A = np.array([[10., -1., 2., 0.],
              [-1., 11., -1., 3.],
              [2., -1., 10., -1.],
              [0., 3., -1., 8.]])
# initialize the RHS vector
b = np.array([6., 25., -11., 15.])

print("System of equations:")
for i in range(A.shape[0]):
    row = ["{0:3g}*x{1}".format(A[i, j], j + 1) for j in range(A.shape[1])]
    print("[{0}] = [{1:3g}]".format(" + ".join(row), b[i]))

x = np.zeros_like(b)

x = gauss_seidel(A, b, x0=x, iters=ITERATION_LIMIT)

print("Solution: {0}".format(x))
error = np.dot(A, x) - b
print("Error: {0}".format(error))

System of equations:
[ 10*x1 +  -1*x2 +   2*x3 +   0*x4] = [  6]
[ -1*x1 +  11*x2 +  -1*x3 +   3*x4] = [ 25]
[  2*x1 +  -1*x2 +  10*x3 +  -1*x4] = [-11]
[  0*x1 +   3*x2 +  -1*x3 +   8*x4] = [ 15]
Iteration 0: x = [ 0.6         2.32727273 -0.98727273  0.87886364], tol = 4.793409090909091
Iteration 1: x = [ 1.03018182  2.03693802 -1.0144562   0.98434122], tol = 0.8531775826446281
Iteration 2: x = [ 1.00658504  2.00355502 -1.00252738  0.99835095], tol = 0.08291831672614614
Iteration 3: x = [ 1.00086098  2.00029825 -1.00030728  0.99984975], tol = 0.012699738431181884
Iteration 4: x = [ 1.00009128  2.00002134 -1.00003115  0.9999881 ], tol = 0.0014610924360077826
Iteration 5: x = [ 1.00000836  2.00000117 -1.00000275  0.99999922], tol = 0.00014260125078868757
Iteration 6: x = [ 1.00000067  2.00000002 -1.00000021  0.99999996], tol = 1.212975945419359e-05
Iteration 7: x = [ 1.00000004  1.99999999 -1.00000001  1.        ], tol = 8.839586538300637e-07
Iteration 8: x = [ 1.  2. -1.  1.], tol = 5.8961691529191285e-08
Solution: [ 1.  2. -1.  1.]
Error: [ 5.57651703e-11 -1.38831879e-09  2.94534175e-10  0.00000000e+00]