Ode - Numerical Solution of Ordinary Differential
Equations
Nick Tufillaro*
Graham A. Ross
Reed College
Portland, Oregon 97202
ABSTRACT
Ode solves the initial value problem for a
family of first-order differential equations when
provided with an explicit expression for each
equation. Ode parses a set of equations, initial
conditions, and control statements and then pro-
vides an efficient numerical solution. Ode makes
the initial value problem easy to express; for
example, the ode program,
# an ode to Euler
y = 1
y' = y
print y from 1
step 0, 1
prints 2.718282.
The Ode User's Manual contains a guide to applying
ode and a discussion of its design and implementa-
tion.
1 January 1982
Copyright Nicholas B. Tufillaro, 1982-1994.
All rights reserved.
_______________________________________________________
* The first author's support from National Science
Foundation CAUSE Grant No. SER 78-06405 is gratefully
acknowledged.
Ode - Numerical Solution of Ordinary Differential
Equations
Nick Tufillaro*
Graham A. Ross
Reed College
Portland, Oregon 97202
1. Introduction
Ode renders a numerical solution to the initial value
problem for many families of first order differential equa-
tions. It is a programming language which resembles the
mathematical language so that the problem posed by a system
is easy to state, thereby making a numerical solution
readily available. Ode solves higher-order systems since a
simple procedure converts an nth' order equation into n
first order equations (see Section 3.2). Three distinct
methods are implemented at present: Runge-Kutta-Fehlberg
(default), Adams-Moulton and Euler. The Adams-Moulton and
Runge-Kutta routines are available with adaptive step size.
The Ode User's Manual falls into two parts. The first
is a tutorial which begins with simple examples and
progresses to show how ode is used with other programs, not-
ably graph(1) and plot(1). Examples also demonstrate that
ode cannot be used blindly! Correct results from ode
presume that the user is aware of certain information about
the particular system of equations, for instance the loca-
tion of any singularities. A discussion of interpreting
results and the rudiments of error analysis are also found.
The second part begins at Section 5; it contains a con-
cise statement of the linguistic and numerical aspects of
ode, and a discussion of its design and implementation.
This is the definitive statement of ode, with the exception
of the code itself. This second part also contains some
information which would be useful in maintaining or modify-
ing ode.
1.1. What is the numerical solution of an ODE?
First a blizzard of definitions. A time honored defin-
ition says a differential equation is an equation involving
an unknown function and its derivatives. A differential
equation is ordinary if the unknown function depends only on
one independent variable, often denoted 't'. The order of
- 2 -
the differential equation is the order of the highest-order
derivative in the equation. One speaks of a family or sys-
tem 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. This "initial condition" specifies a
unique solution. Questions about the existence and unique-
ness of a solution, along with further terminology, are
found in any introductory text (cf. Birkhoff/Rota, Chapter
1).
In practical problems the solution of a differential
equation is usually not known in terms of familiar functions
or it is, as the buzzword goes, intractable. Hence the need
for a numerical solution.
A numerical method for a system of ordinary differen-
tial equations produces approximate solutions at particular
points using only the operations of arithmetic and func-
tional evaluations of the given equations. Ode solves first
order initial value problems of the form:
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,
x(a) = b
y(a) = c
.
.
.
z(a) = d
t = a
where a,b,c,...,d are constants.
In the equations shown above f, g, ..., and h can be
expressed using the usual operators (+,-,*,/,^) and func-
tions listed in section 8.3. In order to solve a system of
equations with ode, it must be possible to solve explicitly
for each derivative. Ode can't solve an equation in which
the first derivative is implicitly defined. For instance,
ode can't solve the system,
y' + y' = 1.
All numerical methods involve the calculation of an
- 3 -
approximate solution at discrete intervals of the indepen-
dent variable 't', where the difference between any two suc-
cessive steps (usually denoted 'h') may be constant or adap-
tive. In general, as 'h' decreases the solution becomes
more accurate. Again, the solution rendered by ode is an
approximation to the actual solution.
Reflection upon the geometric fact expressed by a sin-
gle first order differential equation gives one insight
about how a method works. The set of all solution curves
expresses a vector field or flow. A numerical solution sim-
ply samples the direction of the flow in small discrete
steps and moves accordingly. This is the heart of Euler's
method.
2. Simple examples
Let's begin with simple problems that will illustrate
the essential elements of composing an ode while not getting
bogged down in too many details. The entire specification
for the language may be found in Section 5; look there if
these examples are too elementary.
2.1. Euler
The first ode to compose is for the simplest example
imaginable, to wit,
y'(t) = y(t)
with the initial condition,
y(0) = 1
The solution to this differential equation is a function
with the property that the first derivative of the function
is equal to the function itself. But that's easy, it's
just,
y(t) = e^t,
and in particular,
y(1) = e^1 = 2.718282
which is Euler's constant to seven places.
A similar result is obtained with ode by the sequence
of instructions
ode
y' = y
y = 1
print t, y
step 0, 1
Two columns of numbers will appear. The first is the value
of the independent variable 't', and the second is the value
- 4 -
of the dependent variable 'y'. The last line should be
1 2.718282.
After this output ode will wait for further instructions.
Entering
step 1, 0
should bring you back to where you started, more or less.
To exit ode type a line containing only a '.', and then tap
return. Your shell prompt should appear. Ode also exits
when it sees an end-of-file on its command stream, which you
can send from the terminal by typing ctrl-d.
Exercise 2.1. Run this ode again, except leave out one
or more of the lines. What error messages do you get? If
the program works, what are the default values. For
instance, what does this program do?
ode
y = 1
y' = y
.
Don't forget the period and then tap return.
Each line of the ode is pretty much self-explanatory.
The line
ode
starts up ode, which is now ready to receive information
about the problem. The line
step 0, 1
is the step statement. It sets the beginning and the end of
the domain for the independent variable and causes ode to
perform the numerical step. If it's missing, ode does noth-
ing -- but at least it does it right! The first value of
the step statement (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 = 1 sets
the initial condition for y, which must always be a con-
stant, whereas y' = y defines a function. Whenever a step
statement is encountered the program tries to produce a
solution. The print statement tells what values are to be
printed and in what order; for instance,
print y, y'
- 5 -
prints the value of the function y and its derivative.
Since our solution of this problem is numerical it is
only an approximation. Unlike the exact solution, one must
step from the value of the initial condition of the indepen-
dent variable to the point desired. The actual step, as you
no doubt noticed, takes place by means of a number of
smaller steps; the step size, often denoted 'h', is chosen
to minimize the error. The result is an approximate solu-
tion, an integral curve, on the interval given by the step
statement.
Let's make a file containing the following:
# an ode to Euler
y = 1
y' = y
print t, y, y'
Name this file euler. The '#' is the comment statement,
which may appear anywhere in an ode. Everything from its
appearance to the end of the line is ignored. To execute
the program say
ode -f euler
step 0, 1
and the file euler is read into ode. Ode then waits for any
further statements. Notice that
ode < euler
is not equivalent to
ode -f euler.
The former causes ode to get all its input from the file
euler, while the latter allows subsequent input from the
terminal. (See "Introduction to the UNIX Shell" for more
information on the subject of input/output redirection.)
Exercise 2.2. Ode can take command-line flags to tell
it to do something special. The '-f' flag must always be
followed by a file name; it is a request to read the con-
tents of the file into ode. Use ode with the '-t' flag,
ode -t
What does this do? How about "ode -p 10"? "ode -r 1e-12"?
2.2. Cosine
Consider the differential equation and initial condi-
tions,
- 6 -
y''(t) = -y(t)
y'(0) = 0
y(0) = 1
The solution is
y(t) = cos(t)
In order to solve this problem using ode we must
express this second order equation as two first order equa-
tions. Toward this end perform the following change of
variables. Let,
yprime' = -y.
y' = yprime
These two equations are the same as the single equation
above since,
yprime' = y'' = -y.
The reduction of an 'n-th' order problem to 'n' first order
problems is a standard trick in this game which you will use
often. It may be helpful to think of y' as velocity. Make
a file containing
# cosine : y''(t) = -y(t), y(0) = 1, y'(0) = 0
cosine' = -sine
sine' = cosine
cosine = 1
sine = 0
print t, cosine
Name the file cosine. Now type
ode -f cosine | graph | plot
step 0, 2*PI
.
And sit tight for a few seconds. A graph of the solution
curve will soon appear. Documentation for graph(1) and
plot(1) is found in the UNIX Programmer's Manual, Volume 1.
The '(1)' refers to the section number.
Exercise 2.3. Use ode to solve y''(t) = -y(t) with
y(0) = 0, y'(0) = 1.
Exercise 2.4. While using a graphics terminal try:
ode -f cosine | graph -m 0 | plot
step 0, 2*PI
.
Look up the '-m' flag in graph. What does it do? Do you
see how the step size changes?
- 7 -
In the next section we will deal a bit more with how to
reduce an 'n-th' order problem and with how to make a graph
of the solution curve.
Ode knows a number of functions, to wit, abs, sqrt,
sin, cos, atan, exp, ln, sinh, cosh, and tanh. As you see,
'PI' is a synonym. The names of the functions listed
above are reserved, so, for instance, we couldn't have used
'cos' and 'sin' for our variables in the previous example.
Exercise 2.5. Investigate the differential equation
y'(t) = cos(y)
with ode for different initial conditions.
What is the actual solution? Do the same for
y'(t) = cos(ty).
Exercise 2.6. Write down eighteen first order dif-
ferential equations for the infamous 'three body problem' of
physics and find 'stable' configurations using ode.
3. Advanced Topics
Section 3 highlights some of the finer points of ode.
However this section along with the previous one does not
cover all of ode's capabilities. Read Section 5 after you
have gained some experience with ode to find out about them.
Use a graphics terminal for the following exercises.
Or "xterm" which is capable of tektronix's 4104 emulation.
3.1. Use of Graph and Plot
Look up and read the documentation for graph(1),
plot(1), and spline(1). Graph and plot are general purpose
programs of great pliancy. Graph indulges in copious flags.
Proper use of them allows one to display multiple graphs
with labels on a single screen. We will focus on a few of
the features that often arise while using graph with ode.
Graph has a nasty peculiarity that you should look out
for: it cannot graph more than about 2000 points. This
defect may be circumvented by a careful use of the print
statement. The general structure of the print statement is
print [every ] [from ].
The bracket notation [], means that the enclosed statements
are optional. Up till now we have not used the every state-
ment and the from statement. The '' is a comma-
separated list of variables that occurs in the ode. For
instance in
- 8 -
print t, y, y' every 5 from 1
the is . The statements 'every 5' and
'from 1' do what you hopefully suspect. In this case every
fifth value is printed and the printing will begin when the
independent variable 't', is equal to 1. The use of the
statements every and from to dodge the 2000 point limit is
obvious. The from statement is also handy where we want to
view only the last part of the solution curve.
Exercise 3.1 Use the cosine file of Section 2. Try:
ode -f cosine | graph | plot
print t, cosine from 2*PI
step 0, 4*PI
.
The default for the every statement is 'every 1'. Next try:
ode -f cosine | graph | plot
print t, cosine every 3 from 2*PI
step 0, 10*PI
.
Try other print statements. How big does the interval have
to be until graph silently stops plotting points? Can you
elude the limit with the every statement?
A valuable graph flag is '-m'. It stands for mode. In
particular,
graph -m 0
means print only points rather than connecting them to form
lines. It is possible to draw multiple graphs with this
option.
Exercise 3.2 Try:
ode -f cosine | graph -m 0 | plot
step 0, 2*PI
print t, sine
step 2*PI, 0
.
The ability to draw multiple graphs is big-league. It can
be used to construct a phase plane portrait and to check on
the stability of the solution. More about this a little
later. An ode may contain any number of print statements so
that you can draw multiple graphs. The most recent print
statement is used with each step statement.
Exercise 3.3. A Strange Attractor. An eloquent plot
- 9 -
is given by the Lorentz system.
# Lorentz model, a third order system
# illustration of a "strange attractor"
# interesting cases are r = 26, x = z = 0,
# y = 1, 2.5 < t < 30
# and r = 17, x = z = 0,
# y = 1, 1 < t < 50, print from 15
# Poincare plot is x vs. y
# use with ode -f lorentz | graph -m 0 | plot
r = 26
x' = -3 * ( x - y)
y' = -x * z + r * x - y
z' = x * y - z
x = 0
y = 1
z = 0
print x, y
In this example 'r' is a symbolic constant. Put this ode in
a file named lorentz and construct the suggested plots.
3.2. Solution of Higher Order Systems
A system of differential equations of various orders
can be transformed into a larger system of first order equa-
tions by systematically introducing new variables. This
process of reduction is not unique, in general it is possi-
ble to achieve the reduction in many ways; nevertheless,
there is a systematic way to go about it. Let's transform
the initial value problem
x''' = cos(x) + sin(x') - ex'' + t^2
x''(0) = 13
x'(0) = 7
x(0) = 3
into a form suitable for solution by ode.
Old Var. New Var. Initial Value Diff. Eq.
t t 0 t'=1
x x 3 x'=y
x' y 7 y'=z
x'' z 13 z'=cosx+siny-e+t^2
Exercise 3.4. Solve this system using ode. The
expression 't^2' is valid for ode so you don't need to say
't*t'; in general ode knows 't^u', to raise t to the u-th
power.
The first order equations are often significant. In
physical problems they may be Hamilton's equations of
motion.
- 10 -
Exercise 3.5. Reduction is not unique. Newton's equa-
tion of motion for a simple pendulum is
m l^2 a'' + m g l sin(a) = 0
where 'l' is the length of the string and 'a' is the angle.
Reduce this to two first order equations by the method
illustrated above. Do this before proceeding to the next
paragraph if you know what's good for you.
The solution is
(NOTE: DUE TO FORMATING DIFFICULITIES IN ASCII CONVERSION
THERE MAY BE SOME ERRORS IN THE MATHEMATICAL SECTIONS OF
THE TEXT. NBT: 4 SEPTEMBER 1994. A PROPERLY TYPESET MANUAL
MAY BE SENT TO YOU BY DIRECTLY CONTACTING NICK: nbt@reed.edu).
a' = w
w' = l - g sin(a)
What is the physi-
cal meaning of both sets of equations? How are they related
(what's C)? Chose reasonable initial conditions and solve
both sets of equations with ode.
3.3. Phase Plane Construction
The construction of a phase plane illustrates how ode
may be used in an interactive fashion. The idea is to draw
a set of integral curves from different initial conditions
in the quickest manner possible.
If we have two first order equations then the family
of all solution curves is called the 'phase portrait'; it is
of paramount interest in the qualitative theory of differen-
tial equations and also possesses an aesthetic appeal.
In practice the phase portrait is drawn from a few
well-chosen solution curves. Picking a good set of solution
curves may take some experimentation. The construction of
a phase portrait is illustrated in an example from popula-
tion biology.
Exercise 3.6. A predator-prey problem based on
Volterra's model. 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. A crude model
for the interaction of A and B yields the equations of popu-
lation:
x' = x(a-by)
y' = y(cx-d)
for a, b, c, d > 0 and (0,0) = (b/d,a/c).
- 11 -
It can be shown that if the solution curves are closed, the
fluctuations of x(t) and y(t), starting from any initial
population (greater than zero), are periodic. To draw a
phase portrait try:
ode | graph -m 0 | plot
# an example from population biology
# construction of a phase portrait for
# Predator-prey equations of Volterra and Lotka
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
.
The plot will not begin until the terminating dot is
entered.
3.4. Obscurata
There are a few frills that can come in handy when
things aren't going quite right. For example, The examine
statement can be used to discover pertinent information
about any variable in a system. It's described in section
8.2
The print statement can be used to print out more than
just the value of a variable. If the name of the variable
is followed by a "'", the derivative of the variable is
printed instead. Following the variable name with "~", "!",
or "?" prints the accumulated error (not implemented in the
current version of ode), the absolute single-step error, or
the relative single-step error, respectively.
4. Error
Be warned, this section is incomplete. Entire books
exist on any single subject mentioned (e.g., floating point
error). The goals of this section are modest; first, to
introduce the most basic notions of error analysis as they
relate to ode; second, to steer you around a few of the more
obvious pitfalls. Look through a numerical analysis book
- 12 -
before beginning this section (cf. Atkinson).
4.1. Basics
A cornucopia of definitions whose relevance will soon
be apparent. The error of greatest concern is the differ-
ence between the actual solution and the approximate solu-
tion; this is termed the global or accumulated error, since
the error is built up during each step. Unfortunately, an
approximation of this error is usually not available without
knowledge of the actual solution or some hard core numerical
analysis. Thus there arise a number of more workable
notions of error. The local or single step error is the
difference between the actual solution and the approximate
solution at any single step, assuming the value at the
beginning of the step is correct.
The relative error is the error divided by the approxi-
mate value. Why not the actual value? Well, this is not
known, so it would be a useless definition. Ode, if left to
itself, will track the relative single step error in order
to pick the step size. By default, it maintains an accuracy
of eight digits on any single step. This accuracy can be
adjusted with the '-r' flag.
The sources of error are generally twofold. The first
is due to the limits of machine computation. All computers
work with floating-point numbers, which are not real
numbers, but only an approximation to real numbers. You can
see the difference between the two by using the '-p 17' flag
which will print 17 significant digits. This floating point
error is propagated through each step. However, all calcu-
lations made by ode are done to double precision, so this
tends to be a small source of error.
The second source of error is often called the theoret-
ical truncation error; it is the difference between the
actual solution and the approximate solution due solely to
the numerical method. At the root of many numerical methods
is an infinite series; for ordinary differential equations,
it is a Taylor expansion. Since the computer can not take
all the terms in an infinite series, a numerical method uses
a truncated series; hence the origin of the term. The accu-
mulated error is a sum of the theoretical truncation error
and the floating point error, although in practice the
floating point error is seldom included. Some authors also
term this the actual truncation error. The single step
error for ode consists only of the theoretical truncation
error.
No term in numerical analysis is more democratic than
stability. We will say a numerical method is stable (when
applied to a particular system) if the relative error
decreases as the step size decreases. Both the '-R' and the
- 13 -
'-A' algorithms are stable (based more on empirical results
than theoretical) for a wide class of systems.
4.2. Pitfalls
Here are some common sources of trouble to watch out
for.
4.2.1. Singularities
Look up the definitions of ordinary point, singular
point, regular singular point, and irregular singular point
in your favorite differential equations text. If you don't
have a favorite, try mine, Birkhoff and Rota, Chapter 9.
Ode should not be asked to step over singularities. It
should not be applied to singularities except possibly at an
endpoint. Always attempt to locate and classify all the
singularities for a system.
4.2.2. Ill-Posed Problems
Be careful when chosing the domain (step interval) for
a problem. You must know that the system and the solution
have a well-defined value on the domain. The solution must
also be real. Ode cannot pole-vault complex singularities
from the real line. Whenever you get weird results, e.g.,
the step size suddenly shrinks to the machine limit or the
solution suddenly blows up, check to make sure the problem
is not ill-posed. Keep your eyes peeled for those complex
poles.
As an example of an undefined solution consider the
innocent looking system:
y' = y^2
y(1) = -1
The solution for t > 0 is
y(t) = t-1
With this problem you must avoid any step statement which is
in the neighborhood of t = 0. Try your luck with the step
statement
step 1, -1
for this system. How does ode react?
As a second example consider a system which is unde-
fined for y = 0:
The general solution is
- 14 -
and if y(2) = 2 then
Clearly any step statement which includes negative values
will produce bogus results. Also ode will not be able to
evaluate the system at any value for which y = 0. Try to
solve this system for some bogus intervals using ode. When
using a constant step size be careful not to to 'step over'
such trouble areas. Ode, with an adaptive step size, will
often spot such problems, but not always. Caveat emptor.
4.2.3. Critical points
An autonomous system is one that does not include the
independent variable explicitly in the derivatives. A crit-
ical point for a system is a point such that each derivative
of the family is exactly zero. For instance the system
y' = 2x
x' = 2y
has only one critical point at (x,y) = (0,0).
A critical point is sometimes referred to as a stagna-
tion point. This is because a system at a critical point
will stay there forever, although a system near a critical
point can undergo more violent action.
Exercise 4.1. Solve the system above using ode with
the initial condition x(0) = y(0) = 0. The solution should
be a single point. Now do the same with points near the
critical point. What happens?
Exercise 4.2. Find all the critical points of the
Lorentz system mentioned in exercise 3.3.
Always locate all the critical points of a system
before attempting a solution with ode.
Critical points are classified (e.g., equilibrium, vor-
tex, unstable, stable, etc.) and this knowledge may be of
use; to find out more about this look in a book dealing with
the qualitative theory of differential equations (cf.
Birkhoff/Rota, Chapter 6).
4.2.4. Unsuitable Numerical Method
Lastly, if the results are bad, it just might be that
the numerical method is not suitable to the problem. For
instance, ode at present has no numerical method which han-
dles the so-called 'stiff' problems very well.
Exercise 4.3. Examine the 'stiff system':
- 15 -
for 0 < t < 10. The exact solution is
Solve with ode using various methods, step sizes, and rela-
tive single step error tolerances. Compare with the actual
solution.
4.3. Checks on Error
This section provides some rough and ready heuristic
checks on the accuracy of the solution. This is what they
really do.
4.3.1. Examine stability of solution curves: do they con-
verge?
Use of ode for any serious study should include a check
on stability. Always check how changing the step size or
the tolerance (with the '-r' flag) affects the solution
curve. As the step size decreases, the curve should con-
verge. If it does not then either the step size is not
small enough or the numerical method is unsuitable with
respect to the problem. So apply the following check.
If using an adaptive step size, superimpose the solu-
tion curves for various tolerances (e.g., -r 1e-9, -r 1e-11,
-r 1e-13, ...). If the curves converge then the problem is
to all appearances stable and your accuracy sufficient.
If employing a constant step size, then perform a simi-
lar analysis by halving the step.
Exercise 4.4. A Q.C.D. equation, courtesy of Ray
Mayer. Make a file named qcd containing:
# an equation arising in quantum theory
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:
- 16 -
: sserr is relative single step error limit
for sserr
do
ode -r $sserr < qcd
done | graph -m 0 | plot
The file stability is a shell command file. It is a program
that will superimpose curves of different tolerances. To
execute the program type:
sh stability 1 1e-1 1e-2 1e-3 1e-4 1e-5
and a graph for each tolerance will be drawn.
4.3.2. Check constants of the system: are they constant?
Many systems will have invariant quantities. For
instance, if the problem is a conservative physical system
then the total energy (according to the latest theory) is
constant. Knowledge about the qualitative behavior of any
variable can be used to check the quality of the solution.
4.3.3. 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 solution in
question. If they diverge sharply - that is, if two solu-
tions start out very close and end far apart - then the
quality of the solution may be dubious. On the other hand,
if the curves do not sharply diverge then any error will in
all likelihood not increase by more than an order of magni-
tude. Such problems are sometimes called well-conditioned.
5. Using Ode
This section describes how to use ode on UNIX. It
describes the diagnostic messages produced and the command-
line arguments.
5.1. Command Line Arguments
The prototypical ode command line is:
ode [-f file] [-ARE [step]] [-h hmin [hmax]] [-r rmax
[rmin]] [-e emax [emin]] [-s] [-p prec] [-t]
The following text makes reference to this model in naming
the arguments on the command line.
Ode reads its program from the standard input. A piece
of an ode can be stored in a file. If the -f flag is given,
then the file is read before the standard input.
Only one of the -A, -R, or -E flags may be given. Each
flag specifies a different numerical method. If an optional
- 17 -
step value is present then a constant step size is employed,
instead of an adaptive one.
If -E is specified, then a 'quick and dirty' Euler
scheme is used. The default step value for Euler runs is
0.1. Not recommended for serious applications.
If -R is specified, then a fifth-order Runge-Kutta-
Fehlberg algorithm is employed with an adaptive step size
(default). When a constant step size and no error analysis
are requested, then a classical fourth-order Runge-Kutta
scheme is used.
If -A is specified with or without a step, then a
fourth-order Adams-Moulton predictor-corrector is used. The
Runge-Kutta-Fehlberg algorithm is used to get past values
when needed.
The -h flag indicates a lower and an optional upper
bound on the step size. The numerical routine will not let
the step size go below the value specified for hmin. The
default is to allow the step to shrink to the machine limit.
The value hmax, if given, specifies the maximum value of the
step size. It is useful in preventing the numerical routine
from skipping quickly over an interesting region.
The -r argument sets a ceiling for the 'relative
single-step error'. The -e flag is for the 'absolute
single-step error'. The relative single-step error may
never exceed rmax (default is 10^-9) in any dependent vari-
able. If one does, the solution is abandoned and an
appropriate error message is printed. In a similar fashion,
an optional lower bound may be given to suggest when the
step size should be doubled (default is rmax/1000). Solu-
tions using an adaptive step size monitor the single-step
error to adjust the step size. Thus, a decrease in the
single-step error will cause smaller step sizes to be
picked.
The -s flag suppresses the ceiling on single-step error
so that the program continues if it reaches these limits.
This can result in large numerical errors.
The -p argument causes numeric results to be printed
with prec significant figures. If this flag is given, all
output is in scientific notation.
The -t argument prints a title, naming the output
columns. With this argument, the default print format is
scientific notation.
5.2. Diagnostics
Most messages from the language processor are self-
- 18 -
explanatory. The biggest exception is 'syntax error', indi-
cating a grammatical error. Language error messages are of
the form
ode: nnn: message...
where 'nnn' is the location of the line containing the
error. When the -f flag is used, the phrase '(file)' fol-
lows the 'nnn' for errors encountered inside the file. Sub-
sequently, when ode begins reading the standard input, line
numbers start over from one.
An arithmetic exception arises when an illegal calcula-
tion is attempted. All versions of ode recognize the fol-
lowing exceptional conditions: square root of a negative
number, logarithm of a non-positive number, and negative
number to a non-integer power. In addition, the PDP-11 and
VAX versions of ode discriminate among the hardware faults,
giving: division by zero; floating point overflow (a finite
number becomes indistinguishable from infinity), and, in the
case of the PDP-11, floating point underflow (a number which
should not vanish becomes indistinguishable from zero) where
simpler machines report simply "arithmetic exception".
Because of a quirk in the VAX trap architecture it turns out
to be impractical to diagnose floating-point underflow.
When a calculation underflows, no diagnostic is produced and
the result is zero.
Run-time errors include conditions in which an error
ceiling is exceeded, as well as arithmetic exceptions. When
one occurs, a message describing the problem is emitted and
the solution is abandoned. If appropriate, an effort is
made to indicate which dynamic variable was being calculated
at the time of the exception. A typical message is
ode: real number too large while calculating y'
which indicates that some number involved with calculating
the derivative of y exceeded the machine's size limit for
floating-point numbers (about 1.7x10^38 on the PDP-11 and the
VAX).
6. Limits
6.1. Space Limitations
The data space used in solving an ode consists of both
fixed and dynamic regions. The fixed space is used effi-
ciently regardless of the complexity of the problem being
solved. The dynamic space contains the symbol table, the
expression lists, and several other data structures whose
sizes vary. The space for these structures comes out of a
common pool; they play against each other.
- 19 -
Every identifier in an ode requires 230 bytes; this
space can never be reclaimed. Each operator in an equation
costs 14 bytes. For the sake of speed, some operators are
expanded into a sequence of simpler operations, each of
which requires 14 bytes. In particular, some common powers
are reduced to sequences of square roots, squares, cubes,
and inversions. The extent of these transformations can be
seen with the help of the examine statement. The operator
space for an equation is reclaimed whenever a new equation
is defined for that particular dynamic variable. Each lexi-
cal token in the most complex statement costs 10 bytes, but
this space is reclaimed during the scan. Thus, the most
complex statement circumscribes this requirement. Each ele-
ment in a print statement costs 6 bytes. Upon execution of
the next print statement this space is reclaimed. The total
space available for these dynamic structures amounts to
about 45,000 bytes on the PDP-11. The limit is much higher
on the VAX.
In addition to these absolute space limitations, the
maximum depth of the operand stack used for evaluating
derivatives is thirty cells. This limit could conceivably
be exceeded for a very complex differential equation. The
solution is to simplify the equation if possible, and where
not, to recompile ode with a larger operand stack.
6.2. Speed Limitations.
The time required to solve an ode depends on a great
many factors. A few of them are: number of equations, com-
plexity 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 step
sizes 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 method being used, the diffi-
culty of solution, and perhaps on the magnitude of the error
ceiling itself. A few carefully planned short runs can 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 -- especially on a terminal -- it is
likely that more wall-clock time is spent printing the
results than solving the problem numerically.
7. Internal Design
7.1. Language Implementation Details
The input language is implemented with a Yacc-generated
- 20 -
parser and a Lex-generated lexical analyzer. The benefits
of this approach are great. The basic recognizer was com-
pleted in about two days. And adding comments to the
language required that exactly one line of the lexical
analyzer source code be changed.
Differential equations are compiled into linked lists
of stack operations, one list per dynamic variable. The
function evaluator (eval) traverses one of these lists,
applying the indicated operations to a small operand stack.
At the completion of the traverse, the new derivative is on
top of the stack. The routine pops and returns this result.
Another routine (field) is called by the numerical routines
whenever a re-evaluation of the gradient is required. It
calls eval for each dynamic variable in turn.
Changes to the language involve modifications to the
yacc grammar, the lex rules, and/or the semantics stored
with the grammar and rules. The steps required to add a new
built-in function, for example, are (aside from writing the
function itself):
(1) Invent a new terminal symbol name for the reserved
word.
(2) Add the new reserved word to the lex rules.
(3) Add the new productions to the right-hand side of
"expr" in the grammar, together with whatever semantic
action is required (using the included macros if possi-
ble).
(4) Define a new value for exper and add the stack code
to eval.
(5) Run make.
7.2. Numerical Methods
The numerical and language procedures were designed in
accordance with the Jesuit maxim. The major elements mani-
pulated by both are the symbol table, and the functions
field and printq. Addition of a new numerical routine
requires that the numerical routine be written, and perhaps
that some fields inside the symbol table entry be expanded
or added. None of these changes affect the performance or
compromise the integrity of the parser, although they might
change the space-cost of a dynamic variable. Similarly,
modifications to the lexical analyzer or the grammar only
affect the parse tables and related data structures which
have no effect on the numerical routines.
At each step the numerical routine fills each symbol
table entry with relevant values. Then a call to field
- 21 -
evaluates every differential equation and puts the result in
the symbol table element syrime. The evaluator gets
current values for the dynamics from the field syalue.
Both syalue and syrime are in some sense temporary vari-
ables. All correct results eventually find their way into
syal[0] and syri[0]; these are never touched by the
language processor. Thus the real interface between the
language and numerical procedures is merely syalue and
syrime. It is important to maintain meaningful values in
syalue and syrime, since a floating-point exception
might disrupt the control flow at any moment. It is recom-
mended that the pointer fsp be used in traversing the symbol
table, so that arithmetic exceptions can be diagnosed as
accurately as possible.
Fourth and fifth order methods were chosen as the best
compromise between speed and accuracy. The adaptive step
size has proved a boon in both areas.
7.3. Details of Error and Fault Recovery
Command-line errors are diagnosed at initialization
time and are generally fatal.
Errors in the input stream are diagnosed during the
lexical scan or at parse time. They elicit an appropriate
message. The error recovery method used in the parser is
the vanilla Yacc scheme. (see S. Johnson, "Yacc: Yet
Another Compiler-Compiler", 1978, pp. 16-18).
Run-time errors are detected by the function evaluator,
the numerical routines, or the floating-point unit hardware.
In each case, a message is printed describing the error and
ode returns to a known state awaiting input. The current
step is abandoned, but all variable values are left reflect-
ing the state of the machine just prior to the fault. All
of the numerical routines scan the symbol table computing
values for each dynamic variable, using a common pointer
known to the error recovery routines named fsp. This
pointer provides a way of identifying the offending dynamic
variable when a trap occurs. A small assembly language
assist allows convenient setup and polling of the floating-
point unit status words on the PDP-11. The same function is
performed on the VAX in C.
8. The Ode Input Language
This section provides a concise summary of the input
language to ode.
8.1. Formal Grammar
The following is a formal description of the grammar of
ode's input language. The notation is BNF, for Backus-Naur
- 22 -
form (cf. Wirth). Non-terminal symbols in the grammar are
enclosed in angle brackets. Terminal tokens are in all cap-
ital letters and are defined below. Bare words and symbols
stand for themselves.
::= ... empty ...
|
::= SEP
| IDENTIFIER = SEP
| IDENTIFIER ' = SEP
| print SEP
| step , , SEP
| step , SEP
| examine IDENTIFIER SEP
::=
| ,
::= IDENTIFIER
| IDENTIFIER '
| IDENTIFIER ?
| IDENTIFIER !
| IDENTIFIER ~
::= ... empty ...
| every
::= ... empty ...
| from
::=
::= ( )
| +
| -
| *
| /
| ^
| FUNCTION ( )
| -
| NUMBER
| IDENTIFIER
- 23 -
8.1.1. Operator Precedences
The grammar above is ambiguous. This table summarizes the
precedences and associativities of operators within expres-
sions. Precedences decrease from top to bottom.
Class Operators Associativity
Exponential ^ right
Multiplicative * / left
Additive + - left
8.2. Statement Semantics
There are five kinds of statements:
IDENTIFIER ' =
This defines a first-order differential equation. The
derivative of the IDENTIFIER is given by the .
If a dynamic variable doesn't appear on the left side of a
statement of this form, its derivative is assumed to be
zero. That is, it is a symbolic constant.
IDENTIFIER =
This sets the initial value of the IDENTIFIER to the current
value of the . Dynamic variables which have not
been initialized are set to zero.
step ,
step , ,
These statements describe the behavior of the independent
variable and cause the numerical routine to be executed.
The first is the independent variable's initial
value. The second is its final value. The third is a step
size; if given, it overrides the one given on the command
line (see Section 5). Usually, the step size isn't speci-
fied and it varies as the computation proceeds. However,
the step always terminates with each variable equal to the
correct final value. Another step statement may follow,
after changing any equation, initial condition, or state-
ment.
print [ every ] [ from ]
This statement controls the content and frequency of the
numerical output. A is a comma-separated list
of IDENTIFIERs, where any of the IDENTIFIERs might be fol-
lowed by an apostrophe, denoting the derivative, or a ques-
tion mark, denoting the relative single step error, or an
exclamation mark, denoting the absolute single step error,
or a tilde, denoting the accumulated error (not implemented
- 24 -
at present). The specified values are printed in the order
they are found. Both the every clause and the from clause
are optional. If the every clause is present, a printing
occurs every iterations of the numerical algorithm.
The default is to print on every iteration (i.e. 'every 1').
The first and last values are always printed. If the from
clause is present it means to begin printing when the
independent variable reaches or exceeds . The
default is to print starting with the initial conditions.
If there is no print statement, then the independent vari-
able and all dependent variables which have differential
equations associated with them are printed. The independent
variable is printed first; the dependent variables follow in
the order their equations were given.
examine IDENTIFIER
A table of interesting information about the named variable
is printed on the standard output when this statement is
executed. For example, if the statement
examine y
were encountered after execution of the example given in the
abstract of this paper, the output would be:
"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 simply means that there is a
differential equation describing the behavior of y. The
numeric fields are:
value Current value of the variable.
prime Current derivative of the variable.
sserr Relative single-step error for the last step taken.
aberr Absolute single-step error for the last step taken.
acerr Total error accumulated during the most recent step
statement. Not implemented at present.
The "code" section lists the stack operations required
to compute the derivative of "y" (somewhat reminiscent of a
reverse Polish calculator). The information may come in
handy in discovering whether the precedences in the dif-
ferential equation statement were interpreted correctly or
- 25 -
in determining the time- or space-expense of a particular
calculation. "Push y" means to load y's value on the stack,
which is all that's required to compute its derivative in
this case.
8.3. Lexical Issues and Definition of Terminal Tokens
In the 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 contin-
ued onto a second line by ending the first line with a
backslash ('\'). The combination backslash-newline is
equivalent to a space. Spaces or tabs are required whenever
they are needed to separate identifiers, numbers, and key-
words from one another. Except as separators, they are
ignored.
The terminal tokens referred to by the grammar in Sub-
Section 1 are described here.
FUNCTION
One of the words: abs, atan, cos, cosh, exp, ln, sin,
sinh, sqrt, tanh. These do just what you probably
expect. All angles are expressed in radians. The atan
function gives a value between -J/2 and J/2.
IDENTIFIER
A sequence of alphanumeric characters starting with an
alphabetic character. The first eight characters are
significant. Upper- and lower-case letters are dis-
tinct. In identifiers, the underscore (') is alpha-
betic. Function names and keywords cannot be used as
identifiers, nor can 'PI'.
NUMBER
A non-empty sequence of digits possibly containing a
decimal point and possibly followed by an exponent. An
exponent is 'e' or 'E', followed by an (optionally
signed) one- or two-digit number. All numbers and all
parts of numbers are radix ten. A number cannot con-
tain any white space. The special word 'PI' is a
number equal to not two.
SEP
A semi-colon or a (non-escaped) new-line.
8.4. The Independent Variable
The independent variable is identifiable as the only
variable which never appears on the left side of a differen-
tial equation or an initial value assignment. If there is
more than one such variable, an error message is printed and
the step is abandoned. If there is none, an independent
variable is invented and given the impossible name
- 26 -
'(indep)'.
9. References
Atkinson, K., An Introduction to Numerical Analysis, John
Wiley, New York, 1978. Pages 380-384 contain a nice review
of literature on the numerical solution of ordinary dif-
ferential equations.
Birkhoff, G., Rota, G., Ordinary Differential Equations,
Ginn and Company, New York, 1962.
Wirth, N., Algorithms + Data Structures = Programs,
Prentice-Hall, Englewood Cliffs, N.J., 1976.
Nicholas B. Tufillaro, Tyler Abbott, Jeremiah Reilly, An ex-
perimental approach to nonlinear dynamics and chaos, Addison-
Wesley, 1992.
10. Appendix: Odies but Goodies
# chemical kinetics simulation
# A + 2B <==> C --> D
# rate constants are:
# kf : A + 2B --> C
# kb : C --> A + 2B
# kd : C --> D
a' = kb*c - kf*a*b^2
b' = kb*c - kf*a*b^2
c' = kf*a*b^2 - kb*c - kd*c
d' = kd*c
c = 0
d = 0
a = 0.1
b = 1
kf = 1
kb = 1
kd = 1
print t,c every 10
step 0,10
- 27 -
# Two similar pendula coupled with a horizontal spring.
# illustration of "beating"
# if the motion is undamped then the equations of motion
# are approximated by
xone' = vxone
xtwo' = vxtwo
vxone' = -g*xone/lone + k/mone * ( xtwo - xone)
vxtwo' = -g*xtwo/ltwo - k/ltwo * ( xtwo - xone)
k = 1
g = 9.8
lone = 5
ltwo = 5
mone = 1
mtwo = 1
xone = 0.0
xtwo = 0.3
vxone = 0
vxtwo = 0
print t, xone
# bead sliding on a smooth circular wire of radius a:
# there is viscous damping so the bead should home in
# on the equilibrium point.
# equation of motion is
g = 9.8 # acceleration due to gravity
a = 1 # radius of hoop
w = 10 # angular velocity of hoop
b = 1 # damping coefficient
the' = vthe
vthe' = (w^2)*sin(the)*cos(the) - (g/a)*sin(the) - b * vthe
the = 0.1
vthe = 0
print t, the
# nonlinear system with a limit cycle
yone' = yone + ytwo - yone*(yone^2+ytwo^2)
ytwo' = -yone + ytwo - ytwo*(yone^2+ytwo^2)
yone = 1
ytwo = 2
print yone, ytwo
- 28 -
# A disk dynamo simulation for the earth's magnetic field.
# an attempt to simulate the dramatic switch in the polarity
# every eon or so.
# cf. M. Steele, B.A. thesis, Reed College, 1981 (Physics)
A = 1
V = 1
Q = 14.625
S = 5
w' = Q - z * y - V * w
y' = S * ( A * z - y)
z' = w * y - z
w = 1
y = 1
z = 1
print t, z
step 0, 10
# A damped driven harmonic oscillator.
# y'' = -k/m * y - R/m * y' + cos(w*t)
# if R^2 > 4km, overdamped
# if R^2 = 4km, critically damped
# if R^2 < 4km, damped
# this equation reduces to two first order equations:
y' = vy
vy' = -k/m * y - R/m * vy + cos(w*t)
y = 1
vy = 0
k = 1
m = 1
R = 0.5
w = 2*PI
print t, y
- 29 -
# coupled oscillator : two pendula simulation
m = 1
M = 1.0625
a = 0.5
adot = 0
l = 10
ldot = 0
ldot' = (m*l*adot^2 - M*9.8 + m*9.8*cos(a))/(M + m)
l' = ldot
adot' = -(9.8*sin(a) + 2*adot*ldot) / l
a' = adot
print l, ldot every so often
# A rumor spreads through a closed population of
# constant size N+1.
# At time t the total population can be classified
# into three categories:
# x persons who are ignorant of the rumor;
# y persons who are actively spreading the rumor;
# z persons who have heard the rumor but have
# stopped spreading it;
# If two persons who are spreading the rumor meet
# then they stop spreading it. The contact rate
# between any two categories is constant, u.
# Show that the equations
# x' = -u * x * y,
# y' = u * ( x*y - y*(y - 1) - y*z)
# give a deterministic model of the problem.
# Find the equation of the phase paths and sketch
# the phase diagram. Show that, when initially
# y = 1 and x = N, the number of people who
# ultimately never hear the rumor is s, where
# 2N + 1 -2s + N log( s/N) = 0
x' = -u * x * y
y' = u * ( x*y - y*(y-1) - y*( 100 + 1 - y - x))
y = 1
x = 100
u = 1
print t, x
- 30 -
# planetary orbit simulation
# two suns situated at (0,0) and (-5,0)
# one planet starting out at (1,0)
# x and y are positions
# vx and vy are velocities
vx' = -x/((x^2+y^2)^(3/2)) -(x+5)/(((x+5)^2+y^2)^(3/2))
vy' = -y/((x^2+y^2)^(3/2)) -y/(((x+5)^2+y^2)^(3/2))
y' = vy
x' = vx
x = 1
y = 0
print x,y every 10
# these values seem to give a nice orbit:
# vx = 0
# vy = 1.142
# a more exciting result can be obtained from:
vx = 0
vy = 1.165