GaussSeidel迭代 Jacobi迭代 收斂性

• 迭代法求解線性方程組
• $GaussSeidel$
• • 公式推導
  • • 矩陣形式
    • 非矩陣形式
  • 代碼示例
• $Jacobi$
• • 公式推導
  • • 矩陣形式
    • 非矩陣形式
  • 代碼示例
• 收斂性

迭代法求解線性方程組

目的是對于一個可逆矩陣$A$，利用迭代法求解

$$Ax=b.$$

$GaussSeidel$ 和 $Jacobi$ 都是將 $A$ 矩陣分成 $D+L+U$ 矩陣的形式，構造出$x=Gx+d$的迭代矩陣，來迭代求解出$x$，理想情況下當 $x==Gx+d$ 時算法收斂，且此時$x=A^{?1}b$。其中$DLU$ 分別為 $A$ 矩陣的對角線矩陣，下三角矩陣和上三角矩陣。

在上述中介紹到兩種方法都是通過定義迭代矩陣 $x=Gx+d$ 來求解 $x$ 的。因此接下來分別介紹，GaussSeidel 和Jacobi分別時怎樣定義迭代矩陣的。

$GaussSeidel$

公式推導

矩陣形式

$$A=D+L+U$$
$$Ax=b$$
$$(D+L+U)x=b$$
$$(D+L)x=?Ux+b$$
$$x=?{(D+L)}^{?1}Ux+{(D+L)}^{?1}b$$
$$?x=Gx+d$$
由上述推導可知，在$GaussSeidel$中 $G=?{(D+L)}^{?1}U$， $d={(D+L)}^{?1}b$

非矩陣形式

$$x_i^k=\frac{1}{a_{ii}}(b?\sum _{0<j<i}{a}_{ij}{x}_j^k?\sum _{j>i}{a}_{ij}{x}_j^{k?1}))$$
其中 $x_i^k$ 代表 $x$ 的第$i$個分量在第$k$次迭代中的結果

代碼示例

使用非矩陣形式實現在稀疏矩陣csr_matrix下實現的$GaussSeidel$。

# A csr_matrix import sys import copy import numpy as np import networkx as nx from scipy.sparse import spdiags, tril, triu, coo_matrix, csr_matrixdef gauss_seidel(A, b, x0 = None, maxiter=20):# x0 = spsolve(A, b)if x0 is None:x0 = np.array( [0.0]*A.shape[0] )for _ in range(maxiter):n, _, indptr, indices, data = A.shape[0], A.nnz, A.indptr, A.indices, A.datafor row in range(n):aii, sumi = 0.0, 0.0for j in range(indptr[row], indptr[row+1]):if row == indices[j]:aii = 1.0/data[j]else:sumi += data[j] * x0[ indices[j] ]x0[row] = (aii * (b[row] - sumi) )return x0 if __name__ == '__main__': A = csr_matrix(([10,-2,-1, -2,10,-1, -1,-2,5], ([0,0,0,1,1,1,2,2,2], [0,1,2,0,1,2,0,1,2]) ) , shape=(3,3), dtype=np.float32)x = np.ones(A.shape[0])b = A.dot(x)xk = gauss_seidel(A, b)print(x)print(xk)

輸出如下

$Jacobi$

公式推導

矩陣形式

$$A=D+L+U$$
$$Ax=b$$
$$(D+L+U)x=b$$
$$Dx=?（L+U)x+b$$
$$x=?D^{?1}(L+U)x+D^{?1}b$$
$$?x=Gx+d$$
由上述推導可知，在Jacobi中 $G=?D^{?1}(L+U)$， $d=D^{?1}b$

非矩陣形式

$$x_i^k=\frac{1}{a_{ii}}(b?\sum _{j!=i}{a}_{ij}{x}_j^{k?1})$$

代碼示例

使用矩陣形式實現在稀疏矩陣csr_matrix下實現的$Jacobi$。

# A csr_matrix import sys import copy import numpy as np from scipy.sparse import spdiags, tril, triu, coo_matrix, csr_matrix def jacobi(A, b, x0 = None, maxiter=20):if x0 is None:x0 = np.array( [0.0]*A.shape[0] )D = spdiags(A.diagonal(), np.array([0]), A.shape[0], A.shape[1]) # default with <class 'scipy.sparse.dia.dia_matrix'> L = tril(A) - DU = triu(A) - D dia = generate_dia_inv(D)B = -( dia ).dot(L+U)f = ( dia ).dot(b)for _ in range(maxiter):x0 = B.dot(x0) + freturn x0 def generate_dia_inv(mat):res_dia = mat.diagonal()res_dia = 1.0/res_diares_mat = csr_matrix((res_dia.shape[0], res_dia.shape[0]))res_mat.setdiag(res_dia)return res_mat if __name__ == '__main__': A = csr_matrix(([10,-2,-1, -2,10,-1, -1,-2,5], ([0,0,0,1,1,1,2,2,2], [0,1,2,0,1,2,0,1,2]) ) , shape=(3,3), dtype=np.float32)x = np.ones(A.shape[0])b = A.dot(x)xk = jacobi(A, b)print(x)print(xk)

輸出如下

收斂性

因為兩種方法的根本思想是一致的，都是通過迭代公式 $x=Gx+d$ 來求解 $x$。所以直接通過此迭代公式來看兩種方法收斂的條件。
令誤差向量 $e^i=x^i?x$，其中 $x$ 是 $Ax=b$的精確解， $e^i$是$GaussSeidel$或者$Jacoci$第$i$次迭代所求解。

收斂性證明

由于 $x^i=G{x}^{i?1}+d$

且 $x=Gx+d$

所以有，$e^i=G{e}^{i?1}$

且 $e^{i?1}=G{e}^{i?2}$

$...$

因此，$e^i=G^i{e}^0$

當迭代法收斂時，就意味著 $\underset{i\to \infty }{\lim }{e}^i\to 0$

由特征值和特征向量的定義可知$G{\zeta }_i={\lambda }_i{\zeta }_i$，且特征向量之間線性無關

因此$e^0=\Sigma {\delta }_k{\zeta }_k$, 可以這樣表示，且 ${\delta }_k$是定值

則$e^1=G{e}^0=G\sum {\delta }_k{\zeta }_k=\sum {\delta }_{k}G{\zeta }_k=\sum {\delta }_k{\lambda }_k{\zeta }_k$，由 $G{\zeta }_i={\lambda }_i{\zeta }_i$得。

則$e^2=\sum {\delta }_k{\lambda }_k^2{\zeta }_k$ , $e^3=\sum {\delta }_k{\lambda }_k^3{\zeta }_k$ …

那么 $\underset{i\to \infty }{\lim }{e}^i=\underset{i\to \infty }{\lim }\sum {\lambda }_k^i{\delta }_k{\zeta }_k$

其中${\delta }_k{\zeta }_k$是定值

則當 $abs({\lambda }_k)>1?\underset{i\to \infty }{\lim }{e}^i\to \infty$

當 $abs({\lambda }_k)<1?\underset{i\to \infty }{\lim }{e}^i\to 0$

綜上，當迭代矩陣$G$的所有特征值 $abs({\lambda }_k)<1$時， 該迭代方法收斂。

總結

以上是生活随笔為你收集整理的GaussSeidel迭代 Jacobi迭代及其收敛性的全部內容，希望文章能夠幫你解決所遇到的問題。

如果覺得生活随笔網站內容還不錯，歡迎將生活随笔推薦給好友。