ode Program
The GNU ode utility can produce a numerical solution to the
initial value problem for many systems of first-order ordinary
differential equations (ODE's). ode can also be used to solve
systems of higher-order ODE's, since a simple procedure converts an
n'th-order equation into n first-order equations. The
output of ode can easily be piped to graph, so that one or
more solution curves may be plotted as they are generated.
Three distinct schemes for numerical solution are implemented: Runge--Kutta--Fehlberg (the default), Adams--Moulton, and Euler. The Runge--Kutta--Fehlberg and Adams--Moulton schemes are available with adaptive stepsize.
We begin with some standard definitions. A differential equation
is an equation involving an unknown function and its derivatives. A
differential equation is ordinary if the unknown function
depends on only one independent variable, often denoted t.
The order of the differential equation is the order of the
highest-order derivative in the equation. One speaks of a family, or
system of equations when more than one equation is involved. If
the equations are dependent on one another, they are said to be
coupled. A solution is any function satisfying the
equations. An initial value problem is present when there exist
subsidiary conditions on the unknown function and its derivatives, all
of which are given at the same value of the independent variable. In
principle, such an `initial condition' specifies a unique solution.
Questions about the existence and uniqueness of a solution, along with
further terminology, are discussed in any introductory text. (See
Chapter 1 of Birkhoff and Rota's Ordinary Differential
Equations. For this and other references relevant to ode, see
section Bibliography on ode and solving differential equations.)
In practical problems, the solution of a differential equation is usually not expressible in terms of elementary functions. Hence the need for a numerical solution.
A numerical scheme for solving an initial value problem produces an
approximate solution, using only functional evaluations and the
operations of arithmetic. ode solves first-order initial value
problems of the form:
@ifnottex
x' = f(t,x,y,...,z) y' = g(t,x,y,...,z) . . . z' = h(t,x,y,...,z)
given the initial values for each dependent variable at the initial value of the independent variable t, i.e.,
x(a) = b
y(a) = c
.
.
.
z(a) = d
t = a
@ifnottex where a,b,c,...,d are constants.
@ifnottex
For ode to be able to solve such a problem numerically, the
functions f,g,...,h must be expressed, using the usual operators
(+, -, *, /, and ^), in terms of
certain basic functions that ode recognizes. These are the same
functions that the plotting program gnuplot recognizes.
Moreover, each of f,g,...,h must be given explicitly. ode
cannot deal with a system in which one or more of the first derivatives
is defined implicitly rather than explicitly.
All schemes for numerical solution involve the calculation of an
approximate solution at discrete values of the independent variable
t, where the `stepsize' (the difference between any two
successive values of t, usually denoted h) may be
constant or chosen adaptively. In general, as the stepsize
decreases the solution becomes more accurate. In ode, the
stepsize can be adjusted by the user.
ode
The following examples should illustrate the procedure of stating an
initial value problem and solving it with ode. If these
examples are too elementary, see section The ode input language formally specified, for a formal
specification of the ode input language. There is also a
directory containing examples of ode input, which is distributed
along with the GNU plotting utilities. On most systems it is
installed as `/usr/share/ode' or `/usr/local/share/ode'.
Our first example is a simple one, namely
y'(t) = y(t)
with the initial condition
y(0) = 1
The solution to this differential equation is
@ifnottex
y(t) = e^t.
In particular
@ifnottex
y(1) = e^1 = 2.718282
to seven digits of accuracy.
You may obtain this result with the aid of ode by typing on the
command line the sequence of commands
ode y' = y y = 1 print t, y step 0, 1
Two columns of numbers will appear. Each line will show the value of the independent variable t, and the value of the variable y, as t is `stepped' from 0 to 1. The last line will be
1 2.718282
as expected. You may use the `-p' option to change the precision. If, for example, you type `ode -p 10' rather than `ode', you will get ten digits of accuracy in the output, rather than seven (the default).
After the above output, ode will wait for further instructions.
Entering for example the line
step 1, 0
should yield two more columns of numbers, containing the values of t and y that are computed when t is stepped back from 1 to 0. You could type instead
step 1, 2
to increase rather than decrease t. To exit ode,
you would type a line containing only `.', i.e. a single period,
and tap `return'. ode will also exit if it sees an end-of-file
indicator in its input stream, which you may send from your terminal by
typing control-D.
Each line of the preceding example should be self-explanatory. A
`step' statement sets the beginning and the end of an interval
over which the independent variable (here, t) will range, and
causes ode to set the numerical scheme in motion. The initial
value appearing in the first `step' statement (i.e., 0) and the
assignment statement
y = 1
are equivalent to the initial condition y(0) = 1. The statements
`y' = y' and `y = 1' are very different: `y' = y'
defines a way of computing the derivative of y, while `y
= 1' sets the initial value of y. Whenever a `step'
statement is encountered, ode tries to step the independent
variable through the interval it specifies. Which values are to be
printed at each step is specified by the most recent `print'
statement. For example,
print t, y, y'
would cause the current value of the independent variable t, the variable y, and its derivative to be printed at each step.
To illustrate ode's ability to take its input or the initial part
of its input from a file, you could prepare a file containing the
following lines:
# an ode to Euler y = 1 y' = y print t, y, y'
Call this file `euler'. (The `#' line is a comment line,
which may appear at any point. Everything from the `#' to the
end of the line on which it appears will be ignored.) To process
this file with ode, you could type on your terminal
ode -f euler step 0, 1
These two lines cause ode to read the file `euler', and the
stepping to take place. You will now get three quantities (t,
y, and y') printed at each of the values of t
between 0 and 1. At the conclusion of the stepping, ode
will wait for any further commands to be input from the terminal. This
example illustrates that
ode -f euler
is not equivalent to
ode < euler
The latter would cause ode to take all its input from the file
`euler', while the former allows subsequent input from the
terminal. For the latter to produce output, you would need to include a
`step' line at the end of the file. You would not need to include
a `.' line, however. `.' is used to terminate input only
when input is being read from a terminal.
A second simple example involves the numerical solution of a second-order differential equation. Consider the initial value problem
y''(t) = -y(t) y(0) = 0 y'(0) = 1
Its solution would be
@ifnottex
y(t) = sin(t)
To solve this problem using ode, you must express this
second-order equation as two first-order equations. Toward this end you
would introduce a new function, called yp say, of the independent
variable t. The pair of equations
y' = yp yp' = -y
would be equivalent to the single equation above. This sort of reduction of an n'th order problem to n first order problems is a standard technique.
To plot the variable y as a function of the variable t, you could create a file containing the lines
# sine : y''(t) = -y(t), y(0) = 0, y'(0) = 1 sine' = cosine cosine' = -sine sine = 0 cosine = 1 print t, sine
(y and yp have been renamed sine and cosine, since that is what they will be.) Call this file `sine'. To display the generated data points on an X Window System display as they are generated, you would type
ode -f sine | graph -T X -x 0 10 -y -1 1 step 0, 2*PI .
After you type the ode line, graph -T X will pop
up a window, and after you type the `step' line, the generated
dataset will be drawn in it. The `-x 0 10' and `-y -1 1'
options, which set the bounds for the two axes, are necessary if you
wish to display points in real time: as they are generated.
If the axis bounds were not specified on the command line,
graph -T X would wait until all points are read from the
input before determining the bounds, and drawing the plot.
A slight modification of this example, showing how ode can
generate several datasets in succession and plot them on the same graph,
would be the following. Suppose that you type on your terminal the
following lines.
ode -f sine | graph -T X -C -x 0 10 -y -1 1 step 0, PI step PI, 2*PI step 2*PI, 3*PI .
Then the sine curve will be traced out in three stages. Since the
output from each `step' statement ends with a blank line,
graph -T X will treat each section of the sine curve as a
different dataset. If you are using a color display, each of the three
sections will be plotted in a different color. This is a feature
provided by graph, which normally changes its linemode after each
dataset it reads. If you do not like this feature, you may turn it off
by using graph -T X -B instead of graph -T X.
In the above examples, you could use any of the other variants of
graph instead of graph -T X. For example, you could use
graph -T ps to obtain a plot in encapsulated Postscript format,
by typing
ode -f sine | graph -T ps > plot.ps step 0, 2*PI .
You should note that of the variants of graph, the seven variants
graph -T pnm, graph -T gif, graph -T ai,
graph -T ps, graph -T fig, graph -T pcl and
graph -T hpgl do not produce output in real time, even when the
axis bounds are specified with the `-x' and `-y' options.
So if any of these seven variants is used, the plot will be produced
only when input from ode is terminated, which will occur when you
type `.'.
In the preceding examples, the derivatives of the dependent variables
were specified by comparatively simple expressions. They are allowed to
be arbitrarily complicated functions of the dependent variables and the
independent variable. They may also involve any of the functions that
are built into ode. ode has a fair number of functions
built in, including abs, sqrt, exp, log, log10,
sin, cos, tan, asin, acos, atan, sinh,
cosh, tanh, asinh, acosh, and atanh. Less familiar
functions which are built into it are besj0, besj1,
besy0, besy1, erf, erfc, inverf, lgamma,
gamma, norm, invnorm, ibeta, and igamma. These have
the same definitions as in the plotting program gnuplot. (All
functions take a single argument, except for ibeta, which takes
three, and igamma, which takes two). ode also knows the
meaning of the constant `PI', as the above examples show. The
names of the preceding functions are reserved, so, e.g., `cos' and
`sin' may not be used as names for variables.
Other than the restriction of avoiding reserved names and keywords, the
names of variables may be chosen arbitrarily. Any sequence of
alphanumeric characters starting with an alphabetic character may be
used; the first 32 characters are significant. It is worth noting
that ode identifies the independent variable by the fact that it
is (or should be) the only variable that has not appeared on the left
side of a differential equation or an initial value assignment. If
there is more than than one such variable then no stepping takes place;
instead, an error message is printed. If there is no such variable,
a dummy independent variable is invented and given the name
`(indep)', internally.
ode
We explain here how to use some additional features of ode.
However, the discussion below does not cover all of its capabilities.
For a complete list of command-line options, see section ode command-line options.
It is easy to use ode to create plots of great beauty. An
example would be a plot of a strange attractor, namely the Lorenz
attractor. Suppose that a file named `lorenz' contains the
following lines.
# The Lorenz model, a system of three coupled ODE's with parameter r. x' = -3*(x-y) y' = -x*z+r*x-y z' = x*y-z r = 26 x = 0; y = 1; z = 0 print x, y step 0, 200
Then executing the command
ode < lorenz | graph -T X -C -x -10 10 -y -10 10
would produce a plot of the Lorenz attractor (strictly speaking, a plot of one of its two-dimensional projections). You may produce a Postscript plot of the Lorenz attractor, and print it, by doing something like
ode < lorenz | graph -T ps -x -10 10 -y -10 10 -W 0 | lpr
The `-W 0' ("zero width") option requests that graph -T ps
use the thinnest line possible, to improve the visual appearance of the
plot on a printer or other Postscript device.
Besides plotting a visually striking object in real time, the Lorenz
attractor example shows how statements may be separated by semicolons,
rather than appearing on different lines. It also shows how to use
symbolic constants. In the description read by ode the
parameter r is a variable like x, y, and
z. But unlike them it is not updated during stepping, since no
formula for its derivative r' is given.
Our second example deals with the interactive construction of a `phase portrait': a set of solution curves with different initial conditions. Phase portraits are of paramount interest in the qualitative theory of differential equations, and also possess @ae{}sthetic appeal.
Since a description read by ode may contain any number of
`step' statements, multiple solution curves may be plotted in a
single run. The most recent `print' statement will be used with
each `step' statement. In practice, a phase portrait would be
drawn from a few well-chosen solution curves. Choosing a good set of
solution curves may require experimentation, which makes interactivity
and real-time plotting all-important.
As an example, consider a so-called Lotka--Volterra predator--prey model. Suppose that in a lake there are two species of fish: A (the prey) who live by eating a plentiful supply of plants, and B (the predator) who eat A. Let x(t) be the population of A and y(t) the population of B at time t. A crude model for the interaction of A and B is given by the equations
x' = x(a-by) y' = y(cx-d)
where a, b, c, d are positive constants. To draw a phase portrait for this system interactively, you could type
ode | graph -T X -C -x 0 5 -y 0 5 x' = (a - b*y) * x y' = (c*x - d) * y a = 1; b = 1; c = 1; d = 1; print x, y x = 1; y = 2 step 0, 10 x = 1; y = 3 step 0, 10 x = 1; y = 4 step 0, 10 x = 1; y = 5 step 0, 10 .
Four curves will be drawn in succession, one per `step' line. They
will be periodic; this periodicity is similar to the fluctuations
between predator and prey populations that occur in real-world
ecosystems. On a color display the curves will appear in different
colors, since by default, graph changes the linemode between
datasets. That feature may be turned off by using graph -T X
-B rather than graph -T X.
It is sometimes useful to use ode and graph to plot
discrete points, which are not joined by line segments to form a curve.
Our third example illustrates this. Suppose the file `atwoods'
contains the lines
m = 1 M = 1.0625 a = 0.5; adot = 0 l = 10; ldot = 0 ldot' = ( m * l * adot * adot - M * 9.8 + m * 9.8 * cos(a) ) / (m + M) l' = ldot adot' = (-1/l) * (9.8 * sin(a) + 2 * adot * ldot) a' = adot print l, ldot step 0, 400
The first few lines describe the functioning of a so-called swinging Atwood's machine. An ordinary Atwood's machine consists of a taut cord draped over a pulley, with a mass attached to the cord at each end. Normally, the heavier mass (M) would win against the lighter mass (m), and draw it upward. A swinging Atwood's machine allows the lighter mass to swing back and forth as well as move vertically.
The `print l, ldot' statement requests that the vertical position and vertical velocity of the lighter mass be printed out at each step. If you run the command
ode < atwoods | graph -T X -x 9 11 -y -1 1 -m 0 -S 1 -X l -Y ldot
you will obtain a real-time plot. The `-m 0' option requests that successive data points not be joined by line segments, and the `-S 1' option requests that plotting symbol #1 (a dot) be plotted at the location of each point. As you will see if you run this command, the heavy mass does not win against the lighter mass. Instead the machine oscillates non-periodically. Since the motion is non-periodic, the plot benefits from being drawn as a sequence of unconnected points.
We conclude by mentioning a few features of ode that may be
useful when things are not going quite right. One of them is the
`examine' statement. It may be used to discover pertinent
information about any variable in a system. For details, see section The ode input language formally specified.
Another useful feature is that the `print' statement may be used to print out more than just the value of a variable. As we have seen, if the name of the variable is followed by `'', the derivative of the variable will be printed instead. In a similar way, following the variable name with `?', `!', or `~' prints respectively the relative single-step error, the absolute single-step error, or the accumulated error (not currently implemented). These quantities are discussed in section Numerical error and how to avoid it.
The `print' statement may be more complicated than was shown in the preceding examples. Its general structure is
print <pr-list> [every <const>] [from <const>]
The bracket notation `[...]' means that the enclosed statements are optional. Until now we have not mentioned the `every' clause or the `from' clause. The <pr-list> is familiar, however; it is simply a comma-separated list of variables. For example, in the statement
print t, y, y' every 5 from 1
the <pr-list> is <t, y, y'>. The clauses `every 5' and `from 1' specify that printing should take place after every fifth step, and that the printing should begin when the independent variable t reaches 1. An `every' clause is useful if you wish to `thin out' the output generated by a `step' statement, and a `from' clause is useful if you wish to view only the final portion of a solution curve.
ode command-line options
The command-line options to ode are listed below. There are
several sorts of option:
The following option affects the way input is read.
The following options affect the output format.
The following options specify the numerical integration scheme. Only one of the three basic option `-R', `-A', and `-E' may be specified. The default is `-R' (Runge--Kutta--Fehlberg).
The following options set the error bounds on the numerical solution scheme.
ode to
continue even if this ceiling is exceeded. This may result in large
numerical errors.
Finally, the following options request information.
ode and the plotting utilities
package, and exit.
ode is always in one of two states:
ode moves from the first to the second state after it sees and
processes a `step' line. It returns to the first state after
the generated output has been printed. Errors may occur in either the
`reading' state or the `solving' state, and may terminate computations
or even cause ode to exit. We now explain the possible sorts of
error.
While reading input, ode may encounter a syntax error: an
ungrammatical line that it is unable to parse. (For a summary of its
input grammar, see section The ode input language formally specified.) If so, it emits the error
message
ode::nnn: syntax error
where `nnn' is the number of the line containing the error. When the `-f filename' option is used to specify an input file, the error message will read
ode:filename:nnn: syntax error
for errors encountered inside the input file. Subsequently, when
ode begins reading the standard input, line numbers will start
over again from 1.
No effort is made to recover from syntax errors in the input. However, there is a meager effort to resynchronize, so that more than one syntax error in a file may be found at the same time.
It is also possible that a fatal arithmetic exception (such as a
division by zero, or a floating point overflow) may occur while
ode is reading input. If such an exception occurs, ode
will print an "Floating point exception" error message and exit.
Arithmetic exceptions are machine-dependent. On some machines, the
line
y = 1/0
would induce an arithmetic exception. Also on some machines (not necessarily the same ones), the lines
y = 1e100 z = y^4
@ifnottex
would induce an arithmetic exception. That is because on most
machines, the double precision quantities that ode uses
internally are limited to a maximum size of approximately 1.8x10^308.
When ode is in the `solving' state, i.e., computing a numerical
solution, similar arithmetic exceptions may occur. If so, the solution
will be interrupted and a message resembling
ode: arithmetic exception while calculating y'
will be printed. However, ode will not exit; the exception will
be `caught'. ode itself recognizes the following exceptional
conditions: square root of a negative number, logarithm of a
non-positive number, and negative number raised to a non-integer power.
ode will catch any of these operations before it is performed,
and print an error message specifying which illegal operation it has
encountered.
ode: square root of a negative number while calculating y'
would be a typical error message.
If the machine on which ode is running supports the
`matherr' facility for reporting errors in the computation of
standard mathematical functions, it will be used. This facility reports
domain errors and range errors (overflows, underflows, and losses of
significance) that could occur when evaluating such functions as
`log', `gamma', etc.; again, before they are performed. If
the matherr facility is present, the error message will be fairly
informative. For example, the error message
ode: range error (overflow) in lgamma while calculating y'
could be generated if the logarithmic gamma function `lgamma' is evaluated at a value of its argument that is too large. The generation of any such message, except a message warning of an underflow, will cause the numerical solution to be interrupted.
There is another sort of error that may occur during numerical solution:
the condition that an error ceiling, which the user may set with the
`-r' option or the `-e' option, is exceeded. This too will
cause the numerical solution to be abandoned, and ode to switch
back to reading input.
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.
@ifnottex
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.
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.
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:
@ifnottex
y' = y^2 y(1) = -1The 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' statement
step 1, -1on this system. How does
ode react?
As another example of ill-posedness, consider the system
y'=1/ywhich is undefined at y=0. The general solution is @ifnottex
y = +/- (2(t-C))^(1/2),@ifnottex 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.
y' = 2x x' = 2yhas 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).
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) = 1on the domain [0,1]. The exact solution is @ifnottex
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.
# 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, 5Next 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 -CThis 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 .001and a plot of the solutions with the specified error bounds will be drawn. The convergence, showing stability, should be quite illuminating.
The time required for ode to solve numerically a system of
ordinary differential equations depends on a great many factors. A
few of them are: number of equations, complexity of equations (number
of operators and nature of the operators), and number of steps taken
(a very complicated function of the difficulty of solution, unless
constant stepsizes are used). The most effective way to gauge the time
required for solution of a system is to clock a short or imprecise run
of the problem, and reason as follows: the time required to take two
steps is roughly twice that required for one; and there is a
relationship between the number of steps required and the relative error
ceiling chosen. That relationship depends on the numerical scheme being
used, the difficulty of solution, and perhaps on the magnitude of the
error ceiling itself. A few carefully planned short runs may be
used to determine this relationship, enabling a long but imprecise run
to be used as an aid in projecting the cost of a more precise run over
the same region. Lastly, if a great deal of data is printed, it is
likely that more time is spent in printing the results than in computing
the numerical solution.
ode input language formally specified
The following is a formal specification of the grammar for ode's
input language, in Backus--Naur form. Nonterminal symbols in the
grammar are enclosed in angle brackets. Terminal tokens are in all
capitals. Bare words and symbols stand for themselves.
<program> ::= ... empty ...
| <program> <statement>
<statement> ::= SEP
| IDENTIFIER = <const> SEP
| IDENTIFIER ' = <expression> SEP
| print <printlist> <optevery> <optfrom> SEP
| step <const> , <const> , <const> SEP
| step <const> , <const> SEP
| examine IDENTIFIER SEP
<printlist> ::= <printitem>
| <printlist> , <printitem>
<printitem> ::= IDENTIFIER
| IDENTIFIER '
| IDENTIFIER ?
| IDENTIFIER !
| IDENTIFIER ~
<optevery> ::= ... empty ...
| every <const>
<optfrom> ::= ... empty ...
| from <const>
<const> ::= <expression>
<expression> ::= ( <expression> )
| <expression> + <expression>
| <expression> - <expression>
| <expression> * <expression>
| <expression> / <expression>
| <expression> ^ <expression>
| FUNCTION ( <expression> )
| - <expression>
| NUMBER
| IDENTIFIER
Since this grammar is ambiguous, the following table summarizes the precedences and associativities of operators within expressions. Precedences decrease from top to bottom.
Class Operators Associativity Exponential ^ right Multiplicative * / left Additive + - left
As noted in the grammar, there are six types of nontrivial statement. We now explain the effects (the `semantics') of each type, in turn.
"y" is a dynamic variable value:2.718282 prime:2.718282 sserr:1.121662e-09 aberr:3.245638e-09 acerr:0 code: push "y"The phrase `dynamic variable' means that there is a differential equation describing the behavior of y. The numeric items in the table are:
The grammar for the ode input language contains four types of
terminal token: FUNCTION, IDENTIFIER, NUMBER, and
SEP. They have the following meanings.
gnuplot. All functions take a
single argument, except for ibeta, which takes three, and
igamma, which takes two. For trigonometric functions, all arguments
are expressed in radians. The atan function is defined to give a
value between -PI/2 and PI/2 (inclusive).
In the ode input language, upper and lower-case letters are
distinct. Comments begin with the character `#' and continue to
the end of the line. Long lines may be continued onto a second line by
ending the first line with a backslash (`\'). That is because
the combination backslash-newline is equivalent to a space.
Spaces or tabs are required in the input whenever they are needed to separate identifiers, numbers, and keywords from one another. Except as separators, they are ignored.
ode and solving differential equationsode.
ode: A
numerical simulation of ordinary differential equations,"
pp. 480--481 in Proceedings of the Conference on Computers in
Physics Instruction, Addison--Wesley, 1990.
Go to the first, previous, next, last section, table of contents.