5.13 选读：常微分方程与动力系统视角

本章的最后节把前面各节用到的工具集中到一个具体场景：描述「随时间演化」的方程。这类方程在应用中出现得最多：物理中的运动方程、化学中的反应动力学、生物中的种群模型、机器学习中的参数演化（梯度下降本身就是一条常微分方程的离散化）。本节的目标不是把常微分方程理论讲完整，而是把本章前面各节的工具如何在其中起作用这件事说清楚。

四条线索构成本节。存在唯一性用 Picard 迭代证明，它的核心是 5.1 节的压缩映射原理：把初值问题写成积分方程，再证明右端是压缩映射。相平面分析用 5.8 节的线性化与 Jacobian 特征值，把平衡点分类为结点、鞍点、焦点与中心，稳定性由特征值的实部决定。Lyapunov 方法把「稳定性」的证明化为「找一个沿轨道单调的函数」，梯度流是它最直接的例子（$V$ 沿梯度流单调下降，下降速率由 Hessian 的特征值给出）。数值解法处理「解不能解析写出」的情形，其中刚性是主要困难：显式方法的稳定域有限，隐式方法（如隐式 Euler）对任意步长稳定。实验二的数据把这一点量化：刚性比 $1000$ 的问题上，显式 Euler 在 $h=0.01$ 时误差 $3.5\times {10}^{185}$（完全发散），隐式 Euler 在同一个步长下误差 $5.0\times {10}^{-6}$。

5.13.1 初值问题与存在唯一性

问题的提法与积分方程形式

一阶常微分方程组的初值问题指

$$y^{\prime }=f(t,y),y(t_0)=y_0$$

其中 $y\in {\mathbb{R}}^n$、$f:\mathbb{R}\times {\mathbb{R}}^n\to {\mathbb{R}}^n$。高阶方程可以通过引入中间变量化成一阶方程组（$y^{″}=f(t,y,y^{\prime })$ 令 $v=y^{\prime }$ 得 $y^{\prime }=v$、$v^{\prime }=f(t,y,v)$），因此讨论一阶组不失一般性。

把方程写成积分方程（对 $t$ 从 $t_0$ 积分到 $t$）：

$$y(t)=y_0+{\int }_{t_0}^{t}f(s,y(s))ds$$

这两个形式是等价的（$y$ 连续且满足其一则满足另一）。积分形式的好处是右端只含 $y$ 而不再含 $y^{\prime }$，从而可以把「求导」的运算换成「求积分」的运算，而积分对函数的要求更弱（5.5 节与 5.10 节的讨论）。

Picard 迭代

Picard 迭代（逐步逼近法）把积分方程当作不动点方程来解：

$$y_{k+1}(t)=y_0+{\int }_{t_0}^{t}f(s,y_k(s))ds$$

从常函数 $y_0(t)\equiv y_0$ 出发，逐次积分。Picard-Lindelöf 定理：若 $f$ 在包含 $(t_0,y_0)$ 的矩形区域内连续，且关于 $y$ 满足 Lipschitz 条件

$$\Vert f(t,y_1)-f(t,y_2)\Vert \le L\Vert y_1-y_2\Vert$$

则存在 $\delta >0$ 使初值问题在 $|t-t_0|\le \delta$ 上有唯一解，且该解是 Picard 迭代的一致收敛极限。

证明的骨架与 5.1 节的压缩映射原理完全一致：在连续函数空间 $C[t_0-\delta ,t_0+\delta ]$（配以上确界范数）上定义映射 $T$（右端的积分算子），则

$${\Vert T{y}_1-T{y}_2\Vert }_{\mathrm{\infty }}\le L\delta {\Vert y_1-y_2\Vert }_{\mathrm{\infty }}$$

取 $\delta <\frac{1}{L}$ 使 $T$ 成为压缩映射，由压缩映射原理得不动点存在唯一，且迭代收敛。这里出现的三条要素全部来自前面的章节：上确界范数与完备性（5.6 节与 5.1 节）、Lipschitz 估计（5.4 节的中值定理）、压缩映射原理（5.1 节）。

实验一第一组的数据给出收敛的量化形态（$y^{\prime }=y$、$y(0)=1$）：最大误差从 $n=1$ 的 $7.183\times {10}^{-1}$ 降到 $n=8$ 的 $3.002\times {10}^{-6}$，相邻比值从 $3.29$ 增到 $9.26$，与 $n+1$ 一致。这个形态有解析解释：$y^{\prime }=y$ 的 Picard 迭代给出 $y_n(t)=\sum _{k=0}^n\frac{t^k}{k!}$（指数级数的部分和），误差在 $t=1$ 处是 $\sum _{k>n}\frac{1}{k!}\approx \frac{1}{(n+1)!}$，相邻比值 $\approx n+1$。迭代的收敛速度由级数尾部的量级决定，这一点与其他迭代方法（5.4 节的求根、5.1 节的压缩映射）是同一个道理。

Lipschitz 条件的必要性与解的爆破

Lipschitz 条件是定理的关键假设，去掉它之后存在性仍可能成立，但唯一性会失效。标准反例是 $y^{\prime }=\sqrt{|y|}$、$y(0)=0$：$y\equiv 0$ 是解，$y=\frac{t^2}{4}$（$t\ge 0$）也是解，两条轨道从同一个初值出发分叉。原因是 $f=\sqrt{|y|}$ 在 $y=0$ 附近不满足 Lipschitz 条件（$\frac{\sqrt{ϵ}}{ϵ}={ϵ}^{-1/2}\to \mathrm{\infty }$）。

另一类现象是解的爆破：$f$ 的增长太快时解在有限时间内逃到无穷。实验一第二组用 $y^{\prime }=y^2$、$y(0)=1$ 演示，精确解 $y=\frac{1}{1-t}$ 在 $t=1$ 处发散。Lipschitz 常数在区间 $[0,a]$ 上是 $L=\frac{2}{1-a}$（$f^{\prime }=2y$ 的最大值），压缩条件 $aL<1$ 给出 $a<\frac{1}{3}$。数据给出四条区间上的迭代误差：$a=0.2$ 时 $1.304\times {10}^{-9}$（$aL=0.5$，压缩条件满足）、$a=0.3$ 时 $5.566\times {10}^{-8}$（$aL=0.857$，接近临界）、$a=0.5$ 时 $2.526\times {10}^{-5}$（$aL=2$）、$a=0.8$ 时 $4.203\times {10}^{-2}$（$aL=8$）。

后两行的迭代仍然收敛（误差有限），说明 $aL<1$ 是充分条件而非必要条件：压缩性保证迭代必收敛，不满足压缩性时迭代可能仍然收敛（本例中真实解在 $[0,a]$ 上存在，因此迭代有极限）。误差随 $a$ 增大而上升（从 ${10}^{-9}$ 增到 ${10}^{-2}$），与「区间越长、Lipschitz 常数越大、收敛越慢」一致。这条数据给出的实用判据是：定理的条件是「保证收敛」的充分条件，判断实际收敛性需要看问题本身；而在工程中，条件不满足往往意味着需要缩短时间步或改用隐式方法（5.13.4 节）。

5.13.2 线性系统与相平面

线性系统与矩阵指数

线性常系数系统 $y^{\prime }=Ay$（$A$ 是 $n\times n$ 常矩阵）的解可以直接写出：

$$y(t)=e^{At}y(0),e^{At}=\sum _{k=0}^{\mathrm{\infty }}\frac{(At{)}^k}{k!}$$

矩阵指数的级数与数量指数同形（5.7 节的幂级数结论（对矩阵同样成立）：收敛半径无穷）。若 $A$ 可对角化（$A=P\mathrm{\Lambda }{P}^{-1}$），则 $e^{At}=P{e}^{\mathrm{\Lambda }t}{P}^{-1}$，解是各特征方向上指数项的组合 $e^{{\lambda }_{i}t}$；特征值有重根时出现 $t^k{e}^{\lambda t}$ 型的项。

