Runge - kutta 法

对于初值问题
$$
\begin{cases}u’ = f(t, u)\u(t_0)=u_0\end{cases}
$$

中点公式推导

首先有 $\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 法近似：
$$
\displaystyle u(t_0+\frac{h}{2})\approx u(t_0)+\frac{h}{2}f(t_0,u_0)
$$
代入可得
$$
u(t_1)\approx u(t_0)+hf(t_0+\frac{h}{2},u(t_0)+\frac{h}{2}f(t_0,u_0))
$$
引入记号：
$$
K_1 = f(t_0,u_0)\
K_2 = f(t_0+\frac{h}{2},u_0+\frac{h}{2}K_1)
$$

则数值解
$$
u_1 = u_0+hK_2
$$

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

数值解展开：
$$
u_1 = u_0+hf(t_0,u_0)+\frac{h^2}{2}(f_t + f_uf)(t_0,u_0)+\frac{h^3}{8}(f_{tt}+2f_{tu}f+\cdots)(t_0,u_0)+\cdots
$$
精确解展开：
$$
\begin{align}
u(t_1)&=u(t_0)+hu’(t_0)+\frac{h^2}{2!}u^{(2)}(t_0)+\frac{h^3}{3!}u^{(3)}(t_0)+\cdots\&=u(t_0)+hf(t_0,u_0)+\frac{h^2}{2}(f_t + f_uf)(t_0,u_0)+\frac{h^3}{6}(f_{tt} + 2f_{tu}f +f_{uu}f^2 +f_uf_t+f_u^2f)(t_0,u_0)+\cdots
\end{align}
$$

所以 $|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 法：
$$
\begin{align}
u(t_2)&=u(t_1)+(t_2 - t_1)k_1 \
k_1 &= f(t_1,u(t_1))\
\
u(t_3)&=u(t_2)+(t_3 - t_2)k_2 \
k_2 &= f(t_2,u(t_2))\approx f(t_2,u(t)+(t_2 - t_1)k_1)\
\
u(t_4)&=u(t_3)+(t_4 - t_3)k_3 \
k_3 &= f(t_3,u(t_3))\approx\cdots

\end{align}
$$
故
$$
\begin{align}
u(t_3)&\approx u(t)+(t_2 - t_1)k_1+(t_3 - t_2)k_2 \

u(t_4)&\approx u(t)+(t_2 - t_1)k_1+(t_3 - t_2)k_2+(t_4 - t_3)k_3 \
\vdots\
u(t_m)&\approx u(t)+(t_2 - t_1)k_1+\cdots+(t_m - t_{m - 1})k_{m - 1} \
\end{align}
$$
即通过前一步的近似值和当前点的斜率 $k_i$ 逐步递推后续点的函数值近似

用类似多步法的思想
$$
\begin{align}
u(t_i)&\approx u(t_1)+\rho_{i1}k_1+\cdots+\rho_{i,i - 1}k_{i - 1}\
&=u(t_1)+h(b_{i1}k_1+\cdots+b_{i,i - 1}k_{i - 1})
\end{align}
$$
其中 $\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$ 开始，依次通过
$$
\begin{align}
k_1 &= f(t_1,u_1)\
k_2 &= f(t_2,u(t_2))\approx f(t_1 + a_2h,u_1+hb_{21}k_1)\
k_3 &= f(t_1,u(t_3))\approx f(t_1 + a_3h,u_1+h(b_{31}k_1 + b_{32}k_2))\
\vdots\
k_m &= f(t_1,u(t_m))\approx f(t_1 + a_mh,u_1+h(b_{m1}k_1+\cdots + b_{m,m - 1}k_{m - 1}))
\end{align}
$$
最终数值解
$$
\displaystyle u(t + h)\approx u_1+h(c_1k_1+\cdots + c_mk_m)=\overline{u}
$$
其中满足
$$
h(b_{i1}+\cdots + b_{i,i - 1}) = a_i\cdot h=t_i - t_1
$$

即
$$
b_{i1}+\cdots + b_{i,i - 1}=a_i
$$

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

精确解展开：
$$
\begin{align}
u(t + h)&=u(t)+hu’(t)+\cdots+\frac{h^3}{3!}u^{(3)}(t)+\cdots\
&=u(t)+hf(t,u)+\frac{h^2}{2}(f_t + f_uf)(t,u)+\frac{h^3}{6}(f_{tt}+2f_{tu}f + f_{uu}f^2 + f_uf_t +f_u^2f)(t,u)+\cdots
\end{align}
$$
记
$$
F = f_t + f_uf\
G = f_{tt}+2f_{tu}f + f_{uu}f^2
$$

