Differential Equations on the 49G+ and 50G



#20

That guru of programs for the the HP41, J-M Baillard, recently sent to me some routines to solve first and second-order 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 built-in differential equation solver of my 48G and 49G+. Baillard's programs work really well on the 41-series 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 Bullrisch-Stoer ODE solvers for the 49G series. I don't mind a little challenge, but I don't want to completely reinvent the wheel.

Les


#21

Make sure you have checked the items at the link below.

http://www.hpcalc.org/search.php?query=differential


#22

Thank you! I have done a more specific search of hpcalc.org and haven't come up with anything that specifically refers to Bullrisch-Stoer. 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


#23

Hi everyone,

I have not ( yet? ) written a Bulirsch-Stoer program for the HP48
( it's in the pipeline for the HP-41 to solve systems of ODEs )
but I think that the built-in 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+z-u) z(0)=1
du/dx = xy-zu 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 HP-48
that uses an 8th-order Runge-Kutta method
and it yields almost the same results in 65s
with stepsize = 0.1 but there is no error-estimate.

-Les, thank you for your appreciation
but am I really a guru ? just an amateur.

Regards,
JMB.


#24

This is amazing! Quite ingenious.

The simple system given as an example on pages 19-5 and 19-6 of the HP48G User's Guide led me to believe that the calculator could handle vector-valued equations only when the arrays in question consisted of static numeric constants. Your example demonstrates that the built-in 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+ Runge-Kutta ODE integrator, and I will experiment with it a bit more before embarking on the daunting task of Bulirsch-Stoer for the HP48/49. My RPL skills aren't nearly that good yet.

Les


#25

Quote:
The simple system given as an example on pages 19-5 and 19-6 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


#26

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 16-70 make this clear. One can infer it from the sample problem given on page 16-64 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

#27

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 levels--e.g., the default of .0001, or 1e-8.

Les

#28

Quote:

-Les, thank you for your appreciation
but am I really a guru ? just an amateur.


Actually, you under-value 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

#29

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 built-in RKF function is quite good--provide 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


#30

Hi,
Thank you for all these appreciations!
Here is the 'RK8' program for the HP-48

<< 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.79004561593E-2 * + 4 ROLL .349705286318 * - 4 ROLL
3.30255113145E-2 * - 4 PICK .128986292977 * + 5 PICK + 8
PICK ->NUM 7 PICK * OVER -1.18686838868E-2 * OVER 9 / + 4
PICK 2.00216599311E-3 * + 5 PICK 14 / + 6 PICK + 9 PICK
->NUM 8 PICK * OVER -.632546160696 * OVER .957605344019 *
+ 4 PICK .152777777778 * + 5 PICK 9.08696110082E-3 * - 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 HP-49G or 50G

Regards,
JMB.


#31

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 Gauss-Legendre 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 Runge-Kutta-Fehlberg method. What is your source? Abramowitz and Stegun doesn't go up that high!

Many thanks,

Les


#32

Quote:
I anticipate that accumulated rounding error would readily supersede the O(h^9) global truncation error

Actually, this is incorrect--an eighth order RK method would have O(h^8) global truncation error, by definition, and O(h^9) at each step.

Les

#33

Quote:
-It's probably faster with an HP-49G 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 in--you guessed it--52 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 Dave--a update of the software library is long overdue! ;) Yes, I know you have to squeeze it around 80 hours of week of scanning!

Les


#34

Hi Les,

this 'RK8' program for the HP-48 uses
the same formula as the 'RK8' program for the HP-41
( cf "Runge-Kutta Methods for the HP-41" )
The source is:
J.C. Butcher - "The Numerical Analysis of Ordinary Differential
Equations" - John Wiley & Sons - ISBN 0-471-91046-5

Roundoff errors don't seem very important,
even in the last digits.

Regards,
JMB.


#35

Quote:

Hi Les,

this 'RK8' program for the HP-48 uses
the same formula as the 'RK8' program for the HP-41
( cf "Runge-Kutta Methods for the HP-41" )
The source is:
J.C. Butcher - "The Numerical Analysis of Ordinary Differential
Equations" - John Wiley & Sons - ISBN 0-471-91046-5

Roundoff errors don't seem very important,
even in the last digits.

Regards,
JMB.



I found the HP41 routine and the reference--thank 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 no--doing 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

#36

Quote:

-It seems difficult to do faster in RPL
-I've written an 'RK8' program for the HP-48
that uses an 8th-order Runge-Kutta method
and it yields almost the same results in 65s
with stepsize = 0.1 but there is no error-estimate.


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 ill-posed.

I think your RK8 is in some ways a desirable alternative to the HP48/49 series built-in 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 HP42S--the 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.


#37

Hi,

step doubling is always a good idea to provide an error-estimate.
It can also be used to get a 9th-order 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.


#38

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 better--I 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


Possibly Related Threads...
Thread Author Replies Views Last Post
  [HP Prime] Tips for Solving Differential Equations More Efficiently Chris Pem10 8 1,523 11-21-2013, 08:25 PM
Last Post: Chris Pem10
  HP Prime Solving Nonlinear System of Equations for Complex Results Helge Gabert 11 2,201 09-30-2013, 03:44 AM
Last Post: From Hong Kong
  New article on a new type of neo linear equations Namir 0 651 08-11-2013, 10:27 AM
Last Post: Namir
  Riemann's Zeta Function update (HP-28S, HP-48G/GX/G+, HP-49G/G+/50g) Gerson W. Barbosa 0 551 06-30-2013, 01:01 AM
Last Post: Gerson W. Barbosa
  hp 49g+ and 50g cases William L. Drylie 16 2,419 06-04-2013, 07:28 PM
Last Post: Iqbal
  OT: DIFFERENTIAL EQUATIONS using the HP-48G/GX Ivan Rancati 2 757 06-04-2013, 04:11 AM
Last Post: Frank Boehm (Germany)
  HP-50g How to store many equations - help needed Timo Labrenz 5 1,165 03-03-2013, 12:22 PM
Last Post: Timo Labrenz
  3D Model - HP 50g 49G+ 39G+ 48Gii JJB299 3 803 01-18-2013, 11:51 PM
Last Post: Raymond Del Tondo
  Solving with the built-in equations of the 35s Palmer O. Hanson, Jr. 4 832 10-17-2012, 03:17 PM
Last Post: Dieter
  Equations library on WP34S Alessandro Castellani (Italy) 2 651 08-27-2012, 12:29 PM
Last Post: Alessandro Castellani (Italy)

Forum Jump: