Because of the discrete nature of computers, we can only approximate the
solution to a PDE at certain grid points on the 
plane. We label these
points of approximation 
where 
and 
and employ the abbreviations 
and 
, where
is the numerical approximation at the grid points of interest for
, the actual PDE solution. We are interested, of course, in making
as small as possible. To accomplish this we must consider the
accuracy and the stability of our approximation scheme.
The accuracy of a scheme is really just a matter of minimizing the error term in the Taylor appoximation for the scheme. Because of its simplicity --- and its applicablity to eqns (7)-(9) --- we will use as an example the first order linear traveling wave equation
From Taylor's formula we have that
where 
and also that
where 
, so it follows that
since the 
and the 
terms cancel by applying eq (10). This suggests
the following numerical algorithm to simulate eq (10):
where 
. This algorithm explicitly shows how to
compute values for a new time, 
, given information at the current time,
. To compute the value at a fixed time 
, the algorithm must be run
through
time steps. Since the error from each step is 
, the total error between the
computed solution and the actual solution at 
is 
--- assuming that the errors from each time step can be added
together (which is not the case when the algorithm is unstable as we will soon
discuss). Due to the form of the total error (if the algorithm is stable), we
refer to the algorithm as being ``first order accurate in
x and first order accurate in t". If, instead of using the Taylor expansion in
eq (11), we use
we end up with a similar numerical scheme that also is first order accurate in both x and t:
The scheme in eq (14) is called a downwind scheme; the scheme in eq (15) is called an upwind scheme. We could also subtract eq (15) from eq (11), which yields
which implies that the numerical scheme
is first order accurate in time and second order accurate in space since the
error at a fixed time value is 
. Obviously,
one wishes to employ schemes with higher order errors since these schemes require
less grid points (and therefore less computational time) to obtain a computed
solution that is within a given error of the actual solution.
There is, however, more to numerical simulation of PDEs than just Taylor
expansions and accuracy questions. Consider eq (10) where the domain of
x is infinite and the initial condition 
is specified:
The unique solution to this equation is simply
In other words the value of the solution is constant along each line x=at +
C. These lines are called characteristic lines, and they occur when the
domain of x is finite as well as infinite. In other words the value of the
solution at a point depends only on which characteristic line crosses through
that point. If we consider x to be the abscissa and t the ordinate, then
each characteristic line has a positive slope of
. The constant a is often referred to as the wave speed
since it specifies how fast a disturbance at a point x propagates down to
other values of x. Now consider what happens when we apply the downwind
scheme in eq (14) in an attempt to solve eq (10). The downwind scheme only uses
the points
and
to estimate the value of the point
, but we know from the characteristic lines that the
value at this point only depends on the value of the solution at 
if 
, which is outside the range of the two points being used. This
causes the computer generated solution using the downwind scheme to go haywire
and yield gibberish. One may wonder how can this happen since we know the
scheme is accurate. The answer is that accuracy only considers local
error, whereas the stability problem that is occuring in this case is a global
issue.
The Taylor series analysis of error in determining accuracy assumes that
the solution is correct at 
, and then finds the error that is created when
the solution at 
is approximated. However, if the error accrued in a
step is magnified by each subsequent step, we can no longer just add the
individual errors produced at each step; in fact, when the errors are magnified by
subsequent steps, we observe that the solution is unstable, that is, the global error
quickly blows up. When a scheme attempts to approximate the solution at a point\
whose characteristic line crosses the
line outside the range of points used by the scheme, then the scheme, no
matter how accurate, will always magnify errors in each step and therefore be
unstable. In other words the numerical range of dependence must always contain
the theoretical range of dependence of the solution or instablity will occur.
This principle is called the CFL (Courant-Friedrichs-Lewy) condition (sometimes
referred to simply as the Courant condition) [7]. This guarantees that
the downwind scheme will always be unstable and that the upwind scheme can only
be stable if 
. Even if a scheme satisfies the CFL condition, it may
still be unstable, but we will not be analyzing any schemes that possess this
flaw in this paper.
Before analyzing the FPU nonlinear wave equations, we wish to apply the above
analysis to the linear wave equation, 
, where we use
to denote the constant 
for reasons that will soon be
apparent. If we define 
where 
and 
then we have that the wave equation is
satisfied, 
, and also that the mixed partials are equal,
. This allows us to write the wave equation as the following
first order partial differential system:
Eq (21), of course, is just the matrix form of the travelling wave equation from
eq (10), and its solution contains many of the same properties as eq (10). As
opposed to the solution depending upon one characteristic line, the solution,
, now depends on two characteristic lines whose wave speeds
are c and -c, the eigenvalues of 
. Therefore the slopes of the
two characteristic lines are the reciprocals of the wave speeds, 
and 
. As before, we must satisfy the CFL stability condition, so we
now require that
.
All that remains is to find a sufficiently accurate scheme. To obtain an
approximation for 
we add eq (11) to eq (15), yielding
which we can combine with the analogous expression from the Taylor expansions with respect to time to yield the numerical scheme
which we can reexpress as
We note that this scheme is 2nd order accurate with respect to both time and space.
We are now ready to consider the nonlinear FPU equations. We begin with the quadratic FPU equation
where 
and 
. (Note that 
since 
, like the spring constant k, is always positive in any realistic spring
model.) As with the linear wave equation we can rewrite eq (24) in vector form, but
now the matrix A will depend upon the value of
:
The wave speeds of the two characteristics are still the eigenvalues of
A, which are 
. (Note that we require 
to be
positive since otherwise the wave speeds would be imaginery and the partial
differential equation would no longer be hyperbolic, which leads to equations beyond
our analysis here and causes the scheme to become unstable). Since the speeds are
functions of
, that means the speeds vary, which implies that the characteristics are now
curves whose slopes vary as they progress in time as opposed to being straight lines
as was the case in the linear wave equation. This leads to a stability problem as
the CFL condition now requires that
. If 
is small compared to
, then 
will be positive and the variations in 
may not be
severe enough to violate the CFL condition, but if the nonlinearity becomes more
significant, which means the value of 
becomes larger, eventually the CFL
condition will be violated (or 
will become nonpositive) and
the scheme will become unstable. (Note: in the experiments in section 5, we will set
and use values of 
that are slightly smaller than 1 which will
cause the CFL condition to be violated before 
becomes nonpositive.)
Since by definition
, we see that if
gets large not only does the numerical scheme become unstable, but also
the underlying differential equation becomes less valid since higher order terms
in the Taylor approximation given in eq (4) start to become significant.
The analysis for the cubic FPU equations is similar. We now have that
where 
, which can be expressed in matrix form as
The corresponding CFL condition is 
. Note
that if 
is positive, the wave speeds will be real and we will always have
a hyperbolic differential equation.
There are numerous accurate schemes that one can apply in the region where the algorithm is stable. To simulate the quadratic FPU equation, eq (8), we can combine eq (17) with eq (23) to form
which can be rewritten as a scheme which is second order accurate in time and space:
Eq (28) is a generalization of the FPU quadratic scheme (eq (ii)).
After setting 
and 
, one
can quickly establish that the right hand sides of eq (28) and eq (ii) are
identical.
Similary, we can simulate the cubic FPU equation, eq (9), by again combining eq (17) with eq (23) to obtain the following second order accurate scheme:
The fact that the scheme in eq (30) is not the same as FPU's cubic scheme (eq (iii)) does not pose a problem. Either FPU's scheme or eq (30) can be used to simulate the cubic FPU equation as both schemes are second order accurate and will therefore produce essentially equally valid simulations. (The second order accuracy of the FPU's cubic scheme is quickly established using Taylor series.)
There is another phenomenon that is important in nonlinear hyperbolic partial
differential equations. Because the characteristics are curves, characteristics
within a family can intersect each other. (In the linear case, each family of
characteristic lines has a constant slope: 
for one family,
for the other family.) When the characteristics intersect, a
shock, i.e., a discontinuity in the solution, occurs. While multiple weak
solutions to the nonlinear equation are possible, one usually desires to simulate
the unique observed solution that the phenomena being modeled exhibits. For most
physical phenomena the solution of interest is the Lax entropy solution of the
equation, which can be simulated using monotone schemes [8]. A particularly nice
introduction to shock capturing monotone schemes is presented in [9].
We do not pursue these schemes here, however, as we are interested in other aspects of the transition to instability.