\documentclass [ ] {article}
\usepackage{verbatim}
\usepackage{epsfig}
%
\input defs
\parindent 0in
\topmargin=-0.60in
\hoffset=-0.75in
\textheight=8.75in \vsize\textheight
\textwidth=6.5in \hsize\textwidth
%
\def\IND#1#2{^{{#1}}_{{#2}}}
\def\Lhat{{\hat L}}
\def\uhat{{\hat u}}
\def\ehat{{\hat e}}
\def\fhat{{\hat f}}
\def\tauhat{{\hat \tau}}
\def\uj{u\ind{}{j}}
\def\ujp{u\ind{}{j+1}}
\def\ujph{u\ind{}{\jph}}
\def\ujm{u\ind{}{j-1}}
\def\ubar{\bar u}
\def\XT{{\lb x,t \rb}}
\def\XJT{{\lb x_j,\dt \rb}}
\def\XJN{{\lb x_j,0 \rb}}
\def\DE{\fr{\partial e_2}{\partial t}\XJN}
\def\J{\ind{}{j}}
%
\begin{document}
%
\def\be{\begin{equation}}
\def\ee{\end{equation}}
\def\bea{\begin{eqnarray}}
\def\eea{\end{eqnarray}}
\def\bean{\begin{eqnarray*}}
\def\eean{\end{eqnarray*}}
%
\begin{center}
{\bf
  PHY329~~~~Finite Difference Methods~~~~Oct 28, 1998 \\
  Notes on the 1-d Wave Equation \\
}
\end{center}
%
Recall that we are considering the following problem
(1-d wave equation with unit speed ($c^2=1$) and fixed (Dirichlet)
boundary conditions):

\be
   u_{tt}(x,t) = u_{xx}(x,t) \,, \label{WEQ}
\ee
on the domain
%
\[
   0 \le x \le 1  \quad\quad t  \ge 0 \,,
\]
%
with initial and boundary conditions:
%
\bea
   u(x,0) &=& f(x),  \label{U0} \\
   u_t(x,0) &=& g(x), \label{UDOT0} \\
   u(0,t) &=& u(1,t) \,\,\,=\,\,\, 0 \,. \label{BC}
\eea

\medskip
\begin{center}
{\bf Continuum Solution}
\end{center}
%
Assume for a moment that we are solving~(\ref{WEQ}) on an {\em infinite}
spatial domain:
%
\[
 -\infty < x < +\infty  \quad\quad t > 0 \,.
\]
%
In this case, the general solution can be written as the 
{\em superposition} of an arbitrary {\em left-moving} ``profile''
($v = -c = -1$) and an arbitrary {\em right-moving} ``profile''
($v = +c = +1$), i.e.
%
\be
   u(x,t) = l(x+t) + r(x-t) \,.
\ee
The function $l(x,t) = l(x+t)$ is constant along the ``left-directed''
{\em characteristics} of the wave equation, while the function
$r(x,t) = r(x-t)$ is constant along the ``right-directed'' characteristics
(see Figure 1). 

\begin{figure}[h]
\centerline{\epsfxsize=8cm\epsffile{char.eps}}
\caption{Characteristics of the wave equation: $u_{xx} = u_{tt}$.
         Signals (disturbances) travel along the characteristics 
         (dashed and dotted lines.)}
\end{figure}

\medskip
In fact, it is often convenient to specify the initial conditions for a wave 
equation, such as~(\ref{WEQ}), in terms of the initially left-going
and initially right-going parts of the solution.  In particular, given
the arbitrary profiles (functions), $l(x)$ and $r(x)$, we specify
%
\bea
   u(x,0) &=& l(x) \,+\, r(x) \label{ID1} \,, \\
   u_t(x,0) &=& l'(x) - r'(x) \label{ID2} \,.
\eea
%
where ${}'$ denotes ordinary differentiation (i.e. $l' \equiv
dl/dx$, etc.).