平衡点的分类

$y\equiv 0$ 是所有线性系统的平衡点（$f(0)=0$）。$A$ 的特征值决定轨道的行为，二维情形的分类是标准内容：

两实特征值同号：${\lambda }_1<{\lambda }_2<0$ 时是稳定结点（两条特征方向上的解都衰减，快慢不同），${\lambda }_1>{\lambda }_2>0$ 时是不稳定结点。

两实特征值异号：鞍点，一条方向吸引、一条方向排斥。

共轭复特征值 $\lambda =\alpha \pm i\beta$（$\beta \ne 0$）：$\alpha <0$ 时是稳定焦点（对数螺线向内），$\alpha >0$ 时是不稳定焦点；$\alpha =0$ 时是中心（闭轨道，围绕平衡点周期运动）。

判断稳定性只看实部：$\underset{i}{max}\mathrm{Re}{\lambda }_i<0$ 时渐近稳定，$>0$ 时不稳定，$=0$ 时临界（需要非线性项或更高阶信息才能确定，中心是典型例子）。

实验三用四个线性系统核对分类（初始条件 $(1,1)$，演化到 $T=10$）：

系统	特征值	$|y(T)|$	$y_1$ 的符号变化次数
稳定结点	$-1,-2$	$4.5\times {10}^{-5}$	$0$
鞍点	$+1,-1$	$21916.68$	$0$
稳定焦点	$-0.5\pm i$	$0.009565$	$3$
中心	$\pm i$	$1.421302$	$3$

四条数据的读法：$\parallel y(T)\parallel$ 把稳定性显示出来（稳定结点与稳定焦点衰减，鞍点增长，中心保持在初始的 $\sqrt{2}=1.414214$ 附近）；符号变化次数把「振荡与否」显示出来（结点与鞍点为 $0$，焦点与中心为 $3$）。两类信息合起来就能区分四种类型，这正是「特征值的实部决定稳定性、虚部决定振荡」这条判据的数值形态。

中心一行还有一个细节：$\parallel y(T)\parallel =1.421302$ 比初始值 $1.414214$ 大 $0.5\mathrm{\%}$。这项增长来自格式而非方程：显式 Euler 方法对中心型系统不保持能量，放大因子是 $\sqrt{1+h^2}$（可以验证 ${\Vert (I+hA)y\Vert }^2=(1+h^2)\parallel y{\parallel }^2$ 对 $A=\left(\begin{array}{cc}0& 1\\ -1& 0\end{array}\right)$ 成立），$N={10}^4$ 步后累积因子 $(1+{10}^{-6}{)}^{{10}^4}=1.01005$，与实测的增长一致（$\sqrt{1.01005}=1.00501$，$1.414214\times 1.00501=1.42130$ ✓）。这类现象叫数值耗散或数值反耗散：格式的稳定域包含了虚轴附近的点，但不在虚轴上，因此能量缓慢变化。长时间积分（如分子动力学、天体力学）必须改用辛格式或保能量格式，理由就在这条数据里。

非线性系统与线性化

非线性系统 $y^{\prime }=f(y)$ 的平衡点是 $f(y^{\ast })=0$ 的解。线性化定理（5.8 节的隐函数定理与 Taylor 公式的直接应用）：若 $J=Df(y^{\ast })$ 的所有特征值实部非零（双曲平衡点），则 $y^{\ast }$ 附近的轨道结构与线性系统 $z^{\prime }=Jz$ 的轨道结构拓扑等价，稳定性由 $J$ 的特征值决定。若存在实部为零的特征值（非双曲），线性化不足，需要更高阶的分析（中心的扰动可能变成稳定或不稳定的焦点，取决于高阶项）。

这条定理把非线性问题局部化为线性问题，是相平面分析的基本工具。它也是优化算法的收敛性分析的起点：梯度下降 ${\theta }_{k+1}={\theta }_k-\eta \mathrm{\nabla }V({\theta }_k)$ 在极小点附近的线性化是 ${\theta }_{k+1}-{\theta }^{\ast }\approx (I-\eta H)({\theta }_k-{\theta }^{\ast })$，收敛条件是 $|1-\eta {\lambda }_{max}(H)|<1$ 且 $|1-\eta {\lambda }_{min}(H)|<1$，即 $0<\eta <\frac{2}{{\lambda }_{max}}$。这条条件与第三章的步长选择、5.13.4 节的显式方法稳定条件是同一件事。

5.13.3 Lyapunov 函数与梯度流

Lyapunov 稳定性

Lyapunov 函数是沿轨道单调（不增）的函数，它的存在给出稳定性的判据。设 $V:{\mathbb{R}}^n\to \mathbb{R}$ 在平衡点 $y^{\ast }$ 的邻域内连续可微、$V(y^{\ast })=0$、$V(y)>0$（$y\ne y^{\ast }$），则

$$\dot{V}(y)\text{defeq}\mathrm{\nabla }V(y)\cdot f(y)\le 0\text{（区域内处处成立）}$$

蕴含 $y^{\ast }$ 稳定；若 $\dot{V}<0$（$y\ne y^{\ast }$），则 $y^{\ast }$ 渐近稳定。这条判据的技术要点是不需要解出方程：只需计算 $\dot{V}$ 的符号。

例子：线性系统 $y^{\prime }=Ay$ 取 $V=y^{\mathsf{T}}Py$（$P$ 对称正定），则 $\dot{V}=y^{\mathsf{T}}(A^{\mathsf{T}}P+PA)y$。存在 $P$ 使右端负定，等价于 $A$ 的特征值实部全为负（Lyapunov 方程 $A^{\mathsf{T}}P+PA=-Q$ 对任意正定 $Q$ 有唯一正定解）。这条等价性说明「二次 Lyapunov 函数」对线性系统的稳定性判别是完备的，也给出代数判据的来源。

梯度流

取 $f=-\mathrm{\nabla }V$，方程为 $\dot{y}=-\mathrm{\nabla }V(y)$，称 $V$ 的梯度流。它天然带 Lyapunov 函数：

$$\dot{V}=\mathrm{\nabla }V\cdot (-\mathrm{\nabla }V)=-{\Vert \mathrm{\nabla }V\Vert }^2\le 0$$

等号只在临界点处成立。由 Lyapunov 判据，梯度流的轨道收敛到临界点（在适当条件下收敛到极小点）。

实验四用双井势 $V(\theta )=\frac{({\theta }^2-1{)}^2}{4}$ 与起点 $\theta (0)=1.8$ 演示。$V$ 的极小点在 ${\theta }^{\ast }=\pm 1$，梯度流 $\dot{\theta }=-\theta ({\theta }^2-1)$ 从 $1.8$ 出发单调下降到 ${\theta }^{\ast }=1$：$\theta (T)=1.0000000000$、$V(\theta (T))=3.06\times {10}^{-27}$。单调性检验给出全部步上 $V$ 的增量都不大于零（最大增量 $0.000\times {10}^0$），与 $\dot{V}\le 0$ 的理论结论一致。

收敛速率由线性化给出：在 ${\theta }^{\ast }=1$ 处 $V^{″}(1)=2{\theta }^2\cdot \dots$ 直接计算得 $V^{″}(\theta )=3{\theta }^2-1$，$V^{″}(1)=2$；梯度流的线性化速率是 $V^{″}({\theta }^{\ast })=2$（因为 $\dot{\theta }\approx -V^{″}({\theta }^{\ast })(\theta -{\theta }^{\ast })$）。$V$ 的衰减速率是速率的 $2$ 倍（$V\sim (\theta -{\theta }^{\ast }{)}^2$），因此 $\ln V$ 的斜率应当是 $-4$。实测在区间 $3<t<6$ 上（$V$ 从 $7.224\times {10}^{-7}$ 降到 $4.406\times {10}^{-12}$）斜率是 $-4.004766$，与理论值相符。对照组用单井势 $V=\frac{{\theta }^2}{2}$（$V^{″}=1$），理论斜率 $-2$，实测 $-2.001001$。

