# hp41programs

Runge-Kutta-Nystrom

# Runge-Kutta-Nystrom Methods for the HP-41

Overview

1°)  A 4th-Order Method

a)  Second-order differential equations
b)  Systems of 2 second-order differential equations
c)  Systems of 3 second-order differential equations

2°)  A Runge-Kutta-Nystrom Method with order 6

a)  Second order differential equations
b)  Systems of 2 second-order differential equations

3°)  S-Stage Explicit Runge-Kutta-Nystrom formulae

a)  Second order differential equations
b)  Systems of 2 second-order differential equations

-A differential equation of the type y" = f(x,y) may always be replaced by a system of 2 first order equations:

y' = z
z' = f(x,y)

-However, the same order of accuracy may be obtained with less evaluations of the function if we use formulae specially devised for these problems.

3  evaluations instead of  4  for a  4th-order method
5  evaluations instead of  7  for a  6th-order method
13 evaluations instead of 17 for a 10th-order method ... etc ...

Notes:

-The programs of paragraph 1 are already listed elsewhere in this website
( cf "Explicit Runge-Kutta Methods" & "Systems of Differential Equations" )
-They are just included here for the sake of consistency.

1°)  A 4th-Order Method

a) Second-order differential Equations

-Here, we have to solve:     y" = f(x,y)     with the initial values:     y(x0) = y0     y'(x0) = y'0

-This program uses the 4th-order Runge-Kutta formula:

y(x+h) = y(x) + h ( y'(x) + k1/6 + 2k2/3 )
y'(x+h) = y'(x) + k1/6 + 2k2/3 + k3/6                 where  k1 = h.f (x,y)  ;  k2 = h.f(x+h/2,y+h.y'/2+h.k1/8)  ;  k3 = h.f(x+h,y+h.y'+h.k2/2)

-Note that only 3 evaluations of the function are needed ( for each step ).

Data Registers:    •   R00:  subroutine name                    ( Registers R00 thru R05  are to be initialized before executing "RKS1" )

•   R01 = x0      •   R04 = h/2 = half of the stepsize        R06 to R08: temp
•   R02 = y0      •   R05 = N = number of steps
•   R03 = y'0
Flags: /
Subroutine:  a program which calculates f(x,y) in X-register, assuming  x and y  are in registers X and Y upon entry.

 01  LBL "RKS1"  02  RCL 05  03  STO 08  04  LBL 01  05  RCL 02  06  RCL 01  07  XEQ IND 00  08  RCL 04  09  ST+ 01  10  *  11  STO 06  12  2 13  /  14  RCL 03  15  +  16  RCL 04  17  *  18  RCL 02  19  +  20  RCL 01  21  XEQ IND 00  22  RCL 04  23  ST+ 01  24  * 25  STO 07  26  RCL 03  27  +  28  RCL 04  29  ST+ X  30  *  31  RCL 02  32  +  33  RCL 01  34  XEQ IND 00  35  RCL 04          36  * 37  RCL 06  38  +  39  RCL 07          40  ST+ X  41  ST+ 06  42  ST+ X  43  +  44  3  45  ST/ 06  46  /  47  X<> 03  48  ST+ 03 49  RCL 06   50  +  51  RCL 04          52  ST+ X  53  *  54  ST+ 02  55  DSE 08  56  GTO 01  57  RCL 03  58  RCL 02  59  RCL 01  60  END

( 85 bytes / SIZE 009 )

 STACK INPUTS OUTPUTS Z / y'(x0+N.h) Y / y(x0+N.h) X / x0+N.h

Example:     y(0) = 1  ,   y'(0) = 1  ,   y" = - y ( x2 + y2 ) 1/2

>>> Calculate  y(1) & y'(1)

-Load this short routine:

 01  LBL "T"  02  X^2  03  X<>Y  04  STO Z  05  X^2  06  +  07  SQRT  08  *  09  CHS  10  RTN

"T"  ASTO 00

0  STO 01  STO 03  1  STO 02

-With h = 0.1  and  N = 10  ,  0.05  STO 04  10  STO 05

XEQ "RKS1"  >>>>      x   = 1                                            ---Execution time = 29s---
RDN   y(1) =  0.536630911
RDN   y'(1) = -0.860172085

-With  h = 0.02  &  N = 50 ,  it yields

y(1) =  0.536630617
y'(1) = -0.860171928

Note:

-Simply press R/S to continue with the same h &N

b) Systems of 2 Second-order differential Equations

-"RKS2" solves the system:

y" = f(x,y,z)         y(x0) = y0     y'(x0) = y'0
z" = g(x,y,z)         z(x0) = z0     z'(x0) = z'0

Data Registers:    •   R00:  subroutine name                     ( Registers R00 thru R07  are to be initialized before executing "RKS2" )

•   R01 = x0      •   R04 = h/2 = half of the stepsize       •   R06 = y'0          R08 to R12: temp
•   R02 = y0      •   R05 = N = number of steps            •   R07 = z'0
•   R03 = z0
Flags: /
Subroutine:  a program which calculates f(x;y;z) in Y-register and g(x;y;z) in X-register
assuming x is in X-register , y is in Y-register and z is in Z-register upon entry.

