基于python实现的Jacobi迭代法和Gauss-Seidel迭代法

之前的文章介绍过Gauss消元法、Doolittle和Crout分解等解线性方程组的直接方法。但是直接法不是在任何时候都可用的，它有以下缺点：

• 在有舍入误差的情况下，直接法只能得到方程的近似解。
• 如果待求解问题规模很大，直接法求解的计算量是很大的，更糟的是随着计算量的增大舍入误差积累也会越多，从而得到的解误差更大。所以直接法适合小规模稠密型方程组的求解。
• 另外，直接法的程序设计也比较复杂。

---

迭代法解线性方程组

对于方程组 $Ax=b$ 一般可以将其写成如下形式
$$x=Mx+g$$建立上述方程组的迭代公式：
$$x^{(k+1)}=M{x}^{(k)}+g,k=0,1,2,\dots$$这个迭代公式有什么用呢？

当向量序列 $x^{(k)}$ 收敛于 $x^{\ast }$ 时，有
$$x^{\ast }=M{x}^{\ast }+g$$可以发现 $x^{\ast }$ 即为方程组 $x=Mx+g$ 的解，也就是说，当 $x^{(k)}$ 收敛时，$x^{(k)}$ 即为方程组的解。

以上就是迭代法解方程组的基本原理，只要$x^{(k)}$收敛，就可以通过k步迭代得到方程组的近似解。

下面介绍两种经典的迭代法：Jacobi迭代法和Gauss-Seidel迭代法。

---

Jacobi迭代法

