Runge - kutta 法

对于初值问题

中点公式推导

首先有 $\displaystyle u(t_1)\approx u(t_0)+hf(t_0+\frac{h}{2},u(t_0+\frac{h}{2}))$ ，这里局部截断误差为 $O(h^3)$ 。但 $u(t_0+\frac{h}{2})$ 未知，用 Euler 法近似：

代入可得

引入记号：

则数值解

数值解与精确解展开及误差分析

数值解展开：

精确解展开：

所以 $|u(t_1)-u_1| = O(h^3)$ ，即该 Runge - kutta 方法是二阶的

一般 m 点 Runge - kutta 法思路

已知 $\displaystyle u(t + h)=u(t)+\int_{t}^{t + h}f(\tau,u(\tau))d\tau$ 。在区间 $[t,t + h]$ 上取 $m$ 个点 $t = t_1\leq t_2\leq\cdots\leq t_m = t + h$

类似于多步法，有近似式 $\displaystyle u(t + h)\approx u(t)+h\sum_{i = 1}^{m}c_if(t_i,u(t_i))$ ，其中 $\displaystyle k_i = f(t_i,u(t_i))$ 。但 $u(t_i)$ 未知，需进行近似

利用 Euler 法：

故

即通过前一步的近似值和当前点的斜率 $k_i$ 逐步递推后续点的函数值近似

用类似多步法的思想

其中 $\displaystyle \rho_{i1}+\cdots+\rho_{i,i - 1}=t_i - t_1$

引入记号简化

引入 $\displaystyle t_i = t + a_ih$ ，$i = 2,3,\cdots,m$ ，$t_1 = t$ ，以及 $\displaystyle b_{ij}=\frac{\rho_{ij}}{h}$ ，$i = 2,\cdots,m$ ，$j = 1,\cdots,i - 1$

Runge - kutta 法一般形式（m 级）

格式构建

从 $f(t_1,u_1)=k_1$ 开始，依次通过

最终数值解

其中满足

即

待定系数法（以三级三阶为例）

精确解展开：

记

Runge - kutta 格式展开：

其中 $b_{21}=a_2$

对 $\displaystyle h(c_1k_1 + c_2k_2 + c_3k_3)$ 展开并与精确解 $u(t + h)$ 的展开式对比 $h$ 同次幂系数：

令

展开后得到方程组

这里有 6 个未知数，4 个方程，存在 2 个自由参数，可确定三级三阶方法

经典四阶 Runge - kutta 法格式

注：Runge - kutta 法是单步法，从 $u_n$ 计算 $u_{n + 1}$

与多步法比较及单步法特性

线性单步法：关于 $f$ 线性时，最高为二阶（$O(h^2)$ ）

非线性单步法（如 Runge - kutta 法 ）：因为 $k_1 + k_2 + 2k_3 + k_4$ 不是 $f(t_n,u_n)$ 的线性函数，所以能有更高阶数

Taylor 展开法

其中

舍去余项后，得到 $u_1 = u_0 + h\varphi(t_0, u_0, h)$ ，一般形式为

例如四阶 Runge - kutta 法：

相容性：$\displaystyle \frac{u_{n + 1}-u_n}{h}=\varphi(t_n, u_n, h)$ 逼近微分方程，即

绝对稳定性和绝对稳定域

回顾 Euler 法稳定性：已知 $\displaystyle |e_n| \leq e^{LT}|e_0|$ 。以初值问题 为例，$f(t, u)=-5u$ ，$L = 5$ ，$T = 20$ ，$e^{LT}=e^{100}$ ，精确解 $u(t)=e^{-5t}$ ，$t \to +\infty$ 时，$u(t) \to 0$

当 $h = 0.5$ 时，$u_n=(-1.5)^n \to \pm\infty$ ，与真解不同；当 $h = 0.4$ 时，$u_n=(-1)^n$ ，振荡；当 $h = 0.39$ 时，$u_n=(1 - 5h)^n \to 0$

