三维波传播模拟的多项式特解法及其物理应用

Polynomial Special Solution Method for Three-Dimensional Wave Propagation Simulation and Its Physical Applications

摘要: 本文提出一种基于多项式特解的三维波动方程高精度数值方法，探究在不同传播速率下三维波动方程的数值模拟问题。该方法使用满足波动方程的不同阶的多项式特解作为基函数对波动方程进行数值模拟，本文选取传播速率存在差异的数值算例展开求解并与有限差分法进行对比，结果表明：该算法比传统有限差分法精度更高，易于实现，两个算法的计算误差随波速减小(对应物理波长变短)而增大的规律与波动理论严格一致，验证了模型的物理正确性。本研究可为声波勘探、医学超声、工程波动模拟等领域，提供一种高效高精度的数值工具。

Abstract: This paper proposes a high-accuracy numerical method for the three-dimensional wave equation based on polynomial particular solutions, and investigates the numerical simulation of the three-dimensional wave equation under different propagation velocities. The method employs polynomial particular solutions of various orders that satisfy the wave equation as basis functions to numerically simulate the wave equation. Numerical examples with different propagation velocities are selected for solution and comparison with the traditional finite difference method (FDM). The results show that this algorithm has higher accuracy and easier implementation than the traditional FDM. The rule that the calculation errors of both algorithms increase with the decrease of wave velocity (corresponding to the shortening of physical wavelength) is strictly consistent with the wave theory, which verifies the physical correctness of the model. This work provides an efficient and high-precision numerical tool for applications in seismic exploration, medical ultrasound, and general engineering wave propagation problems.

1. 引言

在物理的领域中，能够对波传播现象的精确模拟在理解声波、弹性波、电磁波等方面具有非常重要的地位。本文针对三维波动方程，提出一种基于多项式特解法的数值模拟算法，重点研究波速变化对波前演化、能量分布及干涉模式的影响。该方法通过构造满足波动方程的不同阶多项式特解作为基函数，实现了对复杂传播速率情形下的高效、稳定模拟。本文验证该方法在模拟波传播、散射及频散问题中的有效性与适用性，为声学设计、地震波分析、光学波导等应用提供数值工具支持。波动方程是一类非常重要的偏微分方程。目前偏微分方程数值解法主要可分为网格类方法与无网格方法两大类。在网格类方法中，有限差分法[1]和有限元法[2]应用广泛；无网格类方法则包括Kansa方法[3]、径向基函数配点法[4]等。其中，基于径向基函数的无网格方法因其实现便捷、精度较高等优势，在各类偏微分方程求解中发展迅速[5]。特解法[6]是在Kansa方法基础上改进而来，研究表明其在数值稳定性方面更具优势[7]。在传统的网格类数值方法中，求解波动方程时我们需要对时间离散格式进行特殊限制以确保稳定性：若采用显式时间格式，则时间步长与空间步长必须满足特定的网格比例条件；若采用隐式格式，则需要通过迭代计算求解代数系统。本文采用满足目标波动方程的多项式特解作为基函数进行数值求解，该方法避免了传统无网格方法中支撑域选取的问题，且所得数值解的稳定性不受时间与空间步长比例关系的限制。为改善系数矩阵的病态性问题，本文所有数值实验均引入多尺度技巧[8]，以确保随着多项式阶数增加，由特解基函数构成的系统矩阵条件数仍处于合理范围。

三维波动方程可描述声波、弹性波或电磁波在均匀介质中的传播过程，本文以如下三维波动方程为例：

