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