说明稳定性概念需要改进，引入绝对稳定性

Euler 法绝对稳定性定义：

将 Euler 法用于形如 $u’=\lambda u$ （$\lambda$ 为实数或复数 ）的方程，若 $|1 + h\lambda| < 1$ ，则称 Euler 法关于 $\overline{h}=h\lambda$ 绝对稳定 。$\overline{D}_A={\overline{h}:|1 + \overline{h}| < 1}$ 称为 Euler 法的绝对稳定域

其绝对稳定域 $\overline{D}_A$ 是以 - 1 为圆心，半径为 1 的圆，即满足 $|1 + h\lambda| < 1$ 时，误差不会无限增长

记 $e_n=u_n-v_n$ ，则

k 步法的绝对稳定性

对于测试方程 $u’=\lambda u$ ，线性 $k$ 步法：

令 $\overline{h}=h\lambda$ ，则

在实际计算中，每步存在舍入误差，设近似解 $\overline{u}_n$ 满足

其中 $\eta_n$ 为舍入误差，且 $|\eta_n| \leq M$ 。记 $\overline{e}_n=\overline{u}_n - u_n$ ，则

记

矩阵

可将误差方程转化为矩阵形式 $E_{n + 1}=\overline{C}E_n+\overline{\eta}_n$ ，后续可通过分析矩阵 $\overline{C}$ 的特征值等性质，确定线性 $k$ 步法的绝对稳定域 ，即确定 $\overline{h}=h\lambda$ 的取值范围，使得在存在舍入误差情况下，误差不会无限增长

将误差方程

写成向量形式

进一步得到

为使误差 $E_n$ 在计算中逐步减小，要求

$\overline{C}$ 的特征方程为

$\rho(\lambda)$ 和 $\sigma(\lambda)$ 分别为 $k$ 步法的第一、第二多项式

根据引理，$\lim_{n \to \infty}\overline{C}^n = 0$ 等价于特征多项式 $\rho(\lambda)-\overline{h}\sigma(\lambda)$ 的根 $|\lambda| < 1$

绝对稳定性及绝对稳定域定义

绝对稳定性定义：如果 $k$ 步法的特征多项式方程 $\rho(\lambda)-\overline{h}\sigma(\lambda)$ 的根 $|\lambda| < 1$ ，则称该 $k$ 步法关于 $\overline{h}=h\lambda$ 绝对稳定

绝对稳定域定义：若存在复数域 $D_A$ ，使 $k$ 步法对 $\forall\overline{h} \in D_A$ 都绝对稳定，则称 $D_A$ 为该 $k$ 步法的绝对稳定域。此时 $\displaystyle \overline{h}=\frac{\rho(\lambda)}{\sigma(\lambda)}$ ，$D_A={\overline{h}:|\lambda(\overline{h})| < 1}$ 。绝对稳定域越大，方法在不同 $\lambda$ 和 $h$ 取值下适应能力越强，越不容易出现误差发散情况

命题 5.1

对于一个相容的稳定多步法，若绝对稳定，则 $\lambda$ 的实部 $Re(\lambda) < 0$ ，即绝对稳定域 $D_A$ 在左半复平面。这是因为对于稳定的数值方法，当应用于稳定的微分方程（一般对应 $Re(\lambda) < 0$ ，此时精确解随时间衰减 ）时，数值方法的绝对稳定域应与之匹配，保证数值解的稳定性

Runge - kutta 法绝对稳定域推导

m 级 m 阶 Runge - kutta 法应用于测试方程

对于测试方程 $u’=\lambda u$ ，其精确解为 $u = e^{\lambda t}$

计算中间斜率值：

这里 $P_i(\lambda h)$ 是关于 $\lambda h$ 的多项式

计算数值解：

其中 $P_m(\lambda h)$ 是关于 $\lambda h$ 的多项式

与精确解 Taylor 展开对比

将精确解 $u(t_{n + 1})$ 在 $t_n$ 处 Taylor 展开：

由于 $m$ 级 $m$ 阶 Runge - kutta 法的精度要求，使得