这两条数据的实用含义：梯度流的收敛速率由极小点处 Hessian 的特征值决定，目标函数的衰减速率是它的两倍。这个关系在优化中直接对应「梯度下降在强凸函数上的线性收敛率」（第三章 3.3 节）：迭代格式的速率是 $1-\eta \lambda$，连续时间梯度流的速率是 $\lambda$，$V$ 的速率是 $2\lambda$，三者的对应关系使连续时间的分析可以用来指导步长的选择。

5.13.4 数值解法与刚性

显式方法与其稳定域

最简单的数值方法是显式 Euler：$y_{k+1}=y_k+hf(t_k,y_k)$，一阶精度。它的稳定性用测试方程 $y^{\prime }=\lambda y$（$\lambda$ 为复数）分析：放大因子是 $R(z)=1+z$（$z=h\lambda$），稳定条件是 $|1+z|\le 1$，即 $z$ 落在以 $-1$ 为圆心、半径 $1$ 的圆盘内。实数 $\lambda <0$ 时给出 $h\le \frac{2}{|\lambda |}$。

四阶 Runge-Kutta（RK4）的放大因子是 $R(z)=1+z+\frac{z^2}{2}+\frac{z^3}{6}+\frac{z^4}{24}$，稳定域在实轴上约 $|z|\le 2.785$，比显式 Euler 宽一些，但仍然是有界的。

刚性

若系统中同时存在时间尺度差异很大的分量（特征值的模相差多个数量级），称问题是刚性的。刚性比定义为

$$\kappa =\frac{\underset{i}{max}\left|\mathrm{Re}{\lambda }_i\right|}{\underset{i}{min}\left|\mathrm{Re}{\lambda }_i\right|}$$

刚性问题的困难在于：解的光滑部分由慢尺度决定（$|\lambda |$ 小的分量），但显式方法的步长被快尺度限制（$h\lesssim \frac{2}{|\lambda {|}_{max}}$），因此需要大量步数才能覆盖感兴趣的区间。

实验二用 $\kappa =1000$ 的问题演示：

$$y^{\prime }=-1000(y-\cos t)-\sin t,y(0)=1,y_{\text{exact}}=\cos t$$

解是光滑的余弦（快分量在瞬态之后消失），但方程的特征值是 $-1000$。显式 Euler 的稳定条件 $h\le \frac{2}{1000}=0.002$。数据：$h=0.01$ 时误差 $3.528\times {10}^{185}$（发散，且发散速度极快——每一步放大 $|1-1000h|=9$）；$h=0.002$（临界）时误差 $2.000\times {10}^{-6}$；$h=0.001$ 时 $5.000\times {10}^{-7}$；$h=0.0005$ 时 $2.500\times {10}^{-7}$；$h=0.0001$ 时 $5.000\times {10}^{-8}$。稳定区间的误差按 $O(h)$ 下降（$h$ 减半误差减半），与此同时也说明即使用到稳定极限，显式方法的精度要求仍然把步长压得很小。

隐式 Euler（$y_{k+1}=y_k+hf(t_{k+1},y_{k+1})$）的放大因子是 $R(z)=\frac{1}{1-z}$，稳定域是「以 $1$ 为圆心的单位圆外部」，包含整个左半平面——A 稳定。数据：同一个 $h=0.01$ 给出的误差是 $4.996\times {10}^{-6}$（与显式方法在 $h=0.002$ 时的精度相当）；$h$ 从 $0.01$ 降到 $0.0001$ 时误差从 $4.996\times {10}^{-6}$ 降到 $5.000\times {10}^{-8}$，始终按 $O(h)$ 下降。这条对照是刚性问题数值解法的核心结论：隐式方法的每步代价更高（需要解方程或线性方程组），但它允许大得多的步长，在刚性问题上总代价通常远低于显式方法。

方法的选取

三条实用原则。

看刚性比：$\kappa$ 在 $10$ 以内时显式方法（RK4 等）通常够用；$\kappa$ 超过 ${10}^3$ 时显式方法的步长限制开始主导代价，应当改用隐式或半隐式方法（隐式 Euler、梯形法、BDF、Rosenbrock 方法）。

看精度与稳定性的平衡：隐式 Euler 只有一阶精度，实际常用二阶的梯形法（A 稳定但不 L 稳定，对极刚分量处理不够好）或 BDF2、Radau 等方法（高阶且 L 稳定）。

检查数值耗散：对需要长时间积分的守恒系统（5.13.2 节的中心型轨道），格式的耗散性质比精度阶更重要，应当用辛格式或保结构格式；对刚性衰减问题，需要格式在 $z\to -\mathrm{\infty }$ 时 $R(z)\to 0$（L 稳定），否则快分量会以虚假的振荡残留下来。实验四的边界层问题（5.12 节）与这里讨论的是同一类困难在不同形式下的表现。

5.13.5 与其他章节的联系

常微分方程与动力系统的视角把本章的许多内容串起来。

与 5.1 节：存在唯一性用压缩映射原理证明，Lipschitz 条件、完备性与不动点迭代是同一条线索。

与 5.4 与 5.8 节：数值解法的精度阶用 Taylor 展开分析（局部截断误差），相平面分类用 Jacobian 的特征值，Lyapunov 函数的构造用二次型与梯度。

与 5.12 节：刚性问题的困难与边界层问题同源（小参数导致的时间或空间尺度分离），处理方法也平行（隐式方法对应指数拟合格式，后者是「在单元内精确求解」的思路）。

与第三章：梯度下降是最优化算法，同时是一条梯度流（连续时间极限 $h\to 0$ 的对偶）；收敛性分析的两条路线（离散的步长条件、连续的 Lyapunov 分析）在这里合流。Adam、动量法等算法的连续时间极限是被修改的动力学系统，稳定性分析的方法与本节相同。

与第二章：随机微分方程是常微分方程加噪声项，其数值解法（Euler-Maruyama、Milstein）的稳定性分析沿用本节的框架；扩散过程与 Fokker-Planck 方程的联系也建立在这里。

5.13.6 数值实验

实验一：Picard 迭代与解的存在唯一性

第一组用线性问题观察迭代的收敛形态，第二组用非线性问题观察压缩条件与解的爆破。

import math
import numpy as np

print("(a) y' = y，y(0) = 1（精确解 e^x），Picard 迭代 y_{n+1} = 1 + int_0^x y_n dt")
xg = np.linspace(0, 1, 2001)
y = np.ones_like(xg)
print(f"{'n':>5} {'最大误差':>16} {'误差比值 (n-1 -> n)':>22}")
prev = None
for n in range(1, 9):
    y = 1.0 + np.array([np.trapezoid(y[:i + 1], xg[:i + 1]) for i in range(len(xg))])
    e = float(np.max(np.abs(y - np.exp(xg))))
    r = "-" if prev is None else f"{prev / e:.2f}"
    print(f"{n:>5} {e:>16.3e} {r:>22}")
    prev = e

print("(b) y' = y^2，y(0) = 1（精确解 1/(1-x)，在 x = 1 处爆破）")
print(f"{'区间右端 a':>10} {'迭代 8 次后的最大误差':>24} {'Lipschitz 常数 2/(1-a)':>24} {'a*L':>10}")
for a in [0.2, 0.3, 0.5, 0.8]:
    xg2 = np.linspace(0, a, 2001)
    y2 = np.ones_like(xg2)
    for n in range(8):
        y2 = 1.0 + np.array([np.trapezoid(y2[:i + 1] ** 2, xg2[:i + 1])
                             for i in range(len(xg2))])
    exact = 1.0 / (1.0 - xg2)
    L = 2.0 / (1.0 - a)
    print(f"{a:>10.1f} {float(np.max(np.abs(y2 - exact))):>24.6e} {L:>24.4f} {a * L:>10.4f}")

