# hp41programs

Runge-Kutta Explicit Runge-Kutta methods for the HP-41

Overview

1°)  The "Classical" 4th-Order Runge-Kutta Method

a)  First order differential equations
b)  System of 2 first-order differential equations
c)  System of 3 first-order differential equations

2°)  A Runge-Kutta Method with order 6

a)  First order differential equations
b)  System of 2 first-order differential equations

3°)  A Runge-Kutta Method with order 8

a)  First order differential equations
b)  System of 2 first-order differential equations

4°)  Explicit Runge-Kutta methods without error-estimate

a)  First order differential equations
b)  System of 2 first-order differential equations

5°)  A Runge-Kutta-Fehlberg Method of order 4 , embedded within 5th order

a)  First order differential equations
b)  System of 2 first-order differential equations

6°)  Explicit Runge-Kutta methods with built-in estimates

a)  First order differential equations
b)  System of 2 first-order differential equations

7°)  Second-Order Equations without dy/dx

-The coefficients of an "optimal" 4th-order method have been added in paragraph 4°) a)
-For more than 2 or 3 equations, cf "Systems of Differential Equations for the HP-41"

1°) The classical 4th order Runge-Kutta method

a) First order differential equations:  dy/dx = f (x;y)

-We solve    dy/dx = f(x,y)   with the initial value:   y(x0) = y0

RK4 uses the following formulae:

k1 = h.f(x;y)
k2 = h.f(x+h/2;y+k1/2)
k3 = h.f(x+h/2;y+k2/2)
k4 = h.f(x+h;y+k3)

and then,   y(x+h) = y(x) + ( k1 + 2.k2 + 2.k3 + k4 ) / 6      This formula duplicates the Taylor series up to the term in h4

Data Registers:          •   R00:  f name                                          ( Registers R00 to R04 are to be initialized before executing "RK4" )

•   R01 = x0      •   R03 = h/2 = half of the stepsize          R05 & R06: temp            ( storing h/2 instead of h saves several bytes )
•   R02 = y0      •   R04 = N = number of steps

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

>>>>>       In other words, this subroutine has to change the stack from

Y = y
X = x        into       X' = f(x;y)

 01  LBL "RK4" 02  RCL 04 03  STO 06 04  LBL 01 05  RCL 02 06  RCL 01 07  XEQ IND 00 08  RCL 03 09  ST+ 01 10  * 11  STO 05 12  RCL 02 13  + 14  RCL 01 15  XEQ IND 00 16  RCL 03 17  * 18  ST+ 05 19  ST+ 05 20  RCL 02 21  + 22  RCL 01 23  XEQ IND 00 24  RCL 03 25  ST+ 01 26  ST+ X 27  * 28  ST+ 05 29  RCL 02 30  + 31  RCL 01 32  XEQ IND 00 33  RCL 03         34  * 35  RCL 05 36  + 37  3 38  / 39  ST+ 02 40  DSE 06 41  GTO 01  42  RCL 02         43  RCL 01 44  END

( 65 bytes / SIZE 007 )

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

Example:     y(0) = 1 and  dy/dx = -2 x.y     Evaluate  y(1)              ( The solution is y = exp ( -x2 ) )

01  LBL "T"
02  *
03  ST+ X
04  CHS
05  RTN

"T"  ASTO 00
0     STO 01
1     STO 02    and if we choose  h = 0.1
0.05  STO 03
10     STO 04    XEQ "RK4"   yields          x = 1                                 ( in 20s )
X<>Y                    y(1) = 0.367881                   ( the exact result is 1/e =  0.367879441 )

-Press R/S to continue with the same h and N.
-Start over with a smaller stepsize to obtain an error estimate ( when h is divided by 2 , errors are approximately divided by 24 = 16 )

b) Systems of 2 first order differential equations:  dy/dx = f (x;y;z)  ;   dz/dx = g(x;y;z)

-The following program solves   dy/dx = f(x,y,z)  ,   dz/dx = g(x,y,z)     with the initial values:    y(x0) = y0  ,   z(x0) = z0

-The preceding formulae are now generalized to a system of 2 differential equations:

k1 = h.f(x;y;z)                                             l1 = h.g(x;y;z)
k2 = h.f(x+h/2;y+k1/2;z+l1/2)                      l2 = h.g(x+h/2;y+k1/2;z+l1/2)
k3 = h.f(x+h/2;y+k2/2;z+l2/2)                      l3 = h.g(x+h/2;y+k2/2;z+l2/2)
k4 = h.f(x+h;y+k3;z+l3)                               l4 = h.g(x+h;y+k3;z+l3)

then,   y(x+h) = y(x) + ( k1 + 2.k2 + 2.k3 + k4 ) / 6
and   z(x+h) = z(x) + ( l1 + 2.l2 + 2.l3 + l4 ) / 6             duplicate the Taylor series through the terms in h4

Data Registers:     ( Registers R00 to R05  are to be initialized before executing "RK4B" )

•   R00:  f name   •   R01 = x0      •   R04 = h/2 = half of the stepsize              R06 to R08: temp            ( storing h/2 instead of h saves several bytes )
•   R02 = y0      •   R05 = N = number of steps
•   R03 = z0
Flags: /
Subroutine:  a program calculating 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.

>>>>>    In other words, this subroutine has to change the stack from