\medskip
\begin{center}
{\bf Generating Initial Data for The FDA}
\end{center}

We first re-introduce the notation wherein a caret (hat)
denotes a quantity which satisfies a difference equation. 
Thus,
$\uhat\IND{n}{j}$ is the quantity 
which satisfies the usual $O(h^2)$ approximation of~(\ref{WEQ}):
\be
   \frac{\uhat \ind\npo j - 2 \uhat\ind n j + \uhat \ind \nmo j}
        {\dt^2} =
   \frac{\uhat \ind n \jpo - 2 \uhat\ind n j + \uhat \ind n \jmo}
        {\dx^2} \,,
\label{DIFFWEQ}
\ee
while $u\IND{n}{j}$
satisfies the {\em differential} (continuum) equation, i.e.
\[
   \lb u_{tt} \rb \ind n j = \lb u_{xx} \rb \ind n j  \,,
    \quad\quad j = 1, 2, \cdots, J \quad\quad n = 0, 1, 2, \cdots \,.
\]
(since the differential equation holds 
{\em everywhere}, it must, in particular, hold at all ``grid points'' 
$(x_j,t^n)$).  We initialize~(\ref{DIFFWEQ}) by specifying
\[
 \uhat\ind{0}{j},  \quad\quad \uhat\ind{1}{j},  \quad\quad j = 1, 2, \cdots, J
\]
{\em freely}, but in a manner consistent with the boundary conditions
\[
  \uhat\ind{0}{1} = \uhat\ind{0}{J} = \uhat\ind{1}{1} = \uhat\ind{1}{J} = 0.
\]

We now assume that our difference solution will admit a 
Richardson expansion (why?)
\be
\uhat\XT = u\XT + h^2 e_2\XT + h^4 e_4\XT + \cdots \,,
\label{RICH}
\ee
and ask what this implies about how accurately we must specify
$\uhat\IND{0}{j}$ and $\uhat\IND{1}{j}$.  First recall that the problem
prescription includes the initial conditions:
\[
u\lb x,0 \rb = f\lb x \rb \,, \quad\quad\quad\quad 
u_t\lb x,0 \rb = g \lb x \rb \,.
\]
Clearly then, we can (and might as well) specify $\uhat\IND{0}{j}$ 
exactly:
\[
\uhat\ind{0}{j} = f_j \,,
\]
which implies (of course) that $e_2(x,0)  =0$.  Then in 
order for (\ref{RICH}) to hold, we must have 
\bean
\uhat\ind{1}{j} &\equiv& \uhat\XJT \\
                &=& u\XJT + h^2 e_2\XJT + O(h^4) \,,
\eean
which implies that 
\[
\uhat\XJT -  u\XJT =  h^2 e_2\XJT + O(h^4) \,.
\]
But
\bean
   e_2\XJT &=& e_2\XJN + \dt \, \DE \, + \, O(\dt^2) \\
           &=& \dt \, \DE + O(\dt^2) = \lambda h \DE + O(\dt^2)
\eean
Thus,
\[
\uhat\XJT -  u\XJT =  h^2 \lb \lambda h \DE \rb \, + \, \cdots = O(h^3) \,,
\]
and {\em we see that we must specify $\uhat\IND{1}{j}$ to $O(h^3)$ accuracy in
order for our solution to be (globally) $O(h^2)$ accurate.}

\medskip
In this example (and elsewhere) we can readily initialize 
$\uhat\IND{1}{j}$ to $O(h^3)$ accuracy simply by 
Taylor-expanding to sufficiently high order, and then
using the equations of motion to eliminate higher time derivatives:
\bean
u\XJT &=& u\XJN + \dt \, u_t\XJN + \fr{1}{2}\dt^2\,u_{tt}\XJN + O(\dt^3) \\
      &=& u\XJN + \dt \, u_t\XJN + \fr{1}{2}\dt^2\,u_{xx}\XJN + O(h^3)\,.
