Picard 迭代与 Newton 迭代：以 Navier–Stokes 方程为例

约 5111 字大约 17 分钟

2026-09-15

求解非线性方程时，常见的思路是：从一个猜测出发，每次解一个较容易的子问题，再逐步修正这个猜测。Picard 迭代和 Newton 迭代都属于这类方法，但它们构造子问题的方式不同：

Picard 迭代：构造不动点映射；在非线性 PDE 中，常表现为冻结一部分系数。
Newton 迭代：对整个残差做一阶线性化，求一个消除当前残差的修正量。

下面先从标量方程说明两者的区别，再推导它们在定常不可压缩 Navier–Stokes 方程中的形式。

1. 从一个标量方程开始

考虑

F(x)=x+x^2-2=0.

它有两个根 $1$ 和 $-2$ 。下面从 $x^0=0$ 出发，观察趋近正根 $x^\ast=1$ 的过程。上标 $k$ 表示迭代次数。

1.1 Picard：把旧值当成系数

将方程写成

(1+x)x=2.

在第 $k$ 步，冻结括号里的 $x$ ，只把后面的 $x$ 作为未知量：

(1+x^k)x^{k+1}=2, \qquad x^{k+1}=\frac{2}{1+x^k}.

于是得到不动点迭代 $x^{k+1}=G(x^k)$ ，其中

G(x)=\frac{2}{1+x}.

在正根处，

G'(1)=-\frac12.

为什么导数决定了误差的变化？定义误差 $e^k=x^k-1$ ，即 $x^k=1+e^k$ 。由于 $1$ 是不动点，满足 $G(1)=1$ ，所以

e^{k+1}=x^{k+1}-1=G(x^k)-G(1)=G(1+e^k)-G(1).

在 $1$ 处对 $G$ 做 Taylor 展开：

G(1+e^k)=G(1)+G'(1)e^k+O\bigl((e^k)^2\bigr).

减去 $G(1)$ ，便得到

e^{k+1}=G'(1)e^k+O\bigl((e^k)^2\bigr) =-\frac12 e^k+O\bigl((e^k)^2\bigr).

这里的 $O\bigl((e^k)^2\bigr)$ 表示余项在误差足够小时，其绝对值不超过某个常数乘以 $(e^k)^2$ 。因此，只有在根附近、忽略二阶及更高阶项时，才有近似关系

e^{k+1}\approx-\frac12 e^k,

而不是对每一步都严格成立的等式。直观上，导数 $G'(1)$ 就是映射 $G$ 在不动点附近对小误差的放大或缩小系数。

本例还可以直接算出精确的误差关系：

e^{k+1}=\frac{2}{1+x^k}-1 =\frac{2}{2+e^k}-1 =-\frac{e^k}{2+e^k}.

当 $e^k$ 趋于零时，分母趋于 $2$ ，所以非零误差的比值 $e^{k+1}/e^k$ 趋于 $-1/2$ 。误差一边改变符号，一边大致缩小为原来的一半，这就是典型的线性收敛。

同一个方程，换一种不动点改写，收敛性质可能完全不同。 例如写成

x^{k+1}=2-(x^k)^2,

则映射在 $x^\ast=1$ 处的导数为 $-2$ ，这个不动点是局部排斥的。从 $x^0=0$ 出发，它实际上会落到另一个根 $-2$ 。所以“写成了等价方程”并不意味着“得到了合适的迭代”。

1.2 Newton：对残差求导

在当前点附近，用切线近似残差：

F(x^k+\delta x) \approx F(x^k)+F'(x^k)\delta x.

令这个线性近似等于零，就得到

(1+2x^k)\delta x^k =-\bigl(x^k+(x^k)^2-2\bigr),

再更新

x^{k+1}=x^k+\delta x^k.

也就是熟悉的 Newton 公式：

x^{k+1} =x^k-\frac{x^k+(x^k)^2-2}{1+2x^k}.

两种方法从同一初值出发的结果如下。表中的误差相对于正根 $1$ 计算，数值经过舍入。

$k$	Picard 的 $x^k$	Picard 误差	Newton 的 $x^k$	Newton 误差
0	0.000000000	$1.000$	0.000000000	$1.000$
1	2.000000000	$1.000$	2.000000000	$1.000$
2	0.666666667	$3.333\times10^{-1}$	1.200000000	$2.000\times10^{-1}$
3	1.200000000	$2.000\times10^{-1}$	1.011764706	$1.176\times10^{-2}$
4	0.909090909	$9.091\times10^{-2}$	1.000045777	$4.578\times10^{-5}$
5	1.047619048	$4.762\times10^{-2}$	1.000000001	$6.985\times10^{-10}$