Z = z
Y = y       to     Y' = f(x;y;z)
X = x               X' = g(x;y;z)

 01  LBL "RK4B" 02  RCL 05 03  STO 08 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 07 14  RCL 03 15  + 16  X<>Y 17  STO 06 18  RCL 02 19  + 20  RCL 01 21  XEQ IND 00 22  RCL 04 23  ST* Z 24  * 25  ST+ 07 26  ST+ 07 27  RCL 03 28  + 29  X<>Y 30  ST+ 06 31  ST+ 06 32  RCL 02 33  + 34  RCL 01 35  XEQ IND 00 36  RCL 04 37  ST+ 01 38  ST+ X 39  ST* Z 40  * 41  ST+ 07 42  RCL 03 43  + 44  X<>Y 45  ST+ 06 46  RCL 02 47  + 48  RCL 01 49  XEQ IND 00 50  RCL 04 51  ST* Z 52  * 53  RCL 07 54  + 55  3 56  / 57  ST+ 03 58  X<>Y 59  RCL 06  60  + 61  3 62  / 63  ST+ 02 64  DSE 08 65  GTO 01 66  RCL 03         67  RCL 02 68  RCL 01 69  END

( 99 bytes / SIZE 009 )

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

Example:     y(0) = 1 ; y'(0) = 0 ;  y'' + 2.x.y' + 2.y = 0   Evaluate y(1) and y'(1)            (  The solution is y = exp ( -x2 ) ;  y' = -2.x exp ( -x2 )  )

-This second order equation is equivalent to the system:             y(0) = 1 ;  z(0) = 0  ;  dy/dx = z , dz/dx = -2 x.z -2 y    where z = y'

01  LBL "U"
02  RCL Z
03  *
04  +
05  ST+ X
06  CHS
07  RTN

"U"  ASTO 00
0      STO 01
1      STO 02
0      STO 03   and if we choose  h = 0.1
0.05  STO 04
10     STO 05    XEQ "RK4B"   yields       x =  1                                        ( in 32s )
RDN                      y(1) =  0.367881                            the exact result is  y(1) = 1/e =  0.367879441
RDN                      z(1) = -0.735762                            the exact result is  z(1) = -2/e = -0.735758883

-Press R/S to continue with the same h and N.

N.B:

-To obtain an error estimate, start over again with a smaller stepsize:
when h is divided by 2 , errors are approximately divided by 16

c) Systems of 3 differential equations: dy/dx = f (x,y,z,u) , dz/dx = g(x,y,z,u) , du/dx = h(x,y,z,u)

-"RK4C" solves   dy/dx = f(x,y,z,u)  ,   dz/dx = g(x,y,z,u)  ,  du/dx = h(x,y,z,u)   with the initial values:    y(x0) = y0  ,   z(x0) = z0  ,  u(x0) = u0

Data Registers:     ( Registers R00 to R06  are to be initialized before executing "RK4C" )

•   R00:  f name   •   R01 = x0      •   R05 = h/2 = half of the stepsize              R07 to R10: temp
•   R02 = y0      •   R06 = N = number of steps
•   R03 = z0
•   R04 = u0
Flags: /
Subroutine:  a program calculating 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.

>>>>>    In other words, this subroutine has to change the stack from

T = u
Z = z                Z' = f(x;y;z;u)
Y = y       to     Y' = g(x;y;z;u)
X = x               X' = h(x;y;z;u)

-Line 89 is a three-byte GTO 01

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

( 133 bytes / SIZE 011 )

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

Example:

dy/dx = -y.z.u                    y(0) = 1         Evaluate  y(1) , z(1) , u(1)  with  h = 0.1   ( whence N = 10 )
dz/dx = x.( y + z - u )         z(0) = 1
du/dx = x.y - z.u                u(0) = 2

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
2     STO 04
0.05  STO 05
10    STO 06    XEQ "RK4C"   >>>>      x   =  1                         ( in 53 seconds )
RDN     y(1) =  0.258209
RDN     z(1) =  1.157620
RDN     u(1) =  0.842179

-Press R/S to continue with the same h and N.

2°) A Runge-Kutta method with order 6

a) First order differential equations:  dy/dx = f (x;y)

-The formula duplicates the Taylor series up to the term in h6
-Methods of this order require 7 stages ( 7 evaluations of the function f within each step)

Data Registers:             ( Registers R00 to R04 are to be initialized before executing "RK6" )

•   R00:  f name   •   R01 = x0      •   R03 =  h = the stepsize                   R05: counter
•   R02 = y0      •   R04 = N = the number of steps      R06 , ...... , R11 = k1 , ..... , k6

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

-Line 164 is a three-byte GTO 01

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

( 219 bytes / SIZE 012 )

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

Example:     y(0) = 1 and  dy/dx = -2 x.y     Evaluate  y(1)

01  LBL "T"
02  *
03  ST+ X
04  CHS
05  RTN

"T"  ASTO 00
0     STO 01
1     STO 02    and if we choose  h = 0.1
0.1    STO 03
10     STO 04    XEQ "RK6"   yields          x = 1                       ( in 84s )
X<>Y                    y(1) = 0.367879436    ( error =  -5 10-9 )

-Press R/S to continue with the same h and N.
-Start over with a smaller stepsize to obtain an error estimate ( when h is divided by 2 , errors are approximately divided by 26 = 64 )

b) Systems of 2 first order differential equations:  dy/dx = f (x;y;z)  ;   dz/dx = g(x;y;z)

Data Registers:                       ( Registers R00 to R05  are to be initialized before executing "RK6B" )

•   R00:  f name   •   R01 = x0      •   R04 = h = stepsize                         R06 to R16: temp
•   R02 = y0      •   R05 = N = number of steps
•   R03 = z0
Flags: /
Subroutine:  a program calculating 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 277 is a three-byte GTO 01

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

( 378 bytes / SIZE 017 )

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