\eean
Recalling (\ref{ID1}-\ref{ID2}), where we specified the initial 
conditions in terms 
of initially left-moving and right-moving profiles, $l(x)$ and $r(x)$,
and their derivatives, $l'(x)$ and $r'(x)$:
\bean
u\lb x,0 \rb &=& l \lb x \rb + r \lb x \rb \,, \\
u_t\lb x,0 \rb &=& l' \lb x \rb - r' \lb x \rb \,,
\eean
and noting that
\[
u_{tt} \lb x,0 \rb = l''\lb x \rb + r'' \lb x \rb  \,,
\]
we have for our initialization:
\bea
 \uhat\ind{0}{j} &=& l\J + r\J                             \,,\\
 \uhat\ind{1}{j} &=& l\J + r\J + \dt \lb l'\J - r'\J \rb +
                     \fr{1}{2} \dt^2 \lb l''\J + r''\J \rb  \,.
\eea

\medskip
\begin{center}
{\bf Implementation Notes}
\end{center}

In treating time-dependent PDEs using FD techniques, it is generally 
not feasible (even in 1-d) to store {\em all} of the data which 
is generated during the evolution---i.e., in the current case, 
all of the $u\IND n j$ for $n = 0, 1, \cdots$.   Indeed, for 
multi-dimensional problems (3-d calculations in particular),
computations are often {\em memory-limited} or {\em memory-bound},
and then it becomes important to implement the FD algorithm using 
as {\em little} storage as possible, using periodic output (usually to 
a file) to save the generated data from intermediate time-steps.  
A natural way to represent (store) $u\IND\nmo j$, $u\IND n j$ and 
$u\IND \npo j$ is to use a two dimensional array:
\begin{verbatim}
     real*8     u(maxj,ntlevs)
\end{verbatim}
where the constant {\tt maxj} is the maximum number of spatial
grid points allowed, and the constant {\tt ntlevs} is the number of 
time-levels of data stored at any instant.   Superficial 
investigation of the difference equation~(\ref{DIFFWEQ}) suggests 
that we need {\tt ntlevs = 3}, but we can actually code~(\ref{DIFFWEQ})
using storage for only {\em two} levels (level $n$ and level $n+1$) since 
$u\IND\nmo j$ is only referenced just before $u\IND\npo j$ is 
computed.  The following code segment also illustrates a technique 
wherein we use integer variables {\tt n}, {\tt nm1} and {\tt np1}
as ``pointers'' to whichever columns of the array {\tt u} currently hold 
the levels $n$, $n-1$ and $n+1$ data respectively.  Using this 
mechanism we can effect (implement) a time-step advance---wherein
the level $n$ data becomes level $n-1$ and level $n+1$ becomes 
level $n$---simply be re-defining the scalars {\tt n}, {\tt nm1} and 
{\tt np1}.  In particular, with this approach there is no need to 
{\em copy} any of the grid function data to advance from
$t=t^n$ to $t=t^{n+1}$.

\begin{verbatim}
      real*8      u(maxj,2)
      integer     n,   nm1,    np1,   nswap

      n   = 2       ! Initialize "pointers"---note 'n-1' and 'n+1'
      nm1 = 1       ! storage is shared; only possible due to 
      np1 = 1       ! explicit nature of FDA

c     Section of code which initializes  u(j,nm1) and u(j,n)
c     goes here.

      do it = 2 , nt       ! Begin time-step loop
         u(1,np1) = 0.0d0
         do j = 2 , nx - 1
            u(j,np1) = 2.0d0 * u(j,n) - u(j,nm1) +
     &          lamsq * (u(j+1,n) - 2.0d0 * u(j,n) + u(j-1,n))
         end do
         u(nx,np1) = 0.0d0

c        Periodic output of u(j,np1) goes here 

c        Swap time level pointers (time step advance)

         nswap = np1
         np1   = n
         n     = nswap
         nm1   = np1
      end do               ! End time step loop
