Backward Euler scheme or Crank-Nicolson method?

[pazdzioch]

Dear Granville,

in order to find a time step suitable for my computations I decided to observe the fastest process appearing in the model I use, namely Alfven waves in plasma physics. The equations, after some derivation may be rewritten in a form of one, wave equation with an analytical solution - simple oscillation has been obtained. The problem now is in the case of numerical computation the amplitude of the oscillation is decaying. I have found out that it happens mainly because of backward Euler scheme. In the case of Crank-Nicolson method the oscillation does not decay. Could you explain it shortly? The most important question is related to a more complex system - is this fast oscillation process damped always then when using the Euler scheme? Is there a difference in efficiency between both methods?

Thank you very much,
Dominik

[sewell]

If you use finite differences or finite elements to solve a problem like
Ut = Uxx (+ Uyy + Uzz)
you will end up having to solve an ODE system of the form:
Ut = AU
where A is a matrix, with negative eigenvalues, some of which are very large and negative. If A = S^{-1}DS, this can be written:
Vt = DV
where V=SU, and D is the diagonal matrix with the eigenvalues of A. Now the ODEs separate out into scalar DEs of the form:
(Vi)t = -lambda*Vi
so let's just look at a simple scalar DE of the form:
Vt = -lambda*V
but remember that some of the lambdas will be very large (the ones which correspond to highly oscillatory (in space) eigenfunctions.)

The backward Euler method is:
V_{n+1} - V_n = -lambda*dt*V_{n+1}
or
V_{n+1) = 1/(1+lambda*dt)*V_n (1)

The Crank-Nicolson method is:
V_{n+1} - V_n = -lambda*dt*(V_n + V_{n+1})/2
or
V_{n+1} = [1 - lambda*dt/2]/[1 + lambda*dt/2]*V_n (2)

You can show that in both cases, as dt --> 0, the solution converges to the correct solution, V_n = V_0*exp(-lambda*t_n).

However, notice that when dt is not yet too small, and lambda is large, corrresponding to large negative eigenvalues of the original system Ut = AU), the corresponding eigenvector is damped out rapidly by the backward Euler method (1) (the factor in front of V_n is small), while the Crank-Nicolson method (2) does not damp it out rapidly (the factor is nearly equal to -1).

This well-know property of Crank-Nicolson, that it does not damp out high frequency components much, is sometimes an advantage (for wave questions such as yours, where you don't want damping), but often is a disadvantage (high frequency "noise" in the initial conditions, if present, just oscillates, see problem 8 of Appendix A of my book).

Granville Sewell

[pazdzioch]

Now it is clear:D

Thanks!:)