-Line 96 is a three-byte  GTO 01

 01  LBL "RKS2"   02  RCL 05   03  STO 12   04  LBL 01   05  RCL 03   06  RCL 02   07  RCL 01   08  XEQ IND 00   09  RCL 04   10  ST+ 01   11  ST* Z   12  *   13  STO 09   14  2   15  /   16  RCL 07   17  +   18  RCL 04   19  *   20  RCL 03 21  +   22  X<>Y   23  STO 08   24  2   25  /   26  RCL 06   27  +   28  RCL 04   29  *   30  RCL 02   31  +   32  RCL 01   33  XEQ IND 00   34  RCL 04   35  ST+ 01   36  ST* Z   37  *   38  STO 11   39  RCL 07   40  + 41  X<>Y   42  STO 10   43  RCL 06   44  +   45  RCL 04   46  ST+ X   47  ST* Z   48  *   49  RCL 03   50  ST+ Z   51  CLX   52  RCL 02   53  +   54  RCL 01   55  XEQ IND 00   56  RCL 04           57  ST* Z   58  *   59  RCL 09   60  + 61  RCL 11   62  ST+ X   63  ST+ 09   64  ST+ X   65  +   66  3   67  ST/ 09   68  /   69  X<> 07   70  ST+ 07   71  RCL 09           72  +   73  X<>Y   74  RCL 08   75  +   76  RCL 10   77  ST+ X   78  ST+ 08   79  ST+ X   80  + 81  3   82  ST/ 08   83  /   84  X<> 06   85  ST+ 06   86  RCL 08           87  +   88  RCL 04   89  ST+ X   90  ST* Z   91  *   92  ST+ 02   93  X<>Y   94  ST+ 03   95  DSE 12   96  GTO 01   97  RCL 03   98  RCL 02   99  RCL 01 100  END

( 139 bytes / SIZE 013 )

 STACK INPUTS OUTPUTS Z / z(x0+N.h) Y / y(x0+N.h) X / x0+N.h

Example:

d2y/dx2 =  -y.z                  y(0) = 2    y'(0) = 1         Evaluate  y(1) , z(1) , y'(1) , z'(1)    with  h = 0.1  &  N = 10
d2z/dx2 =  x.( y + z )         z(0) = 1    z'(0) = 1

 01  LBL "T"  02  RDN  03  STO Z  04  X<>Y  05  CHS  06  ST* Z  07  -  08  R^  09  *  10  RTN

"T"  ASTO 00

0     STO 01              1  STO 06
2     STO 02              1  STO 07
1     STO 03
0.05  STO 04
10    STO 05

XEQ "RKS2"  >>>>     x   =  1                                                                                                  ( in 41 seconds )
RDN   y(1) = 1.531358015     and   R06 = y'(1) = -2.312838895
RDN   z(1) = 2.620254480              R07 = z'(1) =  2.941751649

-With  h = 0.05  it yields   y(1) =  1.531356736      y'(1) = -2.312840085
z(1) =  2.620254295      z'(1) =  2.941748608

-Actually, the exact results - rounded to 9D - are

y(1) =  1.531356646      y'(1) = -2.312840137
z(1) =  2.620254281      z'(1) =  2.941748399

-To continue the calculations, simply press R/S  ( after changing h/2 & N in registers R04 & R05 if needed )

c) Systems of 3 Second-order differential Equations

-"RKS3" solves the system:

y" = f(x,y,z,u)         y(x0) = y0     y'(x0) = y'0
z" = g(x,y,z,u)         z(x0) = z0     z'(x0) = z'0
u" = h(x,y,z,u)        u(x0) = u0     u'(x0) = u'0

Data Registers:    •   R00:  subroutine name                     ( Registers R00 thru R09  are to be initialized before executing "RKS3" )

•   R01 = x0      •   R05 = h/2 = half of the stepsize       •   R07 = y'0
•   R02 = y0      •   R06 = N = number of steps            •   R08 = z'0
•   R03 = z0                   R10 to R16: temp                  •   R09 = u'0
•   R04 = u0
Flags: /
Subroutine:  a program which calculates f(x;y;z;u) in Z-register , g(x;y;z;u) in Y-register and  h(x;y;z;u) in X-register
assuming x is in X-register , y is in Y-register , z is in Z-register and u is in T-register upon entry.

-Line 136 is a three-byte  GTO 01

 01  LBL "RKS3"   02  RCL 06   03  STO 16   04  LBL 01   05  RCL 04   06  RCL 03   07  RCL 02   08  RCL 01   09  XEQ IND 00   10  RCL 05   11  ST+ 01   12  ST* Z   13  ST* T   14  *   15  STO 12   16  2   17  /   18  RCL 09   19  +   20  RCL 05   21  *   22  RCL 04   23  +   24  X<>Y   25  STO 11   26  2   27  /   28  RCL 08   29  + 30  RCL 05   31  *   32  RCL 03   33  +   34  R^   35  STO 10   36  2   37  /   38  RCL 07   39  +   40  RCL 05   41  *   42  RCL 02   43  +   44  RCL 01   45  XEQ IND 00   46  RCL 05   47  ST+ 01   48  ST* Z   49  ST* T   50  *   51  STO 15   52  RCL 09   53  +   54  X<>Y   55  STO 14   56  RCL 08   57  +   58  R^ 59  STO 13   60  RCL 07   61  +   62  RCL 05   63  ST+ X   64  ST* Z   65  ST* T   66  *   67  RCL 04   68  ST+ T   69  CLX   70  RCL 03   71  ST+ Z   72  CLX   73  RCL 02   74  +   75  RCL 01   76  XEQ IND 00   77  RCL 05   78  ST* Z   79  ST* T   80  *   81  RCL 12   82  +   83  RCL 15   84  ST+ X   85  ST+ 12   86  ST+ X   87  + 88  3   89  ST/ 12   90  /   91  X<> 09   92  ST+ 09   93  RCL 12   94  +   95  X<>Y   96  RCL 11   97  +   98  RCL 14   99  ST+ X 100  ST+ 11 101  ST+ X 102  + 103  3 104  ST/ 11 105  / 106  X<> 08  107  ST+ 08 108  RCL 11         109  + 110  R^ 111  RCL 10 112  + 113  RCL 13 114  ST+ X 115  ST+ 10 116  ST+ X 117  + 118  3 119  ST/ 10 120  / 121  X<> 07 122  ST+ 07 123  RCL 10 124  + 125  RCL 05 126  ST+ X 127  ST* Z 128  ST* T 129  * 130  ST+ 02 131  X<>Y 132  ST+ 03 133  R^ 134  ST+ 04 135  DSE 16  136  GTO 01  137  RCL 04         138  RCL 03 139  RCL 02 140  RCL 01 141  END