\end{verbatim}

\medskip
\begin{center}
{\bf Stability Analysis}
\end{center}

One of the most frustrating---yet fascinating---features of FD solutions 
of time dependent problems, is that the discrete solutions often
``blow up''---e.g. floating-point overflows are generated at 
some point in the evolution.  Although ``blow-ups'' 
can sometimes be caused by legitimate (!) ``bugs''---i.e. an 
incorrect implementation---at other times it is simply the 
{\em nature of the FD scheme} which causes problems.   We are 
thus lead to consider the {\em stability} of solutions 
of difference equations (as well as their 
differential-equation progenitors).

\medskip
Let us again consider our prototypical time-dependent 
differential equation~(\ref{WEQ}), and let us now remark
that this is a {\em linear, non-dispersive} wave equation,
a consequence of which is the fact that the ``size'' of 
the waves does {\em not} change with time:
\be
   \Vert u(x,t) \Vert \sim \Vert u(x,0) \Vert \,,
\label{CONSTNORM}
\ee

where $\Vert \cdot \Vert$ is an suitable norm, such as the $L_2$ norm:
\be
   \Vert u(x,t) \Vert \equiv \lb \int_0^1 \, u(x,t)^2 \, dx \rb^{1/2}  \,.
\ee
We will use the property captured by~(\ref{CONSTNORM}) as our working
definition of stability.  In particular, if you believe~(\ref{CONSTNORM}) 
is true for the wave equation, then you believe the wave equation is stable.

\medskip
Fundamentally, if our FDA approximation {\em converges}, then we 
expect the same behaviour for the difference solution (note that in 
this section, we drop the carets 
on solutions of difference equations):
\be
   \Vert u\ind n j \Vert \sim \Vert u \ind 0 j \Vert \,.
   \label{DIFFSTAB}
\ee

Now, we construct our FD solution by {\em iterating in time}, generating
\[
   u\ind 0 j, \,\, u\ind 1 j, \,\, u\ind 2 j, \,\, u\ind 3 j, \,\,
   u\ind 4 j, \,\, \cdots
\]
in succession, using the FD equation
\[
   u\ind{\npo}{j} = 2 u\ind{n}{j} - u\ind{\nmo}{j} + 
     \lambda^2 \lb u\ind{n}{\jpo} - 2 u\ind{n}{j} + u\ind{n}{\jmo} \rb \,.
\]
As it turns out, we are {\em not} guaranteed that~(\ref{DIFFSTAB}) 
holds for all values of $\lambda \equiv \dt/\dx$.  In fact, for certain
$\lambda$ (all $\lambda > 1$, as we shall see), we have
\[
	\Vert u\ind n j \Vert \gg \Vert u \ind 0 j \Vert \,,
\]
and for those $\lambda$, $\Vert u\ind{n}{} \Vert$
{\em diverges} from $u$, even (especially!) as 
$h\to 0$---that is, the difference scheme is {\em unstable}.

\medskip
In fact, for many wave problems (including all linear problems),
given that a FD scheme is {\em consistent} (i.e. so that ${\hat \tau} \to 0$
as $h\to0$), stability is the necessary and sufficient condition for 
convergence (and vice versa).

\medskip
\begin{center}
{\em Heuristic Stability Analysis}
\end{center}
Let us write a general time-dependent FDA in the form
\def\U{{\bf u}}
\def\E{{\bf e}}
\def\G{{\bf G}}
\def\N{\ind{n}{}}
\be
   \U\ind\npo{} = \G [ \U\ind{n}{} ] \,,
