The Crank-Nicolson method

The most popularly used numerical method of solving a stiff system of ODEs such as (11) is the Crank-Nicolson method, chosen because of its unconditional stability and good accuracy. First let us look at the Crank-Nicolson (trapezoidal) method for a single first-order ODE. Applying the Crank-Nicolson method to the model equation

$${{dy}\over {dt}}=\lambda y$$ (21)

yields

$${{y_{j+1}-y_{j}}\over {\Delta t}} = {\lambda \over 2} (y_{j+1} +y_{j})$$ (22)

or

$$y_{j+1} ={ {1+{{\lambda \Delta t}\over 2}} \over {1-{{\lambda \Delta t}\over 2}}}y_{j}$$ (23)

Thus the solution of (21) can be written as

$$y_{j}={\sigma}^{j}y_{0}$$ (24)

where y₀ is the initial condition and

$$\sigma = { {1+{{\lambda \Delta t}\over 2}} \over {1-{{\lambda \Delta t}\over 2}}}$$ (25)

For the migration wavefield extrapolation equation (19), the eigenvalue $\lambda=i\beta$ ($\beta \gt$) is purely imaginary. Thus

$$\sigma = \vert\sigma\vert e^{i\theta}$$ (26)

Because $\vert\sigma\vert=1$, the Crank-Nicolson method as applied to migration wavefield extrapolation does not generate any amplitude error. This result agrees with the one-way migration wavefield extrapolation that is a dampless sinusoidal wave propagation. However, there is phase error. The numerical solution is

$$\theta = 2 \arctan {{ \beta \Delta t }\over 2}$$ (27)

$$y_{j}={\sigma}^{j}y_{0}$$ (28)

$$y_{j}={\sigma}^{j}y_{0}= e^{i\theta j} y_{0}$$ (29)

But the exact solution of the ODE (21) is

$$y(x)= e^{i\beta t} y_{0} = e^{i\beta \Delta t j} y_{0}$$ (30)

The phase error is

$$PE=\beta \Delta t - 2\arctan {{ \beta \Delta t }\over 2}= {{(\beta \Delta t)}^3\over 12} + \ldots$$ (31)

The numerical solution lags behind the exact solution. Thus applying the Crank-Nicolson method to the diffraction equation causes overmigration - the experimental error. And the eigenvalues with larger moduli have larger degrees of overmigration.