输出：

(a) y' = y，y(0) = 1（精确解 e^x），Picard 迭代 y_{n+1} = 1 + int_0^x y_n dt
    n             最大误差        误差比值 (n-1 -> n)
    1        7.183e-01                      -
    2        2.183e-01                   3.29
    3        5.162e-02                   4.23
    4        9.948e-03                   5.19
    5        1.615e-03                   6.16
    6        2.262e-04                   7.14
    7        2.780e-05                   8.14
    8        3.002e-06                   9.26
(b) y' = y^2，y(0) = 1（精确解 1/(1-x)，在 x = 1 处爆破）
    区间右端 a             迭代 8 次后的最大误差     Lipschitz 常数 2/(1-a)        a*L
       0.2             1.304206e-09                   2.5000     0.5000
       0.3             5.565785e-08                   2.8571     0.8571
       0.5             2.526338e-05                   4.0000     2.0000
       0.8             4.203331e-02                  10.0000     8.0000

(a) 迭代的收敛形态。最大误差从 $7.183\times {10}^{-1}$（$n=1$）单调降到 $3.002\times {10}^{-6}$（$n=8$），相邻比值 $3.29$、$4.23$、$5.19$、$6.16$、$7.14$、$8.14$、$9.26$ 逐次增加约 $1$。这条规律与解析解一致：$y^{\prime }=y$ 的迭代给出 $y_n(t)=\sum _{k=0}^n\frac{t^k}{k!}$，误差在 $t=1$ 处约 $\frac{1}{(n+1)!}$，相邻比值约 $n+1$。收敛速度由级数尾部的量级决定，$n=8$ 时误差已在 ${10}^{-6}$ 量级，比线性收敛快得多（阶乘下降）。

(b) 压缩条件与爆破。四行数据把 $aL$ 从 $0.5$ 增到 $8$：$a=0.2$ 时 $aL=0.5$（压缩条件满足），迭代误差 $1.304\times {10}^{-9}$；$a=0.5$ 时 $aL=2$，误差 $2.526\times {10}^{-5}$；$a=0.8$ 时 $aL=8$，误差 $4.203\times {10}^{-2}$。误差随区间变长而增大（三个数量级），与「Lipschitz 常数随 $a$ 增大而增大、收敛变慢」一致。

要注意的是后两行的迭代仍然收敛，尽管 $aL>1$。这解释了 $aL<1$ 的地位：它是压缩映射定理要求的充分条件（保证迭代必收敛且有不依赖具体问题的收敛率），而不是必要条件。本例中真实解在 $[0,a]$（$a<1$）上存在且有界，因此迭代有极限，只是收敛速度不受定理保证。工程上的含义是：定理条件不满足时不能断言发散，但也没有保证，需要通过其他方式（例如数值试验、缩短积分区间、改用隐式方法）来处理。

实验二：刚性方程与显式/隐式方法

用刚性比 $1000$ 的线性问题对比显式与隐式 Euler 的稳定性。

import math
import numpy as np

print("y' = -1000 (y - cos t) - sin t，y(0) = 1，精确解 y = cos t（刚性比 1000）")
print("显式 Euler 的稳定条件 |1 - 1000 h| <= 1，即 h <= 0.002")

def solve_stiff(h, T=2.0, scheme="explicit"):
    N = int(round(T / h))
    ts = np.linspace(0, T, N + 1)
    y = np.zeros(N + 1)
    y[0] = 1.0
    lam = 1000.0
    for i in range(N):
        tn, tn1 = ts[i], ts[i + 1]
        if scheme == "explicit":
            y[i + 1] = y[i] + h * (-lam * (y[i] - math.cos(tn)) - math.sin(tn))
        else:
            y[i + 1] = (y[i] + h * (lam * math.cos(tn1) - math.sin(tn1))) / (1.0 + lam * h)
    return ts, y

print(f"{'h':>10} {'显式 Euler 最大误差':>22} {'显式结果是否发散':>18} "
      f"{'隐式 Euler 最大误差':>22}")
for h in [0.01, 0.002, 0.001, 0.0005, 0.0001]:
    ts, ye = solve_stiff(h, scheme="explicit")
    ts, yi = solve_stiff(h, scheme="implicit")
    ex = np.cos(ts)
    diverge = "是" if float(np.max(np.abs(ye))) > 10 else "否"
    print(f"{h:>10} {float(np.max(np.abs(ye - ex))):>22.3e} {diverge:>18} "
          f"{float(np.max(np.abs(yi - ex))):>22.3e}")

输出：

y' = -1000 (y - cos t) - sin t，y(0) = 1，精确解 y = cos t（刚性比 1000）
显式 Euler 的稳定条件 |1 - 1000 h| <= 1，即 h <= 0.002
         h          显式 Euler 最大误差           显式结果是否发散          隐式 Euler 最大误差
      0.01             3.528e+185                  是              4.996e-06
     0.002              2.000e-06                  否              9.998e-07
     0.001              5.000e-07                  否              4.999e-07
    0.0005              2.500e-07                  否              2.500e-07
    0.0001              5.000e-08                  否              5.000e-08

$h=0.01$ 一行把显式方法的失效写成极端形态：误差 $3.528\times {10}^{185}$，即数值解已经发散（放大因子 $|1-1000\times 0.01|=9$，$200$ 步之后 $9^{200}\approx {10}^{190}$，与实测误差的量级一致）。$h=0.002$ 是稳定条件的临界点，误差 $2.000\times {10}^{-6}$，再缩小步长时误差按 $O(h)$ 下降（$5.000\times {10}^{-7}$、$2.500\times {10}^{-7}$、$5.000\times {10}^{-8}$，每减半误差减半）。

隐式 Euler 的表现完全不同：$h=0.01$（超出显式稳定极限 $5$ 倍）时的误差是 $4.996\times {10}^{-6}$，与显式方法在临界步长处的精度相当；$h$ 从 $0.01$ 降到 $0.0001$ 时误差从 $4.996\times {10}^{-6}$ 降到 $5.000\times {10}^{-8}$，严格按 $O(h)$ 下降，没有任何稳定性限制。两条数据的对照给出刚性问题的核心取舍：隐式方法每步需要求解方程（线性问题是一次除法或一次线性方程组求解，非线性问题需要 Newton 迭代），但允许的步长比显式方法大几个数量级，总代价通常低得多。粗算：显式方法要覆盖 $T=2$ 且满足 $h\le 0.002$ 至少需要 $1000$ 步；隐式方法取 $h=0.01$ 只需 $200$ 步，且误差还略好。

实验三：相平面与线性系统平衡点的分类

对四个线性系统计算特征值、终态范数与振荡指标。

import numpy as np

systems = [
    ("稳定结点", np.array([[-1.0, 0.0], [0.0, -2.0]])),
    ("鞍点", np.array([[1.0, 0.0], [0.0, -1.0]])),
    ("稳定焦点", np.array([[-0.5, 1.0], [-1.0, -0.5]])),
    ("中心", np.array([[0.0, 1.0], [-1.0, 0.0]])),
]
h = 0.001
T = 10.0
N = int(T / h)
print(f"{'类型':>10} {'特征值':>28} {'|y(T)|':>14} {'y1 的符号变化次数':>18}")
for name, A in systems:
    ev = np.linalg.eigvals(A)
    y = np.zeros((N + 1, 2))
    y[0] = [1.0, 1.0]
    for i in range(N):
        y[i + 1] = y[i] + h * (A @ y[i])
    sign_changes = int(np.sum(np.diff(np.sign(y[1:, 0])) != 0))
    evs = ", ".join(f"{v.real:+.4f}{v.imag:+.4f}j" for v in ev)
    print(f"{name:>10} {evs:>28} {float(np.linalg.norm(y[-1])):>14.6f} {sign_changes:>18}")

