von Neumann 稳定性分析

von Neumann 稳定性分析把网格误差分解成不同频率的波，再问：每一种波经过一次更新以后，被放大了多少？ 对常系数、规则网格上的线性格式，平移不会改变波形，只会乘上一个复数，所以一个大型离散算子可以化成逐频率的代数计算。

这种方法尤其适合发现“曲线看起来还很平滑，但微小锯齿误差正在增长”的情形。决定稳定性的往往不是精确解中最明显的低频部分，而是网格允许的最高频率。

形式陈述

格式与放大因子

固定空间步长 $h>0$ 、时间步长 $\tau>0$ ，在无限等距网格或周期网格上考虑标量、常系数的一步格式

U_j^{n+1}=(S_{h,\tau}U^n)_j =\sum_{\ell=-r}^{s}a_\ell(h,\tau)U_{j+\ell}^n.

这里 $S_{h,\tau}$ 是一步传播算子。对模态 $e^{ij\theta}$ ，移位 $\ell$ 个网格点相当于乘 $e^{i\ell\theta}$ ，因此定义

G_{h,\tau}(\theta)=\sum_{\ell=-r}^{s}a_\ell(h,\tau)e^{i\ell\theta}.

周期网格有正整数 $N\ge1$ 个节点时，实际频率集合为 $\Theta_h=\{2\pi m/N:0\le m<N\}$ ；无限网格使用 $[-\pi,\pi]$ 。下文在周期情形使用 $\|U\|_{2,h}^2=h\sum_j|U_j|^2$ ，无限网格则对平方可和序列取相应和式。

从模态计算到范数稳定性

周期网格上的向量 $N^{-1/2}(e^{2\pi i jm/N})_{j=0}^{N-1}$ 构成正交规范基：不同频率的内积由有限等比和给出零，同一频率给出一。因此有限 DFT 保持 Euclidean 范数，乘上共同的 $\sqrt h$ 后也保持这里的网格范数。

无限网格则使用 $\ell^2(\mathbb Z)$ 到圆周 $L^2$ 的对应：令 $\mathcal FU(\theta)=\sum_{j\in\mathbb Z}U_je^{-ij\theta}$ ，级数按 $L^2$ 极限解释。Fourier 级数的 Parseval 等式给出

\frac{h}{2\pi}\int_{-\pi}^{\pi}|\mathcal FU(\theta)|^2\,d\theta =h\sum_{j\in\mathbb Z}|U_j|^2.

移位 $U_j\mapsto U_{j+\ell}$ 在这套对应下成为乘 $e^{i\ell\theta}$ 。于是，两种网格上的一步算子都化成乘相应的 $G_{h,\tau}$ ；这里用的是有限 DFT 或圆周级数对应，而非把序列直接代入 $\mathbb R$ 上的积分变换。对于这里的标量平移不变格式，传播算子的范数因此满足

\|S_{h,\tau}^{\,n}\|_{2,h} =\max_{\theta\in\Theta_h}|G_{h,\tau}(\theta)|^n

（无限网格用本质上确界）。这不是只检查了一组碰巧好算的输入：Fourier 模态完整地分解了整个空间，所以控制所有模态就控制了所有离散误差。[1]

一个最常见的充分条件是

|G_{h,\tau}(\theta)|\le1\qquad\text{对所有允许频率成立}.

它给出每步不增的 $\ell^2$ 范数。对于会真实增长的方程，这个要求可能过强。固定时间 $0\le n\tau\le T$ 上，更一般的稳定性条件是

\sup_{h,\tau}\sup_{n\tau\le T} \sup_{\theta\in\Theta_h}|G_{h,\tau}(\theta)|^n<\infty.

例如 $|G|\le1+C\tau$ 就足够，因为 $(1+C\tau)^n\le e^{CT}$ 。稳定的含义是放大率不随网格加密失控，而不是禁止方程在物理时间内有任何增长。

反过来，设允许的网格族中有一列时间步长 $\tau\to0$ ，并固定 $T>0$ 。若每个周期网格都有实际频率满足 $|G|\ge1+c$ ，其中 $c>0$ 与网格无关，那么取 $n=\lfloor T/\tau\rfloor$ ，相应单位模态的放大量至少为 $(1+c)^n\to\infty$ ，破坏有限时间的一致稳定性。

无限网格上的纯指数序列不属于 $\ell^2(\mathbb Z)$ ，但同样的下界可以用 Fourier 波包实现。有限模板的符号 $G$ 连续；若某点满足 $|G|\ge1+c$ ，它的一个正测度邻域上便有 $|G|\ge1+c/2$ 。取 Fourier 变换支撑在该邻域、网格范数归一化的 $\ell^2$ 输入，则经过 $n$ 步的范数至少为 $(1+c/2)^n$ 。邻域和输入可以随网格改变，仍足以证明传播算子范数无界。

热方程 FTCS 的放大因子 G(θ)=1−4r sin²(θ/2)。当 r=0.6 时，θ=π 处 G=−1.4，模态既翻号又增大。

直觉

一种波为什么可以代表一个计算方向

