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 ( n² + 3.n - 2 ) / 2 coefficients a_i,j , c_i , b_i

    k₁ = h.f(x;y)
    k₂ = h.f(x+c₂.h;y+a_2,1.k₁)
    k₃ = h.f(x+c₃.h;y+a_3,1.k₁+a_3,2.k₂)
    .............................................................                                         (  actually,     a_i,1 + a_i,2 + .......... + a_i,i-1 = c_i  )

    k_s = h.f(x+c_s.h;y+a_s,1.k₁+a_s,2.k₂+ .......... +a_s,s-1.k_s-1)

then,     y(x+h) = y(x) + b₁.k₁+ ................ + b_s.k_s

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

    k₁ = h.f(x , y , y')
    k₂ = h.f(x+c₂.h , y+c₂.h.y'+h A_2,1.k₁ , y'+a_2,1.k₁)
    k₃ = h.f(x+c₃.h , y+c₃.h.y'+h(A_3,1.k₁+A_3,2.k₂) , y'+a_3,1.k₁+a_3,2.k₂)
    .................................................................................

    k_s = h.f(x+c_s.h , y+c_s.h.y'+h (A_s,1.k₁+A_s,2.k₂+ .......... +A_s,s-1.k_s-1) , y'+a_s,1.k₁+a_s,2.k₂+ .......... +a_s,s-1.k_s-1 )

then,     y(x+h) = y(x) + h y'(x) + h [ B₁.k₁+ ................ + B_s.k_s ]   and   y'(x+h) = y'(x) + [ b₁.k₁+ ................ + b_s.k_s ]

-The program listed hereunder takes the coefficients   a_i,j , c_i , b_i  which will be used to calculate y'(x+h)
  and replace them by the coefficients A_i,j , c_i , B_i    to calculate  y(x+h)  ( c_i are unchanged )

-The method is described in reference [1]:
 

  A_j,k = ( c_j - c_k ) a_j,k   k = 2 , ..... , j-1

   B_j  = ( 1 - c_j ) b_j       j = 1 , ..... , s

  A_j,1 = (1/2) c_j² - Sum_k=2,...,j-1  A_j,k    j = 2 , ...... , s       ;      A_2,1 = (1/2) c₂²
 

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 = c₂  |    Rss   = a_2,1
                                      •  R02 = c₃  |  Rss+1 = a_3,1   Rss+2 = a_3,2
                                      •  R03 = c₄  |  Rss+3 = a_4,1   Rss+4 = a_4,2  Rss+5 = a_4,3
                                         ..................................................................................................

                                     •  Rss-1 = c_s | R(s²-s+2)/2 = a_s,1 ........................................................   R(s²+s-2)/2 = a_s,s-1
                                        ____________________________________________________________________________________

                                                          |  R(s²+s)/2 = b₁ .............................................................................................   R(s²+3s-2)/2 = b_s
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"  >>>>  A_2,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] ):

   k₁ = h.f (x,y,y')
   k₂ = h.f( x+h/2 , y+h.y'/2+h.k₁/8 , y'+k₁/2 )
   k₃ = h.f( x+h/2 , y+h.y'/2+h.k₁/8 , y'+k₂/2 )
   k₄ = h.f( x+h , y+h.y'+h.k₃/2 , y'+k₃ )

 y(x+h) = y(x) + h [ y'(x) + ( k₁ + k₂ + k₃ )/6 ]
 y'(x+h) = y'(x) + ( k₁ + 2.k₂ + 2.k₃ + k₄ ) / 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"   >>>>    A_2,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] !