Example:          y(0) = 1 ;  z(0) = 0  ;  dy/dx = z , dz/dx = -2 x.z -2 y    ;    calculate  y(1) and z(1)

01  LBL "U"
02  RCL Z
03  *
04  +
05  ST+ X
06  CHS
07  RTN

"U"  ASTO 00
0      STO 01
1      STO 02
0      STO 03    and if we choose  h = 0.1
0.1    STO 04
10     STO 05    XEQ "RK6B"   yields      x =  1                                ( in 2mn30s )
RDN                      y(1) =  0.367879433             ( y-error = -9 10-9 )
RDN                      z(1) = -0.735758865             ( z-error = 17 10-9 )

-Press R/S to continue with the same h and N.
-Start over with a smaller stepsize to obtain an error estimate ( when h is divided by 2 , errors are approximately divided by 26 = 64 )

3°) A Runge-Kutta method with order 8

a) First order differential equations:  dy/dx = f (x;y)

-Methods of this order require 11 stages. The following one was discovered in 1972 by Cooper and Verner.
-The formula duplicates the Taylor series up to the term in h8
-The following program uses 41 coefficients which are of the form ( a + b 211/2 )/c  where a , b , c are integers.

-They are to be stored in registers R11 thru R51  ( all the digits are correct )

Data Registers:                     ( Registers R00 to R04  and R11 to R51 are to be initialized before executing "RK8" )

•   R00:  f name   •   R01 = x0      •   R03 =  h = the stepsize                      R05  to R10 and R52: temp
•   R02 = y0      •   R04 = N = the number of steps      •  R11 to R51:  the 41 coefficients listed below

R11 = 0.8273268354                   R21 = 0.5766714727                   R31 = 0.1289862930                           R41 = 2.031083139
R12 = 0.5                                     R22 = 0.1855068535                   R32 = -0.01186868389                       R42 = -0.6379313502
R13 = 0.1726731646                   R23 = 0.3865246267                   R33 = 0.002002165993                      R43 = 0.9401094520
R14 = 9                                        R24 = -0.4634553896                 R34 = 0.9576053440                           R44 = -2.229158210
R15= 14                                       R25 = 0.3772937693                   R35 = -0.6325461607                         R45 = 7.553840442
R16 = 0.25                                   R26 = 0.1996369936                   R36 = 0.1527777778                          R46 = -7.164951553
R17 = 0.8961811934                   R27 = 0.09790045616                 R37 = -0.009086961101                     R47 = 2.451380432
R18 = -0.2117115009                  R28 = 0.3285172131                  R38 = 0.03125                                     R48 = -0.5512205631
R19 = 7                                        R29 = -0.3497052863                 R39 = 1.062498447                             R49 = 0.05
R20 = 0.06514850915                 R30 = -0.03302551131               R40 = -1.810863083                           R50 = 0.2722222222
R51 = 2.8125

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

-Line 251 is a three-byte GTO 01

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

( 329 bytes / SIZE 053 )

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

Example:     y(0) = 1 and  dy/dx = -2 x.y     Evaluate  y(1)

-Initialize registers R11 thru R51.

01  LBL "T"
02  *
03  ST+ X
04  CHS
05  RTN

"T" ASTO 00
0    STO 01
1    STO 02    and if we choose  h = 0.1
0.1  STO 03
10   STO 04    XEQ "RK8"   yields            x = 1                        ( in 2mn )
X<>Y                     y(1) = 0.3678794412   correct to 10 D!

-Press R/S to continue with the same h and N.
-Start over with a smaller stepsize to obtain an error estimate ( when h is divided by 2 , errors are approximately divided by 28 = 256 )

b) Systems of 2 first order differential equations:  dy/dx = f (x;y;z)  ;   dz/dx = g(x;y;z)

Data Registers:                      ( Registers R00 to R05 and R18 to R58 are to be initialized before executing "RK8B" )