为看清前面传播算子的作用，取一个离散 Fourier 模态

U_j^n=\widehat U^n(\theta)e^{ij\theta}.

这里 $\theta$ 是每跨过一个网格间距增加的相位。 $\theta\approx0$ 对应缓慢起伏， $\theta=\pi$ 对应正负交替的锯齿序列。物理波数 $k$ 与离散频率的关系是 $\theta=kh$ ；离散频率可限制在一个长度为 $2\pi$ 的基本区间。

由于 $U_{j+\ell}^n=e^{i\ell\theta}U_j^n$ ，代入格式便得到

\widehat U^{n+1}(\theta)=G(\theta)\widehat U^n(\theta), \qquad G(\theta)=\sum_{\ell=-r}^{s}a_\ell e^{i\ell\theta}.

$G$ 称为放大因子或离散 Fourier 符号。经过 $n$ 步，同一模态被乘上 $G(\theta)^n$ 。复数的模 $|G|$ 决定振幅变化，辐角 $\arg G$ 决定位移和相位变化，二者有不同的物理含义。

分析整个网格族时，经常取对全部 $\theta$ 的上确界，以获得不依赖节点数的结论；固定周期网格的精确判定则只需检查 $\Theta_h$ 中实际存在的模态。

补充示意图

热方程 FTCS 的 Fourier 模态、放大因子曲线与最高频翻号 — von Neumann 放大因子

例子与边界

热方程：最短的锯齿决定步长

对热扩散方程 $u_t=\kappa u_{xx}$ ，设 $\kappa>0$ ，使用显式 Euler 与中心差分：

U_j^{n+1}=U_j^n+r(U_{j-1}^n-2U_j^n+U_{j+1}^n), \qquad r=\frac{\kappa\tau}{h^2}.

代入 $e^{ij\theta}$ ，利用 $e^{i\theta}+e^{-i\theta}=2\cos\theta$ ，得到

G(\theta)=1-2r+2r\cos\theta =1-4r\sin^2(\theta/2).

当 $r\ge0$ 时， $G$ 在 $[1-4r,1]$ 内变化。因此 $0\le r\le1/2$ 保证所有频率都满足 $|G|\le1$ 。这给出熟悉的限制

\tau\le\frac{h^2}{2\kappa}.

在包含 $\theta=\pi$ 的偶数节点周期网格上，最高频率正好给出这个阈值；对整个任意加密的周期网格族，它也是自然的一致阈值。个别固定的奇数节点网格因没有恰好的 $\pi$ 模态，精确极限可以略有不同。

取 $r=0.6$ ，最高频模态的 $G(\pi)=-1.4$ 。负号使它每步翻转，模长使它每步增长 40%；十步后，原始误差的振幅变成 $1.4^{10}\approx28.9$ 倍。即使初始数据没有明显锯齿，舍入或边界扰动也可能激发这种方向。

在临界值 $r=1/2$ 时， $G(\pi)=-1$ 。这个模态不增长，却也不衰减，只是不断翻转。因而“稳定”与“充分耗散”是两种不同要求：扩散方程本来应该迅速抹平高频，临界格式却把最高频锯齿保留了下来。

平流方程：同时看振幅和相位

对 $u_t+a u_x=0$ ，先取 $a>0$ 。一阶迎风格式为

U_j^{n+1}=U_j^n-\nu(U_j^n-U_{j-1}^n), \qquad \nu=\frac{a\tau}{h}.

其放大因子为

G(\theta)=1-\nu+\nu e^{-i\theta},

进而

|G(\theta)|^2 =1-4\nu(1-\nu)\sin^2(\theta/2).

因此 $0\le\nu\le1$ 保证不放大。若 $a<0$ ，迎风方向也应反转，不能只把同一左侧模板中的 $a$ 改成负数。

精确平流在一步内只把波平移，放大因子为 $e^{-i\nu\theta}$ ，模始终等于 1。迎风格式在 $0<\nu<1$ 时对非零频率通常有 $|G|<1$ ，于是引入了数值耗散：波形会逐渐变钝。对小 $\theta$ ，

\log G(\theta) =-i\nu\theta-\frac{\nu(1-\nu)}2\theta^2+O(\theta^3).

第一项对应所需平移，第二项对应额外扩散。换成物理波数 $k=\theta/h$ ，这项扩散的主导系数为 $a h(1-\nu)/2$ 。这个展开只描述低频行为；接近网格分辨率极限时，应回到精确的 $G(\theta)$ 。

当 $\nu=1$ 时，更新变成 $U_j^{n+1}=U_{j-1}^n$ ，恰好每步平移一个网格点， $G=e^{-i\theta}$ 。这个特例揭示了数值耗散为何依赖 Courant 数，而不是只由“迎风”这个方法名称决定。[2]

中心空间差分为什么会出问题

若空间导数换成中心差分，仍使用显式 Euler，则

G(\theta)=1-i\nu\sin\theta, \qquad |G(\theta)|^2=1+\nu^2\sin^2\theta.

