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) . 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 . A particularly nice introduction to shock capturing monotone schemes is presented in .
We do not pursue these schemes here, however, as we are interested in other aspects of the transition to instability.