算法设计

用数值方法求解ODE基本上可以分为下述三种方法

• 基于数值积分的方法
• 基于多项式配置或Taylor展开的方法
• 基于有限差分的方法

基于数值积分的方法

求解常微分方程（ODE）初值问题

$$\frac{dy}{dt}=f(t,y),y(t_0)=y_0$$

的一种核心途径是将其写成等价的积分方程：

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

在离散点 $t_n$ 上，递推关系为

$$y(t_{n+1})=y(t_n)+{\int }_{t_n}^{t_{n+1}}f(s,y(s))ds$$

用不同的数值积分公式近似其中的积分，本质上是在构造不同的ODE求解器，就自然导出欧拉法、梯形法、中点法以及亚当斯（Adams）多步法等。

左矩形法 → 欧拉法（显式）

我们用两种视角来看待这种积分器

视角一：Taylor展开

证明：用 Taylor 展开推导显式欧拉法

求解初值问题：

$$\frac{dy}{dx}=f(x,y),y(x_0)=y_0$$

设步长为 $h=x_{k+1}-x_k$，且 $f(x,y)$ 充分光滑。

第一步：Taylor 展开

将精确解 $y(x_{k+1})$ 在 $x_k$ 处作 Taylor 展开：

$$y(x_{k+1})=y(x_k+h)=y(x_k)+h{y}^{\prime }(x_k)+\frac{h^2}{2}{y}^{\prime \prime }(\xi )$$

其中 $\xi \in (x_k,x_{k+1})$（拉格朗日余项）。

第二步：用微分方程替换导数

由原方程 $y^{\prime }(x)=f(x,y(x))$，知 $y^{\prime }(x_k)=f(x_k,y(x_k))$。

引入记号 $y_k\approx y(x_k)$，则上式可写为：

$$y(x_{k+1})=y(x_k)+hf(x_k,y(x_k))+\frac{h^2}{2}{y}^{\prime \prime }(\xi )$$

第三步：截断高阶项，得到数值格式

舍去 $O(h^2)$ 的余项，取近似：

$$\boxed{y_{k+1}=y_k+hf(x_k,y_k)}$$

这就是显式欧拉法。

局部截断误差分析

从 Taylor 展开过程可以直接读出局部截断误差（单步误差）：

$$T_{k+1}=y(x_{k+1})-[y(x_k)+hf(x_k,y(x_k))]=\frac{h^2}{2}{y}^{\prime \prime }(\xi )=O(h^2)$$

所以显式欧拉法的局部截断误差是 $O(h^2)$，而经过 $N\propto 1/h$ 步累积后，全局误差是 $O(h)$，即一阶精度。

视角二：数值积分

从积分方程视角推导显式欧拉法

1. 积分方程是出发点

将微分方程 $\frac{dy}{dt}=f(t,y)$ 在区间 $[t_n,t_{n+1}]$ 上积分，得到等价的积分方程：

$$y(t_{n+1})=y(t_n)+{\int }_{t_n}^{t_{n+1}}f(s,y(s))ds$$

这里，$y(t_n)$ 是已知的当前值，而 $y(t_{n+1})$ 是未知的未来值。要计算未来值，关键就在于怎么处理那个积分项。

2. 左矩形公式：用"旧信息"代替"新积分"

积分项 ${\int }_{t_n}^{t_{n+1}}f(s,y(s))ds$ 的几何意义是曲线 $f(t,y(t))$ 在区间 $[t_n,t_{n+1}]$ 下的面积。

显式欧拉法做了一个极其暴力的近似：它假设在整个区间 $[t_n,t_{n+1}]$ 上，$f$ 的值是不变的，始终等于它在左端点 $t_n$ 处的值 $f(t_n,y_n)$。

这在几何上，就是用一个以左端点函数值为高的矩形面积，去替代曲线下的真实面积。

所以积分近似为：

$${\int }_{t_n}^{t_{n+1}}f(s,y(s))ds\approx (t_{n+1}-t_n)\cdot f(t_n,y_n)=h\cdot f(t_n,y_n)$$

这本质上就是左矩形积分公式（Left Rectangle Rule）。

3. 代回原式，欧拉法诞生

把这个近似代回积分方程：

$$\boxed{y_{n+1}=y_n+h\cdot f(t_n,y_n)}$$

它之所以是 “显式” 的，正是因为这个积分近似只用到了时刻 $t_n$ 的旧信息 $f(t_n,y_n)$，完全没有涉及到我们还不知道的 $y_{n+1}$。

用数值积分的视角更适合建立几何直觉，而用Taylor展开更适合做误差分析

代码实现

我们假定迭代是均匀的

显式欧拉法

#python

import numpy as np
def explicit_euler(f, t0, y0, h, n):
    t = [t0]
    y = [y0]
    for i in range(n):
        y_new = y[-1] + h * f(t[-1], y[-1])
        t_new = t[-1] + h
        t.append(t_new)
        y.append(y_new)
    return t, y
    
#输出:t = [t0, t1, t2, ..., t_n];y = [y0, y1, y2, ..., y_n]

t_arr = np.array(t)
y_arr = np.array(y)

梯形法

几何视角：梯形数值积分

梯形法的几何直觉非常直接：用梯形面积代替曲线下的真实面积。

在区间 $[t_n,t_{n+1}]$ 上，被积函数 $f(t,y(t))$ 在两个端点的值分别为
$f(t_n,y_n)$ 和 $f(t_{n+1},y_{n+1})$。
梯形的面积公式为