Runge - kutta 格式展开：

$$
\begin{align}
\displaystyle \overline{u}&=u(t)+h(c_1k_1+\cdots + c_mk_m)\
k_1 &= f(t_1,u(t_1)) = f\
k_2 &= f(t_1 + a_2h,u_1+hb_{21}k_1)\
&=f+ha_2(f_t + f_uf)+\frac{h^2a_2^2}{2}(f_{tt}+2f_{tu}f + f_{uu}f^2)+o(h^3)\
&=f+ha_2F+\frac{h^2a_2^2}{2}G+o(h^3)\
k_3 &= f+ha_3F+h^2(a_3b_{32}f_uF+\frac{1}{2}a_3^2G)+o(h^3)

\end{align}
$$

其中 $b_{21}=a_2$

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

令
$$
F = f_t + f_uf\
G = f_{tt}+2f_{tu}f + f_{uu}f^2
$$
展开后得到方程组
$$
\begin{cases}c_1 + c_2 + c_3 = 1\c_2a_2 + c_3a_3=\frac{1}{2}\c_2a_2^2 + c_3a_3^2=\frac{1}{3}\c_3a_3b_{32}=\frac{1}{6}\end{cases}
$$
这里有 6 个未知数，4 个方程，存在 2 个自由参数，可确定三级三阶方法

经典四阶 Runge - kutta 法格式
$$
\begin{cases}u_{n + 1}=u_n+\frac{h}{6}(k_1 + 2k_2 + 2k_3 + k_4)\k_1 = f(t_n,u_n)\k_2 = f(t_n+\frac{h}{2},u_n+\frac{h}{2}k_1)\k_3 = f(t_n+\frac{h}{2},u_n+\frac{h}{2}k_2)\k_4 = f(t_n + h,u_n+hk_3)\end{cases}
$$
注：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 展开法
$$
\begin{align}
u(t_1)&=u(t_0)+u’(t_0)h+\frac{u’‘(t_0)}{2!}h^2+\cdots+\frac{u^{(p)}(t_0)}{p!}h^p+O(h^{p + 1})\
&=u(t_0)+h[u’(t_0)+\frac{u’‘(t_0)}{2!}h+\cdots+\frac{u^{(p)}(t_0)}{p!}h^{p - 1}]+O(h^{p + 1})\
&=u(t_0)+h\varphi(t_0,h)+O(h^{p + 1})
\end{align}
$$
其中
$$
\begin{align}
u’(t)&=f(t,u(t))\
u’‘(t)&=\frac{d}{dt}(f(t,u(t)))=f_t + f_uf\
u’‘’(t)&=\frac{d^2}{dt^2}(f(t,u(t)))=\cdots
\end{align}
$$
舍去余项后，得到 $u_1 = u_0 + h\varphi(t_0, u_0, h)$ ，一般形式为
$$
u_{n + 1} = u_n + h\varphi(t_n, u_n, h)
$$
例如四阶 Runge - kutta 法：
$$
\varphi(t_n, u_n, h)=\frac{1}{6}k_1+\frac{1}{3}k_2+\frac{1}{3}k_3+\frac{1}{6}k_4
$$
相容性：$\displaystyle \frac{u_{n + 1}-u_n}{h}=\varphi(t_n, u_n, h)$ 逼近微分方程，即
$$
\lim_{h \to 0}\varphi(t, u(t), h)=f(t, u(t))
$$

绝对稳定性和绝对稳定域

回顾 Euler 法稳定性：已知 $\displaystyle |e_n| \leq e^{LT}|e_0|$ 。以初值问题 $$\displaystyle \begin{cases}u’=-5u \ u(0) = 1\end{cases}$$ 为例，$f(t, u)=-5u$ ，$L = 5$ ，$T = 20$ ，$e^{LT}=e^{100}$ ，精确解 $u(t)=e^{-5t}$ ，$t \to +\infty$ 时，$u(t) \to 0$
$$
\begin{align}
u_{n + 1}&=u_n+hf(t_n, u_n)\
&=u_n - 5hu_n\
&=(1 - 5h)u_n\
&=(1 - 5h)^{n + 1}u_0
\end{align}
$$
当 $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$ 时，误差不会无限增长
$$
u{n+1}=u_n+h\mu u_n=(1+\mu h)u_n\
v_{n+1}=v_n+h\mu v_n=(1+\mu h)v_n
$$
记 $e_n=u_n-v_n$ ，则
$$
\begin{align}
e_n&=(1+\mu h)e_{n-1}\
&=(1+\mu h)^ne_0
\end{align}
$$

