雅可比法

在数值线性代数中，雅可比法（）是一种解对角元素几乎都是各行和各列的绝对值最大的值的线性方程组的算法。求解出每个对角元素并插入近似值。不断迭代直至收敛。这个算法是雅可比矩阵的精简版。方法的名字来源于德国数学家卡尔·雅可比。

描述
给定一个n×n的线性方程组

:A\mathbf x = \mathbf b

其中:

: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}, \qquad \mathbf{x} = \begin{bmatrix} x_{1} \\ x_2 \\ \vdots \\ x_n \end{bmatrix} , \qquad \mathbf{b} = \begin{bmatrix} b_{1} \\ b_2 \\ \vdots \\ b_n \end{bmatrix}.

A 可以分解成对角部分D和剩余部分R:

:A=D+R \qquad \text{其 中} \qquad D = \begin{bmatrix} a_{11} & 0 & \cdots & 0 \\ 0 & a_{22} & \cdots & 0 \\ \vdots & \vdots & \ddots & \vdots \\0 & 0 & \cdots & a_{nn} \end{bmatrix}, \qquad R = \begin{bmatrix} 0 & a_{12} & \cdots & a_{1n} \\ a_{21} & 0 & \cdots & a_{2n} \\ \vdots & \vdots & \ddots & \vdots \\a_{n1} & a_{n2} & \cdots & 0 \end{bmatrix}

线性方程组可以重写为:

:D \mathbf{x} = \mathbf{b} - R \mathbf{x}

雅可比法是一种迭代方法。在每一次迭代中，上一次算出的x被用在右侧，用来算出左侧的新的x。这个过程可以如下表示：

: \mathbf{x}^{(k+1)} = D^{-1} (\mathbf{b} - R \mathbf{x}^{(k)}).

对每个元素可以用以下公式：

: x^{(k+1)}_i = \frac{1}{a_{ii}} \left(b_i -\sum_{j\ne i}a_{ij}x^{(k)}_j\right),\quad i=1,2,\ldots,n.

注意计算 xi(k+1)需要x(k)中除了自己之外的每个元素。 不像高斯-赛德尔迭代，我们不能用 xi(k+1)覆盖xi(k)，因为在接下来的计算中还要用到这些值。这是雅可比和高斯-塞德尔方法最显著的差别，也是为什么前者可以用并行算法而后者不能的原因。最小需要的存储空间是两个长度为n的向量。

算法
选择一个初始猜想值x^{0}

: while 未收敛

:: for i := 1 step until n do

::: \sigma = 0

::: for j := 1 step until n do

:::: if j != i then
::::: \sigma = \sigma + a_{ij} x_j^{(k-1)}
:::: end if
::: end (j-loop)

::: x_i^{(k)} = {{\left( {b_i - \sigma } \right)} \over {a_{ii} }}
:: end (i-loop)
:: 检查是否收敛
: end (未收敛时继续循环)

收敛
标准的收敛情况是当迭代矩阵的谱半径

:\rho(D^{-1}R)

保证收敛的条件是矩阵A为严格或不可约地对角占优矩阵。严格的行对角占优矩阵指对于每一行，对角线上的元素之绝对值大于其余元素绝对值的和，即

:\left | a_{ii} \right | > \sum_{j \ne i} {\left | a_{ij} \right |}.

有时即使不满足此条件，雅可比法仍可收敛。

举例
一个形如 Ax=b 的线性方程，估计初始x^{(0)}：

: A=
\begin{bmatrix}
2 & 1 \\
5 & 7 \\
\end{bmatrix},
\ b=
\begin{bmatrix}
11 \\
13 \\
\end{bmatrix}
, \quad x^{(0)} =
\begin{bmatrix}
1 \\
1 \\
\end{bmatrix} .

我们用上文描述的方程 x^{(k+1)}=D^{-1}(b - Rx^{(k)})来估计x。首先，将等式写为D^{-1}(b - Rx^{(k)}) = Tx^{(k)} + C以方便计算，其中T=-D^{-1}R和C = D^{-1}b。注意 R=L+U中的 L和U是A的严格
递增和递减部分。变成如下数值：