•   R00:  f name   •   R01 = x0      •   R04 = h = stepsize                         R06 to R17: temp    R59: counter
•   R02 = y0      •   R05 = N = number of steps       •   R18 thru  R58 = the same 41 numbers used in the previous program.
•   R03 = z0
Flags: /
Subroutine:  a program calculating 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 444 is a three-byte GTO 01

 01  LBL "RK8B"   02  RCL 05   03  STO 59   04  GTO 01   05  LBL 00   06  XEQ IND 00   07  RCL 04   08  ST* Z   09  *   10  RTN   11  LBL 01   12  RCL 03   13  RCL 02   14  RCL 01   15  XEQ 00   16  STO 12   17  RCL 19   18  *   19  RCL 03   20  +   21  X<>Y   22  STO 06   23  RCL 19   24  *   25  RCL 02   26  +   27  RCL 04   28  RCL 19   29  *   30  RCL 01   31  +   32  XEQ 00   33  STO 13   34  RCL 12   35  +   36  RCL 23   37  *   38  RCL 03   39  +   40  X<>Y   41  STO 07   42  RCL 06   43  +   44  RCL 23   45  *   46  RCL 02   47  +   48  RCL 04   49  RCL 19   50  *   51  RCL 01   52  +   53  XEQ 00   54  STO 14   55  RCL 24   56  *   57  RCL 13   58  RCL 25   59  *   60  +   61  RCL 12   62  RCL 26   63  /   64  +   65  RCL 03   66  +   67  X<>Y   68  STO 08   69  RCL 24   70  *   71  RCL 07   72  RCL 25   73  *   74  +   75  RCL 06 76  RCL 26   77  /   78  +   79  RCL 02   80  +   81  RCL 04   82  RCL 18   83  *   84  RCL 01   85  +   86  XEQ 00    87  STO 15           88  RCL 27   89  *   90  RCL 14   91  RCL 28   92  *   93  +   94  RCL 12   95  RCL 29   96  *   97  +   98  RCL 03   99  + 100  X<>Y 101  STO 09 102  RCL 27 103  * 104  RCL 08 105  RCL 28 106  * 107  + 108  RCL 06 109  RCL 29 110  * 111  + 112  RCL 02 113  + 114  RCL 04 115  RCL 18 116  * 117  RCL 01 118  + 119  XEQ 00 120  STO 16 121  RCL 30 122  * 123  RCL 15 124  RCL 31 125  * 126  + 127  RCL 14 128  RCL 32 129  * 130  + 131  RCL 12 132  RCL 33 133  * 134  + 135  RCL 03 136  + 137  X<>Y 138  STO 10 139  RCL 30 140  * 141  RCL 09 142  RCL 31 143  * 144  + 145  RCL 08 146  RCL 32 147  * 148  + 149  RCL 06 150  RCL 33 151  * 152  + 153  RCL 02 154  + 155  RCL 04 156  RCL 19 157  * 158  RCL 01 159  + 160  XEQ 00  161  STO 17         162  RCL 34 163  * 164  RCL 16 165  RCL 35 166  * 167  + 168  RCL 15 169  RCL 36 170  * 171  + 172  RCL 14 173  RCL 37 174  * 175  + 176  RCL 12 177  RCL 38 178  * 179  + 180  RCL 03 181  + 182  X<>Y 183  STO 11 184  RCL 34 185  * 186  RCL 10 187  RCL 35 188  * 189  + 190  RCL 09 191  RCL 36 192  * 193  + 194  RCL 08 195  RCL 37 196  * 197  + 198  RCL 06 199  RCL 38 200  * 201  + 202  RCL 02 203  + 204  RCL 04 205  RCL 20 206  * 207  RCL 01 208  + 209  XEQ 00 210  STO 13 211  RCL 21 212  / 213  RCL 17 214  RCL 39 215  * 216  + 217  RCL 16 218  RCL 40 219  * 220  + 221  RCL 12 222  RCL 22 223  / 224  + 225  RCL 03 226  + 227  X<>Y 228  STO 07 229  RCL 21 230  / 231  RCL 11 232  RCL 39 233  * 234  + 235  RCL 10 236  RCL 40         237  * 238  + 239  RCL 06 240  RCL 22 241  / 242  + 243  RCL 02 244  + 245  RCL 04 246  RCL 20 247  * 248  RCL 01 249  + 250  XEQ 00 251  STO 14 252  RCL 41 253  * 254  RCL 13 255  RCL 42 256  * 257  + 258  RCL 17 259  RCL 43 260  * 261  + 262  RCL 16 263  RCL 44 264  * 265  + 266  RCL 12 267  RCL 45 268  * 269  + 270  RCL 03 271  + 272  X<>Y 273  STO 08 274  RCL 41 275  * 276  RCL 07 277  RCL 42 278  * 279  + 280  RCL 11 281  RCL 43 282  * 283  + 284  RCL 10 285  RCL 44 286  * 287  + 288  RCL 06 289  RCL 45 290  * 291  + 292  RCL 02 293  + 294  RCL 04 295  RCL 19 296  * 297  RCL 01 298  + 299  XEQ 00 300  STO 15 301  RCL 46 302  * 303  RCL 14 304  RCL 47 305  * 306  + 307  RCL 13 308  RCL 48  309  * 310  + 311  RCL 17         312  RCL 49 313  * 314  + 315  RCL 16 316  RCL 21 317  / 318  + 319  RCL 12 320  RCL 22 321  / 322  + 323  RCL 03 324  + 325  X<>Y 326  STO 09 327  RCL 46 328  * 329  RCL 08 330  RCL 47 331  * 332  + 333  RCL 07 334  RCL 48 335  * 336  + 337  RCL 11 338  RCL 49 339  * 340  + 341  RCL 10 342  RCL 21 343  / 344  + 345  RCL 06 346  RCL 22 347  / 348  + 349  RCL 02 350  + 351  RCL 04 352  RCL 18 353  * 354  RCL 01 355  + 356  XEQ 00 357  X<> 16 358  RCL 55 359  * 360  RCL 15 361  RCL 51 362  * 363  + 364  RCL 14 365  RCL 52 366  * 367  + 368  RCL 13 369  RCL 53 370  * 371  + 372  RCL 17 373  RCL 54 374  * 375  + 376  RCL 16 377  RCL 50 378  * 379  + 380  RCL 03 381  + 382  X<>Y 383  X<> 10  384  RCL 55         385  * 386  RCL 09 387  RCL 51 388  * 389  + 390  RCL 08 391  RCL 52 392  * 393  + 394  RCL 07 395  RCL 53 396  * 397  + 398  RCL 11 399  RCL 54 400  * 401  + 402  RCL 10 403  RCL 50 404  * 405  + 406  RCL 02 407  + 408  RCL 04 409  RCL 01 410  + 411  STO 01 412  XEQ 00 413  RCL 12 414  + 415  RCL 56 416  * 417  RCL 16 418  RCL 14 419  + 420  RCL 57 421  * 422  + 423  RCL 15 424  RCL 58 425  ST/ 09 426  / 427  + 428  ST+ 03 429  X<>Y 430  RCL 06 431  + 432  RCL 56 433  * 434  RCL 10 435  RCL 08 436  + 437  RCL 57 438  * 439  + 440  RCL 09 441  + 442  ST+ 02 443  DSE 59 444  GTO 01 445  RCL 03 446  RCL 02 447  RCL 01 448  END