假设有下列线性方程组
$$\{\begin{array}{l}{a}_{11}{x}_1+a_{12}{x}_2+a_{13}{x}_3=b_1\\ a_{21}{x}_1+a_{22}{x}_2+a_{23}{x}_3=b_2\\ a_{31}{x}_1+a_{32}{x}_2+a_{33}{x}_3=b_3\end{array}$$
将其写成下面形式
$$\{\begin{array}{l}{x}_1=\frac{1}{a_{11}}(b_1-a_{12}{x}_2-a_{13}{x}_3)\\ x_2=\frac{1}{a_{22}}(b_2-a_{21}{x}_1-a_{23}{x}_3)\\ x_3=\frac{1}{a_{33}}(b_3-a_{31}{x}_1-a_{32}{x}_2)\end{array}$$
建立迭代公式
$$\{\begin{array}{l}{x}_1^{(k+1)}=\frac{1}{a_{11}}(b_1-a_{12}{x}_2^{(k)}-a_{13}{x}_3^{(k)})\\ x_2^{(k+1)}=\frac{1}{a_{22}}(b_2-a_{21}{x}_1^{(k)}-a_{23}{x}_3^{(k)})\\ x_3^{(k+1)}=\frac{1}{a_{33}}(b_3-a_{31}{x}_1^{(k)}-a_{32}{x}_2^{(k)})\end{array}$$
对于每一个 $x_i^{(k+1)}$ ，有如下公式
$$\begin{gathered} x_i^{(k+1)}=\frac{1}{a_{ii}}(b_i-\sum _{j=1}^{i-1}{a}_{ij}{x}_j^{(k)}-\sum _{j=i+1}^n{a}_{ij}{x}_j^{(k)}) \\ i=1,2,\cdots ,n,k=0,1,\cdots \end{gathered}$$然后我们就可根据这个公式解线性方程组了。

但是为了便于理论分析，我们把它写成矩阵的形式。
假设方程组为 $Ax=b$，将 $A$ 分裂为
$$A=D-L-U$$其中 $D$ 为对角矩阵 $diag(a_{11},a_{22},\dots ,a_{nn})$，$-L$ 和 $-U$ 分别如下
$$\begin{gathered} -L=\left(\begin{array}{ccccc}0& & & & \\ a_{21}& 0& & & \\ a_{31}& a_{32}& 0& & \\ ⋮& ⋮& \ddots & \ddots & \\ a_{n1}& a_{n2}& \dots & a_{nn-1}& 0\end{array}\right) \\ -U=\left(\begin{array}{ccccc}0& a_{12}& a_{13}& \dots & a_{1n}\\ & 0& a_{23}& \dots & a_{2n}\\ & & \ddots & \ddots & ⋮\\ & & & \ddots & a_{n-1n}\\ & & & & 0\end{array}\right) \end{gathered}$$则方程组和迭代公式可写为
$$\begin{gathered} (D-L-U)x=b \\ x=D^{-1}(L+U)x+D^{-1}b \\ {x}^{(k+1)}=D^{-1}(L+U)x^{(k)}+D^{-1}b(迭代公式) \\ k=0,1,\cdots \end{gathered}$$
可以验证矩阵形式的迭代公式和上面大括号方程组表示的迭代公式一致的。

---

下面根据公式 $x^{(k+1)}=D^{-1}(L+U)x^{(k)}+D^{-1}b$ 编写代码

算法代码

def Jfunction(A,b,k):
    n = A.shape[1]
    D = np.eye(n)
    D[np.arange(n),np.arange(n)] = A[np.arange(n),np.arange(n)]
    LU = D - A
    X = np.zeros(n)
    for i in range(k):
        D_inv = np.linalg.inv(D)
        X =  np.dot(np.dot(D_inv,LU),X) + np.dot(D_inv,b)
    return X

实验
解下列方程组

方程组的精确解为 $x=(1,1,1)$

A = np.array([
    [10,3,1],
    [2,-10,3],
    [1,3,10]
])
b = np.array([14,-5,14])
k = 6
X = Jfunction(A,b,k);

实验结果

---

Gauss-Seidel迭代法

在J迭代法中在计算 $x_i^{(k+1)}$时，前 $i-1$ 个分量已经计算完成了，如果可以利用这些信息，可能加速迭代法的收敛。因此我们对迭代公式做如下修改
$$\{\begin{array}{l}{x}_1^{(k+1)}=\frac{1}{a_{11}}(b_1-a_{12}{x}_2^{(k)}-a_{13}{x}_3^{(k)})\\ x_2^{(k+1)}=\frac{1}{a_{22}}(b_2-a_{21}{x}_1^{(k+1)}-a_{23}{x}_3^{(k)})\\ x_3^{(k+1)}=\frac{1}{a_{33}}(b_3-a_{31}{x}_1^{(k+1)}-a_{32}{x}_2^{(k+1)})\end{array}$$这就是Gauss-Seidel迭代法的迭代公式
对于每一个 $x_i^{(k+1)}$ ，有如下公式
$$\begin{gathered} x_i^{(k+1)}=\frac{1}{a_{ii}}(b_i-\sum _{j=1}^{i-1}{a}_{ij}{x}_j^{(k+1)}-\sum _{j=i+1}^n{a}_{ij}{x}_j^{(k)}) \\ i=1,2,\cdots ,n,k=0,1,\cdots \end{gathered}$$上述公式也可以写成矩阵形式，我们可以从J迭代法的矩阵形式推出G-S迭代法的矩阵形式

按照J迭代法的迭代公式有
$$\begin{gathered} x^{(k+1)}=D^{-1}(L+U)x^{(k)}+D^{-1}b \\ D{x}^{(k+1)}=L{x}^{(k)}+U{x}^{(k)}+b \end{gathered}$$不难发现 $L{x}^{(k)}$ 对应了G-S迭代公式中的 ${\sum }_{j=1}^{i-1}{a}_{ij}{x}_j^{(k+1)}$ 部分，因此按照G-S迭代法的迭代公式，上面迭代公式应该写成
$$D{x}^{(k+1)}=L{x}^{(k+1)}+U{x}^{(k)}+b$$整理，得
$$x^{(k+1)}=(D-L{)}^{-1}U{x}^{(k)}+(D-L{)}^{-1}b$$

---

下面根据公式 $x^{(k+1)}=(D-L{)}^{-1}U{x}^{(k)}+(D-L{)}^{-1}b$ 编写代码

算法代码

def G_Sfunction(A,b,k):
    n = A.shape[1]
    D = np.eye(n)
    D[np.arange(n),np.arange(n)] = A[np.arange(n),np.arange(n)]
    LU = D - A
    L = np.tril(LU)
    U = np.triu(LU)
    D_L = D - L
    X = np.zeros(n)
    for i in range(k):
        D_L_inv = np.linalg.inv(D_L)
        X =  np.dot(np.dot(D_L_inv,U),X) + np.dot(D_L_inv,b)
    return X

实验
解下列方程组（一样的例题）

方程组的精确解为 $x=(1,1,1)$

A = np.array([
    [10,3,1],
    [2,-10,3],
    [1,3,10]
])
b = np.array([14,-5,14])
k = 6
X = G_Sfunction(A,b,k);

实验结果

可以发现G-S迭代法在第4次迭代完成的时候就已经达到了J迭代法第六次迭代完成的效果了。

---

迭代法的特点

• 只能求的近似解
• 可以在迭代的过程中对误差进行矫正
• 如果系数矩阵是稀疏的，即使面对大规模的问题，也可能通过较小的计算量得到方程组的解
• 原理简单方便实现