: D^{-1}=
\begin{bmatrix}
1/2 & 0 \\
0 & 1/7 \\
\end{bmatrix},
\ L=
\begin{bmatrix}
0 & 0 \\
5 & 0 \\
\end{bmatrix}
, \quad U =
\begin{bmatrix}
0 & 1 \\
0 & 0 \\
\end{bmatrix} .

令 T=-D^{-1}(L+U) as
: T=
\begin{bmatrix}
1/2 & 0 \\
0 & 1/7 \\
\end{bmatrix}
\left\{
\begin{bmatrix}
0 & 0 \\
-5 & 0 \\
\end{bmatrix}
+
\begin{bmatrix}
0 & -1 \\
0 & 0 \\
\end{bmatrix}\right\}
=
\begin{bmatrix}
0 & -0.5 \\
-0.714 & 0 \\
\end{bmatrix} .

解出C为：

: C =
\begin{bmatrix}
1/2 & 0 \\
0 & 1/7 \\
\end{bmatrix}
\begin{bmatrix}
11 \\
13 \\
\end{bmatrix}
=
\begin{bmatrix}
5.5 \\
1.857 \\
\end{bmatrix} .

用计算出来的T和C，我们估计x为 x^{(1)}= Tx^{(0)}+C

: x^{(1)}=
\begin{bmatrix}
0 & -0.5 \\
-0.714 & 0 \\
\end{bmatrix}
\begin{bmatrix}
1 \\
1 \\
\end{bmatrix}
+
\begin{bmatrix}
5.5 \\
1.857 \\
\end{bmatrix}
=
\begin{bmatrix}
5.0 \\
1.143 \\
\end{bmatrix} .

继续迭代得：
: x^{(2)}=
\begin{bmatrix}
0 & -0.5 \\
-0.714 & 0 \\
\end{bmatrix}

\begin{bmatrix}
5.0 \\
1.143 \\
\end{bmatrix}
+
\begin{bmatrix}
5.5 \\
1.857 \\
\end{bmatrix}
=
\begin{bmatrix}
4.929 \\
-1.713 \\
\end{bmatrix} .
这个过程一直重复直到收敛（直到\|Ax^{(n)} - b\|足够小）。这个例子在25次迭代之后的解是
: x=\begin{bmatrix}
7.111\\
-3.222
\end{bmatrix}
.

一個用Fortran解的例子
subroutine solve(A,b,x,x0,n)
implicit real*8(a-z)
integer::n,imax=200
integer::i,j,k
real*8::tol=1d-15
real*8::A(n,n),b(n),x(n),x0(n),x1(n),x2(n)

write(102,501)
501 format('Jacobi iteration',/)

x1=x0

do k=1,imax
do i=1,n
s=0
do j=1,n
if (j==i) cycle
s=s+A(i,j)*x1(j)
enddo
x2(i)=(b(i)-s)/A(i,i)
enddo
! the following step check if convergence is reached
dx2=0
do i=1,n
dx2=dx2+(x1(i)-x2(i))**2
enddo
dx2=dsqrt(dx2)
if (dx2(1x,1pd25.15))
enddo
x=x2
end subroutine solve

program main
implicit real*8(a-z)
integer,parameter::n=3
real*8 ::A(n,n),b(n),x(n),x0(n)

open(unit=101,file='result.txt')
open(unit=102,file='it_result.txt')

x0=(/0d0,0d0,0d0/) ! initial guess
b=(/9d0,7d0,6d0/)
A=reshape((/10,-1,0,-1,10,-4,0,-2,10/),(/3,3/))

call solve(A,b,x,x0,n) ! solve Ax=b

write(101,501)'x(1)','x(2)','x(3)',x
501 format('jacobi result',//,(1x,a25),/,(1x,1pd25.15))

end program main

参考文献
外部链接
*

• [http://www.math-linux.com/spip.php?article49 Jacobi Method from www.math-linux.com]
• [https://web.archive.org/web/20070515064856/http://math.fullerton.edu/mathews/n2003/GaussSeidelMod.html Module for Jacobi and Gauss–Seidel Iteration]
• [http://pagerank.suchmaschinen-doktor.de/matrix-inversion.html Numerical matrix inversion]