## Stability Analysis

As in the case of ODEs, the numerical method must be not only accurate but also stable in order to converge to the exact solution when t and h tend to zero. The

 m -1 m +1 m - 2 m m + 2

Figure 2 The schematical graph of the real part of the 'fastest' spatial harmonic in . Its imaginary part looks qualitatively the same.

following von Neumann analysis is the standard tool to test stability of a numerical method. First, we note that if Um and Vm are two slightly different (due to, e.g., different initial conditions) solutions of , then the error e"m = Vn — Un satisfies the same eqn  (since the latter is linear). This error can be expanded into a series of spatial Fourier harmonics:

where ci (n) are amplitudes of the harmonics and ki are inversely proportional to the spatial oscillation periods of the harmonics. (Note that while the variable e"m is real, both the exponentials and their amplitudes ci are complex numbers; this is a standard mathematical procedure that does not restrict the validity of the approach.) The 'slowest' harmonics have ki k 0 and change very slowly with x. The 'fastest' harmonic has ki k k / h, which corresponds to the shortest possible oscillation period of 2h (see Figure 2).

Next, one assumes that the amplitudes ci (n) depend exponentially on the number of iterations: ci (n) = pn. If the subsequent analysis finds that \p\< 1, then the method is stable, and it is unstable otherwise. The substitution of  into  yields the stability condition to be \p\ = \ 1 — 2r+ 2r cos(kih)\< 1, which with the above bounds for ki (min ki = 0 and max ki = -/h) gives the stability condition for :

If this condition is violated (due to taking either r 'too large' or h 'too small'), the fastest spatial harmonic begins to grow exponentially, and the numerical solution soon becomes 'drowned' in a high-frequency noise.

From  one has r = O(h ). Then the accuracy of method  can be expressed only in terms of the mesh size h: O (r) + O (h2) = O (h2) + O (h2) = O (h2). This leads to the following observation. As we noted above, method  is the counterpart of the explicit Euler

 ( ^ ( n +1 J V J w n + 1/2 r ^ ( ^ ( n m J ^ -1 m J K m 1

Figure 3 The stencil of the CN method .

method for ODEs. Other explicit methods may have a higher-order accuracy in the time step; for example, the PDE counterpart for  of the Heun method (see Numerical Methods for Local Models) has accuracy 0(r2). However, this would not increase the overall accuracy of such a method. The reason is that the stability condition for this method is also given by . Then, similarly to the above, the overall accuracy of the counterpart of the Heun method for  is O(r2) + O(h2) = O(h4) + O(h2) = O(h2). 