有限差分方法 (Finite Difference Method, FDM) 是计算机数值模拟最早采用的方法，至今仍被广泛运用。

有限差分方法将求解域划分为差分网格，用有限个网格节点代替连续的求解域，以 Taylor 级数展开等方法，把控制方程中的导数用网格节点上的函数值的差商代替进行离散，从而建立以网格节点上的值为未知数的代数方程组。该方法是一种直接将微分问题变为代数问题的近似数值解法，数学概念直观，表达简单，是发展较早且比较成熟的数值方法。

网格剖分

利用有限差分方法求解偏微分方程问题必须把连续问题进行离散化，首先要求对求解区域给出网格剖分。

设由 $(x,y,z,t)$ 定义的四维求解区域为 $D=D^I \cup D^B$，$D^I$ 和 $D^B$ 分别表示求解区域的内部和边界。

对于直角坐标系，在空间上用平行于坐标轴的直线族将区域划分成一系列六面体，这样的直线称为网格线 (grid line)，它们的交点称为网格点或节点 (grid)。

如果设区域原点坐标为 $(x_0,y_0,z_0)$，每个坐标方向上的网格点等间距分布，分别为 $\Delta x$, $\Delta y$, $\Delta z$，称其为空间步长 (grid spacing)。相应地，在时间上也可以用等间距的时间步长 (time step) $\Delta t$ 来划分求解区域。

$$x_i = x_0 + i \Delta x, \quad y_j = y_0 + j \Delta y, \quad z_k = z_0 + k \Delta z, \quad t_n = t_0 + n \Delta t$$

注

对具体地球物理问题，需要选择合理的网格剖分形式。在许多实际应用中，空间步长相等的六面体网格点是自然而合理的选择。如果需要考虑非均匀体的形状空间变化、不连续面的分布，有时需要选择与之适应的剖分形式。有些情况下需要选择能够简化差分方程的剖分形式。

如果两个节点沿坐标轴只相差一个步长时，称其为相邻的两个节点。如果一个节点的所有相邻节点都属于 $D$，则称此节点为内部节点，否则称为边界节点，分别用 $D^I_\Delta$ 和 $D^B_\Delta$ 表示内部节点和边界节点的集合。

差分形式

对地球物理问题的连续求解区域通过网格剖分离散为空间上的一系列网格点，接下来需要利用一定的差分形式对偏微分方程组中的导数用差商进行近似，从而将偏微分方程组离散化为差分方程组。

对于满足一定条件的函数，我们可以利用 Taylor 级数展开来给出微商的近似形式。

一阶精度近似：

• 一阶导数向前差分形式：

  $$\frac{d}{d x} f(x) = \frac{f(x+\Delta x) - f(x)}{\Delta x} + O(\Delta x)$$

• 一阶导数向后差分形式：

  $$\frac{d}{d x} f(x) = \frac{f(x) - f(x-\Delta x)}{\Delta x} + O(\Delta x)$$

二阶精度近似：

• 一阶导数中心差分形式：

  $$\frac{d}{d x} f(x) = \frac{f(x+\Delta x) - f(x-\Delta x)}{2 \Delta x} + O(\Delta x^2)$$

• 二阶导数中心差分形式：

  $$\frac{d^2}{d x^2} f(x) = \frac{f(x+\Delta x) - 2f(x) + f(x-\Delta x)}{\Delta x^2} + O(\Delta x^2)$$

差分方程

对于具体地球物理问题的偏微分方程组，利用差分格式，可以给出偏导数的微商近似，进一步得到差分方程组。

一个一般的偏微分方程可以表示为

$$L(u(P))=f(P), \quad P \in D^I$$

$$B(u(P))=g(P), \quad P \in D^B$$

那么以上方程的有限差分方程可以表示为

$$L_\Delta(u(P))=f(P), \quad P \in D^I_\Delta$$

$$B_\Delta(u(P))=g(P), \quad P \in D^B_\Delta$$

下面以对流方程（双曲型）的初值问题为例，给出其差分方程的不同形式。

$$\frac{\partial u}{\partial t} + a \frac{\partial u}{\partial x} = 0, \quad x \in \mathbb{R},\quad t > 0$$

$$u(x,0) = g(x), \quad x \in \mathbb{R}$$

取空间步长为 $h$，时间步长为 $\tau$，方程中的微商可以有下面这些差分形式：

