Recursive solution in the time domain

Although the eigenvalue decomposition is unconditionally stable and precise, the N³-dependence of the computation time makes it prohibitively expensive. To make the method affordable, I have derived a scheme that uses a time propagator operator to find the solution at time $t + \Delta t$,using only the solution at time t. The basic principle is the same used by Dablain (1986) to derive a fourth order in time, explicit, finite difference scheme. Contrary to the method presented here, Dablain's scheme used two time-step solutions (${\bf u(t - \Delta t)}$and ${\bf u(t)}$) in order to find ${\bf u(t + \Delta t)}$.

A recursive solution in the time domain can be obtained with the use of the following relations:

$${d^2 {\bf u}(t) \over dt^2} = {\bf A} {\bf u}(t), \mbox{\hsp... ...over dt^4} = {\bf A}^2 {\bf u}(t), \mbox{\hspace{0.5cm}}\ldots,$$

in the Taylor expansion of ${\bf u}(t+dt)$ around ${\bf u(t)}$

$$\begin{eqnarray} {\bf u}(t+dt) & = & {\bf u}(t) + {\bf \dot{u}}(t) dt + {{\bf A}... ...{u}}(t) dt^2 + {{\bf A}^2 \over 6}{\bf u}(t) dt^3 + \ldots \;\; .\end{eqnarray}$$ (7)

Combining these two equations, we obtain the time-propagation equation

$$\begin{array} {cccc} \left( \begin{array} {c} {\bf u} \\ {\... ... {c} {\bf u} \\ \dot{\bf u} \end{array} \right)_t, \end{array}$$

where the forward time-propagation operator $\bf P_{+}$ has the following form:  

$${\bf P_{+}} = \left[ \begin{array} {cc} {\bf I} + {{\bf A} \... ...\right) & {\bf I} + {{\bf A} \over 2} dt^2 \end{array} \right]$$ (8)

for the case of a third-order approximation in time.

Figure  shows one frame of the horizontal displacement field, for the same source and media used to generate the synthetic data of Figure  (but with a finer grid spacing, and at a different time frame). This algorithm was implemented in the Connection Machine C2 with parellization in both spatial axes.

 

s2d2
Figure 7 One frame of the horizontal displacement wavefield, for a vertical source using the recursive time-domain elastic modeling scheme. The upper 3/4 of the model is isotropic, while the bottom part is anisotropic. The model size is 129 by 129 cells and the algorithm is fourth order in time. Press the button to see a movie.

Stability and energy conservation associated with the truncation of the Taylor expansion in equation (8) are discussed in Appendix B.