k 步法的绝对稳定性

对于测试方程 $u’=\lambda u$ ，线性 $k$ 步法：
$$
\sum_{j = 0}^{k} \alpha_j u_{n + j}=h\sum_{j = 0}^{k} \beta_j f_{n + j}=h\lambda\sum_{j = 0}^{k} \beta_j u_{n + j}
$$
令 $\overline{h}=h\lambda$ ，则
$$
\sum_{j = 0}^{k} (\alpha_j - \overline{h}\beta_j)u_{n + j}=0
$$
在实际计算中，每步存在舍入误差，设近似解 $\overline{u}n$ 满足
$$
\sum{j = 0}^{k} (\alpha_j - \overline{h}\beta_j)\overline{u}{n + j}=\eta_n
$$
其中 $\eta_n$ 为舍入误差，且 $|\eta_n| \leq M$ 。记 $\overline{e}n=\overline{u}n - u_n$ ，则
$$
\sum{j = 0}^{k} (\alpha_j - \overline{h}\beta_j)\overline{e}{n + j}=\eta_n
$$
记
$$
\begin{align}
E_n &= (\overline{e}{n + k - 1}, \cdots, \overline{e}_n)^T \
a_j &= \alpha_j - \overline{h}\beta_j\
\overline{\eta}_n &= (a_k^{-1}\eta_n, 0, \cdots, 0)^T
\end{align}
$$
矩阵
$$
\overline{C}=
\begin{bmatrix}

-a_k^{-1}a_{k - 1} & -a_k^{-1}a_{k - 2} &\cdots &-a_k^{-1}a_0 \

1&0&\cdots&0\

\vdots &\ddots&\ddots&\vdots\

0&\cdots&1&0

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

将误差方程
$$
\sum_{j = 0}^{k} (\alpha_j - \overline{h}\beta_j)\overline{e}{n + j}=\eta_n
$$
写成向量形式
$$
E_n=\overline{C}E{n - 1}+\overline{\eta}n
$$
进一步得到
$$
E_n=\overline{C}^nE_0+\sum{l = 0}^{n - 1} \overline{C}^l\overline{\eta}{n - 1 - l}
$$
为使误差 $E_n$ 在计算中逐步减小，要求
$$
\lim{n \to \infty}\overline{C}^n = 0
$$
$\overline{C}$ 的特征方程为
$$
\begin{align}
0=\sum_{j = 0}^{k} a_j\lambda^j&=\sum_{j = 0}^{k} (\alpha_j - \overline{h}\beta_j)\lambda^j\
&=\sum_{j = 0}^{k} \alpha_j\lambda^j-\overline{h}\sum_{j = 0}^{k} \beta_j\lambda^j\
&=\rho(\lambda)-\overline{h}\sigma(\lambda)
\end{align}
$$

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

计算中间斜率值：
$$
\begin{align}
k_1 &= f(t_n, u_n)=\lambda u_n=\lambda u_nP_0(\lambda h)\
\
k_2 &= f(t_n + a_2h, u_n+hb_{21}k_1)\
&=\lambda(u_n+hb_{21}k_1)\
&=\lambda(1 + hb_{21}\lambda)u_n\
&=\lambda u_nP_1(\lambda h)\
\
k_3 &= f(t_n + a_3h, u_n+h(b_{31}k_1 + b_{32}k_2))=\lambda u_nP_2(\lambda h)\
\vdots\
k_{m - 1}&=\lambda u_nP_{m - 1}(\lambda h)

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

计算数值解：
$$
\begin{align}
u_{n + 1}&=u_n+h(c_1k_1+\cdots + c_{m - 1}k_{m - 1})\
&=u_n+h(c_1\lambda u_n + c_2P_1(\lambda h)\lambda u_n+\cdots + c_{m - 1}P_{m - 1}(\lambda h)\lambda u_n)\
&=u_n(1 + P_m(\lambda h))
\end{align}
$$
其中 $P_m(\lambda h)$ 是关于 $\lambda h$ 的多项式

与精确解 Taylor 展开对比

将精确解 $u(t_{n + 1})$ 在 $t_n$ 处 Taylor 展开：
$$
\begin{align}
u(t_{n + 1})&=u(t_n)+u’(t_n)h+\cdots+\frac{u^{(m)}(t_n)}{m!}h^m+O(h^{m + 1})\
&=u(t_n)+u(t_n)\lambda h+\cdots+\frac{u(t_n)\lambda^m h^m}{m!}+O(h^{m + 1})\
&=(1+\lambda h+\cdots+\frac{(\lambda h)^m}{m!})u(t_n)+O(h^{m + 1})

