====== 第十五章 数值解法 ======

===== 15.1 数值解的基本概念 =====

**定义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 Euler方法 =====

**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}$$

二阶精度。

===== 15.3 Runge-Kutta方法 =====

**一般形式**：
$$\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)$，四阶精度。

===== 15.4 线性多步法 =====

**一般形式**：
$$\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}$

  - 二阶：$y_{n+2} = y_{n+1} + \frac{h}{2}(3f_{n+1} - f_n)$
  - 四阶：$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}$

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

**15.4.3 预测-校正方法**

用显式方法预测，隐式方法校正：
  - PECE模式：预测(P) → 估值(E) → 校正(C) → 估值(E)
  - PEC模式：预测 → 估值 → 校正

===== 15.5 稳定性与刚性问题 =====

**定义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 各方法的稳定区域**

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

**15.5.2 刚性问题**

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

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

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

===== 15.6 方程组与高阶方程 =====

**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 边值问题的数值解 =====

**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.8 例题详解 =====

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

*解*：
  - $k_1 = f(0, 0.5) = 0.5 - 0 + 1 = 1.5$
  - $k_2 = f(0.1, 0.5 + 0.1 \times 1.5) = f(0.1, 0.65) = 0.65 - 0.01 + 1 = 1.64$
  - $k_3 = f(0.1, 0.5 + 0.1 \times 1.64) = f(0.1, 0.664) = 0.664 - 0.01 + 1 = 1.654$
  - $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.9 习题 =====

**习题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}$，分析为什么隐式方法更适合。