( 593 bytes / SIZE 060 )

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

Example:     y(0) = 1 ;  z(0) = 0  ;  dy/dx = z , dz/dx = -2 x.z -2 y   ;    calculate  y(1) and z(1)

-Initialize registers R18 thru R58.

01  LBL "U"
02  RCL Z
03  *
04  +
05  ST+ X
06  CHS
07  RTN

"U" ASTO 00
0    STO 01
1    STO 02
0    STO 03   and if we choose  h = 0.1
0.1  STO 04
10   STO 05    XEQ "RK8B"   yields        x = 1                         in 3mn10s
RDN                      y(1) =  0.3678794412   correct to 10D.
RDN                      z(1) = -0.7357588824   z-error = -10-10

-Press R/S to continue with the same h and N.
-Start over with a smaller stepsize to obtain an error estimate ( when h is divided by 2 , errors are approximately divided by 28 = 256 )

4°) Explicit Runge-Kutta Methods without error-estimate

a) First order differential equations:  dy/dx = f (x;y)

-All n stage explicit Runge-Kutta formulae without error estimate are determined by ( n2 + 3.n - 2 ) / 2 coefficients ai,j ; bi ; ci

k1 = h.f(x;y)
k2 = h.f(x+b2.h;y+a2,1.k1)
k3 = h.f(x+b3.h;y+a3,1.k1+a3,2.k2)
.............................................................                                         (  actually,     ai,1 + ai,2 + .......... + ai,i-1 = bi  )

kn = h.f(x+bn.h;y+an,1.k1+an,2.k2+ .......... +an,n-1.kn-1)

then,     y(x+h) = y(x) + c1.k1+ ................ + cn.kn

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

Data Registers:             ( Registers R00 to R05 and    Rn+10 to R(n2+5.n+16)/2     are to be initialized before executing "ERK" )

•   R00:  f name   •   R01 = x0                                      R06 to R09:  temp             •  Rn+10 to R(n2+5.n+16)/2 =  the ( n2 + 3.n - 2 ) / 2 coefficients ai,j ; bi ; ci
•   R02 = y0                                      R10 = k1                                                                            which are to be stored as shown below:
•   R03 = h = stepsize                        R11 = k2
•   R04 = N = number of steps          ............
•   R05 = n  = number of stages        Rn+9 = kn

•  Rn+10 = a2,1                                                                                •  Rn+11 = b2
•  Rn+12 = a3,1        •  Rn+13 = a3,2                                                •  Rn+14 = b3

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

•  R(n2+n+18)/2 = an,1  ................   •  R(n2+3n+14)/2 = an,n-1       •  R(n2+3n+16)/2 = bn

•  R(n2+3n+18)/2 = c1 ............................................................       •  R(n2+5n+16)/2 = cn

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

 01  LBL "ERK" 02  RCL 04 03  STO 06 04  LBL 01 05  RCL 05  06  10.9 07  STO 07         08  + 09  STO 08 10  RCL 05 11  9 12  + 13  .1 14  % 15  11 16  + 17  STO 09 18  RCL 02 19  RCL 01 20  XEQ IND 00 21  RCL 03 22  * 23  STO 10 24  LBL 02 25  RCL 07 26  INT 27  .1 28  % 29  LASTX 30  1/X 31  + 32  STO 07 33  CLX 34  LBL 03 35  RCL IND 07 36  RCL IND 08 37  * 38  + 39  ISG 08 40  ISG 07 41  GTO 03 42  RCL 02 43  + 44  RCL 03 45  RCL IND 08 46  * 47  RCL 01 48  + 49  XEQ IND 00 50  RCL 03 51  * 52  STO IND 09 53  ISG 08 54  ISG 09 55  GTO 02 56  RCL 05 57  ST- 09 58  CLX 59  LBL 04 60  RCL IND 08 61  RCL IND 09  62  * 63  + 64  ISG 08 65  ISG 09 66  GTO 04 67  ST+ 02 68  RCL 03 69  ST+ 01 70  DSE 06 71  GTO 01 72  RCL 02 73  RCL 01 74  END

( 110 bytes / SIZE ( n2+5.n+18 )/2 )

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

-For example, if you are using the classical four stage method:      4   STO 05   and the 13 coefficients:

1/2                 1/2                                             R14                   R15
0    1/2          1/2       are to be stored into        R16  R17          R18        respectively
0     0     1      1                                               R19  R20  R21  R22
1/6  1/3  1/3  1/6                                              R23  R24  R25  R26

------------------------------------------------------------------------------------------------------------------------------------------------------------

-Here are the coefficients of an "optimal" 4th-order 4-stage method:    4   STO 05

R14 =  0.3716151060                                                                                            R15 = 0.3716151060
R16 = -0.1180444797       R17 =  0.7180444797                                                  R18 = 0.6
R19 =  0.5173871366        R20 = -0.5608902997        R21 = 1.043503163         R22 = 1
R23 =  0.1474734369        R24 =  0.3125088197         R25 = 0.3903768538       R26 = 0.1496408895

Example:     y(0) = 1 and  dy/dx = -2 x.y     Evaluate  y(1)

-Initialize registers R14 thru R26.

01  LBL "T"
02  *
03  ST+ X
04  CHS
05  RTN

"T" ASTO 00
0    STO 01
1    STO 02    and with  h = 0.1
0.1  STO 03
10   STO 04
4    STO 05   XEQ "ERK"   yields            x = 1                        ( in 79s )
X<>Y                     y(1) = 0.367879270     the error is only  -1.7 E-7   ( the error was 1.6 E-6  with the "classical" Runge-Kutta method )