对这个例子，可以直接算出 Newton 的误差关系：

e^{k+1} =\frac{(e^k)^2}{3+2e^k}.

在根附近，分母接近 $3$ ，误差主要按平方缩小。这就是二次收敛。它描述的是接近解以后的行为，并不保证前几步就快速下降。

2. 一般非线性方程中的两种方法

考虑

F(U)=0,

其中 $U$ 可以是有限维向量，也可以是函数空间中的未知量。

2.1 Picard 是不动点迭代

先将问题改写成

U=G(U),

再依次计算

U^{k+1}=G(U^k).

一个常用的充分条件是： $G$ 将某个非空完备闭集映到自身，并且在这个集合内满足压缩性质

\|G(U)-G(V)\|\le q\|U-V\|, \qquad 0\le q<1.

此时，不动点存在且唯一，从该集合中的初值出发有

\|U^{k+1}-U^\ast\| \le q\|U^k-U^\ast\|.

对于具有

A(U)U=b

结构的方程，一种常见的 Picard 迭代是

A(U^k)U^{k+1}=b.

这里的“不动点映射”包含了一次线性求解。Picard 并不意味着只需把旧值代入一个简单的显式公式。

2.2 Newton 是完整的一阶线性化

Newton 方法求解

F'(U^k)\delta U^k=-F(U^k),

然后更新

U^{k+1}=U^k+\delta U^k.

在有限维情形， $F'(U^k)$ 就是 Jacobian 矩阵；在函数空间中，它是相应的导数算子。关于这层联系，可以参看 Frechet 导数和 Gateaux 导数。

对于 $F(U)=A(U)U-b$ ，乘积求导给出

F'(U)[\delta U] =A(U)\delta U+A'(U)[\delta U]U.

Picard 冻结 $A(U)$ ，Newton 则还考虑 $U$ 的变化会使算子 $A(U)$ 本身发生怎样的变化。后面的 Navier–Stokes 推导就是这个公式的具体例子。

当解处的导数可逆、导数在附近满足适当的 Lipschitz 连续性、初值足够接近解，而且线性子问题精确求解时，Newton 具有局部二次收敛：

\|U^{k+1}-U^\ast\| \le C\|U^k-U^\ast\|^2.

因此，“Newton 二次收敛”有明确的局部条件，不能理解为任意初值下都比 Picard 好。

3. Navier–Stokes 的非线性来自哪里

考虑常密度、常黏度的定常不可压缩 Navier–Stokes 方程：

\begin{cases} -\nu\Delta u+(u\cdot\nabla)u+\nabla p=f, & \text{在 }\Omega\text{ 内},\\ \nabla\cdot u=0, & \text{在 }\Omega\text{ 内},\\ u=g, & \text{在 }\partial\Omega\text{ 上}. \end{cases}

这里：

$u$ 是速度；
$p$ 是除以常密度后的压力；
$\nu>0$ 是运动黏度；
$f$ 是单位质量上的体力。

为集中讨论线性化，本文在整个边界给定速度，并假定数据足够规则、满足通量相容条件

\int_{\partial\Omega}g\cdot n\,ds=0.

此时压力只确定到一个加法常数，可以加上

\int_\Omega p\,dx=0

固定其唯一表示。其他边界条件需要相应调整弱形式和压力处理。

常黏度下，黏性项、压力梯度和不可压缩约束都是线性的。非线性来自对流项

N(u)=(u\cdot\nabla)u.

按分量写，

N_i(u)=\sum_{j=1}^d u_j\frac{\partial u_i}{\partial x_j}.

这里的 $u$ 有两个角色：前面的 $u_j$ 决定沿哪个方向、以多快的速度输运；后面的 $u_i$ 是被输运的速度分量。两者都未知，所以形成二次非线性。

4. Picard：冻结输运速度，得到 Oseen 问题

已知 $u^k$ 后，将对流项近似为

(u\cdot\nabla)u \quad\longrightarrow\quad (u^k\cdot\nabla)u^{k+1}.

于是每一步求解

\boxed{ \begin{aligned} -\nu\Delta u^{k+1} +(u^k\cdot\nabla)u^{k+1} +\nabla p^{k+1}&=f,\\ \nabla\cdot u^{k+1}&=0. \end{aligned} }

同时施加 $u^{k+1}=g$ 和压力零均值条件。

因为 $u^k$ 已知，这对 $(u^{k+1},p^{k+1})$ 是线性的，通常称为 Oseen 问题。它保留了“已知速度场输运未知速度”的结构。