\label{DIFFFO}
\ee
where $\G$ is some {\em update operator} (linear in our example 
problem), and $\U$ is a column vector containing sufficient 
unknowns to write the problem in first-order-in-time form.  For 
example, if we introduce a new, auxiliary set of unknowns, $v\ind{n}{j}$,
defined by
\[
   v\ind{n}{j} = u\ind{\nmo}{j} \,,
\]
then we can rewrite the differenced-wave-equation~(\ref{DIFFWEQ}) as
\bea
   u\ind{\npo}{j} &=& 2 u\ind{n}{j} - v\ind{n}{j} + 
      \lambda^2 \lb u\ind{n}{\jpo} - 2u\ind{n}{j} + u\ind{n}{\jmo} \rb \,,
  \label{DIFFFO1}  \\
   v\ind{\npo}{j} &=& u\ind n j \,,
  \label{DIFFFO2}
\eea
so with 
\[
   \U\ind{n}{} = [ u\ind n 1, v\ind n 1, \,\, u\ind n 2, v\ind n 2, \, \cdots 
                   \, u\ind n J, v\ind n J] \,,
\]
(for example), (\ref{DIFFFO1}-\ref{DIFFFO2}) is clearly of the 
form~(\ref{DIFFFO}).  Equation (\ref{DIFFFO}) provides us with a 
compact way of describing the solution of the FDA.  Given initial 
data, $\U\ind{0}{}$, the solution after $n$ time-steps is
\be
   \U\N = \G\N \U\ind{0}{} \,,
\label{GITER}
\ee
where $\G\N$ is the $n$-th power of the matrix $\G$.
Now, assume that $\G$ has a complete set of orthonormal eigenvectors
\[
   \E_k, \quad k = 1, 2, \, \cdots \, J \,,
\]
and corresponding eigenvalues
\[
   \mu_k, \quad k = 1, 2, \, \cdots \, J  \,,
\]
so that
\[
   \G \, \E_k = \mu_k \, \E_k, \quad k = 1, 2, \, \cdots \, J  \,.
\]
We can then write the initial data as (spectral decomposition):
\[
   \U\ind{0}{} = \sum_{k=1}^{J} \, c\ind{0}{k} \, \E_k \,,
\]
where the $c\ind{0}{k}$ are coefficients. 
Using~(\ref{GITER}), the solution at time-step $n$ is then
\bea
   \U\N &=& \G\ind{n}{} \lb \sum_{k=1}^{J} \, c\ind{0}{k} \, \E_k \rb     \\
        &=& \sum_{k=1}^{J} \, c\ind{0}{k} \, \lb \mu_k \rb^n \, \E_k \,.
\eea
Clearly, if the difference scheme is to be stable, we must have
\be
\fbox{$
   \vert \mu_k \vert \le 1 \quad\quad k = 1, 2, \, \cdots \, J
$}
\ee
(Note: $\mu_k$ will be complex in general, so $\vert \mu \vert$ 
denotes complex modulus, $\vert \mu \vert \equiv \sqrt{\mu \mu^\star}$).

\smallskip
Geometrically, then, the eigenvalues of the update matrix must lie 
on or within the unit circle (see Figure 2).

\begin{figure}[h]
\centerline{\epsfxsize=8cm\epsffile{ev.eps}}
\caption{Schematic illustration of location in complex plane of 
         eigenvalues of update matrix G.   In this case, all
         eigenvalues (dots) lie on or within the unit circle,  
         indicating that the corresponding finite difference 
         scheme is stable.
        }
\end{figure}

\def\UT{{\tilde {\bf u}}}
\def\GT{{\tilde {\bf G}}}
\medskip
\begin{center}
{\em Von-Neumann (Fourier) Stability Analysis}
\end{center}
Von-Neumann stability analysis is based on the ideas sketched 
above, but additionally assumes that the difference 
equation is linear with constant coefficients, and that the 
boundary conditions are periodic.  We can then use Fourier 
analysis, which has the same benefits in the discrete domain---difference 
operators in real-space variable $x \longrightarrow$ algebraic 
operations in Fourier-space variable $k$---as it does in the continuum
Schematically, instead of writing 
\[
   \U\ind{\npo}{}(x) = \G [ \U\ind{n}{}(x) ] \,,
\]
we consider the Fourier-domain equivalent:
\[
   \UT\ind{\npo}{}(k) = \GT [ \UT\ind{n}{}(k) ] \,,
\]
where $k$ is the wave-number (Fourier-space variable) and 
$\UT$ and $\GT$ are the Fourier-transforms of $\U$ and 
$\G$, respectively.  Specifically, we define the Fourier-transformed
grid function via
\be
   \UT\ind{n}{}(k) = \frac{1}{\sqrt{2\pi}} 
      \int_{-\infty}^{+\infty} \, e^{-ikx} \, \U\ind{n}{}(x) \, dx \, .