-The formula minimizes a bound on the truncation error term.
-However, this "optimal" method is not always so good, and there are examples where it is less accurate than the "classical" RK4

-------------------------------------------------------------------------------------------------------------------------------------------------------------------

-With the 7 stage method used in "RK6":    7  STO 05    and the 34 coefficients:

1/3                                                                        1/3                                                R17                                            R18
0           2/3                                                          2/3                                                R19  R20                                   R21
1/12       1/3    -1/12                                              1/3                                               R22  R23  R24                           R25
25/48   -55/24  35/48       15/8                                5/6      are to be stored into          R26  R27  R28  R29                   R30
3/20   -11/24   -1/8          1/2       1/10                   1/6                                              R31  R32  R33  R34  R35          R36
-261/260  33/13  43/156  -118/39  32/195   80/39      1                                                R37  R38  R39  R40  R41  R42  R43
13/200     0       11/40      11/40     4/25      4/25    13/200                                          R44  R45  R46  R47  R48  R49  R50

................ and so on ..........

-For the other instructions, cf "RK6"
-With the same differential equation, execution time = 2mn50s for a 7 stage method ( instead of 89s )
( the number of bytes is decreased, but execution time is increased ... )

b) Systems of 2 first order differential equations:  dy/dx = f (x;y;z)  ;   dz/dx = g(x;y;z)

Data Registers:                    ( Registers R00 to R06 and    R2n+12 to R(n2+7.n+20)/2     are to be initialized before executing "ERKB" )

•   R00:  f name   •   R01 = x0                                      R07 to R11:  temp             •  R2n+12 to R(n2+7n+20)/2 =  the ( n2 + 3.n - 2 ) / 2 coefficients ai,j ; bi ; ci
•   R02 = y0                                      R12 to Rn+11 = kiy                                                            which are to be stored as shown below:
•   R03 = z0                                      Rn+12 to R2n+11 = kiz
•   R04 = h   = stepsize
•   R05 = N  = number of steps
•   R06 = n   = number of stages

•  R2n+12 = a2,1                                                                                •  R2n+13 = b2
•  R2n+14 = a3,1        •  R2n+15 = a3,2                                              •  R2n+16 = b3

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

•  R(n2+3n+22)/2 = an,1  ................   •  R(n2+5n+18)/2 = an,n-1       •  R(n2+5n+20)/2 = bn

•  R(n2+5n+22)/2 = c1 ............................................................         •  R(n2+7n+20)/2 = cn

Flags: /
Subroutine:  a program calculating 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.

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

( 141 bytes / SIZE ( n2+ 7.n + 22 )/2 )

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

-For example, if you are using the classical four stage method:      4   STO 06   and the 13 coefficients:

1/2                 1/2                                             R20                   R21
0    1/2          1/2       are to be stored into        R22  R23          R24        respectively
0     0     1      1                                               R25  R26  R27  R28
1/6  1/3  1/3  1/6                                              R29  R30  R31  R32

-With the 7-stage method used in "RK6B":         7  STO 06    and the 34 coefficients:

1/3                                                                        1/3                                                R26                                            R27
0           2/3                                                          2/3                                                R28  R29                                   R30
1/12       1/3    -1/12                                              1/3                                               R31  R32  R33                           R34
25/48   -55/24  35/48       15/8                                5/6      are to be stored into          R35  R36  R37  R38                   R39
3/20   -11/24   -1/8          1/2       1/10                   1/6                                              R40  R41  R42  R43  R44          R45
-261/260  33/13  43/156  -118/39  32/195   80/39      1                                                R46  R47  R48  R49  R50  R51  R52
13/200     0       11/40      11/40     4/25      4/25    13/200                                          R53  R54  R55  R56  R57  R58  R59

................ and so on ..........

-For the other instructions, cf "RK6B"
-With the same differential system, execution time = 4mn36s  for a 7-stage method ( instead of 2mn44s )

5°) A Runge-Kutta-Fehlberg Method of order 4 , embedded within 5th order

-In this Runge-Kutta-Fehlberg method, a 4th-order formula is used to compute the solution
and the difference between this 4th-order formula and a 5th-order formula provides an error estimate.

a) First order differential equations:  dy/dx = f (x;y)

Data Registers:                      ( Registers R00 to R04 are to be initialized before executing "RKF" )

•   R00:  f name   •   R01 = x0      •   R03 = h = stepsize                          R05 = error estimate
•   R02 = y0      •   R04 = N = number of steps             R06 thru R10: temp

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

-If you replace line 135 by   ABS   ST+ 05   you'll get a less optimistic error estimate ( in absolute value )
-Lines 150 & 154 are three-byte GTOs

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

( 204 bytes / SIZE 011 )

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

Example:     y(0) = 1 and  dy/dx = -2 x.y     Evaluate  y(1)

01  LBL "T"
02  *
03  ST+ X
04  CHS
05  RTN

"T" ASTO 00
0    STO 01
1    STO 02    and if we choose  h = 0.1
0.1  STO 03
10   STO 04    XEQ "RKF"   yields          x = 1
X<>Y                    y(1) = 0.367879263
and the error estimate    RCL 05      >>>       -9.7 10-8           ( or  5.4 10-7  with   ABS  ST+ 05  at line 135  )
-The error is actually  -1.8 10-7

-Press R/S to continue with the same h and N.

-If you want to use the 5th-order formula to compute y(x) , replace line 135 by  ST+ 02  ST+ 05
( the error will be probably overestimated )
-In the example above, it yields  y(1) = 0.367879453   ( R05 = 9.7 10-8  whereas the error is only 1.2 10-8 )

