2
$\begingroup$

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$

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

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.

$\endgroup$
  • $\begingroup$ 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. $\endgroup$ – VorKir May 5 at 4:52
1
$\begingroup$

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

$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.