\label{FFT}
\ee
For a general difference scheme, we will find that 
\[
   \UT\ind{\npo}{}(k) = \GT(\xi) \, \UT\ind{n}{}(k) \,,
\]
where $\xi \equiv k h$, and we will have to show that 
$\GT(\xi)$'s eigenvalues lie within or on the unit circle for 
all conceivable $\xi$.  The appropriate range for $\xi$ is 
\[
   -\pi \le \xi \le \pi \,,
\]
since the shortest wavelength representable on a uniform mesh with spacing 
$h$ is $\lambda = 2h$  (Nyquist limit), corresponding to a maximum wave number
$k = (2\pi)/\lambda = \pm \pi / h$.

\medskip
Let us consider the application of the Von-Neumann stability 
analysis to our current model problem.  We first define 
a (non-divided) difference operator $D^2$ as follows:
\[
   D^2 \, u(x) = u(x+h) - 2 u(x) + u(x-h) \,.
\]
Then, suppressing the spatial grid index, we can write the 
first-order form of the difference equation~(\ref{DIFFFO1}-\ref{DIFFFO2}) 
as
\bean
   u\ind\npo{} &=& 2 u\ind n{} - v\ind n{} + 
    \lambda^2 \, D^2 \, u\ind n{}\,, \\
   v\ind\npo{} &=& u\ind n{} \,,
\eean
or
\be
   \left[
   \begin{array}{c}
      u \\
      v 
   \end{array}
   \right]^{n+1} =
   \left[
   \begin{array}{cc}
      2 + \lambda^2 \, D^2  & -1 \\
       1                    & 0
   \end{array}
   \right]
   \left[
   \begin{array}{c}
      u \\
      v 
   \end{array}
   \right]^n
    \,.
\label{UPDATE}
\ee

In order to perform the Fourier transform, we need to know the 
action of $D^2$ in Fourier-space.  Using the transform 
inverse to~(\ref{FFT}) we have
\[
   u(x) = \frac{1}{\sqrt{2\pi}} 
     \int_{-\infty}^{+\infty} \, e^{ikx} \, \ut(k) \, dk \,,
\]
so
\bean
   D^2 \, u(x) = u(x+h) - 2 u(x) + u(x-h) 
     &=& \int_{-\infty}^{+\infty}  
       \lb e^{ikh} - 2 + e^{-ikh} \rb e^{ikx} \, \ut(k) \, dk \\
     &=& \int_{-\infty}^{+\infty} 
        \lb e^{i\xi} - 2 + e^{-i\xi} \rb e^{ikx} \, \ut(k) \, dk \,.
\eean

Now consider the quantity $-4\sin^2(\xi/2)$:
\bean
   -4\sin^2 \frac{\xi}{2}&=& -4 \lb \frac{e^{i\xi/2}-e^{-i\xi/2}}{2i}\rb^2 \\
    &=& \lb e^{i\xi/2}-e^{-i\xi/2} \rb^2 =
    e^{i\xi} - 2 + e^{-i\xi} \,,
\eean
so
\[
   D^2 \, u(x) = \frac{1}{\sqrt{2\pi}} 
     \int_{-\infty}^{+\infty} \lb -4\sin^2 \frac{\xi}{2} \rb e^{ikx} \, \ut(k) \, dk \,.
\]