$$\frac{u(x_j,t_{n+1}) - u(x_j,t_n)}{\tau} = \left[\frac{\partial u}{\partial t}\right]_j^n + O(\tau) \tag{1}$$

$$\frac{u(x_j,t_{n})-u(x_j,t_{n-1})}{\tau} = \left[\frac{\partial u}{\partial t}\right]_j^n + O(\tau) \tag{2}$$

$$\frac{u(x_{j},t_{n+1}) - u(x_j,t_{n-1})}{2\tau} = \left[\frac{\partial u}{\partial t}\right]_j^n + O(\tau^2) \tag{3}$$

$$\frac{u(x_{j+1},t_{n}) - u(x_j,t_{n})}{h} = \left[\frac{\partial u}{\partial x}\right]_j^n + O(h) \tag{4}$$

$$\frac{u(x_{j},t_{n}) - u(x_{j-1},t_{n})}{h} = \left[\frac{\partial u}{\partial x}\right]_j^n + O(h) \tag{5}$$

$$\frac{u(x_{j+1},t_{n}) - u(x_{j-1},t_{n})}{2h} = \left[\frac{\partial u}{\partial x}\right]_j^n + O(h^2) \tag{6}$$

利用 (1)(4) 可以得到两层显式格式的差分方程：

$$u_j^{n+1} = u_j^n - a \lambda (u_{j+1}^n - u_{j}^n)$$

其中 $\lambda = \frac{\tau}{h}$ 称为网格比。

        |         |
n+1 ----o---------+----
        |         |
        |         |
  n ----o---------o----
        |         |
        j        j+1

利用 (1)(5) 可以得到另一种两层显式格式的差分方程：

$$u_j^{n+1} = u_j^n - a \lambda (u_{j}^{n} - u_{j-1}^{n})$$

        |         |
n+1 ----o---------+----
        |         |
        |         |
  n ----o---------o----
        |         |
       j-1        j

上面的两种差分方程都是偏心差分格式，利用 (1)(6) 可以得到两层中心差分格式的差分方程：

$$u_j^{n+1} = u_j^n - \frac{a\lambda}{2} (u_{j+1}^n - u_{j-1}^{n})$$

        |         |         |
n+1 ----+---------o---------+----
        |         |         |
        |         |         |
  n ----o---------o---------o----
        |         |         |
       j-1        j        j+1

需要注意的是，用不同差分形式的精度和误差不同，得到的差分方程的数值解的稳定性和收敛性也不同。

差分方程和初始条件、边界条件的离散构成了差分格式。

差分格式的性质

牛顿冷却定律

考虑一个一阶常微分方程：

$$\frac{d T}{d t} = f(T,t), \quad T=T(t), \quad t > 0$$

初始条件为 $T(0) = T_0$

首先求其数值解，对 $T$ 和 $t$ 进行离散：

$$t_j=t_0+jdt, \quad T_j=T(t_j)$$

利用向前差分形式构造显式差分格式：

$$T_{j+1}=T_j+\tau f(T_j,t_j)$$

我们取 $f(T,t)=-T/\tau$，那么差分方程为：

$$T_{j+1} = T_j - \frac{dt}{\tau}T_j = T_j\left(1-\frac{dt}{\tau}\right)$$

对应的数值解为：

$$T_j = T_0 \left(1-\frac{dt}{\tau}\right)^j$$

而方程的解析解为：

$$T(t) = T_0 e^{-t/\tau}$$

当时间步长 $dt$ 趋于零时，差分格式给出的数值解收敛于解析给出的严格解。

当 $0\leq dt < \tau$ 时，数值解是稳定且收敛的。

当 $\tau \leq dt < 2\tau$ 时，数值解振荡，但最终收敛。

当 $dt \geq 2\tau$ 时，数值解振荡且发散。

相容性

设 PDE 表示定解问题的偏微分方程，FDE 表示相应的差分格式，如果当时间和空间步长分别趋于零时，差分格式的截断误差（PDE 和 FDE 的差）趋于零

$$\lim_{\Delta t, \Delta x \to 0} |PDE - FDE| = 0$$

则称差分格式与偏微分方程是相容的，或称差分格式与偏微分方程具有相容性。

相容性描述了当差分网格不断精细直至趋于零时，差分格式是否能够与偏微分方程充分接近的性质。只有满足了相容性，差分格式对偏微分方程的离散才有意义，否则差分格式求解的问题与偏微分方程不具有等价性。

