Taylor 算子

高阶算子的核心思想：通过引入更宽模板（更长的算子长度），利用更多邻域节点的函数值加权组合来逼近导数，使截断误差降至 $O(h^p)$（$p>2$），从而在相同网格步长下获得更高精度，或在相同精度下放宽空间采样准则。

设函数 $f(x)$ 在 $x$ 处充分光滑，其在相邻节点 $x + j\Delta x$ 处的 Taylor 展开为：

$$f(x + j\Delta x) = \sum_{n=0}^{\infty} \frac{(j\Delta x)^n}{n!} f^{(n)}(x)$$

构造线性组合逼近 $m$ 阶导数：

$$f^{(m)}(x) \approx \frac{1}{(\Delta x)^m} \sum_{j=-L}^{R} w_j^{(m)} f(x + j\Delta x)$$

其中 $w_j^{(m)}$ 为待定权重系数，$L$、$R$ 分别为左侧与右侧模板长度。将展开式代入并比较同阶导数项系数，可建立关于 $w_j^{(m)}$ 的线性方程组。

双侧算子

考虑模板 $\{x-2h, x-h, x+h, x+2h\}$，设逼近形式：

$$a f(x-2h) + b f(x-h) + c f(x+h) + d f(x+2h) \approx f^{(m)}(x)$$

将各点 Taylor 展开至四阶并求和，按 $f, f', f'', f''', f^{(4)}$ 分组整理，要求：

• $f$ 项系数为 0（$m>0$ 时）或 1（$m=0$ 插值时）
• $f^{(m)}$ 项系数为 $(\Delta x)^m$
• 其余低阶项系数为 0

插值（$m=0$）方程组：

$$\begin{cases} a + b + c + d = 1 \\ -2a - b + c + 2d = 0 \\ 4a + b + c + 4d = 0 \\ -8a - b + c + 8d = 0 \end{cases} \quad \Rightarrow \quad \begin{bmatrix} a \\ b \\ c \\ d \end{bmatrix} = \begin{bmatrix} 1/6 \\ -2/3 \\ 2/3 \\ -1/6 \end{bmatrix}$$

一阶导数（$m=1$）方程组：

$$\begin{cases} a + b + c + d = 0 \\ -2a - b + c + 2d = 1/\Delta x \\ 4a + b + c + 4d = 0 \\ -8a - b + c + 8d = 0 \end{cases} \quad \Rightarrow \quad \begin{bmatrix} a \\ b \\ c \\ d \end{bmatrix} = \frac{1}{\Delta x} \begin{bmatrix} 1/12 \\ -2/3 \\ 2/3 \\ -1/12 \end{bmatrix}$$

二阶导数（$m=2$）需 5 点模板 $\{x-2h, x-h, x, x+h, x+2h\}$，类似推导得：

$$f''(x) \approx \frac{1}{12(\Delta x)^2} \left[ -f_{i-2} + 16f_{i-1} - 30f_i + 16f_{i+1} - f_{i+2} \right] + O(h^4)$$

单侧算子

在计算域边界处，无法使用双侧模板，需构造单侧差分格式。以右侧边界为例，模板 $\{x, x+h, x+2h, x+3h\}$ 逼近一阶导数：

方程组：

$$\begin{cases} a + b + c + d = 0 \\ 0a + 1b + 2c + 3d = 1/\Delta x \\ 0a + 1b + 4c + 9d = 0 \\ 0a + 1b + 8c + 27d = 0 \end{cases} \quad \Rightarrow \quad \begin{bmatrix} a \\ b \\ c \\ d \end{bmatrix} = \frac{1}{\Delta x} \begin{bmatrix} -11/6 \\ 3 \\ -3/2 \\ 1/3 \end{bmatrix}$$

单侧算子特性：

• 权重系数绝对值随模板长度增加而迅速增大，易引入数值振荡
• 精度阶数受限于模板长度，$p$ 点单侧模板最高可达 $O(h^{p-1})$
• 实际应用中常采用"内侧高阶 + 边界低阶"的混合策略，或引入虚拟节点延拓

