求解线性方程组

在线性代数中，有一类经典问题，就是求解线性方程组，我们熟知的解法有高斯消元法，但是高斯消元法属于直接求解的方法，不适合编程计算，因此引入更适合计算机求解的迭代法。

以线性方程组：

$ \left[ \begin{matrix} 10 & -1 & -2 \ -1 & 10 & -2 \ -1 & -1 & 5 \end{matrix} \right] x = \left[ \begin{matrix} 7.2 \ 8.3 \ 4.2 \end{matrix} \right] $ 为例

Jacobi迭代法

考虑线性方程组$Ax=b$

用L和U分别表示严格下三角矩阵和严格上三角矩阵

可以利用迭代公式

$x^{k+1}=D^{-1}(b-(L+U)x^{k}))$

进行迭代求解，可以通过精度控制迭代次数，还需要控制迭代是否收敛，不然迭代次数再多也无法求解。

迭代法需要给定解列向量的初值，因此初值的选取也很重要。

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23


import numpy as np

def Jacobi(A, B, N, x):
    # Ax=B N为迭代次数
    # 获得对角矩阵的数组
    D = np.diag(A)
    # 获得L+U矩阵，减去对角矩阵即可
    # diagflat为D数组转为对角矩阵
    R = A - np.diagflat(D)
    # R = L + U
    print(R/D)
    for i in range(N):
        # 迭代过程，可以直接/D表示逆矩阵
        x = (B - np.dot(R, x))/D 
    return x

A = np.array([[10, -1, -2], [-1, 10, -2], [-1, -1, 5]]    )
B = np.array([7.2, 8.3, 4.2])
# 设定迭代初始值
x = np.array([1, 1, 1])
answer = Jacobi(A, B, 30, x)

print(answer)

求解结果为 x = [1.1, 1.2, 1.3]

Gauss Seidel迭代法

Gauss Seidel法使用了另一种迭代格式，获得了快的收敛速度

首先我们利用迭代求解的特性，新得到的值总会比老值更优，因此使用新值来进行迭代，即

$x^{k+1}=D^{-1}(b+Lx^{k+1}+Ux^k)$

化简

$x^{k+1}=(D-L)^{-1}(Ux^k+b)$

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26


import numpy as np

def GaussSeidel(A, B, N, x):
    # 下三角矩阵
    L = np.tril(A)
    # 严格上三角矩阵
    U = A - L
    # 获得对角矩阵
    D = np.diagflat(np.diag(A))
    LL = L - D
    # x^{k+1} = (D + L)^-1 (b - Ux^{k})
    # 输出迭代矩阵
    print(np.dot(np.linalg.inv(D-LL), U))
    for i in range(N):
        # 使用linalg包中的inv函数求出L矩阵的逆矩阵
        # 使用上述迭代公式求解
        x = np.dot(np.linalg.inv(L), B - np.dot(U, x))
    return x

A = np.array([[10, -1, -2], [-1, 10, -2], [-1, -1, 5]]    )
B = np.array([7.2, 8.3, 4.2])
# 设定迭代初始值
x = np.array([1, 1, 1])
answer = GaussSeidel(A, B, 30, x)

print(answer)

解得x = [1.1, 1.2, 1.3]

迭代法求解线性方程组

求解线性方程组

Jacobi迭代法

Gauss Seidel迭代法