通常情况下，时间步⻓和空间步⻓之间是独立的。如果上述条件在时间步⻓和空间步⻓之间满足某种关系是才成立，则称为有条件的相容。

稳定性

如果偏微分方程的严格解析解有界，差分格式给出的解也有界，称该差分格式是稳定的；如果差分格式给出的解是无界的，则称该差分格式是不稳定的。如果差分格式给出的解对于所有的时间和空间步长取值都是有界的，称该差分格式是无条件稳定的；如果只是对时间和空间步长的部分取值有界，称它是有条件稳定的；如果对于所有的时间和空间步长取值都是无界的，称该差分格式是无条件非稳定的。

设 $u_j^0$ 的误差是 $\varepsilon_j^0$，传递到待求的 $u_j^{n+1}$ 的误差为 $\varepsilon_j^{n+1}$，考察 $[0,T]$ 时间段，如果存在正常数 $\Delta x_0, \Delta t_0$，及一个非负常数 $K$，使得

$$||\varepsilon_j^{n+1}|| \leq K ||\varepsilon_j^0||$$

对 $0 < \Delta t < \Delta t_0, 0 < \Delta x < \Delta x_0, (n+1) \Delta t < T$ 均成立，则称该差分格式具有稳定性。

稳定性理论分析只能对线性偏微分方程进行，对于非线性偏微分方程，首先必须线性化，然后对线性化以后的偏微分方程的差分格式进行稳定性分析。

稳定性分析方法主要有 Fourier 技术方法、矩阵分析法、能量法等。最常用的稳定性分析方法是 von Neumann 方法。其基本思想是将某一个时间步在某一个网格节点上的离散解用有限 Fourier 级数展开，讨论每一个 Fourier 分量的稳定性。因此，该方法讨论了局部的稳定性。当且仅当每一个 Fourier 分量都稳定，离散解才是稳定的。

Von Neumann 分析方法使用于常系数的线性偏微分方程。尽管有限 Fourier 级数展开隐含了空间上的周期性条件，但是对于不是周期性变化的情况，它也可以给出有用的稳定性分析结果。

收敛性

设 PDE 和 FDE 的解分别为 $u_{i,j,k}^m,U_{i,j,k}^m$，如果当时间和空间步长趋于零时，FDE 解趋于 PDE 解

$$\lim_{h \to 0, \tau \to 0} |u_{i,j,k}^m - U_{i,j,k}^m| = 0$$

则称该差分格式是收敛的。

对连续形式的偏微分方程进行有限差分离散时，差分格式对最终计算结果有重要影响。收敛性表示当差分网格无限细化时，差分方程的解是否具有无限逼近偏微分方程的解的能力。收敛性的判断决定了一个差分格式能否被用来离散偏微分方程，不收敛的差分格式无疑是无法得到偏微分方程的真实解的。

Lax 等价定理

相容性是差分格式的性质，它将差分格式和原始的偏微分方程联系在一起。稳定性和收敛性则是差分格式给出的数值解的性质。

一般来讲，分析稳定性和相容性比较容易，证明收敛性有时是很困难的数学问题。

Lax 等价定理表明，如果逼近一个给定问题的差分格式是相容的，那么该差分格式的收敛性与稳定性互为充分必要条件。

由于判断差分格式的收敛性相对比较困难，该定理提供了通过判断差分格式稳定性来确定收敛性的方法。因此，一般不再讨论收敛性问题，而是讨论差分格式的稳定性。

一维非均匀介质的弹性波动方程

注

此处的一维指的是介质参数和波场只随一个空间坐标变化，但质点振动的方向可以是任意的。

对于一维弹性非均匀模型，设密度和弹性常数是随坐标变化的连续函数，可以将 $\mathbf{u}(x,t)$ 分解为三个标量函数，对应 P 波、SV 波和 SH 波的标量波动方程，由于这三个方程形式相同，可以统一表示为下面三种形式之一。

$$\rho(x) \frac{\partial^2 u(x,t)}{\partial t^2} = \frac{\partial \tau(x,t)}{\partial x} ,\quad \tau(x,t) = E(x) \frac{\partial u(x,t)}{\partial x} \tag{1}$$