输出：

        类型                          特征值         |y(T)|         y1 的符号变化次数
      稳定结点 -1.0000+0.0000j, -2.0000+0.0000j       0.000045                  0
        鞍点 +1.0000+0.0000j, -1.0000+0.0000j   21916.681339                  0
      稳定焦点 -0.5000+1.0000j, -0.5000-1.0000j       0.009565                  3
        中心 +0.0000+1.0000j, +0.0000-1.0000j       1.421302                  3

四行数据把分类的两条判据都显示出来了。稳定性由 $\parallel y(T)\parallel$ 决定：稳定结点（特征值 $-1,-2$）衰减到 $4.5\times {10}^{-5}$，稳定焦点（$-0.5\pm i$）衰减到 $0.009565$，鞍点（$+1,-1$）增长到 $21916.68$，中心（$\pm i$）保持在初始值 $\sqrt{2}=1.414214$ 的量级（末值 $1.421302$）。振荡与否由符号变化次数决定：结点与鞍点为 $0$（单调接近或远离），焦点与中心为 $3$（轨道绕平衡点转动）。

中心一行的末值比初值大 $0.5\mathrm{\%}$（$1.421302$ 对 $1.414214$），这项增长属于格式的性质：显式 Euler 对 $A=\left(\begin{array}{cc}0& 1\\ -1& 0\end{array}\right)$ 的放大因子是 $\sqrt{1+h^2}$（由 ${\Vert (I+hA)y\Vert }^2=(1+h^2)\parallel y{\parallel }^2$ 得），$N={10}^4$ 步的累积因子是 $(1+{10}^{-6}{)}^{{10}^4}=1.01005$，开方得 $1.00501$，$1.414214\times 1.00501=1.42130$ ✓ 与实测一致。中心型轨道上显式 Euler 会让振幅缓慢增长（数值反耗散），这是长时间积分守恒系统时必须使用保结构格式的原因。

实验四：梯度流与 Lyapunov 函数

对双井势与单井势分别运行梯度流，检验 Lyapunov 单调性与收敛速率。

import numpy as np

print("V(theta) = (theta^2 - 1)^2 / 4，梯度流 dtheta/dt = -theta(theta^2 - 1)")
print("起点 theta(0) = 1.8，向极小点 theta* = 1 收敛（线性化速率 2）")
h = 0.001
T = 20.0
N = int(T / h)
ts = np.linspace(0, T, N + 1)
th = np.zeros(N + 1)
th[0] = 1.8
for i in range(N):
    d = th[i] * (th[i] ** 2 - 1.0)
    th[i + 1] = th[i] - h * d
V = (th ** 2 - 1.0) ** 2 / 4.0
dV = np.diff(V)
print(f"  V(theta(T)) = {V[-1]:.6e}，theta(T) = {th[-1]:.10f}")
print(f"  单调性检验：所有步上 V 的增量 <= 0，最大增量 = {float(np.max(dV)):.3e}")
sel = (ts > 3.0) & (ts < 6.0)
slope = float(np.polyfit(ts[sel], np.log(V[sel]), 1)[0])
print(f"  区间 3 < t < 6 上 ln V 的斜率 = {slope:.6f}（理论值 -4）")
print(f"  该区间内 V 从 {float(V[sel][0]):.3e} 降到 {float(V[sel][-1]):.3e}")

print("V(theta) = theta^2/2（单井），起点 theta(0) = 2")
th2 = np.zeros(N + 1)
th2[0] = 2.0
for i in range(N):
    th2[i + 1] = th2[i] - h * th2[i]
V2 = th2 ** 2 / 2
slope2 = float(np.polyfit(ts[sel], np.log(V2[sel]), 1)[0])
print(f"  区间 3 < t < 6 上 ln V 的斜率 = {slope2:.6f}（理论值 -2）")
print(f"  V(theta(T)) = {V2[-1]:.6e}，theta(T) = {th2[-1]:.12f}")

输出：

V(theta) = (theta^2 - 1)^2 / 4，梯度流 dtheta/dt = -theta(theta^2 - 1)
起点 theta(0) = 1.8，向极小点 theta* = 1 收敛（线性化速率 2）
  V(theta(T)) = 3.056885e-27，theta(T) = 1.0000000000
  单调性检验：所有步上 V 的增量 <= 0，最大增量 = 0.000e+00
  区间 3 < t < 6 上 ln V 的斜率 = -4.004766（理论值 -4）
  该区间内 V 从 7.224e-07 降到 4.406e-12
V(theta) = theta^2/2（单井），起点 theta(0) = 2
  区间 3 < t < 6 上 ln V 的斜率 = -2.001001（理论值 -2）
  V(theta(T)) = 8.328351e-18，theta(T) = 0.000000004081

双井势的轨道从 $1.8$ 单调下降到 $1.0000000000$（$V$ 降到 $3.06\times {10}^{-27}$），没有越过极小点也没有跳到另一个井 $\theta =-1$：梯度流沿最陡下降方向运动，在双井势上会停在最近的井中，去哪个井由初值决定（这条性质在机器学习中对应「不同的初始化收敛到不同的解」）。单调性检验给出所有步上 $V$ 的增量都不大于零（最大增量恰好是 $0.000\times {10}^0$），与 $\dot{V}=-\parallel \mathrm{\nabla }V{\parallel }^2\le 0$ 一致。

收敛速率的两条数据给出连续时间版本的「强凸收敛率」：双井势在 ${\theta }^{\ast }=1$ 处 $V^{″}(1)=3{\theta }^2-1=2$，因此 $\ln V$ 的斜率理论值是 $-2V^{″}({\theta }^{\ast })=-4$，实测 $-4.004766$（$V$ 在该区间从 $7.224\times {10}^{-7}$ 降到 $4.406\times {10}^{-12}$，跨越五个半数量级，斜率估计可靠）；单井势 $V^{″}=1$，理论斜率 $-2$，实测 $-2.001001$。目标函数的衰减速率等于极小点处 Hessian 特征值的两倍，这条关系把连续时间的分析与离散算法的收敛率联系起来（第三章 3.3 节的强凸情形：迭代误差按 $(1-\eta \lambda {)}^k$ 衰减，$\eta$ 小时与 $e^{-\eta \lambda k}$ 一致，而 $\eta k$ 对应连续时间 $t$）。

5.13.7 本节小结

常微分方程的四层内容

