An insight into the Engquist-Osher scheme applied to Van Trier's algorithm

The practical implementation of the above formulation is centered around the problem of advancing the computational front. The Engquist-Osher scheme provides a way to compute $\Delta v$ by imposing a time minimization condition along three points of the finite-difference stencil. The scheme calculates ${ \Delta v} \over {\Delta \theta}$ to approximate the partial derivative (with respect to $\theta$), of the function $v(r,\theta)$, by using the values of $v(r,\theta)$in points of minimum traveltime; $v(r,\theta)$ varies only in $\theta$ along the constant radius computational front. Given three consecutive points on the computational front with constant radius, the values of the functions $u(r,\theta)$ and $v(r,\theta)$ vary only in the variable $\theta$.Each function will have the values u_j-1, u_j, u_j+1 and v_j-1, v_j and v_j+1 respectively, in the three points of the stencil. From equation (6) we can write $v(r,\theta)$ as a function of $u(r,\theta)$:

$$F(u)=v(r,\theta)=\sqrt{s^2(r,\theta)-{{u^2(r,\theta)} \over r^2}}.$$

The Engquist-Osher scheme computes $\Delta v$ as  

$${\Delta v} = F(min(u_{j+1},\bar u))-F(min(u_j,\bar u))+ F(max(u_j,\bar u))-F(max(u_{j-1},\bar u)),$$ (7)

where $\bar u =0$, at the point where ${\partial t \over \partial \theta} = 0$. For this case the value of $F(\bar u)=v=slowness$ from equation (6). Because the values of u_j-1, u_j and u_j+1 are compared against zero ($\bar u =0$), the scheme needs only the sign of the function u in the three points of the stencil. Equation (7) allows for eight cases as a function of the positive or negative values of u_j-1, u_j and u_j+1 in the three points of the stencil. The eight cases are shown in Figure 1.

 

Figure 1: Eight possible cases for the Engquist-Osher scheme. The thick vertical bars represent the sign of $u={\partial t \over \partial \theta}$, the continuous line represents the values of the traveltime over the three points of the stencil, the dashed vertical lines represent the location where ${\bar u}={\partial t \over \partial \theta}=0$; the dots on the time line show the coordinates for the function $v(r,\theta)$ which contribute to the difference $\Delta v$.On the right side are the values for the computed $\Delta v$ for each case.

The calculation of $\Delta v$ is done for locations where the value of the traveltime function $t(r,\theta)$, is minimum. From Figure 1 one can see that only cases 1 and 5 will actually give a first order correct value. In both cases the values for the function v are chosen from the three points of the stencil where the value of the time is minimum. While for cases 1 and 5 the accuracy of the scheme is unquestionable, for the rest of the cases some approximations are introduced. The other six cases also calculate the values of $\Delta v$ using the points where the value of the time is minimum. However, the value of $\Delta v$ is divided by a constant $\Delta \theta$,even though the function $v(r,\theta)$ is estimated over a different interval $\Delta \theta$.A potentially more accurate algorithm would calculate the exact value of $\Delta \theta$ for each intermediate case. The algorithm can be designed to calculate the locations of the minimum travel time in the three point stencil interval and the span over the $\theta$ axis, necessary to divide the value $\Delta v$.