RK-RKNG

# Runge-Kutta -> RKNG Coeff. ( 2nd-order ODEs ) for the HP-41

Overview

-Given a Runge-Kutta s-stage formula to solve differential equations y' = f(x,y) ,
the following routine calculates the coefficients of an s-stage Runge-Kutta-Nystrom formula of the same order to solve  y" = f(x,y,y')

-An s-stage explicit Runge-Kutta formula is defined by ( n2 + 3.n - 2 ) / 2 coefficients ai,j , ci , bi

k1 = h.f(x;y)
k2 = h.f(x+c2.h;y+a2,1.k1)
k3 = h.f(x+c3.h;y+a3,1.k1+a3,2.k2)
.............................................................                                         (  actually,     ai,1 + ai,2 + .......... + ai,i-1 = ci  )

ks = h.f(x+cs.h;y+as,1.k1+as,2.k2+ .......... +as,s-1.ks-1)

then,     y(x+h) = y(x) + b1.k1+ ................ + bs.ks

-We can use this formula to solve y" = f(x,y,y')  after replacing this ODE by a system of two 1st-order differential equations.

>>> But it is better to employ a special formula which will use less data registers:

k1 = h.f(x , y , y')
k2 = h.f(x+c2.h , y+c2.h.y'+h A2,1.k1 , y'+a2,1.k1)
k3 = h.f(x+c3.h , y+c3.h.y'+h(A3,1.k1+A3,2.k2) , y'+a3,1.k1+a3,2.k2)
.................................................................................

ks = h.f(x+cs.h , y+cs.h.y'+h (As,1.k1+As,2.k2+ .......... +As,s-1.ks-1) , y'+as,1.k1+as,2.k2+ .......... +as,s-1.ks-1 )

then,     y(x+h) = y(x) + h y'(x) + h [ B1.k1+ ................ + Bs.ks ]   and   y'(x+h) = y'(x) + [ b1.k1+ ................ + bs.ks ]

-The program listed hereunder takes the coefficients   ai,j , ci , bi  which will be used to calculate y'(x+h)
and replace them by the coefficients Ai,j , ci , Bi    to calculate  y(x+h)  ( ci are unchanged )

-The method is described in reference [1]:

Aj,k = ( cj - ck ) aj,k   k = 2 , ..... , j-1

Bj  = ( 1 - cj ) bj       j = 1 , ..... , s

Aj,1 = (1/2) cj2 - Sumk=2,...,j-1  Aj,k    j = 2 , ...... , s       ;      A2,1 = (1/2) c22

Program Listing

Data Registers:           •  R00 = s = number of stages            ( Registers R00 thru R(s^2+3.s-2)/2 are to be initialized before executing "RKNG" )

•  R01 = c2  |    Rss   = a2,1
..................................................................................................

•  Rss-1 = cs | R(s2-s+2)/2 = as,1 ........................................................   R(s2+s-2)/2 = as,s-1
____________________________________________________________________________________

|  R(s2+s)/2 = b1 .............................................................................................   R(s2+3s-2)/2 = bs
Flags: /
Subroutines: /

 01 LBL "RKNG"  02 3  03 RCL 00  04 STO O  05 ST+ Y  06 *  07 2  08 /  09 1  10 ST- O 11 -  12 STO N  13 LBL 00  14 1  15 RCL IND O   16 -  17 ST* IND N  18 DSE N  19 DSE O  20 GTO 00 21 RCL 00  22 1  23 ST- N  24 -  25 STO O  26 LBL 01  27 RCL O  28 STO M  29 RCL IND O   30 X^2 31 2  32 /  33 LBL 02  34 DSE M  35 FS? 30  36 GTO 03  37 RCL IND O  38 RCL IND M  39 -  40 RCL IND N   41 * 42 STO IND N  43 -  44 DSE N  45 GTO 02  46 LBL 03  47 STO IND N   48 DSE N  49 DSE O  50 GTO 01  51 CLA  52 END

( 88 bytes / SIZE s(s+3)/2 )

 STACK INPUT OUTPUT X / A2,1

Example1:  The classical 4th-order Runge-Kutta method is defined by Butcher tableau:

1/2 |   1/2                                                                                              R01  |  R04
1/2 |    0     1/2                           Store these numbers into registers        R02  |  R05  R06
1   |    0      0     1                                                                                 R03  |  R07  R08  R09
--------------------------                                                                      -----------------------------
1/6   1/3   1/3   1/6                                                                               |  R10  R11  R12  R13

-Since it's a 4-stage formula  4  STO 00  XEQ "RKNG"  >>>>  A2,1 = 0.125                                                 ---Execution time = 5s---
-And registers R01 to R13 now contain:

1/2 |   1/8                                                                                              R01  |  R04
1/2 |   1/8    0                                      are in registers                            R02  |  R05  R06
1   |    0      0    1/2                                                                               R03  |  R07  R08  R09
--------------------------                                                                      -----------------------------
1/6   1/6   1/6     0                                                                                |  R10  R11  R12  R13

-So, we've found the formula ( given in reference [4] ):

k1 = h.f (x,y,y')
k2 = h.f( x+h/2 , y+h.y'/2+h.k1/8 , y'+k1/2 )
k3 = h.f( x+h/2 , y+h.y'/2+h.k1/8 , y'+k2/2 )
k4 = h.f( x+h , y+h.y'+h.k3/2 , y'+k3 )

y(x+h) = y(x) + h [ y'(x) + ( k1 + k2 + k3 )/6 ]
y'(x+h) = y'(x) + ( k1 + 2.k2 + 2.k3 + k4 ) / 6

-The coefficients for y' are those of the original formula
-The coefficients for y are those given by "RKNG"

Example2:    A 7-stage 6th-order Runge-Kutta formula was found by J.C. Butcher in 1964 ( cf reference [2] or [3] )

1/3 |    1/3                                                                                                                             R01 |  R07
2/3 |     0           2/3                                                                                                               R02 |  R08  R09
1/3 |   1/12        1/3     -1/12                                               Store these ccoefficients into       R03 |  R10  R11  R12
5/6 |  25/48    -55/24   35/48       15/8                                                                                  R04 |  R13  R14  R15  R16
1/6 |   3/20     -11/24    -1/8         1/2       1/10                                                                      R05 |  R17  R18  R19  R20  R21
1   |-261/260  33/13  43/156  -118/39  32/195   80/39                                                        R06 |  R22  R23  R24  R25  R26  R27
-----------------------------------------------------------------                                           ---------------------------------------------
|  13/200       0      11/40      11/40      4/25     4/25  13/200                                                   |  R28  R29  R30  R31  R32  R33  R34

-Seven stages, so   7  STO 00   XEQ "RKNG"   >>>>    A2,1 = 0.055555556                                     ---Execution time = 13s---
-And we get in registers R01 to R34::

1/3 |    1/18
2/3 |      0             2/9
1/3 |    1/36           0          1/36
5/6 |  125/288  -55/48    35/288     15/16
1/6 |    1/40      11/144     1/16       -1/12       -1/15
1   |-261/260    22/13    43/468   -236/117  16/585   200/117
--------------------------------------------------------------------
|  13/200        0         11/120     22/120      2/75       2/15       0

Notes:

-Of course, there are sometimes a few roundoff-errors.
-The maximum s-value is 23 stages ( SIZE 299 ). s = 24 would require SIZE 324.

References:

[1]  P.W. Sharp & J.M. Fine - A contrast of direct and transformed Nystrom pairs
[2]  J.C. Butcher  - "The Numerical Analysis of Ordinary Differential Equations" - John Wiley & Sons  -   ISBN  0-471-91046-5
[3]  Metin Hatun , Fahri Vatansever - Differential Equation Solver Simulator for Runge-Kutta Methods
[4]  Abramowitz and Stegun - "Handbook of Mathematical Functions" -  Dover Publications -  ISBN  0-486-61272-4

-Many thanks to professor Philip Sharp for sending me the document [1] !