雅可比迭代法

雅可比迭代法（Jacobi iteration）把线性方程组中每个未知量由对应的一行解出，并在下一轮统一使用上一轮的数值更新。它实现简单，所有分量可以并行计算；收敛却需要检查矩阵结构或迭代矩阵，非零对角元本身不足以保证成功。

同时更新两个未知量

考虑 $4x+y=1$、$2x+3y=2$。由每行解出对应变量，得到 $$x^{k+1}=\frac{1-y^k}{4},y^{k+1}=\frac{2-2x^k}{3}.$$ 从 $(0,0)$ 出发，第一轮为 $(1/4,2/3)$；第二轮必须仍同时读取第一轮的旧值，得到 $(1/12,1/2)$；第三轮为 $(1/8,11/18)$。精确解为 $(1/10,3/5)$。

第二轮若在算出 $x^2=1/12$ 后，马上把它用于同一轮的 $y$，就不再是 Jacobi，而是高斯—赛德尔迭代法。是否保存旧向量是算法定义的一部分，并非一个无关的编程优化。

两条直线的交点为(0.1,0.6)。图中的折线连接迭代状态，不代表某一条微分方程轨迹。

矩阵分裂与收敛判据

设 $A=D+L+U$，其中 $D$ 是对角部分，$L,U$ 分别为严格下、上三角部分，且所有 $A_{ii}\ne 0$。Jacobi 更新写成 $$x^{k+1}=D^{-1}(b-(L+U)x^k)=B_J{x}^k+D^{-1}b,B_J=-D^{-1}(L+U).$$ 若精确解为 $x_{\ast }$，误差满足 $e^{k+1}=B_J{e}^k$。对所有初值收敛的充要条件是谱半径 $\rho (B_J)<1$：此时矩阵幂趋零；若有模至少一的特征值，便存在不衰减的相应模式。非对角化情形仍可用 Jordan 块证明，小于一的指数衰减压过块中有限次数的多项式因子。Netlib 给出 Jacobi 分裂及标准收敛条件。

主例 $B_J=\left(\begin{array}{cc}0& -1/4\\ -2/3& 0\end{array}\right)$，特征值为 $\pm 1/\sqrt{6}$，谱半径约为 $0.408248$，故收敛。各分量可换号或起伏，不要求每个坐标单调接近解。

严格对角占优的直接证明

若逐行满足 $|A_{ii}|>\sum _{j\ne i}|A_{ij}|$，令 $$q=\underset{i}{max}\frac{\sum _{j\ne i}|A_{ij}|}{|A_{ii}|}<1.$$ 任意误差向量满足 $$|e_i^{k+1}|\le \sum _{j\ne i}\frac{|A_{ij}|}{|A_{ii}|}|e_j^k|\le q\Vert e^k{\Vert }_{\mathrm{\infty }}.$$ 取最大值得 $\Vert e^k{\Vert }_{\mathrm{\infty }}\le q^k\Vert e^0{\Vert }_{\mathrm{\infty }}$。严格占优也保证 $A$ 非奇异：若 $Av=0$ 且取绝对值最大的非零分量 $v_i$，该行便要求 $|A_{ii}||v_i|\le \sum _{j\ne i}|A_{ij}||v_i|$，与严格不等式矛盾。于是唯一解和收敛都由同一个条件得到。

还可给出后验停止界。记相邻差 $d_k=x^k-x^{k-1}$，同样有 $\Vert d_{k+1}\Vert \le q\Vert d_k\Vert$。从当前点到极限的尾和满足 $$\Vert x_{\ast }-x^k{\Vert }_{\mathrm{\infty }}\le \sum _{j=1}^{\mathrm{\infty }}{q}^j\Vert d_k{\Vert }_{\mathrm{\infty }}=\frac{q}{1-q}\Vert d_k{\Vert }_{\mathrm{\infty }}.$$ 这比只凭相邻两次数字看起来相同更明确，但前提是有可信的收缩常数。主例可取 $q=2/3$，虽比实际谱半径保守，仍给出有效保证。

可复算反例与实现

有解变式取 $A=\mathrm{diag}(2,5)$、$b=(2,10{)}^T$，每行没有其他未知量，一轮即到 $(1,2{)}^T$。反例取 $A=\left(\begin{array}{cc}1& 2\\ 2& 1\end{array}\right)$、$b=(3,3{)}^T$。矩阵可逆、解为 $(1,1{)}^T$，但从零开始的迭代为 $(3,3),(-3,-3),(9,9),\dots$，因为迭代矩阵有模为二的特征值。

下面程序限定在严格行对角占优的有限实矩阵，以便给出收敛保证；其他矩阵即使可能收敛也不在此接口范围内。它从零开始，以实际残差最大范数验收，返回近似及轮数，不修改输入。稠密每轮 $O(n^2)$ 时间、额外向量 $O(n)$ 空间；稀疏矩阵可只遍历非零元。

Python 3

import math


def jacobi(A, b, tol=1e-12, maxiter=1000):
    A, b = [list(map(float, row)) for row in A], list(map(float, b)); n = len(b)
    if not n or len(A) != n or any(len(row) != n for row in A):
        raise ValueError("square matrix matching b required")
    if not all(math.isfinite(v) for row in A for v in row) or not all(math.isfinite(v) for v in b):
        raise ValueError("finite input required")
    if any(abs(A[i][i]) <= sum(abs(A[i][j]) for j in range(n) if j != i) for i in range(n)):
        raise ValueError("this implementation requires strict row diagonal dominance")
    if not math.isfinite(tol) or tol <= 0 or not isinstance(maxiter, int) or maxiter < 0:
        raise ValueError("positive tolerance and nonnegative integer iteration limit required")
    x = [0.0]*n; target = tol*max(1.0, max(map(abs, b)))
    for k in range(maxiter+1):
        residual = [b[i]-sum(a*v for a, v in zip(A[i], x)) for i in range(n)]
        if max(map(abs, residual)) <= target:
            return x, k
        if k == maxiter:
            break
        new = [(b[i]-sum(A[i][j]*x[j] for j in range(n) if j != i))/A[i][i]
               for i in range(n)]
        if not all(math.isfinite(v) for v in new):
            raise ArithmeticError("nonfinite Jacobi iterate")
        x = new
    raise RuntimeError("true residual did not reach tolerance")


print(jacobi([[4, 1], [2, 3]], [1, 2]))  # near ([0.1, 0.6], iteration count)

加权 Jacobi 另用 $x^{new}=x+\omega D^{-1}(b-Ax)$，在多重网格中常用作平滑器。改变 $\omega$ 会改变谱与稳定性条件，不能只因更新便于并行就不再检查收敛。