FINITE-DIFFERENCE SCHEME

Although sophisticated higher-order finite-difference schemes exist for the solution of the flux-conservative equation (Centrella and Wilson, 1984; Hawley et al., 1984), a first-order upwind finite-difference scheme is generally accurate enough for geophysical applications. I implemented a first-order method described by Engquist and Osher (1980).

Depending on the flow direction, upwind finite-difference schemes use a backward or forward finite-difference approximation to the derivative operator. The idea is that, although centered approximations are mathematically more accurate, the underlying physics may dictate otherwise. In other words, it is better to make a big error that is not transported, than to make a small one that blows up. The first-order forward and backward finite differences in space are, respectively:

$$\Delta_+\u\ =\ \u_{j+1} - \u_j,\ \ \ \ \ \ \ \ \Delta_-\u =\ \u_j - \u_{j-1},$$ (5)

where $\u_j$ is the discrete representation of $\u(x)$ on the grid $x_j = j\Delta x$.

Using the above finite-difference scheme, equation (3) is approximated by:  

$${{\u^{n+1}_j-\u^n_j} \over {\Delta z}} \ =\ \left\{ \begin{... ...\over {\Delta x}}\ & {\rm if}\ \u^n_j\gt 0, \end{array}\right.$$ (6)

with n the discrete-depth index of the grid: $z^n = n\Delta z$.However, the resulting scheme is unstable when $\u$ is discontinuous (when the field is going through a ``shock''). The discontinuity occurs at the source: going from left to right through the source, $\u$ changes sign from -s to +s, i.e., rays change direction from left to right.

Engquist and Osher solve the unstability problem in the shock region with the following scheme:

$$\u^{n+1}_j= \u^n_j\ -\ {{\Delta z}\over{\Delta x}}\ \Bigl(\Delta_+F_-(\u^n_j)\ +\ \Delta_-F_+(\u^n_j) \Bigr),$$ (7)

where  

$$F_+(u)\ =\ F({\rm max}(\u,{\overline{\u}})),\ \ \ \ \ \ \ F_-(u) \ =\ F({\rm min}(\u,{\overline{\u}})),$$ (8)

and ${\overline{\u}}$ is the stagnation point, i.e., $F'({\overline{\u}}) = 0$. ${\overline{\u}}= 0$ for the function F under consideration here. Engquist and Osher's scheme reduces to the standard first-order scheme (6) away from the source, and is nonlinearly stable near the source.