定义15.1.1（初值问题）

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

数值解：在离散点 $x_0 < x_1 < \cdots < x_N$ 上求近似值 $y_n \approx y(x_n)$。

步长：$h_n = x_{n+1} - x_n$，常取等步长 $h$。

定义15.1.2（局部截断误差）

单步法 $y_{n+1} = y_n + h\phi(x_n, y_n, h)$ 的局部截断误差： $$T_{n+1} = y(x_{n+1}) - y(x_n) - h\phi(x_n, y(x_n), h)$$

若 $T_{n+1} = O(h^{p+1})$，称方法具有 $p$ 阶精度。

15.2.1 显式Euler方法

$$y_{n+1} = y_n + hf(x_n, y_n)$$

几何意义：用切线近似曲线。

局部截断误差：$T_{n+1} = \frac{h^2}{2}y''(\xi) = O(h^2)$，一阶精度。

整体误差：$O(h)$

15.2.2 隐式Euler方法（后退Euler）

$$y_{n+1} = y_n + hf(x_{n+1}, y_{n+1})$$

需解方程求 $y_{n+1}$，但稳定性更好。

15.2.3 梯形公式

$$y_{n+1} = y_n + \frac{h}{2}[f(x_n, y_n) + f(x_{n+1}, y_{n+1})]$$

局部截断误差 $O(h^3)$，二阶精度。

15.2.4 改进Euler方法（Heun方法）

预测：$\bar{y}_{n+1} = y_n + hf(x_n, y_n)$

校正：$y_{n+1} = y_n + \frac{h}{2}[f(x_n, y_n) + f(x_{n+1}, \bar{y}_{n+1})]$

等价形式： $$\begin{cases}k_1 = f(x_n, y_n) \\ k_2 = f(x_n + h, y_n + hk_1) \\ y_{n+1} = y_n + \frac{h}{2}(k_1 + k_2)\end{cases}$$

二阶精度。

一般形式： $$\begin{cases}k_1 = f(x_n, y_n) \\ k_2 = f(x_n + c_2h, y_n + ha_{21}k_1) \\ \vdots \\ k_s = f(x_n + c_sh, y_n + h\sum_{j=1}^{s-1}a_{sj}k_j) \\ y_{n+1} = y_n + h\sum_{i=1}^s b_ik_i\end{cases}$$

15.3.1 经典四阶Runge-Kutta（RK4）

$$\begin{cases}k_1 = f(x_n, y_n) \\ k_2 = f(x_n + \frac{h}{2}, y_n + \frac{h}{2}k_1) \\ k_3 = f(x_n + \frac{h}{2}, y_n + \frac{h}{2}k_2) \\ k_4 = f(x_n + h, y_n + hk_3) \\ y_{n+1} = y_n + \frac{h}{6}(k_1 + 2k_2 + 2k_3 + k_4)\end{cases}$$

局部截断误差 $O(h^5)$，四阶精度。

一般形式： $$\sum_{j=0}^k \alpha_j y_{n+j} = h\sum_{j=0}^k \beta_j f_{n+j}$$

其中 $\alpha_k = 1$，$|\alpha_0| + |\beta_0| \neq 0$。

15.4.1 Adams-Bashforth方法（显式）

基于外推：$y_{n+k} = y_{n+k-1} + h\sum_{j=0}^{k-1}\beta_j f_{n+j}$

1. 二阶：$y_{n+2} = y_{n+1} + \frac{h}{2}(3f_{n+1} - f_n)$
2. 四阶：$y_{n+4} = y_{n+3} + \frac{h}{24}(55f_{n+3} - 59f_{n+2} + 37f_{n+1} - 9f_n)$

15.4.2 Adams-Moulton方法（隐式）

基于内插：$y_{n+k} = y_{n+k-1} + h\sum_{j=0}^k \beta_j^* f_{n+j}$

1. 二阶（梯形）：$y_{n+1} = y_n + \frac{h}{2}(f_n + f_{n+1})$
2. 四阶：$y_{n+3} = y_{n+2} + \frac{h}{24}(9f_{n+3} + 19f_{n+2} - 5f_{n+1} + f_n)$

15.4.3 预测-校正方法

用显式方法预测，隐式方法校正：

1. PECE模式：预测(P) → 估值(E) → 校正(C) → 估值(E)
2. PEC模式：预测 → 估值 → 校正

定义15.5.1（绝对稳定性）