在固定非零 $\nu$ 、 $\tau\to0$ 的无限网格族上， $\theta=\pi/2$ 附近的 Fourier 波包持续被放大，从而失稳；包含这一实际频率的周期网格族也有同样结论。这个例子说明，“中心差分更高阶”与“完全离散格式更可靠”不是同一个判断。

这里的限定语也有作用：若采用更强的 $\tau=O(h^2)$ 约束，有限时间幂界可以成立，因为 $|G|^{T/\tau}\le\exp(a^2T\tau/(2h^2))$ 。所以正确结论要连同网格缩放一起陈述，而不能仅凭某一步 $|G|>1$ 就忽略它偏离 1 的速度。

向量系统：特征值还不够

若每个网格点存一个向量，Fourier 变换后得到的是矩阵符号

\widehat U^{n+1}(\theta)=G(\theta)\widehat U^n(\theta).

此时需要控制 $\|G(\theta)^n\|$ 。只有当符号矩阵在所用内积下正规，或存在条件数一致有界的对角化时，才能方便地把问题化成单独控制各特征值。

例如

G=\begin{pmatrix}1&1\\0&1\end{pmatrix}, \qquad G^n=\begin{pmatrix}1&n\\0&1\end{pmatrix}.

两个特征值都为 1，幂却随 $n$ 增长。若一步长度 $\tau\to0$ ，固定物理时间内有 $n\approx T/\tau$ ，这种增长破坏了一致稳定性。

但把右上角改成 $\tau$ 后，

G_\tau^n=\begin{pmatrix}1&n\tau\\0&1\end{pmatrix}

在 $n\tau\le T$ 时有界。于是，连“出现 Jordan 块”也不能脱离时间尺度独立判定失稳。真正需要检查的是矩阵幂的统一界。

对于多步时间格式，还需建立放大根方程或伴随矩阵，检查单位圆附近根的结构，并控制启动值。一步格式中的一个标量 $G$ 不能直接覆盖全部多步情形。[3]

推论与应用

与 ODE 绝对稳定域的接口

若线法先把空间离散为 $\dot U=L_hU$ ，某个 Fourier 模态满足

L_he^{ij\theta}=\lambda_h(\theta)e^{ij\theta}.

再使用稳定函数为 $R(z)$ 的时间积分方法，这个模态的完全离散放大因子就是

G_{h,\tau}(\theta)=R\bigl(\tau\lambda_h(\theta)\bigr).

因此，绝对稳定域把时间积分器与空间谱连接起来。例如显式 Euler 的 $R(z)=1+z$ ，其在负实轴上的稳定区间为 $[-2,0]$ ；热方程空间谱最负端约为 $-4\kappa/h^2$ ，于是再次得到 $\tau\lesssim h^2/(2\kappa)$ 。

后向 Euler 的 $R(z)=1/(1-z)$ 会把很负的实特征值强烈压小。梯形法的 $R(z)=(1+z/2)/(1-z/2)$ 虽然在整个左半平面满足 $|R|\le1$ ，但 $z\to-\infty$ 时趋于 $-1$ 。它能够稳定传播刚性模态，却未必把这些模态迅速消除。详细的隐式求解与阻尼比较见隐式 ODE 方法。

这项分析能看见什么，看不见什么

常系数与平移不变性让不同 Fourier 模态互不耦合。变系数会混合频率；非周期边界也可能产生只在边界附近出现的不稳定方向。冻结局部系数所做的 Fourier 检查仍有诊断价值，但一般不能替代完整的边界稳定性或能量估计。

非线性问题中，围绕某个解线性化可检查小扰动，却不能仅凭线性化稳定就得到任意振幅下的非线性稳定结论。范数同样重要：这里通过 Fourier 变换直接得到的是 $\ell^2$ 控制，不是自动得到最大范数稳定、正性保持或离散最大值原理。

使用时可以沿着一条清晰路线：固定方程、边界、范数和允许的网格族，计算实际符号，再检查全部频率的有限时间幂界。得到这个界以后，才与一致性误差结合，判断数值解是否收敛以及能够达到什么阶。

这些分析为有限差分法提供模态稳定性判据，使用时仍须核对边界和系数假设。

参考资料

[1] Lloyd N. Trefethen, Finite Difference and Spectral Methods for Ordinary and Partial Differential Equations, 1996，§§2.2、2.4 的半离散变换与 DFT，Chapter 4 的稳定性；作者开放教材。

[2] Randall J. LeVeque, Finite Difference Methods for Ordinary and Partial Differential Equations, SIAM, 2007，平流、扩散、数值耗散与色散。

[3] John C. Strikwerda, Finite Difference Schemes and Partial Differential Equations, 2nd ed., SIAM, 2004，标量格式、系统、边界与多步稳定性。

形式陈述 ​

格式与放大因子 ​

从模态计算到范数稳定性 ​

一种波为什么可以代表一个计算方向 ​

热方程：最短的锯齿决定步长 ​

平流方程：同时看振幅和相位 ​

中心空间差分为什么会出问题 ​

向量系统：特征值还不够 ​

与 ODE 绝对稳定域的接口 ​

这项分析能看见什么，看不见什么 ​

反向链接