PDE

The acoustic wave equation for the square slowness m and a source q is given in 3D by :

\begin{cases} &mm \frac{d^2 u(x,t)}{dt^2} - \nabla^2 u(x,t) =q \\ &u(.,0) = 0 \\ &\frac{d u(x,t)}{dt}|_{t=0} = 0 \end{cases}

with the zero initial conditons to guaranty unicity of the solution

4th order discretizaition

We have

$ \frac{d^2 u(x,t)}{dt^2} \simeq \frac{u(x,t+dt) - 2 u(x,t) + u(x,t-dt)}{dt^2} - \frac{dt^2}{12}\frac{d^4u}{dt^4}+O(dt^4) $

and by differentiating the PDE we have

$ m \frac{d^4u}{dt^4} = \nabla^2 \frac{d^2 u(x,t)}{dt^2} + \frac{d^2 q}{dt^2}$

and finally by replacing the discrete forms in the original pde we have

$ m \frac{u(x,t+dt) - 2 u(x,t) + u(x,t-dt)}{dt^2} = \nabla^2u(x,t) + q(t) + \frac{dt^2}{12} \left(m^{-1} \nabla^2 (\nabla^2 u(x,t) + q(t) ) + \frac{d^2 q}{dt^2} \right)$

not that the term $\nabla^2 (\nabla^2 u(x,t))$ contains cross derivatives e.g in 2D :

$ \nabla^2 (\nabla^2 u(x,t)) = \frac{d^4 u(x,t)}{dx^4} + \frac{d^4 u(x,t)}{dz^4} + 2 \frac{d^4 u(x,t)}{dx^2dz^2} $

but it can be computed in cascade by first computing $\nabla^2 u(x,t) + q(t)$ then compute $u(x,t+dt)$