( 194 bytes / SIZE 017 )

 STACK INPUTS OUTPUTS T / u(x0+N.h) Z / z(x0+N.h) Y / y(x0+N.h) X / x0+N.h

Example:

d2y/dx2 = -y.z.u                    y(0) = 1        y'(0) = 1         Evaluate  y(1) , z(1) , u(1)  ,   y'(1) , z'(1) , u'(1)      with  h = 0.1  &  N = 10
d2z/dx2 = x.( y + z - u )         z(0) = 1        z'(0) = 1
d2u/dx2 = x.y - z.u                 u(0) = 2       u'(0) = 1

LBL "T"       R^               RCL Z         X<> Z         LASTX            "T"  ASTO 00
R^               ST* Y         R^                ST+ Y         *
CHS            ST+ L         ST+ L          ST* Z          X<>Y
STO L         CLX           ST* Y          X<> T         RTN

0     STO 01
1     STO 02  STO 03  STO 07  STO 08  STO 09
2     STO 04
0.05  STO 05
10    STO 06    XEQ "RKS3"  >>>>     x   = 1                                                                                             ( in 65 seconds )
RDN   y(1) = 0.439528419           y'(1) = -2.101120400 = R07
RDN   z(1) = 2.070938499           z'(1) =   1.269599239 = R08
RDN   u(1) = 1.744522976           u'(1) = -1.704232092 = R09

-With  h = 0.05  it yields

y(1) = 0.439524393          y'(1) = -2.101122784
z(1) = 2.070940521          z'(1) =   1.269597110
u(1) = 1.744524843          u'(1) = -1.704234567

-Actually, the exact results rounded to 9D are:

y(1) = 0.439524100          y'(1) = -2.101122880
z(1) = 2.070940654          z'(1) =   1.269596950
u(1) = 1.744524964          u'(1) = -1.704234756

-Press R/S  to continue the calculations ( after changing h & N in registers R05 & R06 if needed )

2°)  A Runge-Kutta-Nystrom Method with order 6

a) Second-order differential Equations

-Five stages are sufficient to get a 6-th order formula.
-The following one is due to Albrecht:

f0 = f(x0,y0)
f1 = f(x0+ h/4 , y0+ (h/4) y'0+ (h2/32) f0)
f2 = f(x0+ h/2 , y0+ (h/2) y'0+ h2 (-f0/24 + f1/6 )
f3 = f(x0+ 3h/4 , y0+ (3h/4) y'0+ h2 (3f0/32 + f1/8 + f2/16 )
f4 = f(x0+ h , y0+ h y'0+ h2 (3f1/7 - f2/14 + f3/7 )

-Then,    y = y0 + h y'0+ h2 ( 7 f0/90 + 4 f1/15 + f2/15 + 4 f3/45 ) + O(h7)
and      y' = y'0 + h ( 7 f0/90 + 16 f1/45 + 2 f2/15 + 16 f3/45 + 7 f4/90 ) + O(h7)

Data Registers:    •   R00:  subroutine name                    ( Registers R00 thru R05  are to be initialized before executing "RKN6" )

•   R01 = x0      •   R04 = h  = stepsize        R06 to R10: temp
•   R02 = y0      •   R05 = N = number of steps
•   R03 = y'0
Flags: /
Subroutine:  a program which calculates f(x,y) in X-register, assuming  x and y  are in registers X and Y upon entry.

-Line 144 is a three-byte  GTO 01

 01  LBL "RKN6"   02  RCL 05   03  STO 10   04  LBL 01   05  RCL 02   06  RCL 01   07  XEQ IND 00   08  RCL 04   09  *   10  STO 06   11  8   12  /   13  RCL 03   14  +   15  RCL 04   16  *   17  4   18  /   19  RCL 02   20  +   21  RCL 04   22  4   23  /   24  RCL 01   25  +   26  XEQ IND 00    27  RCL 04   28  *   29  STO 07   30  6 31  /   32  RCL 06   33  24   34  /    35  -   36  RCL 03   37  2   38  /   39  +   40  RCL 04   41  *   42  RCL 02   43  +   44  RCL 04   45  2   46  /   47  RCL 01   48  +   49  XEQ IND 00   50  RCL 04   51  *   52  STO 08   53  ST+ X   54  RCL 07   55  4   56  *   57  +   58  RCL 06   59  3   60  * 61  +   62  8   63  /   64  RCL 03   65  3   66  *   67  +   68  RCL 04   69  *   70  4   71  /   72  RCL 02   73  +   74  RCL 04   75  .75   76  *   77  RCL 01   78  +   79  XEQ IND 00   80  RCL 04   81  *   82  STO 09   83  ST+ X   84  RCL 08   85  -   86  RCL 07   87  6   88  *   89  +   90  14 91  /   92  RCL 03   93  +   94  RCL 04   95  ST+ 01   96  *   97  RCL 02   98  +   99  RCL 01 100  XEQ IND 00 101  RCL 04 102  * 103  RCL 06 104  + 105  7 106  * 107  RCL 07 108  RCL 09 109  + 110  32 111  * 112  + 113  RCL 08 114  12 115  * 116  + 117  RCL 09 118  8 119  * 120  RCL 08 121  6 122  * 123  + 124  RCL 07 125  24 126  * 127  + 128  RCL 06         129  7 130  * 131  + 132  90 133  ST/ Z 134  / 135  RCL 03 136  + 137  RCL 04 138  * 139  ST+ 02 140  X<>Y 141  ST+ 03 142  DSE 10 143  GTO 01 144  RCL 03 145  RCL 02 146  RCL 01 147  END

( 178 bytes / SIZE 011 )

 STACK INPUTS OUTPUTS Z / y'(x0+N.h) Y / y(x0+N.h) X / x0+N.h

Example:     y(0) = 1  ,   y'(0) = 1  ,   y" = - y ( x2 + y2 ) 1/2

>>> Calculate  y(1) & y'(1)

-Load this short routine:

 01  LBL "T"  02  X^2  03  X<>Y  04  STO Z  05  X^2  06  +  07  SQRT  08  *  09  CHS  10  RTN

"T"  ASTO 00

0  STO 01  STO 03  1  STO 02

-With h = 0.1  and  N = 10  ,  0.1  STO 04  10  STO 05

XEQ "RKN6"   >>>>     x   = 1                                            ---Execution time = 76s---
RDN   y(1) =  0.536630617
RDN   y'(1) = -0.860171927

-The exact results, rounded to 10D are  y(1) =  0.5366306164  &   y'(1) = -0.8601719268

Note:

-Press R/S  to continue the calculations ( after changing h & N in registers R04 & R05 if need be )

b) Systems of 2 Second-order differential Equations

-"RKN6B" solves the system:

y" = f(x,y,z)         y(x0) = y0     y'(x0) = y'0
z" = g(x,y,z)         z(x0) = z0     z'(x0) = z'0

Data Registers:    •   R00:  subroutine name                     ( Registers R00 thru R07  are to be initialized before executing "RKN6B" )

•   R01 = x0      •   R04 = h =  stepsize                        •   R06 = y'0          R08 to R16: temp
•   R02 = y0      •   R05 = N = number of steps            •   R07 = z'0
•   R03 = z0
Flags: /
Subroutine:  a program which calculates f(x;y;z) in Y-register and g(x;y;z) in X-register
assuming x is in X-register , y is in Y-register and z is in Z-register upon entry.

-Line 255 is a three-byte  GTO 01 ( or replace it by  GTO 16  and line 04 by  LNL 16 )

 01  LBL "RKN6B"   02  RCL 05   03  STO 16   04  LBL 01   05  RCL 03   06  RCL 02   07  RCL 01   08  XEQ IND 00   09  RCL 04             10  ST* Z   11  *   12  STO 12   13  8   14  /   15  RCL 07   16  +   17  4   18  /   19  RCL 04   20  *   21  RCL 03   22  +   23  X<>Y   24  STO 08   25  8   26  /   27  RCL 06   28  +   29  4   30  /   31  RCL 04   32  *   33  RCL 02   34  +   35  RCL 04   36  4   37  /   38  RCL 01   39  +   40  XEQ IND 00   41  RCL 04   42  ST* Z   43  *   44  STO 13 45  6   46  /   47  RCL 12   48  24   49  /   50  -   51  RCL 07             52  2   53  /   54  +   55  RCL 04   56  *   57  RCL 03   58  +   59  X<>Y   60  STO 09   61  6   62  /   63  RCL 08   64  24   65  /   66  -   67  RCL 06   68  2   69  /   70  +   71  RCL 04   72  *   73  RCL 02   74  +   75  RCL 04   76  2   77  /   78  RCL 01   79  +   80  XEQ IND 00   81  RCL 04   82  ST* Z   83  *   84  STO 14   85  ST+ X   86  RCL 13   87  4   88  * 89  +   90  RCL 12   91  3   92  *   93  +   94  32   95  /   96  RCL 07             97  .75   98  STO 15   99  * 100  + 101  RCL 04 102  * 103  RCL 03 104  + 105  X<>Y 106  STO 10 107  ST+ X 108  RCL 09 109  4 110  * 111  + 112  RCL 08 113  3 114  * 115  + 116  32 117  / 118  RCL 06 119  RCL 15 120  * 121  + 122  RCL 04 123  * 124  RCL 02 125  + 126  RCL 04 127  RCL 15 128  * 129  RCL 01 130  + 131  XEQ IND 00 132  RCL 04 133  ST* Z 134  * 135  STO 15 136  ST+ X 137  RCL 14 138  - 139  RCL 13           140  6 141  * 142  + 143  14 144  / 145  RCL 07 146  + 147  RCL 04 148  * 149  RCL 03 150  + 151  X<>Y 152  STO 11 153  ST+ X 154  RCL 10 155  - 156  RCL 09 157  6 158  * 159  + 160  14 161  / 162  RCL 06 163  + 164  RCL 04 165  ST+ 01 166  * 167  RCL 02 168  + 169  RCL 01 170  XEQ IND 00 171  RCL 04 172  ST* Z 173  * 174  RCL 12 175  + 176  7 177  * 178  RCL 13 179  RCL 15           180  + 181  32 182  * 183  + 184  RCL 14         185  12 186  * 187  + 188  X<>Y 189  RCL 08 190  + 191  7 192  * 193  RCL 09 194  RCL 11 195  + 196  32 197  * 198  + 199  RCL 10 200  12 201  * 202  + 203  X<> 11 204  8 205  * 206  RCL 10 207  6 208  * 209  + 210  RCL 09 211  24 212  * 213  + 214  RCL 08 215  7 216  * 217  + 218  90 219  ST/ Z 220  ST/ 11 221  / 222  RCL 06           223  + 224  RCL 04 225  * 226  ST+ 02 227  X<> 11 228  ST+ 06 229  X<>Y 230  X<> 15 231  8 232  * 233  RCL 14 234  6 235  * 236  + 237  RCL 13 238  24 239  * 240  + 241  RCL 12 242  7 243  * 244  + 245  90 246  / 247  RCL 07 248  + 249  RCL 04 250  * 251  ST+ 03 252  X<> 15 253  ST+ 07 254  DSE 16 255  GTO 01 256  RCL 03 257  RCL 02 258  RCL 01 259  END

( 314 bytes / SIZE 017 )

 STACK INPUTS OUTPUTS Z / z(x0+N.h) Y / y(x0+N.h) X / x0+N.h

Example:

d2y/dx2 =  -y.z                  y(0) = 2    y'(0) = 1         Evaluate  y(1) , z(1) , y'(1) , z'(1)    with  h = 0.1  &  N = 10
d2z/dx2 =  x.( y + z )         z(0) = 1    z'(0) = 1

 01  LBL "T"  02  RDN  03  STO Z  04  X<>Y  05  CHS  06  ST* Z  07  -  08  R^  09  *  10  RTN

"T"  ASTO 00

0     STO 01              1  STO 06
2     STO 02              1  STO 07
1     STO 03
0.1    STO 04
10    STO 05

XEQ "RKN6B"  >>>>     x   =  1                                                                                            ---Execution time = 116s---
RDN   y(1) = 1.531356647     and   R06 = y'(1) = -2.312840139
RDN   z(1) = 2.620254282              R07 = z'(1) =  2.941748401

-The precision is almost perfect, since the exact results - rounded to 9D - are

y(1) =  1.531356646      y'(1) = -2.312840137
z(1) =  2.620254281      z'(1) =  2.941748399

-To continue the calculations, simply press R/S  ( after changing h & N in registers R04 & R05 if need be )

3°)  S-Stage Explicit Runge-Kutta-Nystrom formulae

a) Second-order differential Equations

-All s-stage explicit Runge-Kutta-Nystrom formulae without error estimate are determined by ( s2 + 5.s - 2 ) / 2 coefficients ai,j ; bi ; ci

k1 = h.f(x;y)
k2 = h.f(x+c2.h;y+c2.h.y'+h a2,1.k1)
k3 = h.f(x+c3.h;y+c3.h.y'+h(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) )

then,     y(x+h) = y(x) + h y'(x) + h [ b1.k1+ ................ + bs.ks ]   and   y'(x+h) = y'(x) + [ b'1.k1+ ................ + b's.ks ]

-Therefore, a single program may be used for all explicit Runge-Kutta-Nystrom methods,
provided the coefficients are stored in appropriate registers.

Data Registers:      •   R00:  f name         ( Registers R00 to R06  and  Rs+10 to R(s2+7.s+16)/2  are to be initialized before executing "ERKN" )

•   R01 = x0                                      R07 to R09:  temp          •  Rs+10 to R(s2+7.s+16)/2 =  the ( s2 + 5.s - 2 ) / 2 coefficients ai,j ; bi ; ci
•   R02 = y0                                      R10 = k1                                                                            which are to be stored as shown below:
•   R03 = y'0                                     R11 = k2
•   R04 = h = stepsize                        ............
•   R05 = N = number of steps
•   R06 = s  = number of stages        Rs+9 = ks

•  Rs+10 = a2,1                                                                                •  Rs+11 = c2
•  Rs+12 = a3,1        •  Rs+13 = a3,2                                                •  Rs+14 = c3

.........................................................................                                 ...................

•  R(s2+s+18)/2 = as,1  ................   •  R(s2+3s+14)/2 = as,s-1       •  R(s2+3s+16)/2 = cs

•  R(s2+3s+18)/2 = b1 ............................................................       •  R(s2+5s+16)/2 = bs
•  R(s2+5s+18)/2 = b'1 ...........................................................       •  R(s2+7s+16)/2 = b's

Flags: /
Subroutine:  a program calculating f(x;y) assuming x is in X-register and y is in Y-register upon entry.

-Line 97 is a three-byte GTO 01

 01  LBL "ERKN"   02  RCL 05   03  STO 07   04  LBL 01   05  10   06  STO 08   07  3   08  RCL 06   09  ST+ Y   10  *   11  2   12  /   13  8   14  +   15   E3   16  /   17  +   18  RCL 06   19  +   20  STO 09   21  RCL 02 22  RCL 01   23  XEQ IND 00   24  RCL 04   25  *   26  STO 10   27  LBL 02   28  RCL 08   29  INT   30   E3   31  /   32  10   33  +   34  STO 08   35  CLX   36  LBL 03   37  RCL IND 08   38  RCL IND 09   39  *   40  +   41  ISG 09   42  ISG 08 43  GTO 03   44  RCL 03   45  RCL IND 09   46  *   47  +   48  RCL 04   49  *   50  RCL 02   51  +   52  RCL 04   53  RCL IND 09   54  *   55  RCL 01   56  +   57  XEQ IND 00   58  RCL 04   59  *   60  STO IND 08   61  ISG 09   62  GTO 02   63  RCL 06 64  1.001   65  -   66  ST- 08   67  CLX   68  LBL 04   69  RCL IND 08   70  RCL IND 09   71  *   72  +   73  ISG 09   74  CLX   75  ISG 08   76  GTO 04   77  RCL 03   78  +   79  RCL 04   80  ST+ 01   81  *   82  ST+ 02   83  RCL 06   84  ST- 08 85  CLX   86  LBL 05   87  RCL IND 08   88  RCL IND 09   89  *   90  +   91  ISG 09   92  CLX   93  ISG 08   94  GTO 05   95  ST+ 03   96  DSE 07   97  GTO 01   98  RCL 03   99  RCL 02 100  RCL 01 101  END

( 149 bytes / SIZE ( s^2+7.s+18 )/2 )

 STACK INPUTS OUTPUTS Z / y'(x0+N.h) Y / y(x0+N.h) X / x0+N.h

Example:     y(0) = 1  ,   y'(0) = 1  ,   y" = - y ( x2 + y2 ) 1/2

>>> Calculate  y(1) & y'(1)

 01  LBL "T"  02  X^2  03  X<>Y  04  STO Z  05  X^2  06  +  07  SQRT  08  *  09  CHS  10  RTN

"T"  ASTO 00

0  STO 01  STO 03  1  STO 02

-With h = 0.1  and  N = 10  ,  0.1  STO 04  10  STO 05

-If we use the same ( Albrecht ) 5-stage method, store the following coefficients into R15 thru R38

1/32                              1/4                             R15                            R16
-1/24   1/6                     1/2                             R17   R18                   R19
3/32    1/8   1/16           3/4            =              R20  R21  R22            R23
0       3/7  -1/14   1/7    1                              R24  R25  R26  R27   R28
7/90  4/15   1/15  4/45    0                             R29  R30  R31  R32   R33
7/90  16/45 2/15 16/45 7/90                          R34  R35  R36  R37   R38

-Since it is a 5-stage method,    5  STO 06   and

XEQ "ERKN"   >>>>     x   = 1                                            ---Execution time = 147s---
RDN   y(1) =  0.536630617
RDN   y'(1) = -0.860171927

Notes:

-Press R/S  to continue the calculations ( after changing h & N in registers R04 & R05 if need be )
-LBL 04 & LBL 05 are identical, so they could form a subroutine:
-The program would be shorter but slower.

Control Number of the Array of Coefficients for an S-stage method:

 S bbb.eee 5         6         7         8         9        10        11        12        13        14        15        16        17 15.038    16.047    17.057    18.068    19.080    20.093    21.107    22.122    23.138    24.155    25.173    26.192    27.212

-For instance, we can use the coefficients below - taken from reference [3] -  to get a 13-stage 10th-order method.
-These 116 coefficients are to be stored into R23 thru R138 in the following order:

a( 2, 1) =     6.2550354381669436926370260206784041201048E-04  = R23
c( 2)    =     3.5369578561715839852580662156128003290245E-02

a( 3, 1) =     8.3400472508892582568493680275712054934730E-04
a( 3, 2) =     1.6680094501778516513698736055142410986946E-03
c( 3)    =     7.0739157123431679705161324312256006580491E-02

a( 4, 1) =     1.4377592173651130455265029764307572975032E-02
a( 4, 2) =    -2.5520389586889438033660194426507730127969E-02
a( 4, 3) =     2.9955419091121746433503478890580675956400E-02
c( 4)    =     1.9397227470895648005755022968171580307125E-01

a( 5, 1) =     5.6606010103291956408355747332339448084576E-03
a( 5, 2) =     0.0000000000000000000000000000000000000000E+00
a( 5, 3) =     2.2438147205568485440414320313187753497836E-02
a( 5, 4) =     9.6198950134107937593128496677833016937067E-03
c( 5)    =     2.7465849059990290000000000000000000000000E-01

a( 6, 1) =     4.2497841145700712790685880544080818730288E-03
a( 6, 2) =     0.0000000000000000000000000000000000000000E+00
a( 6, 3) =     1.3737389558120657344878988260209037684744E-02
a( 6, 4) =     2.1090823734166889521084117969028056034572E-03
a( 6, 5) =    -2.1698166033104208350129690333992516122959E-04
c( 6)    =     1.9939545825206940000000000000000000000000E-01

a( 7, 1) =     8.4469918165805880106087924734363680346449E-04
a( 7, 2) =     0.0000000000000000000000000000000000000000E+00
a( 7, 3) =     7.6286170426770036487508428629504122228445E-04
a( 7, 4) =    -4.5274889794326362037916138257410810321144E-03
a( 7, 5) =    -9.4704498758374659029865654611073300453711E-05
a( 7, 6) =     4.3839567862160003018135640613696763068192E-03
c( 7)    =     5.2332097109723180000000000000000000000000E-02

a( 8, 1) =     4.0257188676144292921381323534074832686999E-02
a( 8, 2) =     0.0000000000000000000000000000000000000000E+00
a( 8, 3) =     2.8132558570671430000000000000000000000000E-01
a( 8, 4) =    -9.9086973311055342850364525233745309401480E-02
a( 8, 5) =     3.1558678646031991227831775524256371039974E-02
a( 8, 6) =     7.6263009309052972315471957912384344318273E-02
a( 8, 7) =    -2.4509110177537917817170493451316523864377E-01
c( 8)    =     4.1285926718800690000000000000000000000000E-01

a( 9, 1) =    -5.1852206149158582964016598801731374371401E-01
a( 9, 2) =     0.0000000000000000000000000000000000000000E+00
a( 9, 3) =    -4.0220496748396490000000000000000000000000E+00
a( 9, 4) =     1.3339542387088720000000000000000000000000E+00
a( 9, 5) =    -3.6280137620198678475540384353108923505414E-01
a( 9, 6) =    -4.4662932597538443824830112336373381632573E-01
a( 9, 7) =     4.1093353899979717196794462089910810944472E+00
a( 9, 8) =     8.3371860107714029551447509208700700646694E-02
c( 9)    =     5.9440567007045230000000000000000000000000E-01

a(10, 1) =     4.5526515039304022037674417630889739137541E-01
a(10, 2) =     0.0000000000000000000000000000000000000000E+00
a(10, 3) =     3.4410244296399106992261498105244921193431E+00
a(10, 4) =    -1.1420304634021581044575027059094514599562E+00
a(10, 5) =     3.3122504529344460041163434350588640458224E-01
a(10, 6) =     5.2510814510725835145952044841619177157062E-01
a(10, 7) =    -3.4007204966989148959577640759197726181960E+00
a(10, 8) =     1.4959025547025800806021814497275304454960E-02
a(10, 9) =     1.7714834024190792778113334621606086825758E-02
c(10)    =     6.9648498893199050000000000000000000000000E-01

a(11, 1) =    -5.0738859071291316865190644559724141522647E-02
a(11, 2) =     0.0000000000000000000000000000000000000000E+00
a(11, 3) =    -8.5258315802766612821674340523993103859370E-01
a(11, 4) =     2.9256281989303017562028747765972648095241E-01
a(11, 5) =    -4.2631304538837970000000000000000000000000E-01
a(11, 6) =     3.0845126792446679274386181736916000785280E-01
a(11, 7) =     8.2068063064681728553592877499937152955118E-01
a(11, 8) =     2.6353200561785124318114669620069412596405E-01
a(11, 9) =    -3.8002959536498107924589290701514760515073E-02
a(11,10) =     5.0836949957876193531014830503462796310983E-02
c(11)    =     8.5840043338316930000000000000000000000000E-01

a(12, 1) =    -8.5506341904459305448936612697062868735670E-01
a(12, 2) =     0.0000000000000000000000000000000000000000E+00
a(12, 3) =    -4.4172967833111218167905547111397401957735E+00
a(12, 4) =     1.4620414738250776987869127049027576880060E+00
a(12, 5) =     1.5800603670627463865955782231377536485336E+00
a(12, 6) =    -2.1500473087668603894015427990085825826253E+00
a(12, 7) =     5.2192882951843358632420497699580187954903E+00
a(12, 8) =    -7.0122224598841812446074882425859149461800E-01
a(12, 9) =     4.0674284722476557322266460477065027849355E-01
a(12,10) =    -1.0995016400771572122089410865399735394434E-01
a(12,11) =     2.5498951087297871672803054502111066906939E-02
c(12)    =     9.5922053070763063933434705195204984252787E-01

a(13, 1) =     3.6124205767243278950403294321213352485132E+00
a(13, 2) =     0.0000000000000000000000000000000000000000E+00
a(13, 3) =     1.9614376424164347318861614225759179537283E+01
a(13, 4) =    -6.5540954527470471068838634201654381264077E+00
a(13, 5) =    -5.3634775178894355119823336908885299369771E+00
a(13, 6) =     8.9549200633414078297191926517465420210063E+00
a(13, 7) =    -2.2111999575685298066957085805778892539418E+01
a(13, 8) =     3.0666418328538943745146953495973014350614E+00
a(13, 9) =    -1.2974380337600327021222956536808230517716E+00
a(13,10) =     5.9822684827188431105752038905995447885348E-01
a(13,11) =    -2.4107078892562108246647516780668823095046E-02
a(13,12) =     4.5319136185137669988740390100397569527517E-03
c(13)    =     1.0000000000000000000000000000000000000000E+00

b( 1)    =     1.1445045431083081161076675575316257764110E-02
b( 2)    =     0.0000000000000000000000000000000000000000E+00
b( 3)    =     0.0000000000000000000000000000000000000000E+00
b( 4)    =     0.0000000000000000000000000000000000000000E+00
b( 5)    =     0.0000000000000000000000000000000000000000E+00
b( 6)    =     1.5182228814165001267592100569234113490800E-01
b( 7)    =     9.3833383282371058262139365435233240765713E-02
b( 8)    =     1.3138871401731356660181720771852165705253E-01
b( 9)    =     4.2452793993460570894743314052986770767710E-02
b(10)    =     4.6614359052634087726314462069090004222681E-02
b(11)    =     1.9739340751760337610363200447178885100858E-02
b(12)    =     2.7040753297272850676247690093320494184034E-03
b(13)    =     0.0000000000000000000000000000000000000000E+00

bp( 1)   =     1.1445045431083081161076675575316257764110E-02
bp( 2)   =     0.0000000000000000000000000000000000000000E+00
bp( 3)   =     0.0000000000000000000000000000000000000000E+00
bp( 4)   =     0.0000000000000000000000000000000000000000E+00
bp( 5)   =     0.0000000000000000000000000000000000000000E+00
bp( 6)   =     1.8963455766836141912201425856209518713321E-01
bp( 7)   =     9.9015048411147152930074891966312987547087E-02
bp( 8)   =     2.2377720821386995442668671882695678061627E-01
bp( 9)   =     1.0466811506175315699681616849938105031123E-01
bp(10)   =     1.5358172529459859690396485461784888933400E-01
bp(11)   =     1.3940255061073114260364669651015505753606E-01
bp(12)   =     6.6309723413523447970758037743751771753948E-02
bp(13)   =     1.2166025894932047884961697698182018004080E-02  = R138

-With the same differential equation, h = 0.1 & N = 10 ( don't forget to store 13 into R06 ), it yields:

XEQ "ERKN"   >>>>     x   = 1                                            ---Execution time = 9m22s---
RDN   y(1) =  0.5366306165
RDN   y'(1) = -0.8601719269

-The exact results, rounded to 10D are  y(1) =  0.5366306164  &   y'(1) = -0.8601719268
-So, the error is only 1 unit in the last decimal.

-Of course, it's not very fast with a real HP-41, but with a good emulator or with a 41CL...

b) Systems of 2 Second-order differential Equations

-A similar program may be written to solve the system   y" = f(x,y,z)   z" = g(x,y,z)

Data Registers:      •   R00:  fg name         ( Registers R00 to R08  and  R13+2s to R(s2+9.s+22)/2  are to be initialized before executing "ERKN2" )

•   R01 = x0                                      R09 to R12:  temp          •  R13+2s to R(s2+9.s+22)/2 =  the ( s2 + 5.s - 2 ) / 2 coefficients ai,j ; bi ; ci
•   R02 = y0                                      R13 to R12+2s = kyi &  kzi                                                which are to be stored as shown below:
•   R03 = z0
•   R04 = h = stepsize
•   R05 = N = number of steps
•   R06 = y'0
•   R07 = z'0
•   R08 = s  = number of stages

•  R13+2s = a2,1                                                                             •  R14+2s = c2
•  R15+2s = a3,1        •  R16+2s = a3,2                                           •  R17+2s = c3

.........................................................................                                 ...................

•  R(s2+3s+24)/2 = as,1  ................   •  R(s2+5s+20)/2 = as,s-1       •  R(s2+5s+22)/2 = cs

•  R(s2+5s+24)/2 = b1 ............................................................       •  R(s2+7s+22)/2 = bs
•  R(s2+7s+24)/2 = b'1 ...........................................................       •  R(s2+9s+22)/2 = b's

Flags: /
Subroutine:  a program that calculates f(x;y;z) in Y-register and g(x;y;z) in X-register
assuming x is in X-register , y is in Y-register and z is in Z-register upon entry.

-Line 153 is a three-byte GTO 01

 01  LBL "ERKN2"   02  RCL 05   03  STO 09             04  LBL 01   05  13   06  STO 10   07  RCL 08   08  +   09  STO 11   10  5   11  LASTX   12  ST+ Y   13  *   14  2   15  /   16  11   17  +   18   E3   19  /   20  +   21  RCL 08   22  +   23  STO 12   24  RCL 03   25  RCL 02   26  RCL 01   27  XEQ IND 00   28  RCL 04   29  ST* Z   30  *   31  STO IND 11   32  X<>Y 33  STO 13   34  LBL 02   35  RCL 10             36  INT   37   E3   38  /   39  13   40  +   41  STO 10   42  RCL 08   43  +   44  STO 11   45  CLST   46  LBL 03   47  X<>Y   48  RCL IND 10   49  RCL IND 12   50  *   51  +   52  X<>Y   53  RCL IND 11   54  RCL IND 12   55  *   56  +   57  ISG 12   58  ISG 11   59  CLX   60  ISG 10   61  GTO 03   62  RCL 07   63  RCL IND 12   64  * 65  +   66  RCL 04             67  *   68  RCL 03   69  +   70  X<>Y   71  RCL 06   72  RCL IND 12   73  *   74  +   75  RCL 04   76  *   77  RCL 02   78  +   79  RCL 04   80  RCL IND 12   81  *   82  RCL 01   83  +   84  XEQ IND 00   85  RCL 04   86  ST* Z   87  *   88  STO IND 11   89  X<>Y   90  STO IND 10   91  ISG 12   92  GTO 02   93  RCL 08   94  1.001   95  -   96  ST- 10 97  ST- 11   98  CLST   99  LBL 04 100  X<>Y 101  RCL IND 10 102  RCL IND 12 103  * 104  + 105  X<>Y 106  RCL IND 11 107  RCL IND 12 108  * 109  + 110  ISG 12 111  CLX 112  ISG 11 113  CLX 114  ISG 10 115  GTO 04 116  RCL 07           117  + 118  RCL 04 119  * 120  ST+ 03 121  CLX 122  RCL 06 123  + 124  RCL 04 125  ST+ 01 126  * 127  ST+ 02 128  RCL 08 129  ST- 10 130  ST- 11 131  CLST 132  LBL 05 133  X<>Y 134  RCL IND 10 135  RCL IND 12 136  * 137  + 138  X<>Y 139  RCL IND 11 140  RCL IND 12 141  * 142  + 143  ISG 12 144  CLX 145  ISG 11 146  CLX 147  ISG 10 148  GTO 05 149  ST+ 07 150  X<>Y 151  ST+ 06 152  DSE 09 153  GTO 01 154  RCL 03           155  RCL 02 156  RCL 01 157  END

( 225 bytes / SIZE ( s^2+9.s+24 )/2 )

 STACK INPUTS OUTPUTS Z / z(x0+N.h) Y / y(x0+N.h) X / x0+N.h

Example:

d2y/dx2 =  -y.z                  y(0) = 2    y'(0) = 1         Evaluate  y(1) , z(1) , y'(1) , z'(1)    with  h = 0.1  &  N = 10
d2z/dx2 =  x.( y + z )         z(0) = 1    z'(0) = 1

 01  LBL "T"  02  RDN  03  STO Z  04  X<>Y  05  CHS  06  ST* Z  07  -  08  R^  09  *  10  RTN

"T"  ASTO 00

0     STO 01              1  STO 06
2     STO 02              1  STO 07
1     STO 03
0.1    STO 04
10    STO 05

-If we again use the same ( Albrecht ) 5-stage method, store the following coefficients into R23 thru R46

1/32                              1/4                             R23                            R24
-1/24   1/6                     1/2                             R25  R26                    R27
3/32    1/8   1/16           3/4            =              R28  R29  R30            R31
0       3/7  -1/14   1/7    1                              R32  R33  R34  R35   R36
7/90  4/15   1/15  4/45    0                             R37  R38  R39  R40   R41
7/90  16/45 2/15 16/45 7/90                          R42  R43  R44  R45   R46

-Since it is a 5-stage method,    5  STO 08   and

XEQ "ERKN2"  >>>>     x   =  1                                                                                            ---Execution time = 3m41s---
RDN   y(1) = 1.531356646     and   R06 = y'(1) = -2.312840140
RDN   z(1) = 2.620254282              R07 = z'(1) =  2.941748401

Note:

-To continue the calculations, simply press R/S  ( after changing h & N in registers R04 & R05 if need be )

Control Number of the Array of Coefficients for an S-stage method:

 S bbb.eee 5         6         7         8         9        10        11        12        13        14        15        16        17 23.046    25.056    27.067    29.079    31.092    33.106    35.121    37.137    39.154    41.172    43.191    45.211    47.232

-To use the 13-stage 10th-order method mentionned in the paragraph above,
store the same 116 coefficients into R39 to R154.

-With the same system of 2 differential equations, h = 0.1 & N = 10 ( don't forget to store 13 into R08 ), it yields:

XEQ "ERKN2"   >>>>     x   =  1                                                                                            ---Execution time = 14m58s---
RDN   y(1) = 1.531356645     and   R06 = y'(1) = -2.312840138
RDN   z(1) = 2.620254282              R07 = z'(1) =  2.941748401

-Albrecht's method already produced an almost perfect accuracy with this example,
so the superiority of the 10th-order method is not very apparent !
-In most cases, however, it's really worthwhile to use high-order formulae...

Notes:

-These programs employs formulae without error-estimate, so you'll have to calculate again
the solutions with a smaller stepsize - say  h/2 - to get an idea of the precision.
-But in references [2] & [3], you will find the coefficients of pairs of embedded Runge-Kutta-Nystrom formulas
of various orders to control the stepsize more directly.

References:

[1]  Abramowitz and Stegun - "Handbook of Mathematical Functions" -  Dover Publications -  ISBN  0-486-61272-4
[2]  Erwin Fehlberg - "Classical eighth- and lower-order Runge-Kutta-Nystrom formulas with stepsize control
for special second-order differential equations" - NASA TR R-381 & TR R-410
[3]  Philip Sharp - http://www.math.auckland.ac.nz/~sharp/rkn1.dat