$$\rho(x) \frac{\partial^2 u(x,t)}{\partial t^2} = \frac{\partial}{\partial x}\left[E(x) \frac{\partial u(x,t)}{\partial x}\right] \tag{2}$$

$$\frac{\partial v(x,t)}{\partial t} = \frac{1}{\rho(x)} \frac{\partial \tau(x,t)}{\partial x} ,\quad \frac{\partial \tau(x,t)}{\partial t} = E(x) \frac{\partial v(x,t)}{\partial x} \tag{3}$$

对于 P 波，$u$ 表示质点在 $x$ 方向的位移，$E(x)=\lambda(x)+2\mu(x)$；对于 S 波，$u$ 表示质点在 $y,z$ 方向的位移，$E(x)=\mu(x)$。方程 (1) 被称为应力-位移方程，方程 (2) 被称为位移方程，方程 (3) 被称为速度-应力方程。

时间向前、空间中心差分格式

离散格式
对速度-应力方程组：

$$\frac{\partial v}{\partial t} = \frac{1}{\rho(x)} \frac{\partial \tau}{\partial x}, \quad \frac{\partial \tau}{\partial t} = E(x) \frac{\partial v}{\partial x}$$

在规则网格 $(lh, m\Delta t)$ 上，时间采用一阶向前差分，空间采用二阶中心差分，得到：

$$\frac{v_l^{m+1} - v_l^m}{\Delta t} = \frac{1}{\rho_l} \frac{\tau_{l+1}^m - \tau_{l-1}^m}{2h}, \quad \frac{\tau_l^{m+1} - \tau_l^m}{\Delta t} = E_l \frac{v_{l+1}^m - v_{l-1}^m}{2h}$$

稳定性结论：无条件不稳定
将单频谐波试探解 $v_l^m = V e^{i(klh - \omega m\Delta t)}$、$\tau_l^m = T e^{i(klh - \omega m\Delta t)}$ 代入格式，消去振幅后可得增长因子：

$$e^{-i\omega \Delta t} = 1 \pm i \frac{c \Delta t}{h} \sin(kh), \quad c = \sqrt{E/\rho}$$

其模长恒为 $|e^{-i\omega \Delta t}| = \sqrt{1 + \left( \frac{c \Delta t}{h} \sin(kh) \right)^2} > 1$。增长因子模长大于 1 表明数值解振幅随时间呈指数增长，该格式对任意步长组合均不稳定。其物理根源在于速度与应力变量同点定义，时间单侧差分与空间中心差分的相位失配导致数值能量非物理累积。

时间中心、空间中心差分格式

离散格式
将时间微分改用二阶中心差分近似，空间保持中心差分：

$$\frac{v_l^{m+1} - v_l^{m-1}}{2\Delta t} = \frac{1}{\rho_l} \frac{\tau_{l+1}^m - \tau_{l-1}^m}{2h}, \quad \frac{\tau_l^{m+1} - \tau_l^{m-1}}{2\Delta t} = E_l \frac{v_{l+1}^m - v_{l-1}^m}{2h}$$

稳定性结论：有条件稳定
代入谐波试探解并消去振幅，得到离散频散关系：

$$\sin(\omega \Delta t) = \pm \frac{c \Delta t}{h} \sin(kh)$$

为保证角频率 $\omega$ 为实数（解不指数增长），需满足：

$$\left| \frac{c \Delta t}{h} \sin(kh) \right| \leq 1 \quad \Rightarrow \quad \frac{c \Delta t}{h} \leq 1$$

由此导出CFL（Courant–Friedrichs–Lewy）稳定性条件：

$$\Delta t \leq \frac{h}{c_{\max}}$$

其中 $c_{\max}$ 为计算域内的最大波速。物理意义：时间步长不得超过波在单个空间步长 $h$ 内传播所需的时间，否则差分依赖域无法覆盖物理依赖域，导致数值信息传递断裂与误差发散。

交错网格差分格式

网格布置与离散格式
将速度 $v$ 定义在整数网格点 $(lh, m\Delta t)$，应力 $\tau$ 定义在半网格点 $((l+1/2)h, (m+1/2)\Delta t)$，空间与时间均错开半步长。差分格式写为：