用数值方法求解试验方程 $\frac{dy}{dx} = \lambda y$（$\text{Re}(\lambda) < 0$），若对某步长 $h$ 有 $y_n \to 0$（$n \to \infty$），称方法对该 $h$ 绝对稳定。

稳定区域：复平面上使方法绝对稳定的 $z = \lambda h$ 的集合。

15.5.1 各方法的稳定区域

1. 显式Euler：单位圆 $|1+z| < 1$（圆盘）
2. 隐式Euler：$|z-1| > 1$（圆外，包含整个左半平面）
3. 梯形公式：整个左半平面（A-稳定）
4. RK4：较复杂区域，大致为左半平面的一部分

15.5.2 刚性问题

定义：若系统同时存在快变和慢变模式（特征值大小相差悬殊），称为刚性系统。

刚性比：$S = \frac{\max|\text{Re}\lambda_i|}{\min|\text{Re}\lambda_i|}$

当 $S \gg 1$ 时，显式方法需要极小步长才能稳定，隐式方法更合适。

15.6.1 方程组的数值解

$$\mathbf{y}' = \mathbf{f}(x, \mathbf{y}), \quad \mathbf{y} = (y_1, \ldots, y_m)^T$$

RK4推广：各分量分别计算。

15.6.2 高阶方程降阶

$$y^{(m)} = f(x, y, y', \ldots, y^{(m-1)})$$

令 $z_1 = y, z_2 = y', \ldots, z_m = y^{(m-1)}$，化为方程组： $$\begin{cases}z_1' = z_2 \\ z_2' = z_3 \\ \vdots \\ z_m' = f(x, z_1, \ldots, z_m)\end{cases}$$

15.7.1 打靶法（Shooting Method）

将边值问题 $y'' = f(x,y,y')$，$y(a) = \alpha$，$y(b) = \beta$ 转化为初值问题。

猜测 $y'(a) = s$，数值求解至 $x = b$，调整 $s$ 使得 $y(b;s) = \beta$。

用Newton迭代：$s_{n+1} = s_n - \frac{y(b;s_n)-\beta}{\partial y/\partial s|_{s_n}}$

15.7.2 差分法

在内部点 $x_i$ 上用差分近似导数： $$y''(x_i) \approx \frac{y_{i+1} - 2y_i + y_{i-1}}{h^2}$$

得到线性方程组（三对角），可用追赶法求解。

15.7.3 有限元方法

变分原理 + 分片多项式逼近，详见偏微分方程数值解。

例15.1：用RK4求解 $y' = y - x^2 + 1$，$y(0) = 0.5$，取 $h = 0.2$，求 $y(0.2)$。

*解*：

1. $k_1 = f(0, 0.5) = 0.5 - 0 + 1 = 1.5$
2. $k_2 = f(0.1, 0.5 + 0.1 \times 1.5) = f(0.1, 0.65) = 0.65 - 0.01 + 1 = 1.64$
3. $k_3 = f(0.1, 0.5 + 0.1 \times 1.64) = f(0.1, 0.664) = 0.664 - 0.01 + 1 = 1.654$
4. $k_4 = f(0.2, 0.5 + 0.2 \times 1.654) = f(0.2, 0.8308) = 0.8308 - 0.04 + 1 = 1.7908$

$$y_1 = 0.5 + \frac{0.2}{6}(1.5 + 2 \times 1.64 + 2 \times 1.654 + 1.7908) = 0.8293$$

精确解 $y = (x+1)^2 - 0.5e^x$，$y(0.2) = 0.8292986...$，误差极小。

习题15.1：用显式Euler方法（$h = 0.1$）求解 $y' = -y$，$y(0) = 1$，计算 $y(0.5)$ 并与精确解比较。

习题15.2：推导改进Euler方法的局部截断误差，证明其为二阶方法。

习题15.3：用RK4求解初值问题 $y' = x + y$，$y(0) = 1$，取 $h = 0.1$，求 $y(0.4)$。

习题15.4：分析中点法 $y_{n+1} = y_n + hf(x_n + \frac{h}{2}, y_n + \frac{h}{2}f_n)$ 的稳定区域。

习题15.5：用打靶法求解边值问题 $y'' = -y$，$y(0) = 0$，$y(\frac{\pi}{2}) = 1$。

习题15.6：证明梯形公式是A-稳定的。

习题15.7：对于刚性系统 $\begin{cases}y_1' = -1000y_1 + y_2 \\ y_2' = y_1 - y_2\end{cases}$，分析为什么隐式方法更适合。