本节把常微分方程与动力系统的视角组织为四层。存在唯一性层给出 Picard-Lindelöf 定理：把初值问题写成积分方程 $y(t)=y_0+{\int }_{t_0}^{t}f(s,y(s))ds$，Picard 迭代 $y_{k+1}=y_0+\int f(s,y_k)ds$ 在 $L\delta <1$ 时是压缩映射（$L$ 是 Lipschitz 常数、$\delta$ 是时间区间长度），由压缩映射原理得存在唯一；实验一的数据：$y^{\prime }=y$ 的迭代误差从 $7.183\times {10}^{-1}$ 降到 $3.002\times {10}^{-6}$，相邻比值按 $n+1$ 增长（阶乘级数的尾部），$y^{\prime }=y^2$ 的迭代误差随区间长度 $a$ 从 $1.304\times {10}^{-9}$（$aL=0.5$）增到 $4.203\times {10}^{-2}$（$aL=8$），可见 $aL<1$ 是充分条件而非必要条件。相平面层给出线性系统的分类：$y^{\prime }=Ay$ 的解是 $e^{At}y(0)$，特征值的实部决定稳定性、虚部决定振荡（实验三：稳定结点衰减到 $4.5\times {10}^{-5}$、鞍点增长到 $21916.68$、稳定焦点衰减到 $0.009565$、中心保持 $1.414214$ 量级，符号变化次数把结点与焦点区分开）；非线性系统用 Jacobian 线性化，双曲平衡点（特征值实部非零）的局部结构与线性系统等价。Lyapunov 层给出不需要解方程的稳定性判据（$\dot{V}\le 0$ 且 $V>0$），梯度流 $\dot{y}=-\mathrm{\nabla }V$ 满足 $\dot{V}=-\parallel \mathrm{\nabla }V{\parallel }^2\le 0$（实验四：$V$ 的增量全步非正，双井势与单井势的 $\ln V$ 斜率分别为 $-4.004766$ 与 $-2.001001$，与理论值 $-2V^{″}({\theta }^{\ast })$ 一致）。数值解法层给出刚性问题的核心困难与对策：显式 Euler 的稳定条件是 $|1+h\lambda |\le 1$（实轴上 $h\le 2/|\lambda |$），隐式 Euler 的放大因子 $\frac{1}{1-h\lambda }$ 对任意 $h$ 稳定（A 稳定）；实验二的数据（刚性比 $1000$）：显式方法在 $h=0.01$ 时误差 $3.528\times {10}^{185}$（发散），隐式方法在同一个步长下误差 $4.996\times {10}^{-6}$，且步长减小十倍误差减小十倍（$O(h)$）。中心的显式 Euler 积分还给出了「数值反耗散」的例子：振幅按 $(1+h^2{)}^{N/2}$ 缓慢增长（$T=10$ 时 $+0.5\mathrm{\%}$），说明守恒系统的长时间积分需要保结构格式。

至此第五章完成。本章从实数与极限出发，经过数列与级数、连续与微分、积分与测度、多元微积分与变分、渐近分析，最终回到常微分方程与动力系统。这条路线把「极限」这一个工具在六个方向上的展开讲完了：数、函数、函数列、多元、测度、参数极端。第四章的信息论中的核心量（熵、互信息、KL 散度）都以本章积分的严格理论为基础，其中 KL 散度的定义需要 5.10 节的 Radon-Nikodym 定理才能说清楚；而当信息论与最优化用于分析具体算法与模型时，本章的一致收敛条件（5.6 节）、极限交换定理（5.10 节）与渐近分析（5.12 节）会成为判断结论有效性的依据。第三章的梯度下降、近端方法与对偶理论，分别对应本节的梯度流、非光滑动力系统与 Lagrangian 动力学，三条线索在最优化、动力系统与数值分析之间构成一组可以互相翻译的语言。

练习题

第 1 题 概念推导

(a) 把二阶方程 $y^{″}+{\omega }^{2}y=0$（谐振子）化为一阶方程组，求系数矩阵的特征值并判断平衡点的类型。(b) 说明为什么 $y^{\prime }=\sqrt{|y|}$、$y(0)=0$ 的解不唯一，并指出 Picard-Lindelöf 定理的哪一条假设被违反。(c) 证明梯度流满足 $\dot{V}=-\parallel \mathrm{\nabla }V{\parallel }^2$，并说明由此能推出什么结论。

参考答案

(a) 令 $v=y^{\prime }$，得

$$y^{\prime }=v,v^{\prime }=-{\omega }^{2}y⟹{\left(\begin{array}{c}y\\ v\end{array}\right)}^{\prime }=\left(\begin{array}{cc}0& 1\\ -{\omega }^2& 0\end{array}\right)\left(\begin{array}{c}y\\ v\end{array}\right)$$

特征方程是 ${\lambda }^2+{\omega }^2=0$，特征值 $\lambda =\pm i\omega$ 是纯虚数，实部为零。平衡点是中心，轨道是闭曲线（椭圆）。由守恒量可以验证：$\frac{d}{dt}(\frac{v^2}{2}+\frac{{\omega }^2{y}^2}{2})=v{v}^{\prime }+{\omega }^{2}y{y}^{\prime }=v(-{\omega }^{2}y)+{\omega }^{2}yv=0$，即能量 $\frac{v^2}{2}+\frac{{\omega }^2{y}^2}{2}$ 沿轨道守恒。中心属于非双曲平衡点（特征值实部为零），线性化不能判断稳定性（本例中轨道有界但也不收敛到原点，属于「稳定但不渐近稳定」，即 Lyapunov 稳定）。非线性扰动会改变中心的类型（例如加入阻尼 $y^{″}+c{y}^{\prime }+{\omega }^{2}y=0$ 后特征值是 $-\frac{c}{2}\pm i\sqrt{{\omega }^2-c^2/4}$，变成稳定焦点），这正是非双曲情形需要高阶分析的原因。

(b) $f(y)=\sqrt{|y|}$ 在 $y=0$ 处不满足 Lipschitz 条件：对任意 $L$，取 $y_1=ϵ$、$y_2=0$ 得 $\frac{|f(y_1)-f(y_2)|}{|y_1-y_2|}=\frac{\sqrt{ϵ}}{ϵ}={ϵ}^{-1/2}$，$ϵ\to 0$ 时无界，因此不存在统一的 $L$。失唯一性的机制是：方程的两个解 $y\equiv 0$ 与 $y=\frac{t^2}{4}$（$t\ge 0$）在 $t=0$ 处相切（导数都为 $0$），轨道在平衡点处「分叉」而不违反方程。物理上的图像是：$f$ 在零点的斜率无穷（摆离平衡的速度不受限制），质点可以在零点停留任意长的时间再出发。

被违反的假设是 Lipschitz 条件（$f$ 的连续性仍然成立，因此存在性定理（Peano 定理，用 Arzelà-Ascoli 定理证明，5.6 节）仍然给出解的存在，只有唯一性失效）。检查方法：算出 $f^{\prime }$ 或在可疑点附近估计差商的上界，看是否有界。

(c) $V$ 沿轨道的时间导数（链式法则）：$\dot{V}=\frac{d}{dt}V(y(t))=\mathrm{\nabla }V(y)\cdot y^{\prime }(t)$。代入 $y^{\prime }=-\mathrm{\nabla }V$ 得

$$\dot{V}=\mathrm{\nabla }V\cdot (-\mathrm{\nabla }V)=-{\Vert \mathrm{\nabla }V\Vert }^2\le 0$$

由此可以推出三条结论：一是 $V$ 沿轨道单调不增（Lyapunov 函数的核心性质），二是等号只在临界点（$\mathrm{\nabla }V=0$）处成立，因此只要轨道不停在临界点，$V$ 就严格下降，三是 $V$ 有下界时 $\dot{V}$ 的积分有限，${\int }_0^{\mathrm{\infty }}\parallel \mathrm{\nabla }V{\parallel }^{2}dt<\mathrm{\infty }$，从而 $\parallel \mathrm{\nabla }V\parallel$ 在某种意义下趋于零（梯度趋于零），再由 $V$ 的解析性与轨道的有界性（适当条件下）可以证明轨道收敛到临界点。注意「收敛到临界点」与「收敛到极小点」是两件事：鞍点也是临界点，梯度流从鞍点出发（或沿稳定流形接近它）时不会离开，实际中初值的扰动会把它推向极小点方向。这条区别在优化中对应「鞍点问题的困难」（第三章 3.8 节）。

第 2 题 计算推理

(a) 求 $y^{\prime }=-5y$ 的显式 Euler 与隐式 Euler 的稳定步长范围，并比较两者在 $h=0.5$ 时的行为。(b) 对系统 ${\left(\begin{array}{c}{y}_1\\ y_2\end{array}\right)}^{\prime }=\left(\begin{array}{cc}-1& 1\\ 0& -1000\end{array}\right)\left(\begin{array}{c}{y}_1\\ y_2\end{array}\right)$，求刚性比并判断显式 Euler 的步长限制。(c) 求 $V(\theta )=\frac{{\theta }^4}{4}-\frac{{\theta }^2}{2}$（双井势的另一种形式）的极小点与线性化速率。

