Next: Running Time, Previous: Diagnostics, Up: ode

This discussion is necessarily incomplete. Entire books exist on any
subject mentioned below (e.g., floating point error). Our goals are
modest: first, to introduce the basic notions of error analysis as they
apply to `ode`

; second, to steer you around the more obvious
pitfalls. You should look through a numerical analysis text (e.g.,
Atkinson's Introduction to Numerical Analysis) before beginning
this discussion.

We begin with some key definitions. The error of greatest concern is
the difference between the actual solution and the numerical
approximation to the solution; this is termed the *accumulated
error*, since the error is built up during each numerical step.
Unfortunately, an estimate of this error is usually not available
without knowledge of the actual solution. There are, however, several
more usable notions of error. The *single-step error*, in
particular, is the difference between the actual solution and the
numerical approximation to the solution after any single step, assuming
the value at the beginning of the step is correct.

The *relative single-step error* is the single-step error, divided
by the current value of the numerical approximation to the solution.
Why not divided by the current value of the solution itself? The reason
is that the solution is not exactly known. When free to choose a
stepsize, `ode`

will do so on the basis of the relative single-step
error. By default, it will choose the stepsize so as to maintain an
accuracy of eight significant digits in each step. That is, it will
choose the stepsize so as not to violate an upper bound of
10^(-9) on the relative single-step error. This ceiling may be
adjusted with the ‘`-r`’ option.

Where does numerical error come from? There are two sources. The first
is the finite precision of machine computation. All computers work with
floating point numbers, which are not real numbers, but only an
approximation to real numbers. However, all computations performed by
`ode`

are done to double precision, so floating point error tends
to be relatively small. You may nonetheless detect the difference
between real numbers and floating point numbers by experimenting with
the ‘`-p 17`’ option, which will print seventeen significant digits.
On most machines, that is the precision of a double precision
floating point number.

The second source of numerical error is often called the
*theoretical truncation error*. It is the difference between
the actual solution and the approximate solution due solely to the
numerical scheme. At the root of many numerical schemes is an infinite
series; for ordinary differential equations, it is a Taylor
expansion. Since the computer cannot compute all the terms in an
infinite series, a numerical scheme necessarily uses a truncated
series; hence the term. The single-step error is the sum of the
theoretical truncation error and the floating point error, though in
practice the floating point error is seldom included. The single-step
error estimated by `ode`

consists only of the theoretical
truncation error.

We say that a numerical scheme is *stable*, when applied to a
particular initial value problem, if the error accumulated during the
solution of the problem over a fixed interval decreases as the stepsize
decreases; at least, over a wide range of step sizes. With this
definition both the Runge–Kutta–Fehlberg (‘`-R`’) scheme and the
Adams–Moulton (‘`-A`’) scheme are stable (a statement based more
on experience than on theoretical results) for a wide class of problems.

After these introductory remarks, we list some common sources of accumulated error and instability in any numerical scheme. Usually, problems with large accumulated error and instability are due to the single-step error in the vicinity of a `bad' point being large.

- Singularities.
`ode`

should not be used to generate a numerical solution on any interval containing a singularity. That is,`ode`

should not be asked to step over points at which the system of differential equations is singular or undefined.You will find the definitions of singular point, regular singular point, and irregular singular point in any good differential equations text. If you have no favorite, try Birkhoff and Rota's Ordinary Differential Equations, Chapter 9. Always locate and classify the singularities of a system, if any, before applying

`ode`

. - Ill-posed problems.
For

`ode`

to yield an accurate numerical solution on an interval, the true solution must be defined and well-behaved on that interval. The solution must also be real. Whenever any of these conditions is violated, the problem is said to be*ill-posed*. Ill-posedness may occur even if the system of differential equations is well-behaved on the interval. Strange results, e.g., the stepsize suddenly shrinking to the machine limit or the solution suddenly blowing up, may indicate ill-posedness.As an example of ill-posedness (in fact, an undefined solution) consider the innocent-looking problem:

y' = y^2 y(1) = -1

The solution on the domain t > 0 is

y(t) = -1/t.

With this problem you must not compute a numerical solution on any interval that includes t=0. To convince yourself of this, try to use the ‘

`step`’ statementstep 1, -1

on this system. How does

`ode`

react?As another example of ill-posedness, consider the system

y'=1/y

which is undefined at y=0. The general solution is

y = +/- (2(t-C))^(1/2),

so that if the condition y(2)=2 is imposed, the solution will be (2t)^(1/2). Clearly, if the domain specified in a ‘

`step`’ statement includes negative values of t, the generated solution will be bogus.In general, when using a constant stepsize you should be careful not to `step over' bad points or bad regions. When allowed to choose a stepsize adaptively,

`ode`

will often spot bad points, but not always. - Critical points.
An

*autonomous*system is one that does not include the independent variable explicitly on the right-hand side of any differential equation. A*critical point*for such a system is a point at which all right-hand sides equal zero. For example, the systemy' = 2x x' = 2y

has only one critical point, at (x,y) = (0,0).

A critical point is sometimes referred to as a

*stagnation point*. That is because a system at a critical point will remain there forever, though a system near a critical point may undergo more violent motion. Under some circumstances, passing near a critical point may give rise to a large accumulated error.As an exercise, solve the system above using

`ode`

, with the initial condition x(0) = y(0) = 0. The solution should be constant in time. Now do the same with points near the critical point. What happens?You should always locate the critical points of a system before attempting a solution with

`ode`

. Critical points may be classified (as equilibrium, vortex, unstable, stable, etc.) and this classification may be of use. To find out more about this, consult any book dealing with the qualitative theory of differential equations (e.g., Birkhoff and Rota's Ordinary Differential Equations, Chapter 6). - Unsuitable numerical schemes
If the results produced by

`ode`

are bad in the sense that instability appears to be present, or an unusually small stepsize needs to be chosen needed in order to reduce the single-step error to manageable levels, it may simply be that the numerical scheme being used is not suited to the problem. For example,`ode`

currently has no numerical scheme which handles so-called `stiff' problems very well.As an example, you may wish to examine the stiff problem:

y' = -100 + 100t + 1 y(0) = 1

on the domain [0,1]. The exact solution is

y(t) = e^(-100t) + t.

It is a useful exercise to solve this problem with

`ode`

using various numerical schemes, stepsizes, and relative single-step error bounds, and compare the generated solution curves with the actual solution.

There are several rough and ready heuristic checks you may perform on
the accuracy of any numerical solution produced by `ode`

. We
discuss them in turn.

- Examine the stability of solution curves: do they converge?
That is, check how changing the stepsize affects a solution curve. As the stepsize decreases, the curve should converge. If it does not, then the stepsize is not small enough or the numerical scheme is not suited to the problem. In practice, you would proceed as follows.

- If using an adaptive stepsize, superimpose the solution curves for
successively smaller bounds on the relative single-step error (obtained
with, e.g., ‘
`-r 1e-9`’, ‘`-r 1e-11`’, ‘`-r 1e-13`’, ...). If the curves converge then the solution is to all appearances stable, and your accuracy is sufficient. - If employing a constant stepsize, perform a similar analysis by successively halving the stepsize.

The following example is one that you may wish to experiment with. Make a file named

`qcd`containing:# an equation arising in QCD (quantum chromodynamics) f' = fp fp' = -f*g^2 g' = gp gp' = g*f^2 f = 0; fp = -1; g = 1; gp = -1 print t, f step 0, 5

Next make a file named

`stability`, containing the lines:: sserr is the bound on the relative single-step error for sserr do ode -r $sserr < qcd done | spline -n 500 | graph -T X -C

This is a `shell script', which when run will superimpose numerical solutions with specified bounds on the relative single-step error. To run it, type:

sh stability 1 .1 .01 .001

and a plot of the solutions with the specified error bounds will be drawn. The convergence, showing stability, should be quite illuminating.

- If using an adaptive stepsize, superimpose the solution curves for
successively smaller bounds on the relative single-step error (obtained
with, e.g., ‘
- Check invariants of the system: are they constant?
Many systems have invariant quantities. For example, if the system is a mathematical model of a `conservative' physical system then the `energy' (a particular function of the dependent variables of the system) should be constant in time. In general, knowledge about the qualitative behavior of any dependent variable may be used to check the quality of the solution.

- Check a family of solution curves: do they diverge?
A rough idea of how error is propagated is obtained by viewing a family of solution curves about the numerical solution in question, obtained by varying the initial conditions. If they diverge sharply—that is, if two solutions which start out very close nonetheless end up far apart—then the quality of the numerical solution is dubious. On the other hand, if the curves do not diverge sharply then any error that is present will in all likelihood not increase by more than an order of magnitude or so over the interval. Problems exhibiting no sharp divergence of neighboring solution curves are sometimes called

*well-conditioned*.