▼
Posts: 1,368
Threads: 212
Joined: Dec 2006
That guru of programs for the the HP41, JM Baillard, recently sent to me some routines to solve first and secondorder ODEs by the Bullrisch Stoer method, and I have some interest in adapting them for the 49G series, since, to be frank, I am really pretty underwhelmed with the builtin differential equation solver of my 48G and 49G+. Baillard's programs work really well on the 41series and are lightning fast when I run them (no modification necessary) in Free42 both on PC and my Palm TX, so I am interested to see how the algorithm fares on the 49G+
I am finding that translating RPN to RPL is not for the faint of heart, and I am wondering if anyone out there knows if someone has already devised BullrischStoer ODE solvers for the 49G series. I don't mind a little challenge, but I don't want to completely reinvent the wheel.
Les
▼
Posts: 348
Threads: 74
Joined: Jun 2013
Make sure you have checked the items at the link below.
http://www.hpcalc.org/search.php?query=differential
▼
Posts: 1,368
Threads: 212
Joined: Dec 2006
Thank you! I have done a more specific search of hpcalc.org and haven't come up with anything that specifically refers to BullrischStoer. Perhaps some of the package you have have found implement it as an option, but in that case I expect the author(s) would have made specific mention in the description. Will take a look at some of the documentation of the routines and see what turns up.
Les
▼
Posts: 191
Threads: 41
Joined: Jun 2007
Hi everyone,
I have not ( yet? ) written a BulirschStoer program for the HP48
( it's in the pipeline for the HP41 to solve systems of ODEs )
but I think that the builtin RKF is quite good.
In my RKF menu, there are
'X' for the variable
'Y' for the function(s)
'F' for the routine that must calculate dY/dX
'tol' to specify the tolerance
'DIF' that contains the short program
<< { X Y F } tol ROT RKF DROP2 Y >>
For example, to solve the system
dy/dx = yzu y(0)=1
dz/dx = x(y+zu) z(0)=1
du/dx = xyzu u(0)=2
Store 0 in 'X'
Store the vector [ 1 1 2 ] in 'Y'
 10^(11) in 'tol'
and store the following subroutine in 'F'
<< Y OBJ> DROP > Y Z U
<< Y Z U * * NEG
Y Z + U  X *
X Y * Z U *  3 >ARRY
>>
>>
Then to compute y(1) z(1) u(1) place 1 in the stack ( level 1 )
and press the [DIF] key
It gives in 77 seconds
y(1) = 0.258207906452
z(1) = 1.15762398074
u(1) = 0.842178311747
It seems difficult to do faster in RPL
I've written an 'RK8' program for the HP48
that uses an 8thorder RungeKutta method
and it yields almost the same results in 65s
with stepsize = 0.1 but there is no errorestimate.
Les, thank you for your appreciation
but am I really a guru ? just an amateur.
Regards,
JMB.
▼
Posts: 1,368
Threads: 212
Joined: Dec 2006
This is amazing! Quite ingenious.
The simple system given as an example on pages 195 and 196 of the HP48G User's Guide led me to believe that the calculator could handle vectorvalued equations only when the arrays in question consisted of static numeric constants. Your example demonstrates that the builtin ODE integrator can indeed handle a broader range of problems provided that the equations are set up correctly and the derivative vector is computed "on the fly" as your procedure F does.
Thanks for sharing this with us. It helps me put more faith in the HP49+ RungeKutta ODE integrator, and I will experiment with it a bit more before embarking on the daunting task of BulirschStoer for the HP48/49. My RPL skills aren't nearly that good yet.
Les
▼
Posts: 1,368
Threads: 212
Joined: Dec 2006
Quote:
The simple system given as an example on pages 195 and 196 of the HP48G User's Guide
Just a little more info for those who may not know what I am talking about.
The Guide provides as an example the ODE y'' = .5y' + .5y + .5t + 1 with initial condition y(0) = 0 and y'(0) = 0. The example renders this equation as vector equation w' = fw*w + c*(.5t+1), where w = [y y'], fw = [[0 1] [.5 .5]], and c = [0 1].
JMB's approach to this problem would render the desired F subroutine as merely
<< Y OBJ> DROP > Y Z
<< Z
X Y Z + + 2 / 1 +
2 ARRY>
>>
>>
and store 0 in X and [0 0] in variable Y before executing the little routine DIF. This of course solves the system
dy/dx = z y(0) = 0
dz/dx = (x + y + z)/2 + 1 z(0) = 0
which is equivalent to the original second order ODE.
I find JMB's approach much easier to follow.
Les
▼
Posts: 1,368
Threads: 212
Joined: Dec 2006
Quote:
I find JMB's approach much easier to follow.
Just for the record, I should point out that the entry for RKF in the 49G+ AUR does NOT explicitly mention that the solution can be a vector. Nor does the lengthier description in the User's Guide on page 1670 make this clear. One can infer it from the sample problem given on page 1664 which uses a Solve dialog box rather than a command line solution. I suspect there are a lot of functions in the 48 and 49 series calculators whose full power and scope are not fully and clearly disclosed. I never would have known how to use RKF in this way had JMB not pointed it out.
Les
Posts: 1,368
Threads: 212
Joined: Dec 2006
Quote:
It gives in 77 seconds
y(1) = 0.258207906452
z(1) = 1.15762398074
u(1) = 0.842178311747
My 49G+ returns the exact same result with your very same parameters in about 45 seconds. I would gather that the 50G is just as fast.
The solution appears much more quickly for less demanding tolerance levelse.g., the default of .0001, or 1e8.
Les
Posts: 1,368
Threads: 212
Joined: Dec 2006
Quote:
Les, thank you for your appreciation
but am I really a guru ? just an amateur.
Actually, you undervalue your prodigious contributions.
Your higher math programs in the HP41 software library on this site are remarkable for their sheer number and for the fact that they far surpass anything that has been widely available before, including routines in the original User's Library solution books. Your programs are more concise, quicker, and, from what I can tell, about as accurate as one can be on a 10 digit calculator that rounds results every time they are output to the stack.
Frankly, I haven't seen anything quite like this since the 1980s, when our very own Namir Shammas was probably the most prolific contributor of programs to HP41 User's Library catalogue.
Keep on churning out great programs, and hopefully some of the recent work you have shared with me directly will find its way to this site the next time Dave does an update.
Les
Posts: 1,368
Threads: 212
Joined: Dec 2006
Quote:
and store the following subroutine in 'F'
<< Y OBJ> DROP > Y Z U
<< Y Z U * * NEG
Y Z + U  X *
X Y * Z U *  3 >ARRY
>>
>>
I wonder if decomposing the input vector in one step with V> and recomposing the the output in one step with >V3 would speed things up a little?
I have to agree with you now that the builtin RKF function is quite goodprovide that one is aware how to fully take advantage of it. I continue to be intrigued that the AUR entries for both the 49G+ and 48G are so limited. I guess that the authors assume that anyone wanting to use the ODE solver routines are knowledgeable enough about setting up their problems they would have figured it out eventually.
I would be really interested in seeing your RK8 routine as well.
Les
▼
Posts: 191
Threads: 41
Joined: Jun 2007
Hi,
Thank you for all these appreciations!
Here is the 'RK8' program for the HP48
<< OVER 1 SWAP
START DUP 5 PICK >NUM 4 PICK * DUP2 2 / + 6 PICK >NUM
5 PICK * DUP2 + 4 / 4 PICK + 7 PICK >NUM 6 PICK * DUP
.896181193363 * ROT .211711500866 *  3 PICK 7 / + 4 PICK
+ 7 PICK >NUM 6 PICK * OVER .576671472696 * OVER
.065148509147 * + 4 PICK .185506853511 * + 5 PICK + 8 PICK
>NUM 7 PICK * OVER .463455389641 * OVER .386524626691 *
+ 4 PICK .377293769304 * + 5 PICK .199636993645 * + 6 PICK
+ 9 PICK >NUM 8 PICK * OVER .328517213142 * OVER
9.79004561593E2 * + 4 ROLL .349705286318 *  4 ROLL
3.30255113145E2 *  4 PICK .128986292977 * + 5 PICK + 8
PICK >NUM 7 PICK * OVER 1.18686838868E2 * OVER 9 / + 4
PICK 2.00216599311E3 * + 5 PICK 14 / + 6 PICK + 9 PICK
>NUM 8 PICK * OVER .632546160696 * OVER .957605344019 *
+ 4 PICK .152777777778 * + 5 PICK 9.08696110082E3 *  6
PICK 32 / + 7 PICK + 10 PICK >NUM 9 PICK * OVER
1.81086308294 * OVER 1.06249844677 * + 4 PICK 2.03108313917
* + 5 PICK .637931350185 *  6 PICK 9 / + 7 PICK 14 / + 8
PICK + 11 PICK >NUM 10 PICK * OVER 2.22915821019 * OVER
.940109451962 * + 4 PICK 7.55384044212 * + 5 ROLL
7.16495155323 *  5 ROLL 2.45138043242 * + 5 ROLL
.551220563073 *  6 PICK + 9 PICK >NUM 8 PICK * 5 ROLL +
9 * ROT 64 * + 3 ROLLD + 49 * + 180 / +
NEXT
>> ( #62676d , 1032 bytes )
32 bytes may be saved if you replace
10 by 1 ALOG , 49 by 7 SQ ... and so on ...
4 inputs are needed:
level 4: a program ( or its name ) that computes the derivative
of the (n+1)vector [ x y1 ... yn ]
Since dx/dx = 1 the first component is always 1
level 3: the stepsize h
level 2: the number of steps N
level 1: the initial value [ x y1 ... yn ]
For example, the system of 3 differential equations above
with h = 0.1 , N = 10:
<< OBJ> DROP > X Y Z U
<< 1 Y Z U * * NEG
Y Z + U  X *
X Y * Z U *  4 >ARRY
>>
>> ENTER
0.1 ENTER
10 ENTER
[ 0 1 1 2 ] and press the [RK8] key
It yields ( in 61 seconds ) in level 1
[ 1 0.258207906459 1.1576239808 0.842178311686 ]
With h = 0.05 and N = 20 the result is:
[ 1 0.258207906457 1.1576239808 0.842178311706 ]
It's probably faster with an HP49G or 50G
Regards,
JMB.
▼
Posts: 1,368
Threads: 212
Joined: Dec 2006
This is amazing! It will take a few minutes to move it over to my HP49G+ and edit all of the << and > to their proper symbols, but I look forward to experimenting with it.
I anticipate that the routine will give excellent and quick results, especially if one is satisfied with at most 9 or so significant digits. Since the many rational coefficients that are required are rounded to twelve decimal digits, I anticipate that accumulated rounding error would readily supersede the O(h^9) global truncation error at some point. I have encountered the same sort of problem in porting GaussLegendre integration to the 49G+the calculator is very good at computing weights and sample points, but the accumulated rounding errors tend to obscure the benefits of carrying out quadratures to a high degree.
I anticipate the actual formulae are much more complicated than the classical RK4 method or even the somewhat more complex RKF45 RungeKuttaFehlberg method. What is your source? Abramowitz and Stegun doesn't go up that high!
Many thanks,
Les
▼
Posts: 1,368
Threads: 212
Joined: Dec 2006
Quote:
I anticipate that accumulated rounding error would readily supersede the O(h^9) global truncation error
Actually, this is incorrectan eighth order RK method would have O(h^8) global truncation error, by definition, and O(h^9) at each step.
Les
Posts: 1,368
Threads: 212
Joined: Dec 2006
Quote:
It's probably faster with an HP49G or 50G
Actually it is a lot faster.
On my 49G+, the 10 step version of your example returns the exact same result as you get in 26 seconds. The 20 step version does likewise inyou guessed it52 seconds.
Excellent bit of work. I am fascinated by RPL code that has such an obvious awareness of what is happening on the stack at all times. I am not so confident and use a lot of local variables instead, and I believe my programs are not as small or fast as they could be.
Please keep writing great stuff like this and sharing it. And this is a small nudge to Davea update of the software library is long overdue! ;) Yes, I know you have to squeeze it around 80 hours of week of scanning!
Les
▼
Posts: 191
Threads: 41
Joined: Jun 2007
Hi Les,
this 'RK8' program for the HP48 uses
the same formula as the 'RK8' program for the HP41
( cf "RungeKutta Methods for the HP41" )
The source is:
J.C. Butcher  "The Numerical Analysis of Ordinary Differential
Equations"  John Wiley & Sons  ISBN 0471910465
Roundoff errors don't seem very important,
even in the last digits.
Regards,
JMB.
▼
Posts: 1,368
Threads: 212
Joined: Dec 2006
Quote:
Hi Les,
this 'RK8' program for the HP48 uses
the same formula as the 'RK8' program for the HP41
( cf "RungeKutta Methods for the HP41" )
The source is:
J.C. Butcher  "The Numerical Analysis of Ordinary Differential
Equations"  John Wiley & Sons  ISBN 0471910465
Roundoff errors don't seem very important,
even in the last digits.
Regards,
JMB.
I found the HP41 routine and the referencethank you.
I am happily corrected about the rounding errors. I find that integrating the IVP y' = y^2 + 1, y(0) = 0, gives 11 correct digits with 12 steps, and the 12th digit is off by just 1 ULP. Increasing step number doesn't improve the situation. I think that's pretty impressive!
I wonder if there is any way of improving results with Richardson extrapolation? My preliminary impression is nodoing four steps, then eight steps, then extrapolating on the assumption that the procedure is of order 8 doesn't improve things much and doesn't give nearly as good results as just bashing ahead with 12 steps in the first place!
Les
Posts: 1,368
Threads: 212
Joined: Dec 2006
Quote:
It seems difficult to do faster in RPL
I've written an 'RK8' program for the HP48
that uses an 8thorder RungeKutta method
and it yields almost the same results in 65s
with stepsize = 0.1 but there is no errorestimate.
Actually, I have learned that if I can spare a couple of extra minutes on a problem your RK8 method applied twice can generate rough error estimates.
As an 8th order method, stepsize halving should (in theory) roughly improve the error by a factor of 2^8 = 1/256. This means if I run the problem twice with stepsize 2h and then stepzize h, the difference in results can give a rough measure of the estimated error for each run. I won't regurgitate the algebra, which is actually not too difficult, but it seems to me that the error of the 2h run is of the order of 256/255 times this difference, and the error of the h run is roughly the difference divided by 255.
I have experimented with this with some simple problems and find that the error estimates are actually pretty good, especially in those admittedly ideal situations where the analytical solution of the DE is nice and smooth and continuous with an infinity of well behaved higher derivatives that make the Taylor series upon which the error estimates are based a reasonable approximation of reality. The bottom line is that the difference between two runs can yield some meaningful information, and, if one is willing to take the time, one can get a sense of the result the solution is approaching.
I must confess that I do in some cases like routines such as this one where the number of iterations are set a priori by me, since I know it is going to end some time! Adaptive schemes that converge according to some desired tolerance may actually never converge, or take for ever. I would prefer to know as soon as possible whether a particular problem is stiff or just illposed.
I think your RK8 is in some ways a desirable alternative to the HP48/49 series builtin RKF45 adaptive solver. I am also impressed at how well it ports to these calculators, and how quickly I was able to get it onto my calculator via SD card. The HP41 version is excellent too, but is of course slower and really requires a card reader or extended memory so that one need only enter and store all of those constants once. I also think it would well on the HP42Sthe constants can be stored and saved in a vector and swapped for the REGS variable as required. But I would have to hope like heck not to reset the memory at my next battery swap!
Les
Edited: 3 May 2007, 9:27 p.m.
▼
Posts: 191
Threads: 41
Joined: Jun 2007
Hi,
step doubling is always a good idea to provide an errorestimate.
It can also be used to get a 9thorder method by the formula
( 256 y(h/2)  y(h) )/255 if it is applied at each step.
I've searched higher order formulas on the web,
but they seem to be kept secret...
or sold at a high price!
Regards,
JMB.
▼
Posts: 1,368
Threads: 212
Joined: Dec 2006
Quote:
I've searched higher order formulas on the web,
but they seem to be kept secret...
or sold at a high price!
For practical use, I think that higher order routines are really more of a theoretical curiosity than anything.
I think your RK8 routine is about as complex as anyone using a handheld calculator could reasonably desire. Indeed, the 48/49 series implementation of RKF45 is also pretty powerful too. I really wish that the HP literature documented it betterI had to do a bit of outside reading to learn more about the algorithm and develop an appreciation about how clever it was. It really takes a lot of ingenuity to develop formulae that allow for both 4th and 5th order estimations while using the same set of computed coefficients.
Les