$\{\begin{array}{ll}{u}_{tt}\left(\text{x},t\right)-a^2\Delta u\left(\text{x},t\right)=f\left(\text{x},t\right)\hfill & \left(\text{x},t\right)\in \Omega \times \left(0,T\right]\hfill \\ u\left(\text{x},0\right)=u_0\left(\text{x}\right)\hfill & \text{x}\in \Omega \hfill \\ u_t\left(\text{x},0\right)=u_1\left(\text{x}\right)\hfill & \text{x}\in \Omega \hfill \\ u\left(\text{x},t\right)=v\left(\text{x},t\right)\hfill & \left(\text{x},t\right)\in \partial \Omega \times \left(0,T\right]\hfill \end{array},$ (1)

其中$u_0\left(\text{x}\right),u_1\left(\text{x}\right),v\left(\text{x},t\right)$ 为给定函数，$f\left(\text{x},t\right)$ 在物理上可表示为：声学问题中的体积声源分布、弹性波问题中的体力密度、电磁学问题中的电流源或电荷源，$\Delta$ 是拉普拉斯算子，空间求解区域$\Omega$ 是$R^3$ 中的有界连通域，$a$ 为波的传播速率(简称波速)，本文中$T$ 取1。由于传播速率对波动方程具有重要的物理意义，本文将探讨新提出的算法在求解波动方程时不同的传播速率对算法精度和稳定性的影响。

2. 方法论述

多项式特解法

以三维空间波动方程为例，阐述多项式特解法的求解过程。首先，设波动方程(1)的解$u\left(x,y,z,t\right)$ 具有以下形式：

$u\left(x,y,z,t\right)\approx \widehat{u}\left(x,y,z,t\right)={\displaystyle \sum _{i=0}^s{\displaystyle \sum _{j=0}^i{\displaystyle \sum _{k=0}^j{\displaystyle \sum _{l=0}^k{a}_{ijkl}{\Phi }_{ijkl}}}}}\left(x,y,z,t\right),$ (2)

其中，基函数${\Phi }_{ijkl}\left(x,y,z,t\right)$ 为满足

$L{\Phi }_{ijkl}\left(x,y,z,t\right)=x^{i-j-k-l}{y}^j{z}^k{t}^l$ 的多项式特解，微分算子

$L={\partial }_{tt}-a^2\left({\partial }_{xx}+{\partial }_{yy}+{\partial }_{zz}\right)+I,$

其中$I$ 是单位算子。在基于多项式特解的数值方法中，所采用的基函数${\Phi }_{ijkl}\left(x,y,z,t\right)$ 为满足$L{\Phi }_{ijkl}\left(x,y,z,t\right)=x^{i-j-k-l}{y}^j{z}^k{t}^l$ 的解析多项式特解。这些基函数天然满足所求解的偏微分方程，从而保证了数值格式的内在一致性。

本文选取配点形式如下：内配点集合为${\left\{\left({\text{x}}_i,t_i\right)\right\}}_{i=1}^{{}_{n_i}}\subset \Omega \times \left(0,T\right)$ ，边界配点集合为

${\left\{\left({\text{x}}_b,t_b\right)\right\}}_{b=n_i+1}^N\subset \partial \Omega \times \left(0,T\right)\cup \Omega \times \left\{t=0\right\}.$

则方程(2)的待定系数$\left\{a_{ijkl}\right\}$ 可由下式得到

$\{\begin{array}{l}{\displaystyle \sum _{i=0}^s{\displaystyle \sum _{j=0}^i{\displaystyle \sum _{k=0}^j{\displaystyle \sum _{l=0}^k{a}_{ijkl}{x}_r^{i-j-k-l}{y}_r^j{z}_r^k{t}_r^l}}}}-{\displaystyle \sum _{i=0}^s{\displaystyle \sum _{j=0}^i{\displaystyle \sum _{k=0}^j{\displaystyle \sum _{l=0}^k{a}_{ijkl}{\Phi }_{ijkl}\left(x_r,y_r,z_r,t_r\right)}}}}=f\left(x_r,y_r,z_r,t_r\right),r=1,2,\cdots ,n_i\\ {\displaystyle \sum _{i=0}^s{\displaystyle \sum _{j=0}^i{\displaystyle \sum _{k=0}^j{\displaystyle \sum _{l=0}^k{a}_{ijkl}{\Phi }_{ijkl}\left(x_r,y_r,z_r,0\right)}}}}=u_0\left(x_r,y_r,z_r\right)\\ {\displaystyle \sum _{i=0}^s{\displaystyle \sum _{j=0}^i{\displaystyle \sum _{k=0}^j{\displaystyle \sum _{l=0}^k{a}_{ijkl}{\left(\partial {\Phi }_{ijkl}\right)}_t\left(x_r,y_r,z_r,0\right)}}}}=u_1\left(x_r,y_r,z_r\right)\\ {\displaystyle \sum _{i=0}^s{\displaystyle \sum _{j=0}^i{\displaystyle \sum _{k=0}^j{\displaystyle \sum _{l=0}^k{a}_{ijkl}{\Phi }_{ijkl}\left(x_r,y_r,z_r,t_r\right)}}}}=g\left(x_r,y_r,z_r,t_r\right),r=n_i+1,\cdots ,N\end{array},$ (3)

其中$r>\left(s+1\right)\left(s+2\right)\left(s+3\right)\left(s+4\right)/24$ ，$s$ 为多项式特解基函数的阶，将求出的待定系数$\left\{a_{ijkl}\right\}$ 带入式(2)，得到方程(1)的近似解。对于三维空间中二阶常系数偏微分方程，其多项式特解${\Phi }_{ijkl}\left(x,y,z,t\right)$ 的一般表达式已由引理1给出。类似地，该结论可推广至一维与二维情形，相应空间中同类型方程的多项式特解表达式可通过类似推导得到。

引理1 [8] 令$\left(w_1,w_2,w_3,w_4\right)$ 表示$\left(x,y,z,t\right)$ ，考虑下列偏微分方程：

${\displaystyle \sum _{i,j=1}^4{a}_{ij}\frac{{\partial }^2{u}_p}{\partial w_i\partial w_j}}+{\displaystyle \sum _{i=1}^4{b}_i\frac{\partial u_p}{\partial w_i}}+c{u}_p=w_1^i{w}_2^j{w}_3^k{w}_4^l,$ (4)

其中${\left\{a_{ij}\right\}}_{i,j=1}^4$ ，${\left\{b_i\right\}}_{i=1}^4$ 和$c$ 为实常数， $i,j,k,l$ 为非负整数，且$c\ne 0$ ，则(4)的多项式特解由下式给出：

$u_p=\frac{1}{c}{\displaystyle \sum _{r=0}^d{\left(\frac{-1}{c}\right)}^r}{L}^r\left(w_1^i{w}_2^j{w}_3^k{w}_4^l\right).$ (5)

其中$d=i+j+k+l,L={\displaystyle \sum _{i,j=1}^4{a}_{ij}\frac{{\partial }^2}{\partial w_i\partial w_j}}+{\displaystyle \sum _{i=1}^4{b}_i\frac{\partial }{\partial w_i}}$ 。

本文结合波动传播的物理背景，采用均方根误差(RMSE)对算法精度进行量化评估。该指标能够直接反映数值解与真实波场在能量意义上的整体偏差即RMSE越小，表明算法重构的波前相位、振幅及能量分布与物理实际越吻合。这一评估方式既能从统计角度衡量全局精度，也符合波动问题中能量表征的核心物理思想。本文通过计算均方根误差(RMSE)来检验该算法的精度：

$\text{RMSE}=\sqrt{\frac{1}{n_t}{\displaystyle \sum _{j=1}^{n_t}{\left|{\widehat{u}}_j-u_j\right|}^2}},$ (6)

其中$n_t$ 为检验点数目，${\widehat{u}}_j$ 和$u_j$ 分别代表在第$j$ 个检验点处得到的数值解与对应的解析解。

本文在三维波动方程的数值求解中，分别考虑了规则与不规则计算区域，并相应设置了配点集。采用随机方式生成配点，规则区域内随机配点集合为：

${\left\{\left(x_i,y_i,z_i,t_i\right)\right\}}_{i=1}^{n_i}\subset \Omega \times \left(0,T\right)$ ，$n_i=6561$ 。

其中配点均匀分布于整个求解域。对于不规则区域，配点根据区域几何特征以类似原则进行选取，以保证计算域内的点集覆盖。随机边界配点集合选取为

${\left\{\left(x_b,y_b,z_b,t_b\right)\right\}}_{b=1}^{n_b}\subset \Omega \times \left\{t=0\right\}\cup \partial \Omega \times \left(0,T\right)$ ，$n_b=6029$ ，如图1(a)所示。

在不规则空间区域内，随机内配点的集合可表示为：

${\left\{\left(x_i,y_i,z_i,t_i\right)\right\}}_{i=1}^{n_i}\subset \Omega \times \left(0,T\right)$ ，$n_i=7038$ 。

随机边界配点集合取为${\left\{\left(x_b,y_b,z_b,t_b\right)\right\}}_{b=1}^{n_b}\subset \Omega \times \left\{t=0\right\}\cup \partial \Omega \times \left(0,T\right)$ ，$n_b=4530$ ，三维波动方程数值模拟所选取的不规则空间计算域如图1(b)所示。

(a) Collocation points in regular domain (b) Collocation points in irregular domain

Figure 1. Distribution of collocation points in regular and irregular 3d computational domains

图1. 三维规则和不规则计算区域配点图

3. 数值算例

本研究选取波动方程的精确解析解为：$u\left(x,y,z,t\right)=\frac{1}{a}{x}^3{y}^3{t}^3$ ，相应的初值条件为：$u_0\left(x,y,z\right)=u_1\left(x,y,z\right)=0$ 。空间计算区域分别采用图1(a)所示的规则域与图1(b)所示的不规则域进行模拟，相关数值结果汇总于表1与表2，矩阵条件数随多项式特解基函数阶数的变化曲线见图2。

Table 1. RMSEs of example in regular space domain with different a

表1. 规则空间区域内方程在不同a作用下的均方根误差

a	10^−^3	10^−^2	10^−^1	10^0
2	1.13e+01	1.13e+00	1.13e−01	1.14e−02
4	6.23e+00	6.23e−01	6.31e−02	3.33e−03
6	3.94e+00	3.94e−01	3.98e−02	8.74e−04
8	1.01e−08	7.88e−12	2.64e−14	4.68e−13
10	2.48e−07	3.21e−10	5.39e−13	8.77e−13
12	2.07e−06	6.93e−09	3.26e−12	3.97e−12
14	2.20e−05	4.25e−06	2.88e−10	5.92e−12
16	4.11e−05	6.30e−06	1.54e−07	6.98e−12
18	1.99e−05	2.03e−06	1.19e−08	6.57e−12
20	7.98e−05	3.52e−06	3.62e−07	6.74e−12

Table 2. RMSEs of example in irregular space domain with different a

表2. 非规则空间区域内方程在不同a作用下的均方根误差

a	10^−3	10^−2	10^−1	10^0
2	1.81e+01	1.81e+00	1.81e−01	1.79e−02
4	4.13e+00	4.13e−01	4.17e−02	4.22e−03
6	3.00e−01	3.00e−02	4.75e−03	3.01e−04
8	7.71e−14	1.18e−14	1.29e−15	4.63e−13
10	1.60e−13	1.90e−14	1.42e−15	3.42e−13
12	1.17e−10	4.40e−13	3.52e−15	1.05e−12
14	2.50e−09	5.41e−12	2.10e−14	3.66e−13
16	4.78e−08	2.20e−11	1.71e−13	3.53e−13
18	1.83e−07	1.45e−08	6.96e−13	4.22e−13
20	5.46e−07	1.92e−08	7.28e−12	3.59e−13

Figure 2. Condition number vs. order of polynomial particular solution basis functions

图2. 矩阵条件数随多项式特解基函数阶数的变化曲线图

计算结果表明，采用多项式特解法得到的数值解在规则空间区域和非规则空间区域的精度的变化规律大致相同，随着时空多项式特解的阶的增加，计算误差迅速减小，数值解达到收敛状态。继续增加时空多项式特解的阶，计算误差没有明显变化，略有波动，其原因是达到一定精度，所用配置点个数需要大于$\left(s+1\right)\left(s+2\right)\left(s+3\right)\left(s+4\right)/24$ ，其中$s$ 是多项式特解基函数的阶，相应系数矩阵的阶呈几何级数增长，影响了数值解的精度。

本文通过有限差分法离散化三维波动方程的偏导数，将连续的偏微分方程转化为离散的差分方程，再通过时间步进迭代求解空间各点的波场值。采用二阶中心差分格式(空间方向拉普拉斯算子通过二阶中心差分离散，时间方向采用二阶时间步进格式)，网格大小设置为$Nx\times Ny\times Nz=20\times 20\times 20$ (空间区域为${\left[0,1\right]}^3$ )，对应$x,y,z$ 方向网格步长$dx=dy=dz=1/\left(20-1\right)\approx 0.0526$ ，时间步长$dt=0.001$ ，总模拟时间$T=0.5$ ，波速c与待测试参数a一致(a取值分别为1、0.1、0.01、0.001)，且计算过程中严格验证CFL条件(确保CFL数 < 1，满足算法稳定性要求)。计算的均方根误差结果如表3所示，计算的时间与误差和自由度与误差的收敛曲线如图3所示，部分时刻波场三维等值面图如图4所示。

Table 3. RMSEs of example 3 in irregular space domain with different a

表3. 有限差分法在不同a作用下的均方根误差

a	10^−3	10^−2	10^−1	10^0
RMSE	3.78e−02	3.00e−01	3.00e+00	3.00e+01

Figure 3. Degrees of freedom, computational time and error convergence plot

图3. 自由度和计算时间与误差的收敛曲线图

Figure 4. 3D Isosurface plot of wavefield at t = 0.3, t = 0.4 time instant

图4. t = 0.3、t = 0.4时刻波场三维等值面图

结果表明：多项式特解法比通过泰勒展开离散的有限差分法具有更高的格式精度，其数值误差较传统有限差分法降低了很多。这种精度优势在物理上具有深刻含义。波动方程的数值误差主要源于对短波长分量的色散效应。根据物理关系$a=f\lambda$ 其中$f$ 为频率，$\lambda$ 为波长，在固定频率下，介质波速的降低直接导致物理波长变短。因此，“波速减小，数值误差增加”这一现象，并非单纯的计算缺陷，而是数值方法分辨率逼近其物理极限的直接体现，也就是说任何离散方法在模拟波长接近甚至小于网格尺度的波动时，精度都会下降。本研究的结果严格遵循了这一物理规律，从而验证了数值模型的可靠性。更重要的是，多项式特解法因其高阶特性，在相同的离散尺度下，能比低阶有限差分法更有效地解析更短的物理波长。这意味着，对于低速介质(如地下低速带、生物软组织)或高频激励产生的精细波场结构，本方法能以更低的计算成本实现所需的模拟精度。

4. 结论

本文提出采用多项式特解法求解具有不同传播速率的波动方程，该算法使用满足波动方程的多项式特解作为基函数，对转化后的问题进行数值模拟，该算法精度和物理一致性上均展现出显著优势。该算法模拟波动方程局部区域配置点的数值解，具有较高的精度，可应用于研究电磁波在波导中传播的局部性质，尤其是应用相变材料实现波导的光透率的变化，来实现光通讯与光计算；可应用于晶格热传导的局部解以及量子力学绝大多数波动方程的数值解问题；可应用于地震波数值模拟在地震勘探方法研究、地震观测系统优化设计、地震资料解释等方面；可应用于声场问题在数学上的处理，如研究井中各类模式声波的激发和传播规律以及其局部性质。

数值实验结果表明，本方法的计算复杂度随介质波速的降低而增加，与有限差分方法结果类似，这本质上反映了短波长物理现象对数值分辨率的更高要求。该特性与物理理论一致，并凸显了本方法在模拟地下低速异常体、生物软组织声学等具有精细结构或短波长特征的波动问题中的应用潜力。该方法在长时间模拟中具有更高精度与更好稳定性，累积误差显著低于传统时间步进法，但全域离散会增大系统规模，导致内存消耗更高、计算成本更大。与有限差分法相比，该方法在长时间演化中可有效抑制数值色散与漂移，整体精度优势明显；而时间步进法内存占用更低、求解更高效，更适用于大规模短时问题。综合来看，该方法在精度与长时间稳定性上占优，时间步进法在内存与计算效率上更具优势，实际应用可根据模拟时长、精度需求与硬件条件灵活选用。该算例仅讨论方法的基础收敛特性。针对非多项式波形的逼近能力与数值色散特性，后续将重点开展拓展研究：拟选取高斯脉冲传播、Ricker子波辐射等经典波动基准测试算例，验证方法对工程化非多项式波形的适配性，量化分析基函数阶数对脉冲展宽、子波畸变的控制效果，评估数值色散系数在不同波速下的变化规律，为方法的工程化应用奠定基础。