$$\frac{v_l^{m+1} - v_l^m}{\Delta t} = \frac{1}{\rho_l} \frac{\tau_{l+1/2}^{m+1/2} - \tau_{l-1/2}^{m+1/2}}{h}, \quad \frac{\tau_{l+1/2}^{m+1/2} - \tau_{l-1/2}^{m-1/2}}{\Delta t} = E_{l+1/2} \frac{v_{l+1}^m - v_l^m}{h}$$

稳定性与精度特征

1. 稳定性条件：经相同谐波分析可得 $\sin(\omega \Delta t / 2) = \pm \frac{c \Delta t}{h} \sin(kh / 2)$，稳定性判据与传统网格完全一致，仍为 $\Delta t \leq h / c_{\max}$。
2. 精度提升：空间导数采用半步长中心差分，等效空间步长减半。因截断误差量级为 $O(h^2)$，相同名义步长 $h$ 下，交错格式的近似误差降为传统网格的 $1/4$。
3. 物理一致性：变量交错定义更贴合波动方程的特征线传播结构，数值耗散与频散显著低于传统同点网格。

Von Neumann 稳定性分析

Von Neumann 方法是评估线性差分格式局部稳定性的标准工具。其应用需满足以下前提与假设：

1. 线性化与冻结系数假设：理论严格适用于常系数线性方程。对非均匀介质 $\rho(x), E(x)$，采用“冻结系数法”，即在局部网格邻域内将介质参数视为常数，以局部波速 $c_l$ 进行稳定性检验。
2. 周期性边界隐含假设：Fourier 级数展开隐含求解域空间周期性或无限大，忽略边界反射对局部稳定性的干扰。
3. 解的光滑性假设：要求真解充分可微，保证 Taylor 展开与谐波试探解的数学有效性。

分析步骤与判据
将离散解在节点 $(lh, m\Delta t)$ 展开为有限 Fourier 级数，考察单个谐波分量的演化。设向量 $\mathbf{U}^m = [v_l^m, \tau_l^m]^\mathrm{T}$，差分格式可写为矩阵递推形式：

$$\mathbf{U}^{m+1} = \mathbf{G}(k, \Delta t, h) \, \mathbf{U}^m$$

其中 $\mathbf{G}$ 为增长矩阵（Amplification Matrix），其元素由差分格式与波数 $k$ 决定。经过 $n$ 步时间推进，误差传递为 $\mathbf{U}^{n} = \mathbf{G}^n \mathbf{U}^0$。

稳定性充分必要条件：增长矩阵的所有特征值 $\Lambda$ 必须满足谱半径条件

$$\max |\Lambda(\mathbf{G})| \leq 1$$

• 若 $\max|\Lambda| > 1$ 对某些 $k$ 成立，格式不稳定；
• 若 $\max|\Lambda| \leq 1$ 且仅当 $\Delta t, h$ 满足特定比例时成立，为有条件稳定；
• 若对任意 $\Delta t, h$ 均满足，为无条件稳定。

以时间向前/空间中心格式为例，增长矩阵为 $\mathbf{G} = \begin{bmatrix} 1 & i\frac{c\Delta t}{h}\sin(kh) \\ i\frac{c\Delta t}{h}\sin(kh) & 1 \end{bmatrix}$，特征值 $\Lambda = 1 \pm i\frac{c\Delta t}{h}\sin(kh)$，模长恒大于 1，严格证明其无条件不稳定性。

数值频散与空间采样准则

数值相速度与群速度
即使满足 CFL 稳定性条件，有限差分离散仍会引入数值频散，即不同波数分量的传播速度偏离真实介质波速。由交错网格/时间中心格式的频散关系可导出归一化数值相速度 $c_{\text{num}}$ 与群速度 $c_g$：

$$c_{\text{num}} = \frac{\omega}{k} = \frac{2}{k\Delta t} \arcsin\left( \frac{c \Delta t}{h} \sin\frac{kh}{2} \right)$$

$$c_g = \frac{\partial \omega}{\partial k} = c \cdot \frac{\cos(kh/2)}{\sqrt{1 - \left( \frac{c\Delta t}{h} \sin(kh/2) \right)^2}}$$

当空间采样不足（$kh$ 较大）时，$c_{\text{num}}$ 与 $c_g$ 显著偏离 $c$，导致波形畸变与高频成分滞后/超前。

空间采样准则
为控制数值频散在可接受范围内，需保证每个波长内有足够的网格点数。对二阶空间精度格式，频散曲线分析表明：当 $h/\lambda < 0.1$ 时，数值相速度与群速度接近真实值。由此导出空间采样准则：

