8.3 Additional examples using 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 ode Invocation.

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 æ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 Input Language.

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 Numerical Error.

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.