由 $u_{n+1}=u_n(1 + P_m(\lambda h))$ 其特征方程为 $\lambda - [1 + P_m(\lambda h)] = 0$ ，其中

通过求解该不等式确定 的取值范围，从而得到方法的绝对稳定域，判断在哪些 $\lambda$ 和 $h$ 组合下数值解稳定

一阶方程组和刚性问题

常微分方程组（ODEs）相关

常微分方程组：一般形式为

初值条件为 $u_i(t_0)=u_{i0}$ ，$i = 1,\cdots,m$

高阶常微分初值问题：如

令 $u_1 = u$ ，$u_2 = u’$ ，$\cdots$ ，$u_m = u^{(m - 1)}$ ，可转化为常微分方程组

写成向量形式为 $\vec{u}’=\vec{f}(t,\vec{u})$ ，$\vec{u}(t_0)=\vec{u}_0$ 。并且单个方程的数值解法（如 Runge - kutta 法等 ）可推广应用到方程组的求解中

刚性问题（Stiff ODE）

示例分析：对于方程组 $\vec{u}’ = A\vec{u}$ ，$\vec{u} \in \mathbb{R}^3$ ，其中

$A$ 的特征值为 $\lambda_1=-0.1$ ，$\lambda_2=-50$ ，$\lambda_3=-3\times10^4$ 。当 $t \to \infty$ 时，$u_i(t) \to 0$ ，但收敛速度差异极大。使用 Runge - kutta 法时，为保证绝对稳定，要求 $h \leq 10^{-4}$ ；为使解充分接近定态解 0 ，对于收敛最慢的解（对应 $\lambda_1$ ），需 $e^{-0.1T} \approx 0$ ，即 $T \approx 40$ ，则计算步数 $\displaystyle n=\frac{T}{h} \approx 4\times10^5$ ，计算量巨大

定义：设 $\lambda_i$ 为 $A$ 的特征值，若满足：

则称方程组是刚性的，$R$ 为刚性比

对于刚性方程组，需专门的数值方法求解，以避免因步长限制导致计算量过大等问题

外推法

思想：

以 Euler 法 $\displaystyle u_{n + 1}=u_n+hf(t_n,u_n)$ 为例

外推法的思想是通过不同步长 $h$ 和 $\displaystyle \frac{h}{2}$ 分别计算数值解序列 $u_1,u_2,\cdots,u_n$ 和 $\displaystyle \overline{u}_{\frac{1}{2}},\overline{u}_1,\overline{u}_{\frac{3}{2}},\cdots,\overline{u}_{n - \frac{1}{2}},\overline{u}_n$ ，一般情况下 $\overline{u}_1 \neq u_1$

通过 $u_1$ 和 $\overline{u}_1$ 的组合 $\displaystyle 2\overline{u}_1 - u_1=u_0+hf(t_{\frac{1}{2}},\overline{u}_{\frac{1}{2}})$ ，利用 Runge - kutta 法思想，其精度可达到 $u(t_1)+O(h^3)$ ，相比原数值方法精度有所提高

一般外推法误差分析及线性组合推导

对于某数值方法，记步长为 $h$ 时的数值解序列 $u_1,u_2,\cdots,u_n$ 为 $u(t,h)$ ，步长为 $\frac{h}{2}$ 时的数值解序列 $\displaystyle \overline{u}_{\frac{1}{2}},\overline{u}_1,\overline{u}_{\frac{3}{2}},\cdots,\overline{u}_{n - \frac{1}{2}},\overline{u}_n$ 为 $u(t,\frac{h}{2})$

假设该方法的误差有展开式

则

对 $u(t,h)$ 和 $u(t,\frac{h}{2})$ 进行线性组合：

设组合式为

将误差展开式代入：

分子为

所以

通过这种线性组合，提高了数值解的精度阶数，从 $O(h^p)$ 提升到 $O(h^{p + 1})$ ，这就是外推法提高数值计算精度的基本原理