权重系数通用表达式（Fornberg, 1996）

对 $2p$ 阶精度的双侧中心差分算子，一阶导数权重系数为：

$$w_j^{(1)} = \begin{cases} 0, & j=0 \\ \displaystyle \frac{(-1)^{j+1} (p!)^2}{j \cdot (p-j)! (p+j)!}, & j = \pm 1, \pm 2, \dots, \pm p \end{cases}$$

该闭式表达便于编程实现任意阶数算子自动生成。

时间域高阶外推

空间精度提升后，时间离散误差可能成为主导。本节介绍三类时间高阶外推方法。

高阶 Taylor 外推（以声波方程为例）

考虑二维声波方程：

$$\frac{\partial^2 P}{\partial t^2} = c^2 \left( \frac{\partial^2 P}{\partial x^2} + \frac{\partial^2 P}{\partial z^2} \right)$$

在时间层 $t_n$ 处对 $P(t_n + \Delta t)$ 作 Taylor 展开：

$$P^{n+1} = P^n + \Delta t \dot{P}^n + \frac{(\Delta t)^2}{2!} \ddot{P}^n + \frac{(\Delta t)^3}{3!} \dddot{P}^n + \cdots$$

关键观察：波动方程中时间奇阶导数可由空间导数递推表示：

$$\ddot{P} = c^2 \nabla^2 P, \quad \dddot{P} = c^2 \nabla^2 \dot{P}, \quad P^{(4)} = c^4 \nabla^4 P, \quad \dots$$

因此高阶时间外推格式可写为：

$$P^{n+1} = 2P^n - P^{n-1} + \sum_{k=1}^{N} \frac{2(\Delta t)^{2k}}{(2k)!} c^{2k} \nabla^{2k} P^n$$

实现要点：

• 仅需存储两个时间层数据，内存开销低
• 空间高阶导数 $\nabla^{2k} P$ 可通过重复应用空间差分算子计算
• 稳定性条件随阶数提高而更严格，需联合优化 $\Delta t$ 与 $h$

预估 - 校正法（Predictor-Corrector）

以一阶方程 $\partial T/\partial t = f(T,t)$ 为例，修正的 Euler 方法包含两步：

预估步（显式 Euler）：

$$T^* = T^n + \Delta t \cdot f(T^n, t_n)$$

校正步（隐式 Euler 线性化）：

$$T^{n+1} = T^n + \frac{\Delta t}{2} \left[ f(T^n, t_n) + f(T^*, t_{n+1}) \right]$$

特性分析：

• 截断误差 $O(\Delta t^2)$，精度高于纯显式/隐式 Euler
• 稳定性区域与显式 Euler 相同，但数值耗散更小
• 每步需两次函数评估，计算成本约为显式的 2 倍

Runge-Kutta 方法

龙格-库塔（Runge-Kutta, R-K）方法是一类通过计算多个中间斜率并进行线性组合来构造高阶时间外推格式的显式积分方法。其核心思想是：在不直接计算微分方程右端函数高阶偏导数的前提下，通过对 $f(T,t)$ 在不同中间点上的多次求值，使离散格式的 Taylor 展开与精确解的 Taylor 展开在尽可能高的阶数上匹配，从而提升时间积分精度。R-K 方法可视为通过多次斜率计算实现的“高级预估-校正”过程，完全显式，无需求解隐式方程组。

二阶 Runge-Kutta 方法

构造原理与系数匹配
考虑一阶常微分方程初值问题 $\frac{dT}{dt} = f(T,t)$。设更新格式为两斜率线性组合：

$$T_{j+1} = T_j + a k_1 + b k_2$$

