# On the Stability of the Iterated Crank-Nicholson Method in Numerical Relativity

## I Introduction

There is currently a large worldwide effort underway attempting to solve Einstein’s equations numerically for astrophysically interesting scenarios. The problem is extremely challenging technically. Among the difficulties is that of finding a finite-difference scheme that allows a stable time evolution of the system. It is well-known that implicit differencing schemes tend to be stable. However, the difficulty of solving the resulting implicit algebraic equations, especially in three spatial dimensions, has led most researchers to stay with explicit methods and their potential instabilities.

Several years ago, Choptuik proposed solving the implicit Crank-Nicholson scheme by iteration. This would effectively turn it into an explicit scheme, but hopefully by iterating until some convergence criterion was met one would preserve the good stability properties of Crank-Nicholson. The iterated Crank-Nicholson scheme has subsequently become one of the standard methods used in numerical relativity.

In this note, we point out that when using iterated Crank-Nicholson, one should do exactly two iterations and no more. While the limit of an infinite number of iterations is the implicit Crank-Nicholson method, it can in fact be worse to do more than two iterations, and it never helps.

## II Iterated Crank-Nicholson

To understand this paradoxical result, consider differencing the simple advective equation

$$\frac{u}{t}=\frac{u}{x}.$$
(1)   
(Many equations in numerical relativity are generalizations of this form, and the differencing techniques are similar.) A simple first-order accurate differencing scheme is FTCS (Forward Time Centered Space):

$$\frac{u_j^{n+1}u_j^n}{\mathrm{\Delta }t}=\frac{u_{j+1}^nu_{j1}^n}{2\mathrm{\Delta }x}.$$
(2)   
Here $`n`$ labels the time levels and $`j`$ the spatial grid points.

It is a standard textbook result that this scheme is unconditionally unstable. One sees this with a von Neumann stability analysis: Put

$$u_j^n=\xi ^ne^{ikj\mathrm{\Delta }x}$$
(3)   
and find that the amplification factor $`\xi `$ is

$$\xi =1+i\alpha \mathrm{sin}k\mathrm{\Delta }x,$$
(4)   
where $`\alpha =\mathrm{\Delta }t/\mathrm{\Delta }x`$. Since $`|\xi |^2>1`$ for any choice of $`\alpha `$, the method is unconditionally unstable.

Backwards differencing gives a stable scheme:

$$\frac{u_j^{n+1}u_j^n}{\mathrm{\Delta }t}=\frac{u_{j+1}^{n+1}u_{j1}^{n+1}}{2\mathrm{\Delta }x},$$
(5)   
for which

$$\xi =\frac{1}{1+i\alpha \mathrm{sin}k\mathrm{\Delta }x}.$$
(6)   
Now $`|\xi |^2<1`$ for any choice of $`\alpha `$, and so the method is unconditionally stable.

The Crank-Nicholson scheme is a second-order accurate method obtained by averaging equations (2) and (5). Now one finds

$$\xi =\frac{1+\frac{1}{2}i\alpha \mathrm{sin}k\mathrm{\Delta }x}{1\frac{1}{2}i\alpha \mathrm{sin}k\mathrm{\Delta }x}.$$
(7)   
Since $`|\xi |^2=1`$, the method is stable. It is the presence of the quantities $`u^{n+1}`$ on the right hand side of equation (5) that makes the method implicit.

The first iteration of iterated Crank-Nicholson starts by calculating an intermediate variable $`{}_{}{}^{(1)}\stackrel{~}{u}`$ using equation (2):

$$\frac{{}_{}{}^{(1)}\stackrel{~}{u}_{j}^{n+1}u_j^n}{\mathrm{\Delta }t}=\frac{u_{j+1}^nu_{j1}^n}{2\mathrm{\Delta }x}.$$
(8)   
Then another intermediate variable $`{}_{}{}^{(1)}\overline{u}`$ is formed by averaging:

$${}_{}{}^{(1)}\overline{u}_{j}^{n+1/2}=\frac{1}{2}({}_{}{}^{(1)}\stackrel{~}{u}_{j}^{n+1}+u_j^n).$$
(9)   
Finally the timestep is completed by using equation (2) again with $`\overline{u}`$ on the right-hand side:

