Path: utzoo!utgpu!news-server.csri.toronto.edu!cs.utexas.edu!sun-barr!lll-winken!ncis.tis.llnl.gov!dog.ee.lbl.gov!usenet From: austern@ux5.lbl.gov (Matt Austern) Newsgroups: comp.sys.handhelds Subject: Bulirsch-Stoer integration on the HP 48 (Part 1 of 2) Summary: Documentation, description, and examples. Keywords: Bulirsch-Stoer, differential equations Message-ID: <8668@dog.ee.lbl.gov> Date: 18 Dec 90 18:13:28 GMT Sender: usenet@dog.ee.lbl.gov Reply-To: austern@ux5.lbl.gov (Matt Austern) Organization: Lawrence Berkeley Laboratory (theoretical physics group) Lines: 142 X-Local-Date: Tue, 18 Dec 90 10:13:28 PST 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.