线性方程组迭代算法

一般的线性方程组矩阵形式如下

AX=b

迭代解法的思想 $x^{(0)}$ $x^{(1)}$ ，逐步迭代至于收敛。

How to 设计迭代公式？

1 雅克比迭代法（Jacobi）

\begin{cases} a_{11}x_1+a_{12}x_2+&\dots+a_{1n}x_n=b_1 \\a_{21}x_1+a_{22}x_2+&\dots+a_{2n}x_n=b_2 \\ &\dots \\a_{n1}x_1+a_{n2}x_2+&\dots+a_{nn}x_n=b_n \end{cases}

雅克比迭代法建设系数矩阵的对角元素不为0，因此，由上式可得

\begin{cases} x_1=-\frac{1}{a_{11}}(a_{12}x_2+&\dots+a_{1n}x_n-b_1) \\x_2=-\frac{1}{a_{22}}(a_{21}x_2+&\dots+a_{2n}x_n-b_2) \\ &\dots \\x_n=-\frac{1}{a_{nn}}(a_{n1}x_1+&\dots+a_{n,n-1}x_{n-1}-b_n) \end{cases}

将上式写成矩阵形式为

X=-D^{-1}(L+U)X+D^{-1}b

其中，

D= \begin{bmatrix} a_{11}&&&\\&a_{22}&&&\\ &&\ddots& \\&&&a_{nn} \end{bmatrix} L= \begin{bmatrix} 0&\ &\ &\ &\ \\ a_{21}&0&\ &\ &\ \\ a_{31}&a_{32}&0&\ &\ \\ \vdots&\vdots&\ &\ddots\\ a_{n1}&a_{n2}&a_{n3}&\cdots&0 \end{bmatrix} U= \begin{bmatrix} 0&a_{12}&a_{13}&\ &a_{1n}\\ \ &0&a_{23}&\cdots&a_{2n}\\ \ &\ &\ddots&\ &\vdots\\ \ &\ &\ &0&a_{n-1,n}\\ \ &\ &\ &\ &0 \end{bmatrix}

$B=-D^{-1}(L+U), d=D^{-1}b$ ，则迭代公式变为

X=BX+d

$B$ 是迭代矩阵，在进行迭代计算时，式(3)变为

\begin{cases} x_1^{(k+1)}=-\frac{1}{a_{11}}(a_{12}x_2^{(k)}+&\dots+a_{1n}x_n^{(k)}-b_1) \\x_2^{(k+1)}=-\frac{1}{a_{22}}(a_{21}x_1^{(k)}+&\dots+a_{2n}x_n^{(k)}-b_2) \\ &\dots \\x_n^{(k+1)}=-\frac{1}{a_{nn}}(a_{n1}x_1^{(k)}+&\dots+a_{n,n-1}x_{n-1}^{(k)}-b_n) \end{cases}

Example


x
22
1
import numpy as np
2
3
4
def jacobi(A, b, x0, tol):
5
    D = np.diag(np.diag(A)) #对角矩阵
6
    L = np.tril(A,-1) #下三角矩阵的下一格
7
    U = np.triu(A, 1) #上三角矩阵的上一格
8
    B = -np.dot(np.linalg.inv(D),(L+U))
9
    d = np.dot(np.linalg.inv(D), b)
10
    x = x0
11
    while np.linalg.norm(np.dot(A,x)-b) > tol:
12
        x = np.dot(B, x) + d
13
    return x
14
15
16
if __name__ == '__main__':
17
    A = np.array([[2, 1], [1, -1]], dtype=np.float)
18
    b = np.array([[3], [0]], dtype=np.float)
19
    x0 = np.array([[0], [0]])
20
    tol = 0.01
21
    x = jacobi(A, b, x0, tol)
22
    print(x)

2 高斯-赛德尔迭代法

从雅克比到高斯-赛德尔

$k+1$ $x_i^{(k+1)}$ $x_1^{(k+1)}, x_2^{(k+1)},\cdots,x_{i-1}^{(k+1)}$ $x_i^{(k+1)}$ 时，雅克比迭代法并没有利用这些新计算出来的值，如果这些值能得以使用，则可以提高计算速度，这就是高斯-赛德尔迭代法的来源。

\begin{cases} x_1^{(k+1)}=-\frac{1}{a_{11}}(a_{12}x_2^{(k)}+a_{13}x_3^{(k)}+&\dots+a_{1n}x_n^{(k)}-b_1) \\x_2^{(k+1)}=-\frac{1}{a_{22}}(a_{21}x_1^{(k+1)}+a_{13}x_3^{(k)}+&\dots+a_{2n}x_n^{(k)}-b_2) \\ &\dots \\x_n^{(k+1)}=-\frac{1}{a_{nn}} (a_{n1}x_1^{(k+1)}+a_{13}x_3^{(k+1)}+&\dots+a_{n,n-1}x_{n}^{(k+1)}-b_n) \end{cases}

$k+1$ 部分移项到左边后，用方程组用矩阵的形式表示为

(D+L)X^{(k+1)}=-UX^{(k)}+b

写成迭代方程为

X^{(k+1)}=GX^{(k)}+d

$G=-(D+L)^{-1}U, d=(D+L)^{-1}b$


xxxxxxxxxx
10
1
def gsdddy(A, b, x0, tol):
2
    D = np.diag(np.diag(A)) #对角矩阵
3
    L = np.tril(A,-1) #下三角矩阵的下一格
4
    U = np.triu(A, 1) #上三角矩阵的上一格
5
    G = np.dot(-np.linalg.inv(D+L), U)
6
    d = np.dot(np.linalg.inv(D+L), b)
7
    x = x0
8
    while np.linalg.norm(np.dot(A, x)-b) > tol:
9
        x = np.dot(G,x) + d
10
    return x