hp41programs

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
                                      •  R02 = c3  |  Rss+1 = a3,1   Rss+2 = a3,2
                                      •  R03 = c4  |  Rss+3 = a4,1   Rss+4 = a4,2  Rss+5 = a4,3
                                         ..................................................................................................

                                     •  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] !