Implicit finite difference in time-space domain with the helix transform

Application of spectral factorization and the helix operator

The linear system which must be solved according to Equation 9 has the form:

$$\left( \begin{array}{ccccc} U_0 & U_1 & 0 & 0 \\ U_1 & U_0 &... ..._1 \\ P^{t+1}_2 \\ P^{t+1}_3 \\ P^{t+1}_4 \end{array}\right) = \left( \begin{array}{ccccc} V_0 & V_1 & 0 & 0 \\ V_1 & V_0 &... ...rray}{c} P^t_1 \\ P^t_2 \\ P^t_3 \\ P^t_4 \end{array}\right)$$

$$+ \left( \begin{array}{ccccc} W_0 & W_1 & 0 & 0 \\ W_1 & W_0 &... ...1 \\ P^{t-1}_2 \\ P^{t-1}_3 \\ P^{t-1}_4 \end{array}\right).$$ (10)

For simplicity, we can combine all the constants into one: $\alpha = \frac{C^2 \Delta t^2}{3 \Delta x^2}$ . The matrix coefficients in Equation 10 (the finite-difference weights) are then:
$U_0 = 1 + 2 \alpha, U_1 = -\alpha;$
$V_0 = 2 - 2 \alpha, V_1 = \alpha;$
$W_0 = -1 - 2 \alpha, W_1 = \alpha.$

In shorter notation, Equation 10 reads:

$$U P^{t+1} = V P^t + W P^{t-1}.$$ (11)

The solution of this linear system is:

$$P^{t+1} = U^{-1} \left(V P^t + W P^{t-1}\right).$$ (12)

To solve this system, we must perform polynomial division. The system is tridiagonal (and easily solvable) only for 1 dimension. For multiple dimensions, matrix $U$ is block diagonal. Additional non-zero elements appear at a certain offset from the diagonal, making the solution process more complicated. However, using spectral factorization, the finite-difference weights of matrix $U$ can be factorized into a set of causal filter coefficients $u$ and it's time reverse $u^T$ . Using the helical approach to deconvolution, the system can be recast as:

$$P^{t+1} = (u^T u)^{-1} \left(V P^t + W P^{t-1}\right);$$ (13)

$$P^{t+1} = u^{-1} (u^T)^{-1} \left(V P^t + W P^{t-1}\right).$$ (14)

As stated above, polynomial division is equal to deconvolution. This means that the polynomial division in Equation 12 can be achieved by a set of two deconvolutions of the spectrally factorized coefficients $u$ of matrix $U$ . One deconvolution is done along the data in the reverse direction (application of the adjoint of the filter):

$$y_k = x_k - \sum_{i=1}^{N_a} a_i' y_{k-i},$$ (15)

where $a'$ is the time reversed filter coefficients of $u$ . The other deconvolution is done in the forward direction:

$$x_k = y_k - \sum_{i=1}^{N_a} a_i x_{k-i}.$$ (16)

I used the SEPlib module polydiv, which uses the helical coordinates to do the deconvolutions (the polynomial division) in equations 15 and 16.
The wavefield propogation is done by the following sequence:

1. Spectrally factorize the coefficients of matrix $U$ .
2. Multiply the saved wavefield at time $t-1$ by the coefficients of matrix $W$ .
3. Multiply the saved wavefield at time $t$ by the coefficients of matrix $V$ .
4. Sum the results of the previous 2 steps into a result vector.
5. Uncoil the factorized coefficients $u$ over the result vector (eq. 15).
6. Coil the factorized coefficients $u$ over the result vector (eq. 16).

This sequence is repeated for each time step.

Implicit finite difference in time-space domain with the helix transform