$$h < \frac{\lambda_{\min}}{10} = \frac{c_{\min}}{10 f_{\text{ac}}}$$

其中 $\lambda_{\min}$ 为需精确模拟的最小波长，$c_{\min}$ 为模型最小波速，$f_{\text{ac}}$ 为有效最高频率。若采用四阶或更高阶空间差分格式，该准则可放宽至每波长 5–6 个网格点。

工程实施建议
实际模拟中，时间步长 $\Delta t$ 需同时满足：

1. 稳定性约束：$\Delta t \leq h / c_{\max}$（全局最大波速控制）
2. 频散约束：$h < c_{\min} / (10 f_{\text{ac}})$（最小波长与最高频率控制） 二者需联合优化，在保证计算稳定的前提下兼顾精度与效率。

显式与隐式差分格式

有限差分格式按未知量在不同时间层上的依赖关系，可分为显式格式与隐式格式。显式格式中，下一时间层的解可由当前时间层已知值通过显式代数运算直接得出；隐式格式中，下一时间层的解隐含在耦合的代数方程组中，需联立求解。以下结合典型算例系统阐述两类格式的构造、稳定性特征及适用场景。

1. 常微分方程算例：牛顿冷却定律

考虑一阶常微分方程初值问题：

$$\frac{dT}{dt} = -\frac{T}{\tau}, \quad T(0) = T_0$$

其中 $\tau > 0$ 为冷却时间尺度，解析解为 $T(t) = T_0 e^{-t/\tau}$。

（1）显式格式（向前差分）
利用一阶向前差分离散时间导数，右侧函数值取在 $j$ 层：

$$\frac{T^{j+1} - T^j}{\Delta t} = -\frac{T^j}{\tau} \quad \Rightarrow \quad T^{j+1} = \left(1 - \frac{\Delta t}{\tau}\right) T^j$$

递推可得 $T^j = \left(1 - \frac{\Delta t}{\tau}\right)^j T_0$。

• 稳定性条件：要求增长因子模长 $\left|1 - \frac{\Delta t}{\tau}\right| \leq 1$，即 $0 \leq \Delta t \leq 2\tau$ 时格式稳定。
• 物理一致性：当 $0 \leq \Delta t < \tau$ 时解单调衰减，符合热力学规律；当 $\tau < \Delta t < 2\tau$ 时解振荡衰减；当 $\Delta t > 2\tau$ 时解振荡发散，格式失效。
• 结论：显式格式为有条件稳定，时间步长受物理参数严格限制。

（2）隐式格式（向后差分）
利用一阶向后差分离散时间导数，右侧函数值取在 $j+1$ 层：

$$\frac{T^{j+1} - T^j}{\Delta t} = -\frac{T^{j+1}}{\tau} \quad \Rightarrow \quad T^{j+1} = \frac{T^j}{1 + \Delta t/\tau}$$

• 稳定性：因 $\Delta t, \tau > 0$，分母恒大于 1，故 $|T^{j+1}| < |T^j|$ 对任意 $\Delta t$ 成立。格式为无条件稳定。
• 收敛性：当 $\Delta t \to 0$ 时，$(1 + \Delta t/\tau)^{-j} \to e^{-t/\tau}$，数值解收敛于解析解。
• 精度特征：与显式格式同为时间一阶精度。大步长下虽稳定，但截断误差累积可能导致数值解偏离物理衰减曲线。

2. 显隐混合格式：Crank-Nicolson 格式

为兼顾稳定性与精度，采用时间中点 $t_{j+1/2}$ 处的差分近似，将方程右侧在 $j$ 层与 $j+1$ 层取算术平均：

$$\frac{T^{j+1} - T^j}{\Delta t} = \frac{1}{2} \left[ f(T^j) + f(T^{j+1}) \right]$$

应用于牛顿冷却定律：

$$\frac{T^{j+1} - T^j}{\Delta t} = -\frac{1}{2\tau} (T^j + T^{j+1}) \quad \Rightarrow \quad T^{j+1} = T^j \frac{1 - \Delta t/(2\tau)}{1 + \Delta t/(2\tau)}$$