\end{align}
$$
由于 $m$ 级 $m$ 阶 Runge - kutta 法的精度要求，使得
$$
1 + P_m(\lambda h)=1+\lambda h+\cdots+\frac{(\lambda h)^m}{m!}
$$
由 $u_{n+1}=u_n(1 + P_m(\lambda h))$ 其特征方程为 $\lambda - [1 + P_m(\lambda h)] = 0$ ，其中
$$
\begin{align}
\lambda_1 &= 1 + P_m(\lambda h)=1+\lambda h+\cdots+\frac{(\lambda h)^m}{m!}\
D_A &={\overline{h}:|\lambda_1(\overline{h})| < 1}={\overline{h}:|1+\overline{h}+\cdots+\frac{\overline{h}^m}{m!}| < 1}\
& \quad \quad \quad \quad \quad \overline{h} =\lambda h\
\end{align}
$$
通过求解该不等式确定 $$\overline{h}$$ 的取值范围，从而得到方法的绝对稳定域，判断在哪些 $\lambda$ 和 $h$ 组合下数值解稳定

一阶方程组和刚性问题

常微分方程组（ODEs）相关

常微分方程组：一般形式为
$$
\begin{cases}u_1’ = f_1(u_1,\cdots,u_m)\\vdots\u_m’ = f_m(u_1,\cdots,u_m)\end{cases}
$$
初值条件为 $u_i(t_0)=u_{i0}$ ，$i = 1,\cdots,m$

高阶常微分初值问题：如
$$
\begin{cases}u^{(m)} = f(t,u,u’,\cdots,u^{(m - 1)})\u(t_0)=u_0, u’(t_0)=u_0’,\cdots,u^{(m - 1)}(t_0)=u_0^{(m - 1)}\end{cases}
$$
令 $u_1 = u$ ，$u_2 = u’$ ，$\cdots$ ，$u_m = u^{(m - 1)}$ ，可转化为常微分方程组
$$
\begin{cases}u_1’ = u_2\u_2’ = u_3\\vdots\u_m’ = f(t,u_1,u_2,\cdots,u_m)\end{cases}
$$
写成向量形式为 $\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=\begin{bmatrix}-0.1&-49.9&0\0&-50&0\0&70&-3\times10^4\end{bmatrix}
$$
$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$ 的特征值，若满足：
$$
Re(\lambda_i) < 0 ，i = 1,\cdots,m\
\frac{\max|Re(\lambda_i)|}{\min|Re(\lambda_i)|} = R \gg 1
$$
则称方程组是刚性的，$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$
$$
\begin{align}
u_1 &= u_0+hf(t_0,u_0)\
\overline{u}1&=\overline{u}{\frac{1}{2}}+\frac{h}{2}f(t_{\frac{1}{2}},\overline{u}{\frac{1}{2}})=u_0+\frac{h}{2}f(t_0,u_0)+\frac{h}{2}f(t{\frac{1}{2}},\overline{u}_{\frac{1}{2}})
\end{align}
$$
通过 $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)+\varepsilon(t)h^p+O(h^{p + 1})
$$
则
$$
u(t,\frac{h}{2})=u(t)+\varepsilon(t)(\frac{h}{2})^p+O(h^{p + 1})
$$
对 $u(t,h)$ 和 $u(t,\frac{h}{2})$ 进行线性组合：

设组合式为
$$
\frac{2^p u(t,\frac{h}{2})-u(t,h)}{2^p - 1}
$$
将误差展开式代入：

分子为
$$
\begin{align}
&2^p\left[u(t)+\varepsilon(t)(\frac{h}{2})^p+O(h^{p + 1})\right]-\left[u(t)+\varepsilon(t)h^p+O(h^{p + 1})\right]\
&=2^p u(t)+\varepsilon(t)h^p+2^pO(h^{p + 1})-u(t)-\varepsilon(t)h^p - O(h^{p + 1})\
&=(2^p - 1)u(t)+O(h^{p + 1})
\end{align}
$$
所以
$$
\frac{2^p u(t,\frac{h}{2})-u(t,h)}{2^p - 1}=u(t)+O(h^{p + 1})
$$
通过这种线性组合，提高了数值解的精度阶数，从 $O(h^p)$ 提升到 $O(h^{p + 1})$ ，这就是外推法提高数值计算精度的基本原理