Most implicit schemes are based on a Crank-Nicholson formulation of the
extrapolation problem. This method calculates a centered finite difference
approximation to the *z*-derivative and averages the right hand side of the
equation over the two depths:

This problem is now in the form of a set of linear equations *Ax*=*y*. If
the second-order-in-space finite difference operator is used
to evaluate the Laplacian, the matrix *A* is pentadiagonal.
Direct methods, such as Gaussian elimination, may be used to solve the
system. For the particular structure of the system a
block-Gaussian-elimination scheme may be used. For an grid
of data in *x* and *y* the operation count for the direct solution is
*O*(*n ^{2}*). However, this method does not vectorize or parallelize well
because of recursions in the algorithm. Methods that are more suitable
for vectorization, such as block-cyclic-reduction, have a higher
operation count, . These costs are usually considered
too high for routine application in seismic data processing. If a more
accurate Laplacian operator is used the operator

Because of the high cost, this system is not usually solved in 3-D migration algorithms. Instead, a ``splitting'' approximation is made and the equations are rewritten in an approximate form that produces two sets of tridiagonal equations (Claerbout, 1985). However, the splitting approximation is not an isotropic approximation. If a time slice is made of the impulse response of the migration algorithm the wavefront is not circular. Waves propagating at to the computational grid travel slower than those along the axes of the grid. Hence, we would prefer to find a method that solves the pentadiagonal system.

12/18/1997