Following Steven G. Johnson's excellent PML (perfectly matched layer - an absorbing boundary condition) tutorial, I came-up with the following implementation of a PML for the 1D scalar wave equation using finite differences:
Scalar wave equation: \begin{equation} u_{xx}(x, t) = \frac{u_{tt}(x, t)}{c^2(x)} \tag{1} \end{equation} where subscripts indicate partial differentiation, so $u_{xx}$ is the second partial differential of $u$ in space.
Split (1) into: \begin{equation} \phi_t(x,t) = u_x(x,t)\tag{2a} \end{equation} \begin{equation} u_t(x,t) = c^2(x) \phi_x(x,t)\tag{2b} \end{equation}
Verification: \begin{equation} \mbox{$x$ partial derivative of (2a): } \phi_{xt}(x,t) = u_{xx}(x,t)\tag{3a} \end{equation} \begin{equation} \mbox{$t$ partial derivative of (2b): } u_{tt} = c^2(x) \phi_{xt}(x,t)\tag{3b} \end{equation} substitute (3a) in (3b): $u_{tt}(x,t) = c^2(x) u_{xx}(x,t)$, which is (1).
For the PML region, we use the wavefield at $x' = x + i \frac{f(x)}{\omega}$ (where $\omega$ is angular frequency), as this will exponentially damp the wavefield.
Making this change of coordinates, and with $\sigma(x) = f_x(x)$, $u_x(x,t)$ is replaced by $\frac{u_x(x',t)}{1 + i \frac{\sigma(x)}{\omega}}$.
Applying this to (2): \begin{equation} \phi_t(x',t) = \frac{u_x(x',t)}{1 + i \frac{\sigma(x)}{\omega}}\tag{4a} \end{equation} \begin{equation} u_t(x',t) = c^2(x)\frac{\phi_x(x',t)}{1 + i\frac{\sigma(x)}{\omega}}.\tag{4b} \end{equation} ($c(x')=c(x)$ and $\sigma(x')=\sigma(x)$)
Assuming that $u$ and $\phi$ can be written in the form $A e^{i(kx - \omega t)}$, $u_t = -i\omega u$. Applying this to (4) and multiplying both sides by the denominator of the right hand sides: \begin{equation} (1 + i \frac{\sigma(x)}{\omega}) (-i\omega) \phi(x',t) = u_x\tag{5a} \end{equation} \begin{equation} (1 + i\frac{\sigma(x)}{\omega}) (-i\omega) u = c^2(x)\phi_x(x',t)\tag{5a} \end{equation}
Multiplying, \begin{equation} (1 + i\frac{\sigma(x)}{\omega}) (-i\omega) = -i\omega + \sigma(x). \end{equation} Using this in (5), and noting that $-i\omega u = u_t$, \begin{equation} \phi_t(x',t) = u_x(x',t) - \sigma(x)\phi(x',t)\tag{6a} \end{equation} \begin{equation} u_t(x',t) = c^2(x)\phi_x(x',t) - \sigma(x)u(x',t)\tag{6b} \end{equation}
Taking the $x$ partial derivative of (6a), the $t$ partial derivative of (6b), and substituting, as in (3), we get \begin{equation} u_{tt}(x',t) = c^2(x)(u_{xx}(x',t) - (\sigma(x)\phi(x',t))_x) - \sigma(x)u_t(x',t)\tag{7} \end{equation}
Expanding the derivative $(\sigma(x)\phi(x',t))_x$ using the product rule: \begin{equation} u_{tt}(x',t) = c^2(x)(u_{xx}(x',t) - (\sigma_x(x)\phi(x',t) + \sigma(x)\phi_x(x',t))) - \sigma(x)u_t(x',t)\tag{8} \end{equation}
Using central differences for $u_{tt}$ and $u_t$: \begin{equation} \frac{u(x',t+1) - 2u(x',t) + u(x',t-1)}{\Delta t^2} = c^2(x)(u_{xx}(x',t) - (\sigma_x(x)\phi(x',t) + \sigma(x)\phi_x(x',t))) - \sigma(x)\frac{u(x',t+1) - u(x',t-1)}{2\Delta t} \end{equation}
Rearranging: \begin{equation} u(x',t+1) = \frac{c^2(x)\Delta t^2}{1 + \Delta t\frac{\sigma(x)}{2}}(u_{xx}(x',t) - (\sigma_x(x)\phi(x',t) + \sigma(x)\phi_x(x',t))) + \Delta t \sigma(x)\frac{u(x',t-1)}{2 + \Delta t\sigma(x)} + \frac{2u(x',t) - u(x',t-1)}{1 + \Delta t \frac{\sigma(x)}{2}}\tag{9} \end{equation}
Using a forward difference for $\phi_t$ in (6a): \begin{equation} \frac{\phi(x',t+1) - \phi(x',t)}{\Delta t} = u_x(x',t) - \sigma(x)\phi(t), \end{equation}
Rearranging: \begin{equation} \phi(x',t+1) = \Delta t u_x(x',t) + \phi(x',t) (1 - \Delta t \sigma(x))\tag{10} \end{equation}
(9) shows how to update $u$, and (10) shows how to update $\phi$.
This seems to work fairly well, but I noticed that the implementation in PySIT seems to be a bit different.
In my notation, I believe the PySIT version would be written: \begin{equation} u(x',t+1) = \frac{c^2(x)\Delta t^2}{1 + \Delta t\frac{\sigma(x)}{2}} (u_{xx}(x',t) + \phi_x(x',t)) + \Delta t\sigma(x)\frac{u(x',t-1)}{2 + \Delta t\sigma(x)} + \frac{2u(x',t) - u(x',t-1)}{1 + \Delta t \frac{\sigma(x)}{2}}\tag{11} \end{equation} \begin{equation} \phi(x',t+1) = -\Delta t \sigma(x) u_x(x',t) + \phi(x',t) (1 - \Delta t \sigma(x))\tag{12} \end{equation}
These are close to my own, but there are some discrepancies that trouble me:
- In (11), there is a $+ \phi_x(x',t)$ term instead of $-(\sigma_x(x)\phi(x',t) + \sigma(x)\phi_x(x',t))$ (compared to (9))
- In (12), $u_x$ is multiplied by $-\sigma$ (compared to (10))
The PySIT implementation seems to work as well, and has the advantages of $\phi$ being zero outside the PML region (where $\sigma=0$, so we can just use the regular wave equation there and avoid calculating $\phi$), and also not requiring $\sigma_x$.
How does my derivation need to be altered to arrive at the PySIT equations?