Matrix representation of implicit operators

The diffraction term of the 45$^\circ$ equation Claerbout (1985) can be rewritten as the following matrix equation, by inserting the rational part of the implicit extrapolator () into equation ():

$$\begin{eqnarray} \left({\bf I} + \alpha {\bf D} \right) {\bf q}_{z+\Delta z} & =... ...{z} \\ {\bf A}{\bf q}_{z+\Delta z} & = & {\bf A}^\ast {\bf q}_{z}\end{eqnarray}$$ (38) (39)

where the complex coefficient, $\alpha$ can be calculated, and ${\bf D}$ is a finite-difference representation of the Laplacian, $\nabla^2$.

The right-hand-side of equation () is known. The challenge is to find the vector ${\bf q}_{z+\Delta z}$ by inverting the matrix, ${\bf A}$.Given the wavefield on the surface, this equation provides a way to downward-continue in depth.

The matrices in equation () represent convolution with a scaled finite-difference Laplacian, with its main diagonal stabilized. In the two-dimensional problem, the $\nabla^2$ operator acts only in the x-direction, and can be represented by the three-point convolutional filter, d=(1,-2,1). The matrix, ${\bf A}$,therefore, has a tridiagonal structure, which can be inverted efficiently with a recursive solver.

In three-dimensional wavefield extrapolation, the $\nabla^2$ operator acts in both the x and y-directions. ${\bf A}$ therefore represents a 2-D convolution, and d can be represented by the a simple 5-point filter,

$$d = \left[ \begin{array} {ccc} & 1 & \\ 1 & -4 & 1\\ & 1 & \end{array} \right]$$ (40)

or a more isotropic 9-point filter Iserles (1996),  

$$d = \left[ \begin{array} {ccc} 1/6 & 2/3 & 1/6 \\ 2/3 & -10/3 & 2/3\\ 1/6 & 2/3 & 1/6\end{array} \right]$$ (41)

The vectors ${\bf q}_z$ and ${\bf q}_{z+\Delta z}$ contain the wavefield at every point in the (x,y)-plane. Therefore, the convolution matrices that operate on them are square with dimensions $N_x N_y \times N_x N_y$. As an illustration, for a $4 \times 2$ spatial plane, the structure of matrix ${\bf D}$ with the five-point approximation and transient boundary conditions, will be the blocked-tridiagonal matrix  

$${\bf D} = \left[ \begin{array} {cccc\vert cccc} -4 & 1 & . ... ... & 1 \\ . & . & . & 1 & . & . & 1 & -4 \\ \end{array}\right]$$ (42)

This blocked system cannot be easily inverted, even for the case of constant velocity, since the missing coefficients on the second diagonals break the Toeplitz structure.