$$\text{面积}\approx h\cdot \frac{f(t_n,y_n)+f(t_{n+1},y_{n+1})}{2}.$$

将此近似代入积分递推式，就得到梯形法的格式：

$$\boxed{y_{n+1}=y_n+\frac{h}{2}[f(t_n,y_n)+f(t_{n+1},y_{n+1})]}.$$

几何直觉

• 梯形比只用一个端点值的矩形（显式欧拉法）更贴近真实曲线，因此精度更高
• 梯形的右端点高依赖于未知的 $y_{n+1}$，这意味着公式两端都含有 $y_{n+1}$，因此梯形法是隐式的

代数视角：Taylor 展开与误差分析

代数推导能严格求出局部截断误差，并揭示低阶项完全抵消的数学结构

假设 $f(t,y)$ 充分光滑，记精确解为 $y(t)$，且 $y_n=y(t_n)$

第一步：展开精确解

将 $y(t_{n+1})$ 在 $t_n$ 处 Taylor 展开：

$$\begin{array}{rl}y(t_{n+1})& =y(t_n+h)\\ & =y(t_n)+h{y}^{\prime }(t_n)+\frac{h^2}{2}{y}^{\prime \prime }(t_n)+\frac{h^3}{6}{y}^{\prime \prime \prime }({\xi }_1),{\xi }_1\in (t_n,t_{n+1}).\end{array}$$

第二步：展开导数项 $f(t_{n+1},y_{n+1})$

注意到 $y^{\prime }(t)=f(t,y(t))$，将 $y^{\prime }(t_{n+1})$ 在 $t_n$ 处展开：

$$\begin{array}{rl}{y}^{\prime }(t_{n+1})& =f(t_{n+1},y(t_{n+1}))\\ & =y^{\prime }(t_n)+h{y}^{\prime \prime }(t_n)+\frac{h^2}{2}{y}^{\prime \prime \prime }({\xi }_2),{\xi }_2\in (t_n,t_{n+1}).\end{array}$$

第三步：构造梯形和

计算数值格式等号右侧的表达式：

$$\begin{array}{rl}& y(t_n)+\frac{h}{2}[y^{\prime }(t_n)+y^{\prime }(t_{n+1})]\\ =& y(t_n)+\frac{h}{2}[y^{\prime }(t_n)+(y^{\prime }(t_n)+h{y}^{\prime \prime }(t_n)+\frac{h^2}{2}{y}^{\prime \prime \prime }({\xi }_2)\left)\right]\\ =& y(t_n)+h{y}^{\prime }(t_n)+\frac{h^2}{2}{y}^{\prime \prime }(t_n)+\frac{h^3}{4}{y}^{\prime \prime \prime }({\xi }_2).\end{array}$$

第四步：计算局部截断误差

局部截断误差定义为真实解与数值解（从精确数据出发）的差：
$$T_{n+1}=y(t_{n+1})-\left[y(t_n)+\frac{h}{2}(f(t_n,y(t_n))+f(t_{n+1},y(t_{n+1})))\right].$$
将两个展开式代入，相减：
$$\begin{array}{rl}{T}_{n+1}& =\left[y(t_n)+h{y}^{\prime }(t_n)+\frac{h^2}{2}{y}^{\prime \prime }(t_n)+\frac{h^3}{6}{y}^{\prime \prime \prime }({\xi }_1)\right]\\ & -\left[y(t_n)+h{y}^{\prime }(t_n)+\frac{h^2}{2}{y}^{\prime \prime }(t_n)+\frac{h^3}{4}{y}^{\prime \prime \prime }({\xi }_2)\right]\\ & =\frac{h^3}{6}{y}^{\prime \prime \prime }({\xi }_1)-\frac{h^3}{4}{y}^{\prime \prime \prime }({\xi }_2)\\ & =O(h^3).\end{array}$$
解释

• 常数项、$h$ 项和 $h^2$ 项全部抵消，只剩下 $h^3$ 项
• 因此梯形法的局部截断误差为 $O(h^3)$，比显式欧拉（$O(h^2)$）高一阶
• 经过 $N\propto 1/h$ 步累积后，全局误差为 $O(h^2)$，即方法具有二阶精度

代数视角的优势
精确量化误差，揭示对称性（梯形法取两端点平均）带来的低阶抵消效应。

总结

设步长 $h=t_{n+1}-t_n$，近似解记为 $y_n$

数值积分规则	公式	方法	局部截断误差
左矩形	$${\int }_{t_n}^{t_{n+1}}f(s,y(s))ds\approx hf(t_n,y_n)$$	显式欧拉	$O(h^2)$
右矩形	$${\int }_{t_n}^{t_{n+1}}f(s,y(s))ds\approx hf(t_{n+1},y_{n+1})$$	隐式欧拉	$O(h^2)$
梯形	$${\int }_{t_n}^{t_{n+1}}f(s,y(s))ds\approx \frac{h}{2}\left[f(t_n,y_n)+f(t_{n+1},y_{n+1})\right]$$	梯形法（隐式）	$O(h^3)$
中点	$${\int }_{t_n}^{t_{n+1}}f(s,y(s))ds\approx hf\left(t_n+\frac{h}{2},y_{n+1/2}\right)\text{需预测}$$	显式中点法	$O(h^3)$

补充：数值求解的形式

数值求解常微分方程（ODE）的最终目的，就是得到两组一一对应的离散数据：

• 一组时间点 $t_0,t_1,t_2,\dots$
• 一组对应的数值解 $y_0,y_1,y_2,\dots$