Runge-Kutta-Verner method with error-estimate



#4

Hi!

I've found on the web the coefficients of a Runge-Kutta-Verner
method of order 8 embedded within 9th-order.
Here is a program that uses these formulae for the HP-48:

------------------------- 'RKV8E9' ----------------------------

<< OVER 1 SWAP
START DUP 5 PICK ->NUM 4 PICK * DUP 12 / 3 PICK + 6 PICK ->NUM
5 PICK * 2 * OVER + 27 / 3 PICK + 6 PICK ->NUM 5 PICK * DUP2 3 *
+ 24 / 4 PICK + 7 PICK ->NUM 6 PICK * DUP 2.23331697733 * ROT
-2.39805710715 * + 3 PICK .624672095524 * + 4 PICK + 7 PICK ->NUM
6 PICK * DUP .245675793931 * 3 PICK .273953453873 * + 4 PICK
4.36700683814E-2 * + 5 PICK + 8 PICK ->NUM 7 PICK * DUP
7.02490360993E-3 * ROT 1.34776896111E-2 * - ROT .181531822412 * +
3 PICK 6.16216474034E-2 * + 4 PICK + 7 PICK ->NUM 6 PICK * DUP
.341657217459 * 3 PICK .250935375134 * + 4 PICK 7.40740740741E-2 *
+ 5 PICK + 8 PICK ->NUM 7 PICK * DUP -.03515625 * 3 PICK
.340504422039 * + 4 PICK .120433077961 * + 5 PICK .07421875 * +
6 PICK + 9 PICK ->NUM 8 PICK * DUP -.296296296296 * 3 PICK 16 / -
4 PICK .310705427943 * + 5 PICK .305035312798 * + 6 PICK
7.63888888889E-2 * + 7 PICK + 10 PICK ->NUM 9 PICK * DUP
-.26063186736 * 3 PICK 7.27205419732E-2 * + 4 PICK .011746330035
* - 5 PICK .37852828889 * + 7 PICK 7.11293665317E-2 * + 8 PICK +
11 PICK ->NUM 10 PICK * DUP -328.691358025 * 3 PICK 241.543474144
* - 4 PICK 626.941484898 * + 5 PICK 113.719201869 * + 6 PICK
413.485511011 * + 7 PICK 574.436392562 * - 8 PICK 8.14163971388
* - 9 PICK + 12 PICK ->NUM 11 PICK * DUP 1.80288461538E-3 * 3 PICK
2.49192782526 * + 4 PICK 7.69118880716E-2 * - 5 PICK .690428248367
* - 6 PICK .228863388685 * + 7 PICK 1.9030978898 * - 8 PICK
.69337350173 * + 9 PICK .087803759282 * + 10 PICK + 13 PICK ->NUM
12 PICK * DUP 1.77777777778 * 3 PICK 2.35042735043E-2 * - 4 PICK
4.86229819563 * - 5 PICK 4.88888888889 * - 6 PICK 2.27160493827 *
+ 7 PICK 6.22916666667 * - 8 PICK 7.48550699458 * + 9 PICK
5.5746781906 * + 10 PICK .105709876543 * - 11 PICK + 14 PICK ->NUM
13 PICK * OVER 26.25 / 4 PICK 1920 / - 5 PICK .564373897707 * +
6 PICK .867950955719 * - 7 PICK .701048612636 * + 8 PICK
.125460177737 * - 9 PICK 7.51961665302E-2 * + 10 PICK
.272718881509 * - 11 PICK 5.46035999955E-2 * + 12 PICK + 15 PICK
->NUM 14 PICK * DUP 5.40772532189 * 4 PICK .422317596567 * + 5
PICK 2.87223506108E-3 * - 6 PICK 18.706283702 * - 7 PICK
7.80529021378 * + 8 PICK 3.73079546162 * + 9 PICK 1.38102723196
* + 10 ROLL 6.81524922033 * + 10 ROLL 5.45684741856 * - 10 PICK
.396401690533 * - 11 PICK + 14 PICK ->NUM 13 PICK *
-5.54761904762E-2 * SWAP .36 * - OVER 3.21428571429E-2 * + 3 PICK
.12 * + 4 PICK 5.76923076923E-4 * + 5 PICK 1.05025641026 * + 6
PICK 1.05 * - 7 PICK .56 * + 8 PICK .315 * - 9 PICK .0175 * +
1 0 PUT 'ERROR'
IFERR STO+
THEN STO
END
3.21428571429E-2 * SWAP .342857142857 * + SWAP 4.12087912088E-4 *
+ SWAP .750183150183 * + SWAP .717857142857 * - SWAP .72380952381
* + SWAP .192857142857 * - SWAP 6.13095238095E-2 * + +
NEXT
>> ( 2368.5 bytes / #57545d )


*** PURGE 'ERROR' before the first execution.

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, 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

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 [RK8E9] key

-It yields ( in 137 seconds ) in level 1

[ 1 0.258207906449 1.15762398084 0.842178311722 ]

and 'ERROR' contains [ 0 2.7631E-11 3.27392E-11 7.636E-11 ]

Regards,
JMB.


#5

JM, this is an excellent bit of work as usual.

Your sample problem is somewhat faster on my 49G+, and I would expect similar performance on the 50G.

I find in some examples I have worked with that the error estimate is a little optimistic. This is understandable--users will know that the error estimates of the Runge-Kutta methods are based on the higher derivative terms of Taylor series expansions, and the assumption that these derivative values remain roughly the same across each subinterval is often not a reasonable one. Moreover, the estimated error is typically used in more sophisticated routines to adjust stepsize--do better than a desired tolerance at the step, increase the stepsize, do worse, decrease it.

Your routine raises some interesting questions about the pros and cons using higher order methods vs. repeated application of lower order methods. In which case is the computational work load best balanced by the performance? I believe that the 48 and 49 series calcs use RKF45 as the basis of an adaptive stepsize routine.

One day, perhaps after all of the hubbub regarding the 35s dies down , I will post some thoughts on this, and perhaps issue a challenge on how to best find results to some challenging nonstiff IVPs.

Les


#6

Hi,

the error-estimates are algebraically added after each step.
Adding the magnitudes could be a less optimistic alternative:

Simply add OBJ-> EVAL ->LIST ABS OBJ-> ->ARRY

just before 1 0 PUT 'ERROR' near the end

I was worried when I saw that some coefficients
of this Runge-Kutta-Verner method were of the order of 600:
I thought it was going to produce great roundoff errors,
but fortunately, the corresponding k-value ( namely k12 )
is weighted with a small 3/7280 and indeed,
the results are usually satisfactory, even in the last decimal:

For instance: y'= -2xy , y(0) = 1

h = 0.1 gives y(1) = 0.367879441185
h = 0.05 ----- y(1) = 0.367879441171 ( exact )
h = 0.025 ---- y(1) = 0.367879441171 ( still exact! )

At least, 'RKV8E9' gives an idea of the obtained accuracy
with 16 evaluations of the function per step,
instead of 11+2*11=33 evaluations if one uses 'RK8' with h and then with h/2.
This may be an advantage if the functions are complicated.

In "Numerical Recipes" they don't seem to like high-order Runge-Kutta formulas.
On the other hand, they praise Bulirsch-Stoer methods
( "Runge-Kutta is for ploughing the fields,
Burlisch-Stoer is a high-strung racehorse" )
and it is a little contradictory:
though quite fantastic, these methods require more than 11
evaluations per step to achieve an 8th-order formula!
Moreover, roundoff errors are also greater.

I've always been fascinated by Runge-Kutta methods,
especially if they are compared to the Taylor's method
and the huge amount of calculus they require!
And the size of the non-linear systems that must be solved to find high-order formulae!
It seems like a touch of magic!

Of course, the classical 'RK4' is probably enough for an HP-41,
but we can try more accurate formulae with a 49G or 5OG.

In fact, I wrote 'RK8' and 'RKV8E9' for completeness, perhaps it will be useful.
However, I must say that, even now, I still prefer the HP-41 programming language
with its direct arithmetic in the 4-level stack,
and all the trickeries we can ( must? ) find to create neat programs.

Speed may be a problem, but thanks to Warren Furlows and his excellent V41,
HP-41 programs can run faster than HP-49 programs!

But I stop here my "philosophical' point of view...

Regards,
JMB.


Possibly Related Threads…
Thread Author Replies Views Last Post
  A fast Bernoulli Number method for the HP Prime Namir 16 5,651 11-22-2013, 04:46 PM
Last Post: Namir
  MS advert shows spreadsheet with obvious error BruceH 3 2,248 11-14-2013, 09:50 AM
Last Post: Bill (Smithville, NJ)
  HP Prime: Rounding error in determinant Stephan Matthys 3 1,699 10-25-2013, 09:29 PM
Last Post: Walter B
  Prime Error or Mine? toml_12953 12 3,798 10-22-2013, 10:35 AM
Last Post: toml_12953
  Explaination on How to Reset Caculator in Users guie: error Harold A Climer 5 2,643 10-15-2013, 02:11 AM
Last Post: cyrille de Brébisson
  Repair of HP-34C - Error 0 and Error 9 Jeff Kearns 3 1,753 10-11-2013, 12:29 PM
Last Post: Randy
  A brand new calculator benchmark: "middle square method seed test" Pier Aiello 25 7,195 09-13-2013, 01:58 PM
Last Post: Pier Aiello
  Do You Think a Listing Error Was Made? Jim Johnson 13 3,980 09-04-2013, 09:23 AM
Last Post: Eddie W. Shore
  HP Prime: Error reports : Here Patrice 111 24,001 07-24-2013, 05:52 PM
Last Post: Thomas Klemm
  Downhill Simplex Method (Nelder Mead) for WP-34S Marcel Samek 0 920 07-07-2013, 03:05 AM
Last Post: Marcel Samek

Forum Jump: