低阶差分算子只使用少量相邻节点。它们形式简单，但在波动问题中往往需要非常密集的网格才能控制相位误差。高阶差分算子通过扩大空间模板、匹配更多 Taylor 展开项，在不显著减小网格步长的情况下提高导数近似精度。

高阶算子的分析重点包括两部分：一是如何求出不同节点的权重；二是这些权重对有效波数、数值频散和稳定性的影响。

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
• 仍为条件稳定；对负实轴问题的稳定区间为 $[-2,0]$，与显式 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 次右端函数评估，内存开销极低（仅需存储当前层状态），易于与显式空间差分格式耦合。

交错网格高阶差分算子

传统中心差分在整数网格点 $x_i=ih$ 处利用 $x_i$ 两侧的整数节点逼近导数。交错网格算子则在半网格点 $x_{i+1/2}$ 处计算导数，例如二阶精度格式为

$$f'_{i+1/2}=\frac{f_{i+1}-f_i}{h}+O(h^2).$$

利用更宽的对称模板可以构造四阶精度格式：

$$f'_{i+1/2} =\frac{1}{h}\left[ \frac{9}{8}(f_{i+1}-f_i) -\frac{1}{24}(f_{i+2}-f_{i-1}) \right]+O(h^4).$$

一般的 $2M$ 阶交错网格一阶导数可以写为

$$f'_{i+1/2}\approx \frac{1}{h} \sum_{m=1}^{M}a_m \left(f_{i+m}-f_{i-m+1}\right),$$

其中系数 $a_m$ 由 Taylor 展开匹配得到。与同阶传统网格算子相比，交错网格算子在波数域内通常具有更宽的高精度区间，尤其适合速度—应力一阶方程组。

注

交错网格的优势不仅来自几何位置错开，还来自物理量与微分关系的对应。例如速度位于整数节点、应力位于半整数节点时，应力梯度恰好在速度节点上计算，能够避免多余插值。

边界附近缺少对称节点时，需要采用单侧高阶算子、虚拟节点或降低局部阶数。若直接截断内部算子，通常会破坏整体精度，甚至影响稳定性。

傅里叶微分算子

设周期函数 $f(x)$ 的傅里叶变换为 $\hat f(k)$：

$$f(x)=\frac{1}{2\pi}\int_{-\infty}^{\infty}\hat f(k)e^{ikx}\,dk.$$

对上式求导可得

$$\frac{df}{dx} =\frac{1}{2\pi}\int_{-\infty}^{\infty}ik\hat f(k)e^{ikx}\,dk,$$

因此

$$\mathcal F\left[\frac{df}{dx}\right]=ik\hat f(k), \qquad \mathcal F\left[\frac{d^2f}{dx^2}\right]=-k^2\hat f(k).$$

傅里叶微分算子的计算步骤为：

1. 对离散函数值作快速傅里叶变换，得到 $\hat f(k)$；
2. 一阶导数乘以 $ik$，二阶导数乘以 $-k^2$；
3. 作逆快速傅里叶变换，返回空间域导数。

对于网格间距 $h$、节点数 $N$、区间长度 $L=Nh$ 的离散模型，可解析的波数范围满足

$$|k|\leq k_N=\frac{\pi}{h},$$

其中 $k_N$ 为 Nyquist 波数。傅里叶算子在离散可解析波数范围内直接使用真实波数 $k$，而有限差分算子使用修正波数 $k^*$，两者的偏差决定空间离散误差。

有效波数与差分精度

将平面波

$$f(x)=e^{ikx}$$

代入二阶中心差分一阶导数算子，有

$$\frac{f(x+h)-f(x-h)}{2h} =i\frac{\sin(kh)}{h}e^{ikx}.$$

因此二阶中心差分的有效波数为

$$k^*=\frac{\sin(kh)}{h},$$

而精确微分对应 $k^*=k$。定义无量纲波数 $\theta=kh$，则

$$\frac{k^*}{k}=\frac{\sin\theta}{\theta}.$$

当 $\theta\ll 1$ 时二者接近；接近 Nyquist 波数时偏差迅速增大。高阶中心差分的一般有效波数为

$$k^*h=2\sum_{m=1}^{M}a_m\sin(m\theta),$$

交错网格算子则为

$$k^*h=2\sum_{m=1}^{M}a_m \sin\left[\left(m-\frac12\right)\theta\right].$$

比较 $k^*h$ 与精确值 $\theta$，可以直接判断不同算子的有效波数范围。算子阶数越高，曲线在更大的 $\theta$ 范围内贴近直线 $k^*h=\theta$，因而每个波长所需的节点数可以减少。

注意

高阶格式降低的是空间截断误差，并不会自动消除时间离散误差、边界反射或模型采样误差。空间阶数提高后，时间离散往往成为新的主要误差来源。

不同微分算子的比较

算子	空间模板	主要特点	主要限制
低阶中心差分	局部、窄模板	实现简单，单步成本低	高频误差大，数值频散明显
高阶 Taylor 差分	局部、宽模板	可系统构造，精度随阶数提高	边界处理复杂，通信宽度增加
交错网格高阶差分	局部、错位模板	有效波数范围较宽，适合一阶波动方程	变量布局和边界实现更复杂
傅里叶微分算子	全局	对光滑周期函数具有谱精度	隐含周期条件，全局通信成本高

高阶差分可看作局部多项式近似；傅里叶微分则是全局三角函数近似。当问题具有规则网格、周期性或可接受周期延拓时，傅里叶算子能提供很高精度；当模型复杂、并行区域分解或局部边界处理更重要时，高阶有限差分通常更灵活。

版权所有

版权归属：Staaaaaaaaar

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

编辑此页

贡献者: Staaaaaaaaar

上一页有限差分方法

下一页傅里叶伪谱方法