其中 $k_1 = f(T_j, t_j)\Delta t$，$k_2 = f(T_j + \alpha k_1, t_j + \beta \Delta t)\Delta t$。将 $k_2$ 在 $(T_j, t_j)$ 处作二元 Taylor 展开至一阶，并与精确解 $T(t_{j+1}) = T_j + \Delta t f + \frac{(\Delta t)^2}{2}(f_t + f f_T) + O(\Delta t^3)$ 对比同阶项系数，可得待定系数需满足的方程组：

$$a + b = 1, \quad b\alpha = \frac{1}{2}, \quad b\beta = \frac{1}{2}$$

该方程组存在无穷多组解。采用一组经典解：$a=0, b=1, \alpha=\frac{1}{2}, \beta=\frac{1}{2}$，得到二阶 R-K 格式（中点格式）：

$$\begin{aligned} k_1 &= f(T_j, t_j) \Delta t \\ k_2 &= f\left(T_j + \frac{1}{2}k_1, t_j + \frac{1}{2}\Delta t\right) \Delta t \\ T_{j+1} &= T_j + k_2 \end{aligned}$$

数值特性

• 精度：局部截断误差 $O(\Delta t^3)$，全局误差 $O(\Delta t^2)$，显著优于一阶显式 Euler 格式。
• 稳定性：绝对稳定区域为复平面上满足 $\left|1 + z + \frac{z^2}{2}\right| \leq 1$ 的区域（$z=\lambda \Delta t$），实轴稳定区间为 $[-2, 0]$。对衰减型问题要求 $\Delta t \leq 2/|\lambda|$。
• 物理意义：利用起点斜率预估中点状态，再以中点斜率推进整步长。在波场模拟中表现出优于简单 Euler 格式的相位保持特性，数值耗散较低。

四阶 Runge-Kutta 方法（经典格式）

为进一步提升精度并扩大稳定区域，引入四个中间斜率 $k_1, k_2, k_3, k_4$，通过加权平均构造四阶格式。其标准形式为：

$$\begin{aligned} k_1 &= f(T_j, t_j) \Delta t \\ k_2 &= f\left(T_j + \frac{1}{2}k_1, t_j + \frac{1}{2}\Delta t\right) \Delta t \\ k_3 &= f\left(T_j + \frac{1}{2}k_2, t_j + \frac{1}{2}\Delta t\right) \Delta t \\ k_4 &= f\left(T_j + k_3, t_j + \Delta t\right) \Delta t \\ T_{j+1} &= T_j + \frac{1}{6}(k_1 + 2k_2 + 2k_3 + k_4) \end{aligned}$$

算法实现步骤

1. 计算起点斜率 $k_1$：基于当前层 $(T_j, t_j)$ 的右端函数值。
2. 第一次中点预估 $k_2$：利用 $k_1/2$ 推进半步长，计算中点斜率。
3. 第二次中点修正 $k_3$：利用 $k_2/2$ 重新推进半步长，修正中点斜率，提高中点状态估计的可靠性。
4. 终点预估 $k_4$：利用 $k_3$ 推进全步长，计算终点斜率。
5. 加权更新：以 $1:2:2:1$ 的权重组合四个斜率，得到下一层解。该权重分配等效于 Simpson 积分公式，保证了四阶代数精度。

数值特性

• 精度：局部截断误差 $O(\Delta t^5)$，全局误差 $O(\Delta t^4)$。无需显式计算高阶偏导数即可实现高阶精度，编程实现直观。
• 稳定性：绝对稳定区域显著扩大，实轴稳定区间约为 $[-2.785, 0]$。允许在刚性较弱或中等刚性问题中采用更大时间步长，误差累积速度远低于低阶格式。
• 计算成本：每时间步需 4 次右端函数评估，内存开销极低（仅需存储当前层状态），易于与显式空间差分格式耦合。

版权所有

版权归属：Staaaaaaaaar

许可证：署名-非商业性-相同方式共享 4.0 国际 (CC-BY-NC-SA-4.0)

编辑此页

贡献者: Staaaaaaaaar

上一页有限差分方法