b) Systems of 2 first order differential equations:  dy/dx = f (x;y;z)  ;   dz/dx = g(x;y;z)

Data Registers:                     ( Registers R00 to R05  are to be initialized before executing "RKFB" )

•   R00:  f name   •   R01 = x0      •   R04 = h = stepsize                         R06 = y-error estimate           R08 to R16: temp
•   R02 = y0      •   R05 = N = number of steps           R07 = z-error estimate
•   R03 = z0
Flags: /
Subroutine:  a program calculating 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.

-If you replace line 201 by   ABS   ST+ 07   you'll get a less optimistic error estimate ( in absolute value )
-If you replace line 215 by   ABS   ST+ 06   you'll get a less optimistic error estimate ( in absolute value )
-Lines 241 & 246 are three-byte GTOs

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

( 343 bytes / SIZE 017)

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

Example:     y(0) = 1 ;  z(0) = 0  ;  dy/dx = z , dz/dx = -2 x.z -2 y   ;    calculate  y(1) and z(1)

01  LBL "U"
02  RCL Z
03  *
04  +
05  ST+ X
06  CHS
07  RTN

"U" ASTO 00
0    STO 01
1    STO 02
0    STO 03    and if we choose  h = 0.1
0.1  STO 04
10   STO 05    XEQ "RKFB"   yields       x  =  1                      ( in 2mn16s )
RDN                      y(1) =  0.367879517
RDN                      z(1) = -0.735759034

and the error    RCL 06     >>>    -8.7 10-8    ( or  6.5 10-7  with   ABS  ST+ 05 at line 218 )      actually   y-error =    7.6 10-8
estimates:        RCL 07     >>>    -2.1 10-7    ( or    8  10-7  ----    ABS  ST+ 06  at line 204 )    actually    z-error =  -1.5 10-7

-Press R/S to continue with the same h and N.

-If you want to use the 5th-order formula to compute y(x) & z(x) , replace line 215 by  ST+ 02  ST+ 06  and line 201 by  ST+ 03  ST+ 07
( the errors will be probably overestimated )
-In the example above, it yields:

y(1) =  0.367879439   ( R06 = 8.7 10-8  whereas the error is only -2 10-9 )
z(1) = -0.735758876   ( R07 = 2.1 10-7  whereas the error is only  6 10-9 )

6°) Explicit Runge-Kutta Methods with built-in estimates

a) First order differential equations:  dy/dx = f (x;y)

-All n stage explicit Runge-Kutta formulae with error estimate are determined by ( n2 + 5.n - 2 ) / 2 coefficients ai,j ; bi ; ci ; di

k1 = h.f(x;y)
k2 = h.f(x+b2.h;y+a2,1.k1)
k3 = h.f(x+b3.h;y+a3,1.k1+a3,2.k2)
.............................................................                                         (  actually,     ai,1 + ai,2 + .......... + ai,i-1 = bi  )

kn = h.f(x+bn.h;y+an,1.k1+an,2.k2+ .......... +an,n-1.kn-1)

then,     y(x+h) = y(x) + c1.k1+ ................ + cn.kn    and the difference between the higher-order and the lower-order formulae gives an error estimate:

err-estim =  d1.k1+ ................ + dn.kn

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

Data Registers:             ( Registers R00 to R05 and    Rn+11 to R(n2+7.n+18)/2     are to be initialized before executing "RKE" )

•   R00:  f name   •   R01 = x0                                     R06 = y-error             •  Rn+11 to R(n2+7.n+18)/2 =  the ( n2 + 5.n - 2 ) / 2 coefficients ai,j ; bi ; ci ; di
•   R02 = y0                                      R07 to R10: temp                                                         which are to be stored as shown below:
•   R03 = h = stepsize                        R11 to Rn+10 = ki
•   R04 = N = number of steps
•   R05 = n  = number of stages

•  Rn+11 = a2,1                                                                                •  Rn+12 = b2
•  Rn+13 = a3,1        •  Rn+14 = a3,2                                                •  Rn+15 = b3

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

•  R(n2+n+20)/2 = an,1  ................   •  R(n2+3n+16)/2 = an,n-1       •  R(n2+3n+18)/2 = bn

•  R(n2+3n+20)/2 = c1 ............................................................       •  R(n2+5n+18)/2 = cn

•  R(n2+5n+20)/2 = d1 ............................................................       •  R(n2+7n+18)/2 = dn

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

-Replace line 67 by  ABS  ST+ 06   after replacing line 70 by  RCL IND 10  if you want a less optimistic error-estimate ( in magnitude )
-Line 85 is a three-byte GTO 10

 01  LBL "RKE"  02  CLX 03  STO 06 04  LBL 10 05  RCL 04 06  STO 07 07  LBL 01 08  RCL 05 09  11.9 10  STO 08 11  + 12  STO 09         13  RCL 05 14  10 15  + 16  .1 17  % 18  12 19  + 20  STO 10 21  RCL 02 22  RCL 01 23  XEQ IND 00 24  RCL 03 25  * 26  STO 11 27  LBL 02 28  RCL 08 29  INT 30   E3 31  / 32  11 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 02 45  + 46  RCL 03 47  RCL IND 09 48  * 49  RCL 01 50  + 51  XEQ IND 00 52  RCL 03 53  * 54  STO IND 10 55  ISG 09 56  ISG 10 57  GTO 02 58  RCL 05 59  ST- 10 60  RCL 09 61  + 62  0 63  LBL 04 64  RCL IND Y 65  RCL IND 10 66  * 67  ST- 06 68  CLX 69  RCL IND 09 70  LASTX 71  * 72  + 73  ISG Y 74  ISG 09 75  ISG 10 76  GTO 04 77  ST+ 02 78  RCL 03 79  ST+ 01 80  DSE 07 81  GTO 01  82  RCL 02         83  RCL 01 84  RTN 85  GTO 10 86  END

