波动方程数值求解（二）

波动方程数值解是波动方程正演、逆时偏移和全波形反演的核心技术之一。本文采用高阶有限差分对波动方程进行了离散，进而实现了对波动方程的数值求解，模拟出其在介质中的传播过程。

一、声波波动方程高阶离散

由以下离散格式：

则二维波动方程高阶有限差分离散格式为：

其中am 为差分系数满足如下方程：

取
二、均匀介质二维声波方程求解
求取am 差分系数：

def high_order_finite_difference_coefficient(M, r):
    j0 = np.arange(1, M)  # 行
    m  = np.arange(1, M)  # 列
    j, m = np.meshgrid(j0, m, indexing='ij')  # indexing='ij'返回[行，列]， indexing='xy' 返回[列，行]
    A = m**(2*j) * (np.cos(np.pi / 8) ** (2*j) + np.sin(np.pi / 8) ** (2*j))
    b = r ** (2 * j0 -2)
    a_m = linalg.solve(A, b)
    a0 = - 2 * a_m.sum()
    return np.concatenate(([a0], a_m))

高阶离散迭代格式:

JJ = np.arange(M,Nz+1-M)
II = np.arange(M,Nx+1-M)
II,JJ = np.meshgrid(II,JJ)
M_sum = 0
for m in range(1, M):
    M_sum += a_M[m] * (u[t,II-m,JJ] + u[t,II+m,JJ] + u[t,II,JJ-m] + u[t,II,JJ+m])

u[t+1,II,JJ] = s_t[t]*h_x_z[II,JJ] + r ** 2 * (2 * a_M[0] * u[t,II,JJ] + M_sum) + 2 * u[t,II,JJ] - u[t-1,II,JJ]

求得结果与二阶结果对比如下图所示：快照时间分别为 0.4s 和 0.92s

[1]. Yang Liu, Mrinal K. Sen. A new time–space domain high-order fifinite-difference method for the acoustic wave equation. Journal of Computational Physics, 2009,228(23): 8779-8806.