注意，这不是把整个对流项都替换成 $(u^k\cdot\nabla)u^k$ 。后一种做法相当于把对流项整体移到右端，形成另一种不动点迭代，收敛性质也会不同。本文所说的 Navier–Stokes Picard 迭代，特指上述 Oseen 形式。

如果迭代收敛到 $(u^\ast,p^\ast)$ ，在适当的连续性条件下取极限，就重新得到原来的非线性方程。

5. Newton：对两个速度因子同时求导

5.1 展开对流项

令当前速度受到扰动 $u^k+\delta u$ ，直接展开：

\begin{aligned} N(u^k+\delta u) ={}&(u^k\cdot\nabla)u^k\\ &+(u^k\cdot\nabla)\delta u\\ &+(\delta u\cdot\nabla)u^k\\ &+(\delta u\cdot\nabla)\delta u. \end{aligned}

最后一项关于 $\delta u$ 是二次的。一阶线性化保留前面的常数项和两个线性项，因此

\boxed{ N'(u^k)[\delta u] =(u^k\cdot\nabla)\delta u +(\delta u\cdot\nabla)u^k. }

两项分别表示：

$(u^k\cdot\nabla)\delta u$ ：原来的流场输运速度扰动；
$(\delta u\cdot\nabla)u^k$ ：输运速度的扰动作用于原来的速度梯度。

第二项正是 Picard 冻结输运速度时没有保留的反馈。

5.2 增量形式：先解修正量，再更新

定义动量残差和连续性残差

\begin{aligned} R_m^k&=-\nu\Delta u^k+(u^k\cdot\nabla)u^k+\nabla p^k-f,\\ R_c^k&=\nabla\cdot u^k. \end{aligned}

Newton 修正量满足

\boxed{ \begin{aligned} -\nu\Delta\delta u^k +(u^k\cdot\nabla)\delta u^k +(\delta u^k\cdot\nabla)u^k +\nabla\delta p^k&=-R_m^k,\\ \nabla\cdot\delta u^k&=-R_c^k. \end{aligned} }

求解后更新

u^{k+1}=u^k+\delta u^k, \qquad p^{k+1}=p^k+\delta p^k.

若当前速度已满足边界条件 $u^k=g$ ，则修正量必须满足

\delta u^k=0 \quad\text{在 }\partial\Omega\text{ 上}.

若当前压力为零均值，则也令 $\int_\Omega\delta p^k\,dx=0$ 。

这里不能随意把第二行右端写成零。只有当前速度已经满足相应的不可压缩约束时，才有 $\nabla\cdot\delta u^k=0$ ；一般情况下，Newton 也需要修正当前的散度误差。

5.3 新解形式：右端为什么多出一个对流项

将 $\delta u^k=u^{k+1}-u^k$ 和 $\delta p^k=p^{k+1}-p^k$ 代入，可以得到等价的全步长形式：

\boxed{ \begin{aligned} -\nu\Delta u^{k+1} +(u^k\cdot\nabla)u^{k+1} +(u^{k+1}\cdot\nabla)u^k +\nabla p^{k+1} &=f+(u^k\cdot\nabla)u^k,\\ \nabla\cdot u^{k+1}&=0. \end{aligned} }

右端多出的 $(u^k\cdot\nabla)u^k$ 来自展开并整理旧值项。也可以直接从

N(u^{k+1}) \approx (u^k\cdot\nabla)u^{k+1} +(u^{k+1}\cdot\nabla)u^k -(u^k\cdot\nabla)u^k

看出这个符号。

增量形式的右端是负残差，新解形式的右端是整理后的载荷；两者不能混用。 上式描述的是完整 Newton 步。如果使用阻尼，应再对求得的候选解与当前解做加权更新。

6. 弱形式和离散矩阵是什么样的

有限元通常从弱形式出发。令

V_0=H_0^1(\Omega)^d, \qquad V_g=\{v\in H^1(\Omega)^d:v|_{\partial\Omega}=g\}, \qquad Q=L_0^2(\Omega),

其中 $L_0^2(\Omega)$ 表示零均值平方可积函数空间。

用 $(\cdot,\cdot)$ 表示 $L^2$ 内积，并定义

\begin{aligned} a(u,v)&=\nu(\nabla u,\nabla v),\\ c(w;u,v)&=((w\cdot\nabla)u,v),\\ b(v,q)&=-(q,\nabla\cdot v). \end{aligned}

原问题是求 $(u,p)\in V_g\times Q$ ，使得对任意 $(v,q)\in V_0\times Q$ ，

\begin{aligned} a(u,v)+c(u;u,v)+b(v,p)&=(f,v),\\ b(u,q)&=0. \end{aligned}

Picard 每一步满足

\begin{aligned} a(u^{k+1},v)+c(u^k;u^{k+1},v)+b(v,p^{k+1})&=(f,v),\\ b(u^{k+1},q)&=0. \end{aligned}

定义当前弱残差

\mathcal R^k(v,q) =a(u^k,v)+c(u^k;u^k,v)+b(v,p^k)+b(u^k,q)-(f,v).

Newton 则求 $(\delta u^k,\delta p^k)\in V_0\times Q$ ，使

\begin{aligned} &a(\delta u^k,v) +c(u^k;\delta u^k,v) +c(\delta u^k;u^k,v)\\ &\qquad+b(v,\delta p^k)+b(\delta u^k,q) =-\mathcal R^k(v,q). \end{aligned}

选取有限元基函数，并处理速度边界自由度和压力常数模态后，Picard 线性子问题与 Newton 修正方程的矩阵分别具有结构

K_{\mathrm P}^k= \begin{bmatrix} A+N(U^k)&B^T\\ B&0 \end{bmatrix}, \qquad K_{\mathrm N}^k= \begin{bmatrix} A+N(U^k)+W(U^k)&B^T\\ B&0 \end{bmatrix}.

这里 $U^k$ 表示离散速度系数向量，且：

$A$ 来自黏性项 $a$ ；
$N(U^k)$ 来自 $c(u^k;\delta u,v)$ ；
$W(U^k)$ 来自 $c(\delta u;u^k,v)$ ；
$B$ 来自所定义的散度耦合 $b$ 。

Newton 多出的正是 $W(U^k)$ ，未知量的数量没有因此增加。

对于这里采用的标准、未加稳定化的混合有限元，速度和压力空间应满足离散 inf-sup 条件。线性化方法本身不能修复不稳定的空间配对；使用稳定化方法时，则需要对完整的离散残差一致地线性化。

实际实现中，可以用 Oseen 形式构造 Picard 矩阵，也可以对弱残差自动求导得到 Newton Jacobian；FENaPack 的 Navier–Stokes 示例展示了这两种构造方式。

7. 为什么一个通常线性收敛，一个局部二次收敛

前面的对流项展开已经揭示了关键：Newton 保留了所有一阶项，遗漏项是

(\delta u\cdot\nabla)\delta u,

它对扰动是二次的。Picard 则省略了一项一阶反馈。

还可以直接比较误差方程。令 $u^\ast$ 为真解，

e^k=u^k-u^\ast, \qquad \pi^{k+1}=p^{k+1}-p^\ast.

因为各步和真解具有相同的速度边界值，误差在边界上为零。以下比较针对精确线性求解和全步长更新。

7.1 Picard 的误差右端是一阶的

将真解方程与 Picard 方程相减并整理，得到

-\nu\Delta e^{k+1} +(u^k\cdot\nabla)e^{k+1} +\nabla\pi^{k+1} =-(e^k\cdot\nabla)u^\ast.

右端关于旧误差 $e^k$ 是一阶的。若左端 Oseen 问题有适当的稳定性估计，就可以得到形如

\|e^{k+1}\|_V\le q\|e^k\|_V

的界；还需要 $q<1$ 才能据此推出收敛。小数据或黏性足够强，是获得这类压缩估计的典型条件。

7.2 Newton 的误差右端是二次的

对 Newton 的新解形式做同样的操作，可得

\begin{aligned} -\nu\Delta e^{k+1} +(u^k\cdot\nabla)e^{k+1} +(e^{k+1}\cdot\nabla)u^k +\nabla\pi^{k+1} =(e^k\cdot\nabla)e^k. \end{aligned}

右端已经是旧误差的二次项。在解附近，如果线性化算子的逆一致有界，并且对流项满足相应的连续性估计，就能得到

\|e^{k+1}\|_V\le C\|e^k\|_V^2.

这说明 Newton 的快速收敛来自完整的一阶抵消，也说明了为什么 Jacobian 接近奇异时，这个结论可能失效。

比较项	Picard / Oseen	Newton
线性化方式	冻结输运速度	对整个残差求导
对流线性项	保留一个	保留两个
典型收敛性质	压缩条件下通常线性收敛	适当条件下局部二次收敛
初值要求	压缩区域内可收敛，不保证任意初值	局部理论要求初值足够接近解
每步工作	更新 Oseen 矩阵并求解	更新完整 Jacobian 并求解
总体成本	可能需要较多外层迭代	步数可能更少，每步求解难度也可能更高