求解二维声波方程的高精度Runge-Kutta 方法

陈 丽，朱兴文，张朝元*

（1.大理大学工程学院，云南大理 671003；2.大理大学数学与计算机学院，云南大理 671003）

过去的地震波计算方法〔1-5〕因存在数值频散的不良缺点满足不了现在大尺度数值模拟的特性，杨顶辉教授团队〔6〕在2003 年将近似解析离散化（nearly analytic discrete，NAD）算子引入到模拟地震波中，已发展了系列空间离散精度为四阶的有关NAD 算子的数值算法〔7-10〕，这些算法能有效抑制数值频散现象的发生。

为了在地震波模拟时能更好地压制数值频散现象的发生，基于二维声波方程，选择八阶精度的NAD 算子〔11〕去离散常微分方程中的高阶空间偏导数，选择三阶Runge-Kutta 方法〔12〕去离散常微分方程中的时间导数，从而获得八阶NAD-RK 算法。针对该方法，分析了其理论误差，结合八阶Lax-Wendroff（LWC）格式和八阶交错网格（SG）格式两种高阶方法进行数值误差比较研究，详细推导其稳定性条件。

1 算法推导

设二维声波方程为：

此处，d 表示位移，c 表示二维声波速度。

再令U＝（D，P）T，则方程（2）式变为下式：

式（3）转化为标准的常微分方程模式，这样就可用常微分方程模式的求解程序来转化为对二维声波方程式（1）的对应式（3）的求解。过程如下：

第一，由于位移d 和粒子速度p 的二阶和三阶偏导数存在高阶偏微分算子矩阵L 中，先利用高阶偏导数的计算公式〔11〕对式（3）右边的L 进行空间高阶偏导数离散化。

第二，在离散化高阶空间偏导数后，式（3）成为一个关于时间t 的半离散化常微分方程模式，选用三阶Runge-Kutta 方法〔12〕进行求解。三阶Runge-Kutta 方法的计算公式为：

其中，Δt 为时间步长，U（1）、U（2）为中间变量，Un＝U（nΔt）。

消去式（4）中的变量U（1）、U（2），可以加快计算速度和减少存储量，得计算公式：

把U＝（D，P）T代入式（5）中，得到如下分量的计算公式：

这里A2＝A·A。

上面计算公式（5）或式（6a）和式（6b）称为求解二维声波方程式（1）的八阶NAD-RK 算法。

2 理论误差

对于求解二维声波方程的八阶NAD-RK 算法，由泰勒级数展开公式知式（5）右边的二阶和三阶空间偏导数计算公式的截断误差为O（Δx8+Δz8）〔11〕，即八阶NAD-RK 算法在空间上是一种八阶精度格式。采用三阶Runge-Kutta 方法求解常微分方程（5）式时，时间导数误差为O（Δt3）〔12〕。因此，求解二维声波方程的八阶NAD-RK 算法的理论误差为O（Δt3+Δx8+Δz8）。

3 数值误差

下面研究八阶NAD-RK 算法的数值误差。考虑以下二维初值问题：

其中，频率f0＝25 Hz，t=0 s 时刻声波的入射角度θ0=π/4。容易求得该初值问题的解析解为：

定义数值解的相对误差〔11〕为：

其中d（tn，xi，zj）为精确解，为数值解。

在进行数值计算时，计算参数选择为：计算域0＜x，z≤10 km、频率f0=20 Hz、波速c=4 km/s、网格尺寸Δx=Δz=50 m、时间步长Δt=0.001 s。图1 给出了式（9）中相对误差Er（%）的计算结果。图1 中的3条曲线分别对应八阶NAD-RK 算法、八阶LWC 格式和八阶SG 格式。从图1 可以看出，在3 种数值方法中，八阶NAD-RK 算法的数值误差最小。

图1 八阶NAD-RK 算法、八阶LWC 格式和八阶SG格式在求解初值问题（7）时在空间步长和时间步长分别为Δx=Δz=50 m 和Δt=0.001 s 情况下的相对误差

4 稳定性条件

为了保持数值迭代的稳定性，需要考虑如何选择合适的时间和空间步长，即Δt，Δx 和Δz。在条件h=Δx=Δz 下，所定义的库朗数α=cΔt/h 在数学上给出了速度c 与时间和空间两个网格步长h 之间的关系。需要确定α 的范围，使提出的方法稳定。在本部分，根据傅里叶分析和Yang 等〔13〕提出的分析过程，推导二维声波情况下八阶NAD-RK 算法的稳定性条件。

设

然后将谐波解：

带入式（5），结合式（6a）和（6b）两个高阶偏导数的离散公式，可得到如下方程：

其中H 为增长矩阵。

由于增长矩阵H 中部分元素hij的表达式太复杂，在这里只展示了H 的前两行元素的表达式，如下：

这里，αmax是最大库朗数。

5 结论

为有效抑制地震波数值模拟过程中的数值频散而获得更高的模拟精度和更快的计算速度，本文在杨顶辉教授团队的研究基础上，基于二维声波方程，结合八阶NAD 算子和三阶Runge-Kutta 方法，推导出八阶NAD-RK 算法。针对该方法，分析了其理论误差，结合八阶LWC 格式和八阶SG 格式两种高阶方法进行了数值误差比较研究，详细推导了其稳定性条件。研究结果揭示，八阶NAD-RK 算法的理论误差和稳定性条件分别为O（Δt3+Δx8+Δz8）和α≤αmax≈0.541 6；同八阶LWC 格式和八阶SG 格式相比，八阶NAD-RK 算法具有最小的数值误差。