austern@ux5.lbl.gov (Matt Austern) (12/19/90)
I have written a program for the hp 48 which solves systems of
differential equations using the Bulirsch-Stoer algorithm. This
posting is a description of the program, and the next posting is the code
itself.
SUMMARY OF THE METHOD
If you look at the code, it will immediately be apparent that this
algorithm is immensely more complicated than, say, Runge-Kutta.
Although this complexity makes a single step rather slow, the
advantage of Bulirsch-Stoer is that it is often possible to choose
extremely large step sizes while maintaining reasonable accuracy.
The Bulirsch-Stoer algorithm is described in Numerical Recipes, by
Press, Flannery, Teukolsky, and Vettering. (If you do any numerical
work and you don't already have a copy of that book, by the way, you
should get one.) I'll give a brief summary, just so that this
program won't be completely a black box.
If you have a system of equations of the form y' = f(x, y), where y
and y' are in general vectors, and the initial conditions (x0, y0),
then the problem is to find y(x1) for some x1. (Normally, we require
x1 - x0 to be small; here, we make no such requirement.) Define H =
x1 - x0.
We divide H into n subintervals, and then use n iterations of a
simple first-order method (specifically, the modified midpoint
method) to find y(x1). Define this value to be y(x1, h), where
h = H/n.
The Bulirsch-Stoer algorithm is to compute y(x1, h) for several
different values of h, and then extrapolate down to h=0. The
extrapolation algorithm provides an estimate of its accuracy; when
this accuracy is good enough, we return the answer. Convergence is
somewhat quicker than might be expected; it turns out that
y(x1, h) contains only even powers of h, so the extrapolation is
actually in h^2.
This algorithm can occasionally fail. The rational function
extrapolation might encounter a pole, in which case there's nothing
to do but quit. The program checks for this failure, so that at
least there won't be a mysterious blowup. (This is rare, but it
happened to me once. The solution is to use a different
extrapolation procedure in those cases. If I get energetic, I may
implement polynomial extrapolation for the 48sx as a fallback.) It
is also conceivable that the step size will have to be reduced so far
that x+h will be indistinguishable from x. I have never seen this
happen (my guess is that it would only happen when you're close to
some singular point in the solution of the differential equation, or
when you have a system of equations which vary on very different
scales), but there is code to check for that too.
THE PROGRAM:
There are actually four programs here: EXTRAP, NDESTEP, DESTEP1, and
DE. (There are also two objects that are intended for internal use.
They are named in lowercase.)
EXTRAP is the extrapolation routine. It takes two arguments, both
lists. Level 2 contains a list of x values, and level 1 a list of y
values. EXTRAP returns two values: in level 2 it returns y(0), and
in level 1 an estimate of the error in that prediction. (EXTRAP,
incidentally, is where most of the time gets spent. I've tried to
make the code in its inner loop efficient, but I'd welcome any
suggestions for speeding it up.)
DESTEP1 solves a single first-order differential equation. It takes
five arguments. In order, level 5 to level 1, they are: ydot,
tolerance, stepsize, x, y. ydot is a function: it takes x and y (on
the stack, in that order) and returns y'(x, y). Tolerance is the
required accuracy (i.e., if Delta-y is the error, then we demand
Delta-y / y < tolerance), and stepsize, x, and y are
self-explanatory. DESTEP1 returns the same information, to make it
convenient to take another step. ydot and tolerance are left
unchanged on the stack; the new step size is that recommended by the
algorithm (you don't have to follow the recommendation, of course!);
and x and y are the new values.
(Note: if you are too greedy with your step size and your tolerance,
DESTEP1 might have to choose a smaller step than you asked for. If
so, x will be incremented by the step actually used. This is
slightly unusual, though: Bulirsch-Stoer can handle most reasonable
requests.)
NDESTEP is just the same, except that y is a vector, and ydot must
take a vectorial y as an argument and return a vectorial ydot as a
result. (Most of the code in NDESTEP, actually, is identical to that
in DESTEP1. Merging them would be quite easy if you want to conserve
memory.)
DE is intended for interactive use; it's essentially just a
cosmetic shell for NDESTEP and DESTEP1. It takes and returns only
four arguments: tolerance is taken from the display format. (In STD
format, it uses a tolerance of 0.0001.) It chooses to call NDESTEP
or DESTEP1, depending on the type of argument given. It returns x
and y as tagged objects; if the user provided tags then those tags
are preserved, otherwise it uses the rather unimaginative labels "x"
and "y". Finally, NDESTEP and DESTEP1 use user flags; DE saves and
restores the old values.
EXAMPLE PROBLEM:
Solve the equation x^2 y''(x) + x y'(x) + x^2 y(x) = 0, subject
to the initial conditions y(0) = 1, y'(0) = 0. Specifically, find y(5).
Solution:
This is equivalent to the system
y1' = - y1/x - y
y' = y1,
with y(0) = 1, y'(0) = 0.
Put the following four entries on the stack:
\<< OBJ\-> DROP \-> x y1 y
\<< IF x 0 == THEN y NEG ELSE y1 x / y + NEG END
y1 2 \->ARRY
\>>
\>>
5
0
[0 1]
Then, with the display mode set to 3 FIX, hit DE. After 67 seconds,
I get the result that y(5) = -0.178, and y'(5) = 0.328.
In fact, the solution to this equation is a Bessel function, J0(x).
The tabulated answer is J0(5) = -0.177597, and J0'(5) = -J1(5) = 0.327579.
67 seconds is a long time, but we were able to span the entire
interval from 0 to 5 in a single step, and still get three digits of
accuracy.
Enjoy! There are obviously better tools than the HP 48 for serious
numerical calculations, but it's very convenient to have a quick way
of playing with a differential equation.
--
Matthew Austern austern@lbl.bitnet Proverbs for paranoids, 3: If
(415) 644-2618 austern@ux5.lbl.gov they can get you asking the wrong
austern@lbl.gov questions, they don't have to worry
about answers.austern@ux5.lbl.gov (Matt Austern) (12/19/90)
The values returned by BYTES are #60761d and 3026.5.
************************* CUT BELOW THIS LINE *************************
%%HP: T(3)A(R)F(.);
DIR
@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
@ DE takes 4 args: ydot, step, x0, y0. ydot is a function that computes
@ the derivative of y; it must take x and y on the stack (x in level 2, y
@ in level 1), and return y in level 1 of the stack. Step is the step size
@ that the user requests, and x0 and y0 are the initial conditions.
@ It returns the same four pieces of information, in the same positions
@ on the stack: the function ydot (unchanged), a recommendation for the
@ size of the next step, and the new values of x and y. Normally, the
@ new x will be x0 + step, but occasionally it will be necessary to
@ take a smaller step than the user requested.
@ x and y will always be returned as tagged objects. The user may provide
@ tags; otherwise, the tags will be "x" and "y".
@
DE
\<<
IF -49 FC? -50 FC? AND @ Start by getting precision from display.
THEN 4 @ Use 4 digits as default, if in STD mode.
ELSE
0
0 3 FOR I
2 * -48 I + FS? +
NEXT
END
NEG ALOG @ Convert 3 digits, for example, to 10^-3.
RCLF @ Store flags, since program uses them.
\-> ydot step x0 y0 tol flags
\<<
IF x0 TYPE 12 == @ Result will be displayed as a tagged
THEN x0 OBJ\-> SWAP @ object. If user has provided a tag,
ELSE "x" x0 @ use it; otherwise, use "x" as default.
END
IF y0 TYPE 12 == @ Now do the same for y.
THEN y0 OBJ\-> SWAP
ELSE "y" y0
END
\-> xt x yt y
\<< ydot tol step x y
IF y TYPE DUP 0 == SWAP 1 == OR
THEN DESTEP1 @ Call DESTEP1 if y is a real or complex
ELSE NDESTEP @ number, or NDESTEP otherwise. It is
END @ assumed that if y isn't a number it is
@ a vector; if user puts in something
@ weird, program will bomb.
yt \->TAG SWAP @ Tag the new y value.
xt \->TAG SWAP @ Tag the new x value.
4 ROLL DROP @ Get rid of the tolerance that is returned.
flags STOF @ Restore old flag settings.
\>>
\>>
\>>
@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
@
@ DESTEP1 takes 5 arguments: ydot, tolerance, stepsize, x, y. These
@ arguments mean the same thing here as with DE; we're just adding
@ tolerance. This is the maximum fractional error permitted in the
@ result; i.e., we must have (Delta y) / y < tolerance. (That's not
@ quite the way I write the test, though; y == 0.0 may be rare, but it has
@ been known to happen.)
DESTEP1
\<< { } { } 0
\-> ydot tol step x y htab ytab iter
@ ydot, tol, step, x, and y are the arguments,
@ and are documented above. htab is a table
@ of the subintervals used, and ytab is a
@ table of the result for each choice of
@ number of subdivisions.
\<< 20 CF @ Clear flag 20 to indicate that no
@ extrapolation has converged yet.
WHILE 20 FC? iter 15 < AND @ Loop until an extrapolation converges, or
REPEAT @ until we have done too many iterations.
divisions 'iter' INCR GET @ Find number of subdivisions for this
DUP step SWAP / @ iterations, and then find the size h of
IF DUP x + x == @ each subdivision. Check to make sure
THEN @ that it isn't so small that x+h = h.
DROP2 @ If it is, then clean up the stack and
"Stepsize underflow" @ signal an error.
DOERR
END
DUP x y ydot EVAL * y + @ If N is the number of subdivisions,
y x 4 PICK + 3 ROLL 4 ROLL @ we will now do N iterations of the
\-> h @ modified midpoint method. Most of
\<< 4 ROLL 1 - @ this is done on the stack, but what's
1 SWAP FOR n @ going on is the recursion relation
DUP2 ydot EVAL h * @ <next y> = <last y> + 2*h*y'(x, y)
2 * 4 PICK + 4 ROLL @ <next x> = x + h,
DROP 3 ROLL h + @ where <last y> is saved from before
SWAP @ the current iteration and where y
@ represents y at the current iteration.
NEXT @ The last iteration is a bit different:
DUP2 ydot EVAL h * @ <new y> = 0.5*(y + <last y>
@ + h * y'(x + h, y)).
+ SWAP DROP + .5 * @ Store this is table of y values by
1 \->LIST 'ytab' STO+ @ prepending it to the list.
h SQ 1 \->LIST 'htab' STO+
@ Store current value of h^2 the same way.
@ (Isn't it neat that STO+ works with lists?)
\>>
IF iter 1 \=/ @ Now do an extrapolation with our values of
THEN @ h^2 and y(x+step, h^2). We can't
@ extrapolate with a single point, though!
htab 1 6 SUB @ Only use the most recent 6 values; adding
ytab 1 6 SUB @ older ones isn't very helpful. (Why 6?
EXTRAP @ Folklore. You can try other values.)
IF OVER ABS tol * < @ Is the error in this extrapolation small
THEN 20 SF @ enough? If so, set flg 20 to mark success.
ELSE DROP @ If not, get rid of the extrapolated y value
END @ on the stack, and go to another iteration.
END
END
IF 20 FS? @ Did we get convergence? If so, choose a
THEN @ new step size and return. The criterion
@ for choosing a new step size is that we
@ want the 6th extrapolation to converge.
@ If it took more, then shrink the step size;
@ if it took fewer, then expand it. If we're
@ at exactly 6 or 7, expand or shrink
@ slightly. (Why 6? Folklore. Experiment
@ with different values if you like.)
CASE iter 6 == THEN 1.2 END
iter 7 == THEN .95 END 16
divisions iter GET / END
step *
SWAP x step + @ Now put the stack in the right order
SWAP tol 4 ROLLD @ so that arguments are returned the
ydot 5 ROLLD @ same way as they were given.
ELSE @ If control got to here, then the step
@ size that was requested was too large.
@ Shrink it by a large factor and try again
@ recursively.
ydot tol step 250 / x y
DESTEP1
END
\>>
\>>
@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
@
@ Arguments for NDESTEP are exactly the same as for DESTEP1, except that y
@ is a vector instead of a number. Similarly, of course, the function ydot
@ must accept arguments of x and y, with y a vector, and must return a
@ vector. There is no error checking for correct dimensionality.
@
@ The meaning of tolerance is also slightly different. This is the
@ maximum fractional error for *any* of the components of y. This
@ could therefore introduce large inefficiencies if you're trying to
@ solve a system of equations that vary on very different length
@ scales.
@
@ I'm not going to bother to comment this. The code here is identical
@ to DESTEP1, since + and * can take vectorial arguments, except for
@ the part that does the extrapolation and checks it for convergence.
@ All the hard work is done in the function lextrap, which I did comment; if
@ you read that, it should be completely clear what's happening in NDESTEP.
@
NDESTEP
\<< { } { } 0
\-> ydot tol step x y htab ytab iter
\<< 20 CF
WHILE 20 FC? iter 15 < AND
REPEAT
divisions 'iter' INCR GET
DUP step SWAP /
IF DUP x + x ==
THEN
DROP2
"Stepsize underflow"
DOERR
END
DUP x y ydot EVAL * y +
y x 4 PICK + 3 ROLL 4 ROLL
\-> h
\<< 4 ROLL 1 -
1 SWAP FOR n
DUP2 ydot EVAL h *
2 * 4 PICK + 4 ROLL
DROP 3 ROLL h +
SWAP
NEXT
DUP2 ydot EVAL h *
+ SWAP DROP + .5 *
1 \->LIST 'ytab' STO+
h SQ 1 \->LIST 'htab' STO+
\>>
IF iter 1 \=/
THEN
htab 1 6 SUB ytab 1 6 SUB tol
IF lextrap
THEN 20 SF
END
END
END
IF 20 FS?
THEN
CASE iter 6 == THEN 1.2 END
iter 7 == THEN .95 END 16
divisions iter GET / END
step *
SWAP x step +
SWAP tol 4 ROLLD
ydot 5 ROLLD
ELSE
ydot tol step 250 / x y
NDESTEP
END
\>>
\>>
@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
@
@ EXTRAP takes two arguments: x, y. Both are lists. Returns two
@ numbers: y(0) and y_err. y(0) is the value of y extrapolated to x=0,
@ and y_err is an estimate of the error.
@
EXTRAP
\<< DUP DUP2 SIZE DUP 3 PICK SWAP GET
\-> X Y C D N RES @ C and D are temporary variables (lists),
@ used in a recursion relation. N is the
@ number of entries, and RES will be the
@ final result.
\<< 1 N 1 - FOR COL @ Loop over "columns" in a 2-d tableux.
1 N COL - FOR I @ Loop over entries in the current column.
D I GET C I 1 + GET @ What we are doing here is computing the
\-> DI CI1 @ values C and D will have in the next
@ column. We only use one element of C and
@ one of D at a time, so squirrel them away
@ to save time.
\<< X I GET @ The recursion relation is
X I COL + GET / @ (x(i)/x(i+col))*D(i)*(C(i+1)-D(i))
DUP DI * CI1 - @ C'(i) = ----------------------------------
CI1 DI - DUP @ (x(i)/x(i+col)) * D(i) - C(i+1)
CI1 * SWAP DI * @
4 ROLL * ROT @ C(i+1) * (C(i+1) - D(i))
@ D'(i) = -------------------------------- .
@ (x(i)/x(i+col)) * D(i) - C(i+1)
@
IF DUP 0 == @ This algorithm occasionally fails; check
THEN @ for division by 0 which would signal that.
3 DROPN X Y
"Extrapolation failed"
DOERR
END
DUP ROT SWAP
/ 'C' I ROT PUT @ Store new values of C and D.
/ 'D' I ROT PUT
\>>
NEXT
'RES' @ The final result is given by
D N COL - GET @ the sum of D(N - col) for all col.
STO+ @ This in, in fact, only one possible choice;
@ any of several sums over C or D will work.
NEXT
RES D 1 GET ABS @ The last correction to the result is our
\>> @ error estimate.
\>>
@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
@
@ divisions is a vector containing the number of subintervals to use
@ in successive iterations. There are 15 elements here; this is the
@ maximum number of iterations that we will use. The numbers here are
@ from Numerical Recipes, and have no theoretical basis. Once again, they
@ are folklore; changing them will reduce or increase efficiency, but
@ won't give you wrong answers.
divisions
[ 2 4 6 8 12 16 24 32 48 64 96 128 192 256 384 ]
@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
@
@ lextrap is an internal function used by NDESTEP. Arguments are x, y,
@ tol. x is a list of numbers, y is a list of vectors, and tol is the
@ maximum error permitted. Extrapolates list of vectors to x=0, and
@ returns 1 in level 1 of the stack if extrapolation has sufficiently
@ small error, 0 if not. If the extrapolation succeeded, returns y(0)
@ in level 2, otherwise returns nothing but the 0 to signify failure.
@
lextrap
\<< OVER DUP SIZE SWAP
1 GET SIZE 1 GET 0
\-> xl yl tol N n comp @ xl is list of x values, yl is list of y
@ values, each of which is a vector.
@ N is number of elements in the lists, and
@ n is number of elements in each y vector.
@ comp is an iteration variable: which
@ component of y vector are we looking at?
\<< 21 SF @ No extrapolations have failed.
WHILE 21 FS? @ Quit when an extrapolation fails, or when
comp n < @ all extrapolations are done successfully.
AND
REPEAT
xl 'comp' INCR @ Look at the next component.
1 N FOR i
yl i GET @ Put the current component of all the
OVER GET SWAP @ y values on the stack,
NEXT @ ...
DROP @ ...
N \->LIST @ and assemble them into a list.
EXTRAP @ Now do an extrapolation with this list.
IF OVER ABS tol * >
THEN @ If this extrapolation didn't converge, then
comp DROPN @ clean up the stack, and clear flag 21 to
21 CF @ signal failure.
0
END
END
IF 21 FS? @ If we're done with this loop, and flag 21
THEN n \->ARRY 1 @ is still set, then all extrapolations have
END @ succeeded. The results are on the stack;
\>> @ turn them into a list.
\>>
END
************************* CUT ABOVE THIS LINE *************************
--
Matthew Austern austern@lbl.bitnet Proverbs for paranoids, 3: If
(415) 644-2618 austern@ux5.lbl.gov they can get you asking the wrong
austern@lbl.gov questions, they don't have to worry
about answers.