参考答案

(a) 显式 Euler 的放大因子 $R(z)=1+z$（$z=-5h$），稳定条件 $|1-5h|\le 1$ 给出 $0\le h\le 0.4$。隐式 Euler 的放大因子 $R(z)=\frac{1}{1+5h}$，对任意 $h>0$ 都有 $|R|<1$，无条件稳定。

在 $h=0.5$ 时：显式 Euler 的放大因子是 $1-5\times 0.5=-1.5$，$|R|=1.5>1$，数值解按 ${1.5}^n$ 发散（并且交替变号）；隐式 Euler 的放大因子是 $\frac{1}{1+2.5}=0.2857$，数值解按 ${0.2857}^n$ 收敛到零。两种方法在同一个步长下的结论相反，这是 A 稳定性的具体含义：隐式方法能处理显式方法无法处理的步长。代价是隐式方法每步需要求解方程（本例中是一次除法），精度都是一阶。

(b) 矩阵是下三角，特征值是 $-1$ 与 $-1000$（对角元）。刚性比 $\kappa =1000$。显式 Euler 的稳定条件由模最大的特征值决定：$|1-1000h|\le 1$，即 $h\le 0.002$。解的光滑分量（$e^{-t}$）在 $t\in [0,5]$ 上变化很小，但显式方法仍必须用 $h\le 0.002$，需要 $2500$ 步以上才能覆盖这个区间。刚性比的实用含义是「步长限制与解的变化尺度之间的失配倍数」：这里失配 $1000$ 倍，即显式方法需要多花 $1000$ 倍的步数来完成同等精度的工作。若改用隐式方法（如隐式 Euler 或 BDF2），步长可以由精度要求决定（例如 $h=0.05$），总步数降到 $100$ 量级。

(c) $V^{\prime }(\theta )={\theta }^3-\theta =\theta ({\theta }^2-1)$，临界点是 $\theta =0,\pm 1$。$V^{″}(\theta )=3{\theta }^2-1$，在 $\theta =\pm 1$ 处 $V^{″}=2>0$（极小），在 $\theta =0$ 处 $V^{″}=-1<0$（极大）。因此极小点是 ${\theta }^{\ast }=\pm 1$，$V(\pm 1)=\frac{1}{4}-\frac{1}{2}=-\frac{1}{4}$；线性化速率是 $V^{″}({\theta }^{\ast })=2$，$V-V^{\ast }$ 的衰减速率是 $2V^{″}=4$，与本节实验四的双井势结论相同（两种形式的势函数只差一个平移与比例，动力学的速率由二阶导数决定）。这条计算说明收敛速率只依赖极小点处的二阶信息，势函数的整体形状（井的深浅、势垒高低）只影响「收敛到哪个井」与「瞬态阶段的长短」。

第 3 题 代码验证

用本节代码做三组实验并解释：(a) 在实验三里把中心系统的 $h$ 从 $0.001$ 增到 $0.01$，观察振幅的增长比例是否与 $(1+h^2{)}^{N/2}$ 相符；(b) 在实验二里把刚性比从 $1000$ 改成 $10000$，重新确定显式 Euler 的临界步长，并说明隐式方法的误差是否改变；(c) 在实验四里把起点改成 $\theta (0)=0.5$，观察轨道收敛到哪个极小点，并说明为什么。

参考答案

(a) 放大因子是 $\sqrt{1+h^2}$，$N=T/h=10/h$ 步的累积放大是 $(1+h^2{)}^{T/(2h)}\approx e^{hT/2}$。$h=0.01$ 时 $e^{0.01\times 10/2}=e^{0.05}=1.0513$，即振幅增长 $5.1\mathrm{\%}$（$h=0.001$ 时是 $e^{0.005}=1.0050$，增长 $0.5\mathrm{\%}$，与实验数据相符）。增大步长会加速数值反耗散，因为每一步的相对增长 $h^2$ 与步数的乘积随 $h$ 线性增长。这条对照的实用含义是：判断格式的长期行为要看「$N$ 步累积的放大」，而不是单步放大因子与 $1$ 的差（后者是 $O(h^2)$，看起来很小，但乘以 $N=T/h$ 后变成 $O(h)$）。

(b) 刚性比 $10000$ 时显式 Euler 的稳定条件是 $|1-10000h|\le 1$，即 $h\le 0.0002$，比原来小 $10$ 倍。隐式 Euler 的误差不变：它的放大因子 $\frac{1}{1+10000h}$ 仍然对任意 $h$ 稳定，而误差的主项由方法的截断误差决定（一阶方法误差约 $Ch$，$C$ 由解的二阶导数决定），与刚性比无关。因此数据表会显示：显式列在 $h\le 0.0002$ 时正常（误差按 $O(h)$ 下降），$h>0.0002$ 时发散；隐式列在所有 $h$ 上误差都是 $0.5h$ 量级。这条练习给出刚性问题的判断依据：刚性比只影响显式方法的步长限制，不影响隐式方法的精度，因此刚性比越大，隐式方法的相对优势越明显。

(c) $\theta (0)=0.5$ 位于 $(0,1)$ 之间，$V^{\prime }(\theta )=\theta ({\theta }^2-1)<0$ 在该区间上（$\theta >0$、${\theta }^2<1$），梯度流 $\dot{\theta }=-V^{\prime }(\theta )>0$，因此 $\theta$ 单调增加到 ${\theta }^{\ast }=1$。轨道收敛到 $1$ 而不是 $-1$，尽管起点 $0.5$ 到两个极小点的距离分别是 $0.5$ 与 $1.5$（近的那个被选中）。要收敛到 $-1$ 需要 $\theta (0)<0$。这条性质是「梯度流不跨越势垒」的表现：势垒在 $\theta =0$ 处（$V(0)=0$ 而 $V(\pm 1)=-1/4$，$\theta =0$ 是局部极大），梯度流沿最陡下降方向运动，无法越过势垒，因此收敛结果由初值所在的「吸引域」决定。离散算法（梯度下降）在大步长时可能因为数值误差跳过势垒，这与连续时间的结论不同，也是随机初始化的理论依据之一。

第 4 题 综合应用

某数值程序需要求解一个由 $n=50$ 个状态变量组成的常微分方程组（来自某物理模型的半离散化），时间区间 $[0,100]$。已知系统的 Jacobian 特征值分布在 $[-{10}^5,-{10}^{-3}]$。请设计求解方案：给出刚性的量化判断、方法选择（显式或隐式、具体格式）、步长策略（固定或自适应）、以及精度与代价的估计。

参考答案

刚性的量化判断。刚性比 $\kappa =\frac{{10}^5}{{10}^{-3}}={10}^8$，属于典型的刚性问题。显式方法的步长限制由模最大的特征值决定：对显式 Euler 是 $h\le \frac{2}{{10}^5}=2\times {10}^{-5}$，对 RK4 是 $h\le \frac{2.785}{{10}^5}=2.8\times {10}^{-5}$。要覆盖 $[0,100]$ 需要 $5\times {10}^6$ 步（显式 Euler）或 $3.6\times {10}^6$ 步（RK4），每步 $50$ 次函数求值（$n=50$），总求值次数在 $2\times {10}^8$ 量级——代价过高。结论：不能使用显式方法（除非只关心 ${10}^{-3}$ 量级以前的时间）。

