Runge-Kutta-Verner method with error-estimate

[Jean-Michel]

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.

[Les Wright]

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

[Jean-Michel]

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.