摘要: 本贴从零开始学习正演的数值模拟方法.

1. 偏微分基础

引例: 物体从一维坐标的原点开始移动, 在 $t$ 时刻, 它在坐标轴的位置由函数 $s(t)$ 确定, 则速度为位置变化量与时间的比值:
$$v(t)=\frac{\mathrm{d}s(t)}{\mathrm{d}t}=\underset{\Delta t\to 0}{\lim }\frac{s(t+\Delta t)-s(t)}{\Delta t}\text{\hspace{1em}(1)}$$
加速度为速度变化量与时间的比值:
$$a(t)=\frac{\mathrm{d}v(t)}{\mathrm{d}t}=\underset{\Delta t\to 0}{\lim }\frac{v(t)-v(t-\Delta t)}{\Delta t}=\underset{\Delta t\to 0}{\lim }\frac{s(t+\Delta t)-2s(t)+s(t-\Delta t)}{\Delta t^2}\text{\hspace{1em}(2)}$$

推广 1: 给定一个单变量函数
$$y=f(x)\text{\hspace{1em}(3)}$$
其一阶导数记为
$$y^{\prime }=\frac{\mathrm{d}f(x)}{\mathrm{d}x}\text{\hspace{1em}(4)}$$
二阶导数记为
$$y^{\prime \prime }=\frac{{\mathrm{d}}^{2}f(x)}{\mathrm{d}{x}^2}\text{\hspace{1em}(5)}$$

推广 2: 给定一个二变量函数
$$z=f(x,y)\text{\hspace{1em}(6)}$$
其针对 $x$ 偏导的为
$$\frac{\partial z}{\partial x}=\underset{\Delta x\to 0}{\lim }\frac{f(x+\Delta x,y)-f(x,y)}{\Delta x}\text{\hspace{1em}(7)}$$
即 $x$ 发生了变化, 而 $y$ 并没变化. 二阶偏导为
$$\frac{{\partial }^{2}z}{\partial x^2}=\underset{\Delta x\to 0}{\lim }\frac{f(x+\Delta x,y)-2f(x,y)+f(x-\Delta x,y)}{\Delta x^2}\text{\hspace{1em}(8)}$$

另外有:
$$\frac{{\partial }^{2}z}{\partial x\partial y}=\underset{\Delta x\to 0,\Delta y\to 0}{\lim }\frac{f(x+\Delta x,y+\Delta y)-f(x,y+\Delta y)-f(x+\Delta x,y)+f(x,y)}{\Delta x\Delta y}\text{\hspace{1em}(9)}$$
$$\frac{{\partial }^{2}z}{\partial y\partial x}=\frac{{\partial }^{2}z}{\partial x\partial y}\text{\hspace{1em}(10)}$$
在进行数值模拟的时候, 不可能取 $\Delta x\to 0$, 因此 (8) 式简化为
$$\frac{{\partial }^{2}z}{\partial x^2}\approx \frac{f(x+\Delta x,y)-2f(x,y)+f(x-\Delta x,y)}{\Delta x^2}\text{\hspace{1em}(11)}$$
其中 $\Delta x$ 越小越准确, 但涉及的计算量越大, 我们只能取一个折中.

注 1: 为统一起见, 即使一元函数, 以后也常使用 $\partial$ 而不是 $\mathrm{d}$.

2. 波动方程

2.1 弦振动 (横波) 方程

参见全波形反演的深度学习方法: 第 2 章 正演, 根据牛顿第二定律
$$F=ma\text{\hspace{1em}(12)}$$
弦振动方程为
$$\frac{{\partial }^{2}u(x,t)}{\partial t^2}=c^2\frac{{\partial }^{2}u(x,t)}{\partial x^2}+f(x,t)\text{\hspace{1em}(13)}$$
其中 $c^2=T/\rho$, $f(x,t)=F(x,t)/\rho$, 左式的物理意义是瞬时加速度 $a$, 右式第一项的物理意义是 单位质量所受的力 $F$, $c$ 的物理意义是速度.

进一步忽略重力 $F(x,t)$ 的作用, 可以推出一维齐次波动方程的解:
$$\frac{{\partial }^{2}u(x,t)}{\partial x^2}=\frac{1}{c^2}\frac{{\partial }^{2}u(x,t)}{\partial t^2}\text{\hspace{1em}(14)}$$

2.2 声波 (纵波) 方程

声波仅有纵波. 考虑二维的情况, 它满足
$$\frac{1}{v^2}\frac{{\partial }^{2}U}{\partial t^2}=\frac{{\partial }^{2}U}{\partial x^2}+\frac{{\partial }^{2}U}{\partial z^2}\text{\hspace{1em}(15)}$$
其中 $U$ 指压力.

图 1 矩阵网格剖分

为了便于数值模拟, 将平面进行离散化, 仅考虑某些网格交叉点, 质量、压力等仅存在于这些点 (称为质点, 不知是否专业). 这样, 我们只考察第 $i$ 行第 $j$ 列的质点在时间 $k$ 的压力
$$U_{i,j}^k\text{\hspace{1em}(16)}$$
将 (11) 式按照变量名改造后代入 (15) 式可得
$$\frac{1}{v^2}\frac{U_{i,j}^{k+1}-2U_{i,j}^k+U_{i,j}^{k-1}}{\Delta t^2}=\frac{U_{i+1,j}^k-2U_{i,j}^k+U_{i-1,j}^k}{\Delta x^2}+\frac{U_{i,j+1}^k-2U_{i,j}^k+U_{i,j-1}^k}{\Delta y^2}\text{\hspace{1em}(17)}$$
其中 $k+1$ 表示下一个时间点, $i+1$ 表示下一个质点.