$$\frac{u_j^{n+1}u_j^n}{\mathrm{\Delta }t}=\frac{{}_{}{}^{(1)}\overline{u}_{j+1}^{n+1/2}{}_{}{}^{(1)}\overline{u}_{j1}^{n+1/2}}{2\mathrm{\Delta }x}.$$
(10)   
(Iterated Crank-Nicholson can alternatively be implemented by averaging the right-hand side of equation (1). For linear equations, this is completely equivalent.)

Iterated Crank-Nicholson with two iterations is carried out in the same way. After steps (8) and (9), we calculate

$`{\displaystyle \frac{{}_{}{}^{(2)}\stackrel{~}{u}_{j}^{n+1}u_j^n}{\mathrm{\Delta }t}}`$ $`=`$ $`{\displaystyle \frac{{}_{}{}^{(1)}\overline{u}_{j+1}^{n+1/2}{}_{}{}^{(1)}\overline{u}_{j1}^{n+1/2}}{2\mathrm{\Delta }x}},`$ (11)   
$`{}_{}{}^{(2)}\overline{u}_{j}^{n+1/2}`$ $`=`$ $`\frac{1}{2}({}_{}{}^{(2)}\stackrel{~}{u}_{j}^{n+1}+u_j^n).`$ (12)   
Then the final step is computed analogously to equation (10):

$$\frac{u_j^{n+1}u_j^n}{\mathrm{\Delta }t}=\frac{{}_{}{}^{(2)}\overline{u}_{j+1}^{n+1/2}{}_{}{}^{(2)}\overline{u}_{j1}^{n+1/2}}{2\mathrm{\Delta }x}.$$
(13)   
Any number of iterations can be carried out in the same way.

Now consider the stability of these iterated schemes. If we define $`\beta =(\alpha /2)\mathrm{sin}k\mathrm{\Delta }x`$, and call the FTCS scheme (2) the zeroth-order method, then direct calculation shows that the amplification factors are

$`{}_{}{}^{(0)}\xi `$ $`=`$ $`1+2i\beta ,`$ (14)   
$`{}_{}{}^{(1)}\xi `$ $`=`$ $`1+2i\beta 2\beta ^2,`$ (15)   
$`{}_{}{}^{(2)}\xi `$ $`=`$ $`1+2i\beta 2\beta ^22i\beta ^3,`$ (16)   
$`{}_{}{}^{(3)}\xi `$ $`=`$ $`1+2i\beta 2\beta ^22i\beta ^3+2\beta ^4,`$ (17)   
and so on. As one would expect, these are exactly the same values one gets by expanding equation (7) in powers of $`\beta `$ and truncating at the appropriate point.

To check stability, compute $`|\xi |^2`$ for each of these expressions. You find an alternating pattern. Levels 0 and 1 are unstable; levels 2 and 3 are stable provided $`\beta ^21`$; levels 4 and 5 are unstable; levels 6 and 7 are stable provided $`\beta ^21`$; and so on. Since the stability requirement must hold for all wave numbers $`k`$, it translates into $`\alpha ^2/41`$, or $`\mathrm{\Delta }t2\mathrm{\Delta }x`$. This is just the Courant condition (the 2 occurs because of the 2 in eqn. ).

Now we see the resolution of the paradox: While the magnitude of the amplification factor for iterated Crank-Nicholson does approach 1 as the number of iterations becomes infinite, the convergence is not monotonic. The magnitude oscillates above and below 1 with ever decreasing oscillations. All the cases above 1 are unstable, although the instability might be very slowly growing for a large number of iterations.

The accuracy of the scheme is determined by the truncation error. This remains second order in $`\mathrm{\Delta }t`$ and $`\mathrm{\Delta }x`$ from the first iteration on. Doing more iterations changes the stability behavior, but not the accuracy. Since the smallest number of iterations for which the method is stable is two, there is no point in carrying out more iterations than this.

Note that there was nothing special about using the advective equation (1) for this analysis. Similar behavior is found for the wave equation, written in first-order form

$`{\displaystyle \frac{u}{t}}`$ $`=`$ $`v,`$ (18)   
$`{\displaystyle \frac{v}{t}}`$ $`=`$ $`{\displaystyle \frac{^2u}{x^2}},`$ (19)   
with the standard centered difference formula for the second derivative term. One recovers the usual Courant condition (without the factor of 2) for the stable cases.

## III Acknowledgements

This work was supported in part by NSF grant PHY 99-00672 and NASA Grant NAG5-7264 to Cornell University.