查看“︁雅可比法”︁的源代码

{{NoteTA
|1=zh-hans:迭代;zh-hant:疊代;
|G1 = Math
}}

在数值线性代数中，'''雅可比法'''（{{lang|en|Jacobi Method}}）是一种解对角元素几乎都是各行和各列的绝对值最大的值的线性方程组的算法。求解出每个对角元素并插入近似值。不断迭代直至收敛<ref name="ppt">[http://www.bb.ustc.edu.cn/jpkc/xiaoji/szjsff/kj2/chap6.ppt] {{Wayback|url=http://www.bb.ustc.edu.cn/jpkc/xiaoji/szjsff/kj2/chap6.ppt |date=20220724043038 }},解线性方程组的迭代法</ref>。这个算法是[[雅可比矩阵]]的精简版。方法的名字来源于[[德国]]数学家[[卡尔·雅可比]]。

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

:<math>A\mathbf x = \mathbf b</math>

其中:

:<math>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}.</math>

''A'' 可以分解成[[对角矩阵|对角]]部分''D''和剩余部分''R'':

:<math>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} </math>

线性方程组可以重写为:

:<math>D \mathbf{x} = \mathbf{b} - R \mathbf{x} </math>

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

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


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

:<math> 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. </math>

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

==算法==
选择一个初始猜想值<math>x^{0}</math> <br>
:    while 未收敛 <br>
::   for i := 1 step until n do <br>
:::  <math> \sigma = 0 </math> <br>
:::   for j := 1 step until n do <br>
::::  if j != i then
:::::      <math> \sigma  = \sigma  + a_{ij} x_j^{(k-1)} </math>
::::  end if
:::    end (j-loop) <br>
:::    <math>  x_i^{(k)}  = {{\left( {b_i  - \sigma } \right)} \over {a_{ii} }} </math>
::   end (i-loop)
::   检查是否收敛
:    end (未收敛时继续循环)

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

:<math>\rho(D^{-1}R) < 1. </math><ref name="ppt" />

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

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

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

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

:<math> 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} .</math>

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

:<math> 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} .</math>

令 <math> T=-D^{-1}(L+U) </math> as
:<math> 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}  .</math>

解出C为：

:<math> 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} .</math>

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

:<math> 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}  .</math>

继续迭代得：
:<math> 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}  .</math>
这个过程一直重复直到收敛（直到<math>\|Ax^{(n)} - b\|</math>足够小）。这个例子在25次迭代之后的解是
:<math> x=\begin{bmatrix}
7.111\\
-3.222
\end{bmatrix}
.</math>

==一個用Fortran解的例子==
<syntaxhighlight lang="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<tol) exit
    x1=x2
    write(102,502)k,x1 ! record the iteration process
    502 format(1x,i3,<n>(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',//,<n>(1x,a25),/,<n>(1x,1pd25.15))

end program main
</syntaxhighlight>

==参考文献==
<references />

== 外部链接 ==
* {{MathWorld|urlname=JacobiMethod|title=Jacobi method|author=Black, Noel; Moore, Shirley; and Weisstein, Eric W.}}
* [http://www.math-linux.com/spip.php?article49 Jacobi Method from www.math-linux.com] {{Wayback|url=http://www.math-linux.com/spip.php?article49 |date=20130305230852 }}
* [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] {{Wayback|url=http://pagerank.suchmaschinen-doktor.de/matrix-inversion.html |date=20140901133130 }}

[[Category:数值线性代数]]
[[Category:带有伪代码示例的条目]]
[[Category:带有代码示例的条目]]