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

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$$

• It will be completely stable if you discretize the term $2iu$ with a time average: $i(\epsilon^{n+1}+\epsilon^n)$. – David Ketcheson Apr 28 '19 at 11:51

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.
• It is not a proof. You can have the largest eigenvalue behaving as $1 + C \tau$, and the scheme will be stable, because $(1 + C\tau)^{\frac{T}{\tau}}$ will stay bounded. – VorKir May 5 '19 at 4:52
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$$.