方法选择。刚性问题的标准选择是隐式多步或隐式 Runge-Kutta（IRK）。三条候选：BDF（向后差分公式）方法的 BDF2 到 BDF5（L 稳定，阶数 2 到 5，每步需要解一次隐式方程，用 Newton 迭代，Jacobian 可以复用多个步）；Radau IIA（L 稳定的 IRK，3 阶或 5 阶，每步需要解 $3n$ 或 $5n$ 维的非线性方程组）；梯形法（A 稳定但不 L 稳定，对 $\lambda \to -\mathrm{\infty }$ 放大因子趋于 $-1$，快分量会以虚假振荡残留）。对本问题，推荐 BDF2 或 BDF3（配合自适应步长）：它们 L 稳定（快分量被压制而不是振荡）、阶数中等（精度合适）、每步代价是 $n^3$ 级别的线性代数（$n=50$ 时是 $1.25\times {10}^5$ 次运算，可以接受）。若需要更高精度或更强稳定性，改用 Radau IIA 的 5 阶版本。

步长策略。用自适应步长：步长的选择由局部误差估计（嵌入方法或步长加倍比较）与稳定性两方面决定。L 稳定方法没有稳定性限制，因此步长只受精度与 Newton 迭代收敛性的约束，初始阶段（快分量尚未衰减）需要小步长，之后可以迅速放大到与慢尺度匹配的步长（$\frac{1}{{10}^{-3}}={10}^3$ 量级的上界，实际取到局部误差允许的范围）。实现细节：Newton 迭代用解析或数值 Jacobian（$n=50$ 时解析 Jacobian 的代价约 $n^2$，可以接受），Jacobian 在若干个步内复用（每次 Newton 迭代的收敛速度随复用而变差，需要平衡）。

精度与代价的估计。BDF2 的局部截断误差是 $O(h^3)$、整体误差 $O(h^2)$；要达到相对误差 ${10}^{-6}$，步长约 ${10}^{-3}$（由误差公式外推），总步数约 ${10}^5$，每步的代价是 Newton 迭代（$2$ 到 $3$ 次迭代，每次求解 $n\times n$ 线性方程组，$O(n^3)=1.25\times {10}^5$ 次运算）加上一次 Jacobian 求值与分解，总代价约 ${10}^5\times 4\times 1.25\times {10}^5=5\times {10}^{10}$ 次浮点运算——在现代机器上约分钟级。与显式方法的 $2\times {10}^8$ 次函数求值（每次 $50$ 次运算，即 ${10}^{10}$ 次运算）相比，隐式方法的代价更高，但隐式方法能覆盖整个区间 $[0,100]$，而显式方法在同样的计算量下只能覆盖到 $t\approx 0.05$（约 $2500$ 步 $\times 2\times {10}^{-5}$）。这条对比是刚性问题的核心判断依据：当需要覆盖比刚性时间尺度长几个数量级的时间区间时，隐式方法的优势是决定性的；若只关心快尺度附近的短时间行为，显式方法更简单也更快。

验证方案：用两种不同阶数的方法（如 BDF2 与 BDF3）或两种不同步长序列计算，检查精度阶是否与理论相符（步长减半误差降到四分之一）；用守恒量或已知的稳态（若模型有稳态解，检查长期极限是否一致）做物理合理性检验；检查 Newton 迭代的收敛性（残差是否按预期下降）与步长序列的合理性（步长变化不应剧烈振荡，否则说明误差估计器有问题）。报告中给出刚性比、所用格式、稳定性质（A/L 稳定）、步长策略与总计算代价，这些信息合起来才能让结果可复核。

常见错误

错误 1 · 把稳定性与精度混为一谈

原因：稳定性（误差不放大）与精度（误差的大小）是两个独立的问题。显式 Euler 在 $h=0.002$ 处（稳定条件的临界点）误差是 $2.000\times {10}^{-6}$，而隐式 Euler 在 $h=0.01$ 处误差是 $4.996\times {10}^{-6}$——精度量级相当，但步长相差 $5$ 倍。反过来，一个稳定的格式可能精度很差（隐式 Euler 只有一阶），一个高精度的格式可能在刚性问题上不稳定（RK4 稳定域有限）。本节实验二的数据把这两条性质放在同一张表里，容易读成「隐式方法更精确」——准确的读法是「隐式方法在刚性问题上可以用更大的步长达到同等精度」。

解决：分别检查两条性质。稳定性用稳定域分析（把 $h\lambda$ 代入放大因子，检查模是否小于 $1$），精度用局部截断误差的阶数（与步长的幂次）。选择方法时先确定问题的刚性，再在「允许的步长」与「需要的精度」之间做取舍。报告中分开给出稳定性的依据与精度的依据，不要用「误差小」代替「稳定」。

错误 2 · 用线性化的结论判断非双曲平衡点

原因：线性化定理只对双曲平衡点（Jacobian 特征值实部全非零）成立。非双曲情形（存在实部为零的特征值）下，高阶项决定稳定性，线性化给出的结论可能完全错误。典型例子：$y^{\prime }=-y^3$ 的线性化给出 $\lambda =0$（非双曲），但实际解按 $t^{-1/2}$ 衰减（渐近稳定，但衰减是多项式速率而不是指数速率）；$y^{\prime }=y^3$ 的解在有限时间爆破（不稳定）；$y^{\prime }=ϵy-y^3$ 的稳定性随 $ϵ$ 的符号变化（分岔）。中心型平衡点（如本节的谐振子）是典型的非双曲情形：线性化给出闭轨道，但任意小的扰动可能变成稳定或不稳定的焦点。

解决：先检查特征值的实部是否全非零（计算 $\underset{i}{min}|\mathrm{Re}{\lambda }_i|$ 并与矩阵的尺度比较）。非双曲时用 Lyapunov 函数、中心流形约化或正规型方法做高阶分析，或者用数值方法做参数扫描（观察轨道的长期行为随参数的变化）。报告中说明平衡点的类型（双曲或非双曲）与所用判据，不要对非双曲情形直接引用线性化结论。

错误 3 · 在刚性问题上沿用显式方法的步长

原因：显式方法的步长限制由刚性尺度决定，与解的光滑程度无关。本节实验二的例子中，解是光滑的余弦（变化尺度是 $O(1)$），但步长被限制在 $0.002$（由特征值 $-1000$ 决定）。$h=0.01$ 时误差是 $3.528\times {10}^{185}$——看起来「步长只超出了 5 倍」却导致完全发散，因为误差按 $9^n$ 增长。这类失效在调试时容易被误认为是「模型错误」或「参数错误」，实际上格式与问题的刚性不匹配。

解决：先估计 Jacobian 的特征值范围（解析求导、数值差分或用 Gershgorin 圆盘估计），算出刚性比与显式方法的步长限制，再决定方法。刚性比超过 ${10}^3$ 时改用隐式方法（BDF、Radau、隐式 Euler 起步），或采用隐式-显式（IMEX）分区处理（把快分量隐式、慢分量显式）。实现后用「步长加倍/减半时结果是否稳定」做检验，并检查是否存在步长的硬性上限（若存在，说明问题刚性而格式不合适）。

错误 4 · 用数值解的长时行为代替理论分析

原因：数值格式本身可能改变系统的定性行为，因此长期积分的结果不能直接当作方程的结论。本节实验三给出例子：中心型系统（真实轨道闭合、能量守恒）用显式 Euler 积分后振幅按 $(1+h^2{)}^{N/2}$ 缓慢增长（$T=10$ 时 $+0.5\mathrm{\%}$），定性行为被改变（从中心变成缓慢外扩的螺旋）。同类问题：耗散性格式会把守恒系统「磨平」、辛格式可能保持结构但引入其他误差、显式方法的长期能量漂移会导致混沌系统的统计性质变化。把这类数值效应当作物理或模型的结论会得出错误判断。

解决：长时积分前先确定系统的定性性质（守恒量、平衡点的类型、是否耗散），再选择与之一致的格式（守恒系统用辛格式或保能量格式，耗散系统用 L 稳定格式）。用守恒量或 Lyapunov 函数监控数值解（本节实验四的单调性检验就是这种监控）：守恒量的漂移量应当与格式的阶数和步长相符，出现系统性漂移说明格式不合适。报告中说明格式的性质（是否保结构）与监控量的漂移量，不要只报告终态的数值。