4. 数值求解基础与误差来源

动力系统初值问题通常写作 y' = f(t, y)，y(t₀) = y₀。除少数特殊形式外，绝大多数非线性方程没有解析解，必须借助数值方法逐步推进。数值求解的本质，是把连续时间区间上的微分方程转化为有限个离散节点上的代数递推关系——这就是离散化思想。理解离散化带来的两类误差，以及方法本身是否可靠（稳定、收敛），是正确使用数值解、判断结果可信度的前提。

一、离散化：从连续到离散的转化

以最基础的前向欧拉法为例。将区间 [t₀, T] 划分为步长 h 的网格 tₙ = t₀ + nh，用一阶差商近似导数：

y'(tₙ) ≈ (y(tₙ₊₁) − y(tₙ)) / h
⇒ yₙ₊₁ = yₙ + h·f(tₙ, yₙ)

这里 yₙ 是数值近似，y(tₙ) 是精确解。离散化把"求一条连续曲线"变成"从初值出发逐点递推"，代价是引入了近似。不同方法对应不同的离散格式：后向欧拉用 f(tₙ₊₁, yₙ₊₁)（隐式，需解方程），梯形法取两端平均，四阶龙格–库塔法在一步内取四个斜率加权组合。步长 h 越小，离散点越密，但计算量越大，且舍入误差会累积——这正是误差分析要回答的问题。

二、截断误差：方法本身的"先天误差"

截断误差来源于用差分代替导数时丢弃的高阶项。把精确解在 tₙ 处做泰勒展开：

y(tₙ₊₁) = y(tₙ) + h·y'(tₙ) + (h²/2)·y''(tₙ) + O(h³)

前向欧拉只保留了到 h·y' 项，因此局部截断误差为 (h²/2)·y''(ξ)，量级 O(h²)；经过 N = (T−t₀)/h 步累积，全局截断误差为 O(h)，故称一阶方法。四阶龙格–库塔的局部截断误差为 O(h⁵)，全局为 O(h⁴)。这解释了为什么提高阶数能在相同步长下获得更高精度，但高阶方法每步计算量也更大，需权衡。

三、舍入误差：计算机有限字长的"后天误差"

计算机用有限位浮点数表示实数，每次运算都会产生舍入。单步舍入误差约为机器精度 ε ≈ 2.2×10⁻¹⁶（双精度），看似微不足道，但它在递推中会被不断放大或累积。总误差大致由两部分竞争：截断误差随 h 减小而下降，舍入误差累积却随步数 N = (T−t₀)/h 增大而上升，量级约 ε/h。二者叠加存在一个"最优步长"：步长过大会被截断误差主导，步长过小则舍入误差抬头，总误差反而增大。实际计算中应避免盲目追求极小步长。

四、稳定性与收敛性：方法可靠性的两个维度

收敛性回答"步长趋于零时，数值解是否趋于精确解"。一个相容（截断误差随 h→0 趋于零）且稳定的方法必然收敛（Lax 等价定理）。稳定性回答"计算过程中微小扰动（初值扰动、舍入扰动）是否被控制住而不被无限放大"。

稳定性常用模型方程 y' = λy（λ 为复数）分析。前向欧拉的递推放大因子为 R(z) = 1 + z，z = λh；要求 |R(z)| ≤ 1，得到稳定域 |1 + λh| ≤ 1。对实负特征值 λ < 0，需 h ≤ 2/|λ|。若系统含快速衰减模态（|λ| 很大），前向欧拉被迫取极小步长，否则数值解会振荡发散——这就是刚性问题的典型表现。后向欧拉的放大因子为 1/(1−z)，对 Re(z) < 0 恒有 |R| ≤ 1，是 A-稳定的，适合刚性问题，但每步需解隐式方程。

五、实践中的判断与操作

面对一个具体初值问题，可按以下思路组织数值求解与验证。以 y' = −50y，y(0) = 1，求解到 T = 0.5 为例（λ = −50，精确解 y = e^{−50t}，属刚性倾向问题）：

1. 选方法：若用前向欧拉，稳定要求 h ≤ 2/50 = 0.04。取 h = 0.05 会越界，数值解出现符号交替放大；取 h = 0.01 则稳定但需 50 步。
2. 估误差：前向欧拉全局误差 O(h)，h = 0.01 时误差量级约 10⁻²，对精度要求高时不够。
3. 换方法：改用后向欧拉或四阶龙格–库塔。后向欧拉无条件稳定，可用 h = 0.05 甚至更大；四阶 RK 在 h = 0.01 时全局误差约 O(h⁴) ≈ 10⁻⁸，但稳定域有限，仍需满足 |1 + λh + (λh)²/2 + …| ≤ 1 的约束。
4. 验证：将步长减半重算，比较两次结果之差；若差值随 h 按方法阶数下降，说明截断误差主导、结果可信；若差值不降反升，说明已进入舍入误差主导区或方法不稳定。

总结来看，离散化是数值求解的起点，截断误差刻画方法的"先天精度"，舍入误差反映计算的"后天限制"，二者共同决定总误差；稳定性保证扰动不被放大，收敛性保证方法在极限意义下正确。只有同时把握这四个概念，才能在给定问题上选对方法、取对步长，并对数值结果给出有依据的可信度判断。