Crank-Nicolson shouldn't work - but it does. Von neumann / Fourier analysis

146 Views Asked by At

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}(\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}) +2i\epsilon_j^n$

After rearranging it a bit, we can arrive at:

$\displaystyle -\frac{ki}{2h^2}(\epsilon_{j-1}^{n+1}+\epsilon_{j+1}^{n+1}-2\epsilon_j^{n+1})+\epsilon_{j}^{n+1} = \frac{ki}{2h^2}(\epsilon_{j-1}^{n}+\epsilon_{j+1}^{n}-2\epsilon_j^n)+(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}(2r\cos(\omega h)-2r)+r=\frac{ki}{2h^2}(2\cos(\omega h)-2)+1+2ki$

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

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

So

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

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.

1

There are 1 best solutions below

0
On

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.