波动方程的有限差分解（地震数据正演）

有限差分法

• 有限差分法是一种求解微分方程的数值方法，对于微分方程来说微分项很难处理，有限差分的思路便是用近似的方法处理微分项，例如对于一个一阶微分方程

  $$\begin{aligned} u^{'}(x) + c(x)u(x) &= f(x)\\ u(a) &= d \end{aligned}\tag{1}$$

  由导数定义有

  $$u^{'}(x) = \lim\limits_{\Delta x \rightarrow 0}\frac{u(x+\Delta x) - u(x)}{\Delta x} \approx \frac{u(x + \Delta x) - u(x)}{\Delta x} \tag{2}$$

  这里的几何意义便是使用割线斜率近似替代切线斜率，显然$\Delta x$越小则越精确，在后续求解中，我们将$\Delta x$记为步长。

  由此在求解时我们要对求解区域进行离散化，即将区间按步长分为一个个离散点。设对于区间$[a,b]$共分为$n$个区间，步长为$h$，坐标为${(x_0, u_0), (x_1, u_1),\dots, (x_n,u_n)}$，则我们可以得到各点处导数值得近似表示

  $$u_i^{'} = \frac{1}{h}(u_{i+1}- u_i)\tag{3}$$

  带入原方程得

  $$(u_{i+1} - u_i)/h + c(x_i)u_i = f(x_i)$$

  其中$x_i = a+ih, i=0,1,2,\dots$，则我们有$n$个方程$(i=0,1,\dots,n-1)$，包含$n$个未知数$(u_1,u_2,\dots,u_n)$。此时我们对方程变形得到

  $$u_{i+1} = (1-hc(x_i))u_i + hf(x_i)\tag{4}$$

  则上述方程组可写为矩阵的形式

  $$\begin{bmatrix} 1&&&\\ C_1&1&&\\ &\ddots&\ddots&\\ &&c_{n-1}&1 \end{bmatrix}\begin{pmatrix} u_1\\ u_2\\ \vdots\\ n_n \end{pmatrix} = \begin{pmatrix} hf(x_0)\\ hf(x_1)\\ \vdots\\ hf(x_{n-1}) \end{pmatrix} - \begin{pmatrix} C_0u_0\\ 0\\ \vdots\\ 0 \end{pmatrix}\tag{5}$$

  其中$C_i = hc(x_i)-1$，上述矩阵方程记为$\mathbf{A}\boldsymbol{u} = \mathbf{F}$，则解为$\boldsymbol{u} = \mathbf{A}^{-1}\mathbf{F}$（这里显然$r(\mathbf{A}) = n$即可逆）。