( 129 bytes / SIZE ( n2+ 7.n + 20 )/2 )

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

Example:     Here are the 51 coefficients of an 8-stage method with order 5 with a 6th-order error-estimating formula.
Since n = 8 , these numbers are to be stored from R19 to R69

R19 =  1/18                                                                                                                                                                                          R20 = 1/18
R21 = -1/12          R22 = 1/4                                                                                                                                                                R23 = 1/6
R24 = -2/81          R25 = 4/27        R26 =  8/81                                                                                                                                   R27 = 2/9
R28 =  40/33        R29 = -4/11       R30 = -56/11        R31 = 54/11                                                                                                     R32 = 2/3
R33 = -369/73      R34 = 72/73      R35 = 5380/219    R36 = -12285/584  R37 = 2695/1752                                                             R38 = 1
R39 = -8716/891  R40 = 656/297  R41 = 39520/891  R42 = -416/11        R43 = 52/27          R44 = 0                                              R45 = 8/9
R46 =  3015/256  R47 = -9/4         R48 = -4219/78     R49 = 5985/128     R50 = -539/384     R51 = 0            R52 = 693/3328        R53 = 1

R54 = 3/80           R55 = 0              R56 = 4/25            R57 = 243/1120     R58 = 77/160         R59 = 73/700    R60 = 0                    R61 = 0
R62 = 33/640       R63 = 0              R64 = -132/325     R65 = 891/2240     R66 = -33/320       R67 = -73/700  R68 = 891/8320        R69 = 2/35

With our usual equation:   y(0) = 1 and  dy/dx = -2 x.y     Evaluate  y(1)  with  h = 0.1  and  N = 10

01  LBL "T"
02  *
03  ST+ X
04  CHS
05  RTN

"T" ASTO 00
0    STO 01
1    STO 02
0.1  STO 03
10   STO 04
8    STO 05    XEQ "RKE"   yields          x = 1                                    ( in 3mn56s )
X<>Y                    y(1) = 0.367879457
and the error estimate                         R06 = -13 10-9
-The error is actually  16 10-9

-Press R/S to continue ( after changing  h and N in registers R03 and R04 if needed ).

b) Systems of 2 first order differential equations:  dy/dx = f (x;y;z)  ;   dz/dx = g(x;y;z)

Data Registers:             ( Registers R00 to R06 and    R2n+14 to R(n2+9.n+24)/2     are to be initialized before executing "RKEB" )

•   R00:  f name   •   R01 = x0                                   R07 = y-error             •  R2n+14 to R(n2+9.n+24)/2 =  the ( n2 + 5.n - 2 ) / 2 coefficients ai,j ; bi ; ci ; di
•   R02 = y0                                   R08 = z-error                                                                which are to be stored as shown below:
•   R03 = y0
•   R04 = h =  stepsize                        R09  to  R13: temp
•   R05 = N = number of steps           R14    to Rn+13 = ki(y)
•   R06 = n  = number of stages         R14+n to R2n+13 = ki(z)

•  R2n+14 = a2,1                                                                                •  R2n+15 = b2
•  R2n+16 = a3,1        •  R2n+17 = a3,2                                              •  R2n+18 = b3

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

•  R(n2+3n+26)/2 = an,1  ................   •  R(n2+5n+22)/2 = an,n-1       •  R(n2+5n+24)/2 = bn

•  R(n2+5n+26)/2 = c1 ............................................................         •  R(n2+7n+24)/2 = cn

•  R(n2+7n+26)/2 = d1 ............................................................         •  R(n2+9n+24)/2 = dn

Flags: /
Subroutine:    a program calculating 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.

-Replace line 79  by  ABS  ST+ 07  and line 77 by  ABS  ST+ 08   if you want a less optimistic error-estimate ( in magnitude )
-Lines 83 & 88 are three-byte  GTOs

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

( 164 bytes / SIZE ( n2+ 9.n + 26 )/2 )

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

Example:     y(0) = 1 ;  z(0) = 0  ;  dy/dx = z , dz/dx = -2 x.z -2 y   ;    calculate  y(1) and z(1)  with  h = 0.1  and  N = 10

-If you want to use the same 8-stage 5th-order method embedded within 6th-order, store the same 51 coefficients into  R30 thru R80

01  LBL "U"
02  RCL Z
03  *
04  +
05  ST+ X
06  CHS
07  RTN

"U" ASTO 00
0    STO 01
1    STO 02
0    STO 03
0.1  STO 04
10   STO 05
8    STO 06   XEQ "RKEB"   yields       x  =  1                      ( in 6mn34s )
RDN                      y(1) =  0.367879378
RDN                      z(1) = -0.735758757

and the error estimates    RCL 07  =    -85 10-9       actually    y-error =   -63 10-9
RCL 08  =    153 10-9       actually    z-error =   126 10-9

-Press R/S to continue ( after changing  h and N in registers R04 and R05 if needed ).

7°)  Second-order Equations without dy/dx

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

 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

Notes:

-Simply press R/S to continue with the same h &N
-For more than 1 equation, cf "Systems of Differential Equations for the HP-41"

Reference:

[1]  J.C. Butcher  - "The Numerical Analysis of Ordinary Differential Equations" - John Wiley & Sons  -   ISBN  0-471-91046-5