Path: utzoo!utgpu!news-server.csri.toronto.edu!cs.utexas.edu!swrinde!zaphod.mps.ohio-state.edu!sdd.hp.com!ucsd!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 2 of 2) Summary: Code for integration of differential equations. Keywords: Bulirsch-Stoer, differential equations Message-ID: <8670@dog.ee.lbl.gov> Date: 18 Dec 90 18:15:24 GMT Sender: usenet@dog.ee.lbl.gov Reply-To: austern@ux5.lbl.gov (Matt Austern) Organization: Lawrence Berkeley Laboratory (theoretical physics group) Lines: 326 X-Local-Date: Tue, 18 Dec 90 10:15:24 PST 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 * @ = + 2*h*y'(x, y) 2 * 4 PICK + 4 ROLL @ = x + h, DROP 3 ROLL h + @ where 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 * @ = 0.5*(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.