• 精度提升：Taylor 展开分析表明，该格式时间截断误差为 $O(\Delta t^2)$，精度高于纯显式/隐式格式的一阶精度。
• 稳定性：增长因子模长恒 $\leq 1$，对任意 $\Delta t > 0$ 均稳定。当 $\Delta t \leq 2\tau$ 时数值解单调衰减，符合物理预期。
• 计算特征：属隐式格式，需每步求解代数方程，但精度与稳定性综合表现最优。

3. 偏微分方程算例：扩散方程

考虑一维热传导/扩散方程：

$$\frac{\partial u}{\partial t} = \kappa \frac{\partial^2 u}{\partial x^2}$$

其中 $\kappa$ 为扩散系数。

（1）显式格式（FTCS）
时间向前差分，空间二阶中心差分：

$$u_l^{m+1} = u_l^m + s \left( u_{l+1}^m - 2u_l^m + u_{l-1}^m \right), \quad s = \frac{\kappa \Delta t}{h^2}$$

• Von Neumann 稳定性分析：代入试探解 $u_l^m = e^{i(klh - \omega m\Delta t)}$，得增长因子 $G = 1 - 4s \sin^2(kh/2)$。稳定性要求 $|G| \leq 1$，导出：

  $$s \leq \frac{1}{2} \quad \Rightarrow \quad \Delta t \leq \frac{h^2}{2\kappa}$$

• 局限性：时间步长与空间步长的平方成正比。为提高空间精度将 $h$ 减半，$\Delta t$ 需缩减至 $1/4$，计算步数剧增，实际应用中效率受限。

（2）Crank-Nicolson 格式
对空间二阶导数在 $m$ 与 $m+1$ 层取平均：

$$\frac{u_l^{m+1} - u_l^m}{\Delta t} = \frac{\kappa}{2h^2} \left[ (u_{l+1}^m - 2u_l^m + u_{l-1}^m) + (u_{l+1}^{m+1} - 2u_l^{m+1} + u_{l-1}^{m+1}) \right]$$

整理并引入 $s = \kappa \Delta t / h^2$，得隐式方程组：

$$-\frac{s}{2} u_{l-1}^{m+1} + (1+s) u_l^{m+1} - \frac{s}{2} u_{l+1}^{m+1} = \frac{s}{2} u_{l-1}^m + (1-s) u_l^m + \frac{s}{2} u_{l+1}^m$$

写为矩阵形式 $\mathbf{A} \mathbf{U}^{m+1} = \mathbf{B} \mathbf{U}^m$，其中 $\mathbf{A}, \mathbf{B}$ 为三对角矩阵。

• 稳定性：增长因子 $G = \frac{1 - 2s \sin^2(kh/2)}{1 + 2s \sin^2(kh/2)}$，恒满足 $|G| \leq 1$，格式无条件稳定。
• 求解策略：每时间步需求解三对角线性方程组，可采用 Thomas 算法（追赶法）高效求解，时间复杂度 $O(N)$。虽单步成本高于显式，但允许取更大 $\Delta t$，在精细网格下整体计算效率更优。

4. 椭圆型方程求解：松弛方法（Relaxation Method）

对于与时间无关的 Poisson 方程 $\nabla^2 \Phi = f$，采用空间中心差分离散：

$$\frac{\Phi_{i+1,j} - 2\Phi_{i,j} + \Phi_{i-1,j}}{h^2} + \frac{\Phi_{i,j+1} - 2\Phi_{i,j} + \Phi_{i,j-1}}{h^2} = f_{i,j}$$

整理得 Jacobi 迭代格式：

$$\Phi_{i,j}^{(m+1)} = \frac{1}{4} \left( \Phi_{i+1,j}^{(m)} + \Phi_{i-1,j}^{(m)} + \Phi_{i,j+1}^{(m)} + \Phi_{i,j-1}^{(m)} \right) - \frac{h^2}{4} f_{i,j}$$

该格式将节点值更新为四周邻点值的算术平均与源项修正之和。迭代从初始猜测（如全零场）开始，循环更新直至相邻迭代步的残差满足收敛准则：

$$\max_{i,j} \left| \Phi_{i,j}^{(m+1)} - \Phi_{i,j}^{(m)} \right| < \epsilon$$

该方法本质上是隐式求解椭圆边值问题的迭代实现，适用于稳态场模拟。

版权所有

版权归属：Staaaaaaaaar

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

编辑此页

贡献者: Staaaaaaaaar

下一页高阶差分算子