In summary, under Fourier transformation, we have
\bean
   \U(x) &\longrightarrow& \UT(k) \,, \\ 
   D^2 \U(x) &\longrightarrow& -4\sin^2 \frac{\xi}{2} \UT(k)  \,.
\eean

Using this result in the Fourier transform of~(\ref{UPDATE}), we see that we 
need to compute the eigenvalues of 
\[
   \left[
   \begin{array}{cc}
      2 - 4\lambda^2\sin^2 (\xi / 2) & -1  \\
      1                         & 0
   \end{array}
   \right] \,,
\]
and determine the conditions under which the eigenvalues lie on
or within the unit circle.  The characteristic equation (whose 
roots are the eigenvalues) is 

\[
   \Bigg\vert
   \begin{array}{cc}
      2 - 4\lambda^2\sin^2 (\xi/2)  - \mu & -1 \\
      1                         & -\mu \,,
   \end{array}
   \Bigg\vert
\,\,=\,\, 0
\]
or
\[
   \mu^2 + \lb 4 \lambda^2 \sin^2 \frac{\xi}{2} - 2\rb \mu + 1 = 0 \, .
\]

This equation has roots
\[
   \mu(\xi) = \lb 1 - 2\lambda^2 \sin^2 \frac{\xi}{2} \rb
              \pm \lb \lb 1 - 2\lambda^2 \sin^2 \frac{\xi}{2}\rb - 1\rb^{1/2}
\,.
\]
We now need to find sufficient conditions for
\[
   \vert \mu(\xi) \vert \le 1,
\]
or equivalently
\[
   \vert \mu(\xi) \vert^2 \le 1.
\]
To this end, we note that we can write
\[
   \mu(\xi) = (1 - Q) \,\pm\, ((1 - Q)^2 - 1)^{1/2} \,,
\]
where the quantity, $Q$ 
\[
   Q \equiv 2 \lambda \sin^2 \frac{\xi}{2} \,,
\]
is {\em real} and {\em non-negative} ($Q \ge 0$).  There are now 
two cases to consider:
\begin{enumerate}
   \item $(1 - Q)^2 - 1 \le 0$ \,,
   \item $(1 - Q)^2 - 1 > 0$ \,.
\end{enumerate}
In the first case, $((1 - Q)^2 - 1)^{1/2}$ is purely imaginary, so 
we have
\[
   \vert \mu(\xi) \vert^2 = (1 - Q)^2 + (1 - (1 - Q)^2) = 1 \,.
\]
In the second case,  $(1 - Q)^2 - 1 > 0 \,\longrightarrow\, (1 - Q)^2 > 1 
\,\longrightarrow\, Q  > 2$, and then we have 
\[
   1 - Q - ((1 - Q^2) - 1)^{1/2} \, < \, -1 \,,
\]
so, in this case, our stability criterion will {\em always} be violated.  We 
thus conclude that a necessary condition for Von-Neumann stability is 
\[
(1 - Q)^2 - 1 \le 0  \,\longrightarrow\,
(1 - Q)^2 \le 1 \,\longrightarrow\,  Q \le 2 \,.
\]
Since $Q \equiv 2 \lambda \sin^2 (\xi/2)$ and 
$\sin^2 (\xi/2) \le 1$, we must therefore have
\[
\lambda \equiv \frac{\dt}{\dx} \le 1 \,,
\]
for stability of our scheme~(\ref{DIFFWEQ}).  This condition is often called 
the CFL condition---after Courant, Friedrichs and Lewy who 
derived it in 1928 (the ratio $\lambda=\dx/\dt$ is also frequently called 
the {\em Courant number}).   In practical terms, we must limit
time-discretization scale , $\dt$, to values no larger than 
the space-discretization scale, $\dx$.   Furthermore, this 
type of instability has a ``physical'' interpretation, often
summarized by the statement {\em the numerical domain of 
dependence of an explicit difference scheme must contain the 
physical domain of dependence}.  
\end{document}
