The Museum of HP Calculators

Explicit Runge-Kutta methods for the HP-41

This program is by Jean-Marc Baillard and is used here by permission.

This program is supplied without representation or warranty of any kind. Jean-Marc Baillard and The Museum of HP Calculators therefore assume no responsibility and shall have no liability, consequential or otherwise, of any kind arising from the use of this program material or any part thereof.

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

 2°)  A Runge-Kutta Method with order 6

   a)  First order differential equations ( improved )
   b)  System of 2 first order differential equations   ( improved )

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

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

   a)  First order differential equations  ( improved )
   b)  System of 2 first order differential equations   ( improved )

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

   a)  First order differential equations    ( new )
   b)  System of 2 first order differential equations   ( new )
 
 

-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:     ( Registers R00 to R04 are to be initialized before executing "RK4" )

    R00:  f name      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.

                                                                                                                  Z = z
      >>>>>    In other words, this subroutine has to change the stack from   Y = y       into     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                            -----------------  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.

                                                                                                                  T = u
                                                                                                                  Z = z                   Z' = f(x;y;z;u)
      >>>>>    In other words, this subroutine has to change the stack from   Y = y       into     Y' = g(x;y;z;u)
                                                                                                                  X = x                  X' = h(x;y;z;u)
 

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                 a three-byte GTO
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.
 
 

  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              a three-byte GTO
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.
 
 

  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         a three-byte GTO
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.
 
 

  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            a three-byte GTO
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.
 

  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             a three-byte GTO
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

-------------------------------------------------------------------- ( new )  ------------------------------------------------------------------------------------

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

  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               if you replace this line by   ABS   ST+ 05   you'll get a less optimistic error estimate ( in absolute value )
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            a three-byte GTO
151  RCL 02
152  RCL 01
153  RTN
154  GTO 10            a three-byte GTO
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.
 
 

  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               if you replace this line by   ABS   ST+ 07   you'll get a less optimistic error estimate ( in absolute value )
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               if you replace this line by   ABS   ST+ 06   you'll get a less optimistic error estimate ( in absolute value )
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            a three-byte GTO
242  RCL 03
243  RCL 02
244  RCL 01
245  RTN
246  GTO 10            a three-byte GTO
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 estimates    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
                                      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 )
                                          and  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.
 

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                 replace this line by  ABS  ST+ 06   after replacing line 70 by  RCL IND 10  if you want a less optimistic error-estimate ( in magnitude )
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            a three-byte GTO
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.
 

  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                 replace line 79  by  ABS  ST+ 07  and line 77 by  ABS  ST+ 08   if you want a less optimistic error-estimate ( in magnitude )
  78  X<>Y
  79  ST- 07
  80  RCL 04
  81  ST+ 01
  82  DSE 09
  83  GTO 01               a three-byte GTO
  84  RCL 03
  85  RCL 02
  86  RCL 01
  87  RTN
  88  GTO 10               a three-byte GTO
  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 ).
 
 

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

Go back to the HP-41 software library
Go back to the general software library
Go back to the main exhibit hall