Stability of Crank-Nicolson for $u_t = iu_{xx}+2iu$

Question

I want to use the Crank-Nicolson scheme to solve the equation $$u_t = iu_{xx}+2iu$$

Here's the analysis: Suppose we make a grid, with $k = dt$ and $h = dx$, the usual notation, and also $u_j^n = u(x_j,t_n)$ . The round-off error equation for the crank-nicolson scheme for this equation is:

$$\displaystyle\frac{\epsilon_{j}^{n+1}-\epsilon_j^n}{k} = \frac{i}{2h^2}\left(\epsilon_{j-1}^n+\epsilon_{j+1}^n-2\epsilon_j^n+\epsilon_{j-1}^{n+1}+\epsilon_{j+1}^{n+1}-2\epsilon_{j}^{n+1}\right) +2i\epsilon_j^n$$

After rearranging it a bit, we can arrive at:

$$\displaystyle -\frac{ki}{2h^2}\left(\epsilon_{j-1}^{n+1}+\epsilon_{j+1}^{n+1}-2\epsilon_j^{n+1}\right)+\epsilon_{j}^{n+1} = \frac{ki}{2h^2}\left(\epsilon_{j-1}^{n}+\epsilon_{j+1}^{n}-2\epsilon_j^n\right)+(1+2ki)\epsilon_j^n$$

By writing the Fourier series of epsilon, $\epsilon_j^n = r^ne^{i\omega jh}$ and dividing by $\epsilon_j^n$ we have

$$\displaystyle -\frac{ki}{2h^2}\big(2r\cos(\omega h)-2r\big)+r=\frac{ki}{2h^2}\big(2\cos(\omega h)-2\big)+1+2ki$$

Using the identity $1-\cos(\theta) = 2\sin^2\big(\frac{\theta}{2}\big)$ we finally have

$$\displaystyle \left(1+\frac{2ki}{h^2}\sin^2\bigg(\frac{\omega h}{2}\bigg)\right)r = 1 + 2ki -\frac{2ki}{h^2}\sin^2\bigg(\frac{\omega h}{2}\bigg)$$

So

$$r = \frac{1 + 2ki -\frac{2ki}{h^2}\sin^2\big(\frac{\omega h}{2}\big)}{1+\frac{2ki}{h^2}\sin^2\big(\frac{\omega h}{2}\big)}$$

Why is this a problem? - Consider the Fourier mode $\omega \approx \frac{2\pi}{h}$. For that mode, the $\sin$ is nearly zero. For that fourier mode, we would have roughly speaking $r = \frac{1+2ki}{1}$, and that's a problem because $|r|^2 = \frac{1+4k^2}{1}= 1 + 4k^2 > 1$, so I expect the error to explode, and the scheme to be unstable.

But it isn't. Infact, the code I wrote that implements this scheme works. Why? I expected it to be unstable.

another thing to note is that the problem occurs also at $\omega = 0$

Answer 1

The scheme is indeed unstable. It explodes - but very very slowly. By printing the maximum eigenvalue of the operator i confirmed the instability. It's greater than 1. Then why does it work? because it's $1.000053263$ and my t_final is small.

Answer 2

Von Neumann (Fourier modes) stability analysis gives you only a sufficient condition for stability if you compare the amplifying coefficient $r$ with 1. If you have amplifying coefficient bounded by $1 + C\tau$, then after making $\frac{T}{\tau}$ time steps your error will be bounded by $(1 + C \tau)^{\frac{T}{\tau}}$ which has a bounded limit when $\tau \rightarrow 0$.