hp41programs

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

Overview

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

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

2°) A Runge-Kutta Method with order 6

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

3°) A Runge-Kutta Method with order 8

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

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

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

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

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

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

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

7°) Second-Order Equations without dy/dx

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

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

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

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

RK4 uses the following formulae:

         k₁ = h.f(x;y)
         k₂ = h.f(x+h/2;y+k₁/2)
         k₃ = h.f(x+h/2;y+k₂/2)
         k₄ = h.f(x+h;y+k₃)

and then, y(x+h) = y(x) + ( k₁ + 2.k₂ + 2.k₃ + k₄ ) / 6 This formula duplicates the Taylor series up to the term in h⁴

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

• R01 = x₀ • R03 = h/2 = half of the stepsize R05 & R06: temp ( storing h/2 instead of h saves several bytes )
• R02 = y₀ • 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(x₀+N.h)
X	/	x₀+N.h

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

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 2⁴ = 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(x₀) = y₀ , z(x₀) = z₀

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

                         k₁ = h.f(x;y;z)                                             l₁ = h.g(x;y;z)
                         k₂ = h.f(x+h/2;y+k₁/2;z+l₁/2)                      l₂ = h.g(x+h/2;y+k₁/2;z+l₁/2)
                         k₃ = h.f(x+h/2;y+k₂/2;z+l₂/2)                      l₃ = h.g(x+h/2;y+k₂/2;z+l₂/2)
                         k₄ = h.f(x+h;y+k₃;z+l₃)                               l₄ = h.g(x+h;y+k₃;z+l₃)

then, y(x+h) = y(x) + ( k₁ + 2.k₂ + 2.k₃ + k₄ ) / 6
and z(x+h) = z(x) + ( l₁ + 2.l₂ + 2.l₃ + l₄ ) / 6 duplicate the Taylor series through the terms in h⁴

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

•   R00: f name   •   R01 = x₀      •   R04 = h/2 = half of the stepsize              R06 to R08: temp            ( storing h/2 instead of h saves several bytes )
                            •   R02 = y₀      •   R05 = N = number of steps
                            •   R03 = z₀
Flags: /
Subroutine: a program calculating f(x;y;z) in Y-register and g(x;y;z) in X-register
                      assuming x is in X-register , y is in Y-register and z is in Z-register upon entry.

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

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

01 LBL "RK4B"
02 RCL 05
03 STO 08
04 LBL 01
05 RCL 03
06 RCL 02
07 RCL 01
08 XEQ IND 00
09 RCL 04
10 ST+ 01
11 ST* Z
12 *
13 STO 07
14 RCL 03

15 +
16 X<>Y
17 STO 06
18 RCL 02
19 +
20 RCL 01
21 XEQ IND 00
22 RCL 04
23 ST* Z
24 *
25 ST+ 07
26 ST+ 07
27 RCL 03
28 +

29 X<>Y
30 ST+ 06
31 ST+ 06
32 RCL 02
33 +
34 RCL 01
35 XEQ IND 00
36 RCL 04
37 ST+ 01
38 ST+ X
39 ST* Z
40 *
41 ST+ 07
42 RCL 03

43 +
44 X<>Y
45 ST+ 06
46 RCL 02
47 +
48 RCL 01
49 XEQ IND 00
50 RCL 04
51 ST* Z
52 *
53 RCL 07
54 +
55 3
56 /

57 ST+ 03
58 X<>Y
59 RCL 06
60 +
61 3
62 /
63 ST+ 02
64 DSE 08
65 GTO 01
66 RCL 03
67 RCL 02
68 RCL 01
69 END

( 99 bytes / SIZE 009 )

STACK	INPUTS	OUTPUTS
Z	/	z(x₀+N.h)
Y	/	y(x₀+N.h)
X	/	x₀+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 ( -x² ) ; y' = -2.x exp ( -x² ) )

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

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

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

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

N.B:

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

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

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

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

•   R00: f name   •   R01 = x₀      •   R05 = h/2 = half of the stepsize              R07 to R10: temp
                            •   R02 = y₀      •   R06 = N = number of steps
                            •   R03 = z₀
                            •   R04 = u₀
Flags: /
Subroutine: a program calculating f(x;y;z;u) in Z-register , g(x;y;z;u) in Y-register and h(x;y;z;u) in X-register
                      assuming x is in X-register , y is in Y-register , z is in Z-register and u is in T-register upon entry.

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

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

-Line 89 is a three-byte GTO 01

01 LBL "RK4C"
02 RCL 06
03 STO 10
04 LBL 01
05 RCL 04
06 RCL 03
07 RCL 02
08 RCL 01
09 XEQ IND 00
10 RCL 05
11 ST+ 01
12 ST* Z
13 ST* T
14 *
15 STO 09
16 RCL 04
17 +
18 X<>Y
19 STO 08

20 RCL 03
21 +
22 R^
23 STO 07
24 RCL 02
25 +
26 RCL 01
27 XEQ IND 00
28 RCL 05
29 ST* Z
30 ST* T
31 *
32 ST+ 09
33 ST+ 09
34 RCL 04
35 +
36 X<>Y
37 ST+ 08
38 ST+ 08

39 RCL 03
40 +
41 R^
42 ST+ 07
43 ST+ 07
44 RCL 02
45 +
46 RCL 01
47 XEQ IND 00
48 RCL 05
49 ST+ 01
50 ST+ X
51 ST* Z
52 ST* T
53 *
54 ST+ 09
55 RCL 04
56 +
57 X<>Y

58 ST+ 08
59 RCL 03
60 +
61 R^
62 ST+ 07
63 RCL 02
64 +
65 RCL 01
66 XEQ IND 00
67 RCL 05
68 ST* Z
69 ST* T
70 *
71 RCL 09
72 +
73 3
74 /
75 ST+ 04
76 X<>Y

77 RCL 08
78 +
79 3
80 /
81 ST+ 03
82 R^
83 RCL 07
84 +
85 3
86 /
87 ST+ 02
88 DSE 10
89 GTO 01
90 RCL 04
91 RCL 03
92 RCL 02
93 RCL 01
94 END

( 133 bytes / SIZE 011 )

STACK	INPUTS	OUTPUTS
T	/	u(x₀+N.h)
Z	/	z(x₀+N.h)
Y	/	y(x₀+N.h)
X	/	x₀+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 h⁶
-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 = x₀ • R03 = h = the stepsize R05: counter
• R02 = y₀ • R04 = N = the number of steps R06 , ...... , R11 = k₁ , ..... , k₆

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

-Line 164 is a three-byte GTO 01

01 LBL "RK6"
02 RCL 04
03 STO 05
04 LBL 01
05 RCL 02
06 RCL 01
07 XEQ IND 00
08 RCL 03
09 *
10 STO 06
11 3
12 /
13 RCL 02
14 +
15 RCL 03
16 3
17 /
18 RCL 01
19 +
20 STO 10
21 XEQ IND 00
22 RCL 03
23 *
24 STO 07
25 1.5
26 /
27 RCL 02
28 ST+ Y
29 CLX
30 RCL 03
31 LASTX
32 /
33 RCL 01
34 +

35 XEQ IND 00
36 RCL 03
37 *
38 STO 08
39 RCL 07
40 4
41 *
42 X<>Y
43 -
44 RCL 06
45 +
46 12
47 /
48 RCL 02
49 +
50 RCL 10
51 XEQ IND 00
52 RCL 03
53 *
54 STO 09
55 18
56 *
57 RCL 08
58 7
59 *
60 +
61 RCL 07
62 22
63 *
64 -
65 RCL 06
66 5
67 *
68 +

69 9.6
70 /
71 RCL 02
72 +
73 RCL 03
74 1.2
75 /
76 RCL 01
77 +
78 XEQ IND 00
79 RCL 03
80 *
81 STO 10
82 5
83 /
84 RCL 09
85 +
86 RCL 08
87 4
88 /
89 -
90 RCL 06
91 18
92 *
93 RCL 07
94 55
95 *
96 -
97 60
98 /
99 +
100 2
101 /
102 RCL 02

103 +
104 RCL 03
105 6
106 /
107 RCL 01
108 +
109 XEQ IND 00
110 RCL 03
111 ST+ 01
112 *
113 STO 11
114 80
115 *
116 RCL 07
117 99
118 *
119 +
120 RCL 09
121 118
122 *
123 -
124 20
125 *
126 RCL 10
127 ST+ 11
128 128
129 *
130 +
131 RCL 08
132 ST+ 09
133 215
134 *
135 +
136 RCL 06

137 783
138 *
139 -
140 780
141 /
142 RCL 02
143 +
144 RCL 01
145 XEQ IND 00
146 RCL 03
147 *
148 RCL 06
149 +
150 13
151 *
152 RCL 11
153 32
154 *
155 +
156 RCL 09
157 55
158 *
159 +
160 .5
161 %
162 ST+ 02
163 DSE 05
164 GTO 01
165 RCL 02
166 RCL 01
167 END

( 219 bytes / SIZE 012 )

STACK	INPUTS	OUTPUTS
Y	/	y(x₀+N.h)
X	/	x₀+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 2⁶ = 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 = x₀      •   R04 = h = stepsize                         R06 to R16: temp
                            •   R02 = y₀      •   R05 = N = number of steps
                            •   R03 = z₀
Flags: /
Subroutine: a program calculating f(x;y;z) in Y-register and g(x;y;z) in X-register
                      assuming x is in X-register , y is in Y-register and z is in Z-register upon entry.

-Line 277 is a three-byte GTO 01

01 LBL "RK6B"
02 RCL 05
03 STO 16
04 GTO 01
05 LBL 00
06 XEQ IND 00
07 RCL 04
08 ST* Z
09 *
10 RTN
11 LBL 01
12 RCL 03
13 RCL 02
14 RCL 01
15 XEQ 00
16 STO 11
17 3
18 /
19 RCL 03
20 +
21 X<>Y
22 STO 06
23 3
24 /
25 RCL 02
26 +
27 RCL 04
28 3
29 /
30 RCL 01
31 +
32 STO 10
33 XEQ 00
34 STO 12
35 1.5
36 /
37 RCL 03
38 ST+ Y
39 X<> Z
40 STO 07
41 LASTX
42 /
43 RCL 02
44 ST+ Y
45 CLX
46 RCL 04
47 LASTX

48 /
49 RCL 01
50 +
51 XEQ 00
52 STO 13
53 RCL 12
54 4
55 *
56 X<>Y
57 -
58 RCL 11
59 +
60 12
61 /
62 RCL 03
63 +
64 RCL 07
65 4
66 *
67 R^
68 STO 08
69 -
70 RCL 06
71 +
72 12
73 /
74 RCL 02
75 +
76 RCL 10
77 XEQ 00
78 STO 14
79 18
80 *
81 RCL 13
82 7
83 *
84 +
85 RCL 12
86 22
87 *
88 -
89 RCL 11
90 5
91 *
92 +
93 9.6
94 /

95 RCL 03
96 +
97 X<>Y
98 STO 09
99 18
100 *
101 RCL 08
102 7
103 *
104 +
105 RCL 07
106 22
107 *
108 -
109 RCL 06
110 5
111 *
112 +
113 9.6
114 /
115 RCL 02
116 +
117 RCL 04
118 1.2
119 /
120 RCL 01
121 +
122 XEQ 00
123 STO 15
124 10
125 /
126 RCL 14
127 2
128 /
129 +
130 RCL 13
131 8
132 /
133 -
134 RCL 12
135 11
136 *
137 24
138 /
139 -
140 RCL 11
141 .15

142 *
143 +
144 RCL 03
145 +
146 X<>Y
147 STO 10
148 10
149 /
150 RCL 09
151 2
152 /
153 +
154 RCL 08
155 8
156 /
157 -
158 RCL 07
159 11
160 *
161 24
162 /
163 -
164 RCL 06
165 .15
166 *
167 +
168 RCL 02
169 +
170 RCL 04
171 6
172 /
173 RCL 01
174 +
175 XEQ 00
176 RCL 12
177 99
178 *
179 X<>Y
180 STO 12
181 80
182 *
183 +
184 RCL 14
185 118
186 *
187 -
188 20

189 *
190 RCL 15
191 ST+ 12
192 128
193 *
194 +
195 RCL 13
196 ST+ 14
197 215
198 *
199 +
200 RCL 11
201 783
202 *
203 -
204 780
205 /
206 RCL 03
207 +
208 RCL 07
209 99
210 *
211 R^
212 STO 07
213 80
214 *
215 +
216 RCL 09
217 118
218 *
219 -
220 20
221 *
222 RCL 10
223 ST+ 07
224 128
225 *
226 +
227 RCL 08
228 ST+ 09
229 215
230 *
231 +
232 RCL 06
233 783
234 *
235 -

236 780
237 /
238 RCL 02
239 +
240 RCL 01
241 RCL 04
242 +
243 STO 01
244 XEQ 00
245 RCL 11
246 +
247 13
248 *
249 RCL 12
250 32
251 *
252 +
253 RCL 14
254 55
255 *
256 +
257 .5
258 %
259 ST+ 03
260 R^
261 RCL 06
262 +
263 13
264 *
265 RCL 07
266 32
267 *
268 +
269 RCL 09
270 55
271 *
272 +
273 .5
274 %
275 ST+ 02
276 DSE 16
277 GTO 01
278 RCL 03
279 RCL 02
280 RCL 01
281 END

( 378 bytes / SIZE 017 )

STACK	INPUTS	OUTPUTS
Z	/	z(x₀+N.h)
Y	/	y(x₀+N.h)
X	/	x₀+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 2⁶ = 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 h⁸
-The following program uses 41 coefficients which are of the form ( a + b 21^1/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 = x₀ • R03 = h = the stepsize R05 to R10 and R52: temp
• R02 = y₀ • R04 = N = the number of steps • R11 to R51: the 41 coefficients listed below

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

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

-Line 251 is a three-byte GTO 01

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

44 +
45 RCL 05
46 RCL 19
47 /
48 +
49 RCL 02
50 +
51 RCL 03
52 RCL 11
53 *
54 RCL 01
55 +
56 XEQ 00
57 STO 08
58 RCL 20
59 *
60 RCL 07
61 RCL 21
62 *
63 +
64 RCL 05
65 RCL 22
66 *
67 +
68 RCL 02
69 +
70 RCL 03
71 RCL 11
72 *
73 RCL 01
74 +
75 XEQ 00
76 STO 09
77 RCL 23
78 *
79 RCL 08
80 RCL 24
81 *
82 +
83 RCL 07
84 RCL 25
85 *
86 +

87 RCL 05
88 RCL 26
89 *
90 +
91 RCL 02
92 +
93 RCL 03
94 RCL 12
95 *
96 RCL 01
97 +
98 XEQ 00
99 STO 10
100 RCL 27
101 *
102 RCL 09
103 RCL 28
104 *
105 +
106 RCL 08
107 RCL 29
108 *
109 +
110 RCL 07
111 RCL 30
112 *
113 +
114 RCL 05
115 RCL 31
116 *
117 +
118 RCL 02
119 +
120 RCL 03
121 RCL 13
122 *
123 RCL 01
124 +
125 XEQ 00
126 STO 06
127 RCL 14
128 /
129 RCL 10

130 RCL 32
131 *
132 +
133 RCL 09
134 RCL 33
135 *
136 +
137 RCL 05
138 RCL 15
139 /
140 +
141 RCL 02
142 +
143 RCL 03
144 RCL 13
145 *
146 RCL 01
147 +
148 XEQ 00
149 STO 07
150 RCL 34
151 *
152 RCL 06
153 RCL 35
154 *
155 +
156 RCL 10
157 RCL 36
158 *
159 +
160 RCL 09
161 RCL 37
162 *
163 +
164 RCL 05
165 RCL 38
166 *
167 +
168 RCL 02
169 +
170 RCL 03
171 RCL 12
172 *

173 RCL 01
174 +
175 XEQ 00
176 STO 08
177 RCL 39
178 *
179 RCL 07
180 RCL 40
181 *
182 +
183 RCL 06
184 RCL 41
185 *
186 +
187 RCL 10
188 RCL 42
189 *
190 +
191 RCL 09
192 RCL 14
193 /
194 +
195 RCL 05
196 RCL 15
197 /
198 +
199 RCL 02
200 +
201 RCL 01
202 RCL 03
203 ST+ 01
204 RCL 11
205 *
206 +
207 XEQ 00
208 X<> 09
209 RCL 48
210 *
211 RCL 08
212 RCL 44
213 *
214 +
215 RCL 07

216 RCL 45
217 *
218 +
219 RCL 06
220 RCL 46
221 *
222 +
223 RCL 10
224 RCL 47
225 *
226 +
227 RCL 09
228 RCL 43
229 *
230 +
231 RCL 02
232 +
233 RCL 01
234 XEQ 00
235 RCL 05
236 +
237 RCL 49
238 *
239 RCL 09
240 RCL 07
241 +
242 RCL 50
243 *
244 +
245 RCL 08
246 RCL 51
247 /
248 +
249 ST+ 02
250 DSE 52
251 GTO 01
252 RCL 02
253 RCL 01
254 END

( 329 bytes / SIZE 053 )

STACK	INPUTS	OUTPUTS
Y	/	y(x₀+N.h)
X	/	x₀+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 2⁸ = 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 = x₀      •   R04 = h = stepsize                         R06 to R17: temp    R59: counter
                            •   R02 = y₀      •   R05 = N = number of steps       •   R18 thru R58 = the same 41 numbers used in the previous program.
                            •   R03 = z₀
Flags: /
Subroutine: a program calculating f(x;y;z) in Y-register and g(x;y;z) in X-register
                       assuming x is in X-register , y is in Y-register and z is in Z-register upon entry.

-Line 444 is a three-byte GTO 01

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

76 RCL 26
77 /
78 +
79 RCL 02
80 +
81 RCL 04
82 RCL 18
83 *
84 RCL 01
85 +
86 XEQ 00
87 STO 15
88 RCL 27
89 *
90 RCL 14
91 RCL 28
92 *
93 +
94 RCL 12
95 RCL 29
96 *
97 +
98 RCL 03
99 +
100 X<>Y
101 STO 09
102 RCL 27
103 *
104 RCL 08
105 RCL 28
106 *
107 +
108 RCL 06
109 RCL 29
110 *
111 +
112 RCL 02
113 +
114 RCL 04
115 RCL 18
116 *
117 RCL 01
118 +
119 XEQ 00
120 STO 16
121 RCL 30
122 *
123 RCL 15
124 RCL 31
125 *
126 +
127 RCL 14
128 RCL 32
129 *
130 +
131 RCL 12
132 RCL 33
133 *
134 +
135 RCL 03
136 +
137 X<>Y
138 STO 10
139 RCL 30
140 *
141 RCL 09
142 RCL 31
143 *
144 +
145 RCL 08
146 RCL 32
147 *
148 +
149 RCL 06
150 RCL 33

151 *
152 +
153 RCL 02
154 +
155 RCL 04
156 RCL 19
157 *
158 RCL 01
159 +
160 XEQ 00
161 STO 17
162 RCL 34
163 *
164 RCL 16
165 RCL 35
166 *
167 +
168 RCL 15
169 RCL 36
170 *
171 +
172 RCL 14
173 RCL 37
174 *
175 +
176 RCL 12
177 RCL 38
178 *
179 +
180 RCL 03
181 +
182 X<>Y
183 STO 11
184 RCL 34
185 *
186 RCL 10
187 RCL 35
188 *
189 +
190 RCL 09
191 RCL 36
192 *
193 +
194 RCL 08
195 RCL 37
196 *
197 +
198 RCL 06
199 RCL 38
200 *
201 +
202 RCL 02
203 +
204 RCL 04
205 RCL 20
206 *
207 RCL 01
208 +
209 XEQ 00
210 STO 13
211 RCL 21
212 /
213 RCL 17
214 RCL 39
215 *
216 +
217 RCL 16
218 RCL 40
219 *
220 +
221 RCL 12
222 RCL 22
223 /
224 +
225 RCL 03

226 +
227 X<>Y
228 STO 07
229 RCL 21
230 /
231 RCL 11
232 RCL 39
233 *
234 +
235 RCL 10
236 RCL 40
237 *
238 +
239 RCL 06
240 RCL 22
241 /
242 +
243 RCL 02
244 +
245 RCL 04
246 RCL 20
247 *
248 RCL 01
249 +
250 XEQ 00
251 STO 14
252 RCL 41
253 *
254 RCL 13
255 RCL 42
256 *
257 +
258 RCL 17
259 RCL 43
260 *
261 +
262 RCL 16
263 RCL 44
264 *
265 +
266 RCL 12
267 RCL 45
268 *
269 +
270 RCL 03
271 +
272 X<>Y
273 STO 08
274 RCL 41
275 *
276 RCL 07
277 RCL 42
278 *
279 +
280 RCL 11
281 RCL 43
282 *
283 +
284 RCL 10
285 RCL 44
286 *
287 +
288 RCL 06
289 RCL 45
290 *
291 +
292 RCL 02
293 +
294 RCL 04
295 RCL 19
296 *
297 RCL 01
298 +
299 XEQ 00
300 STO 15

301 RCL 46
302 *
303 RCL 14
304 RCL 47
305 *
306 +
307 RCL 13
308 RCL 48
309 *
310 +
311 RCL 17
312 RCL 49
313 *
314 +
315 RCL 16
316 RCL 21
317 /
318 +
319 RCL 12
320 RCL 22
321 /
322 +
323 RCL 03
324 +
325 X<>Y
326 STO 09
327 RCL 46
328 *
329 RCL 08
330 RCL 47
331 *
332 +
333 RCL 07
334 RCL 48
335 *
336 +
337 RCL 11
338 RCL 49
339 *
340 +
341 RCL 10
342 RCL 21
343 /
344 +
345 RCL 06
346 RCL 22
347 /
348 +
349 RCL 02
350 +
351 RCL 04
352 RCL 18
353 *
354 RCL 01
355 +
356 XEQ 00
357 X<> 16
358 RCL 55
359 *
360 RCL 15
361 RCL 51
362 *
363 +
364 RCL 14
365 RCL 52
366 *
367 +
368 RCL 13
369 RCL 53
370 *
371 +
372 RCL 17
373 RCL 54
374 *
375 +

376 RCL 16
377 RCL 50
378 *
379 +
380 RCL 03
381 +
382 X<>Y
383 X<> 10
384 RCL 55
385 *
386 RCL 09
387 RCL 51
388 *
389 +
390 RCL 08
391 RCL 52
392 *
393 +
394 RCL 07
395 RCL 53
396 *
397 +
398 RCL 11
399 RCL 54
400 *
401 +
402 RCL 10
403 RCL 50
404 *
405 +
406 RCL 02
407 +
408 RCL 04
409 RCL 01
410 +
411 STO 01
412 XEQ 00
413 RCL 12
414 +
415 RCL 56
416 *
417 RCL 16
418 RCL 14
419 +
420 RCL 57
421 *
422 +
423 RCL 15
424 RCL 58
425 ST/ 09
426 /
427 +
428 ST+ 03
429 X<>Y
430 RCL 06
431 +
432 RCL 56
433 *
434 RCL 10
435 RCL 08
436 +
437 RCL 57
438 *
439 +
440 RCL 09
441 +
442 ST+ 02
443 DSE 59
444 GTO 01
445 RCL 03
446 RCL 02
447 RCL 01
448 END

( 593 bytes / SIZE 060 )

STACK	INPUTS	OUTPUTS
Z	/	z(x₀+N.h)
Y	/	y(x₀+N.h)
X	/	x₀+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 2⁸ = 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 ( n² + 3.n - 2 ) / 2 coefficients a_i,j ; b_i ; c_i

    k₁ = h.f(x;y)
    k₂ = h.f(x+b₂.h;y+a_2,1.k₁)
    k₃ = h.f(x+b₃.h;y+a_3,1.k₁+a_3,2.k₂)
    .............................................................                                         ( actually,     a_i,1 + a_i,2 + .......... + a_i,i-1 = b_i )

k_n = h.f(x+b_n.h;y+a_n,1.k₁+a_n,2.k₂+ .......... +a_n,n-1.k_n-1)

then, y(x+h) = y(x) + c₁.k₁+ ................ + c_n.k_n

-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(n²+5.n+16)/2 are to be initialized before executing "ERK" )

•   R00: f name   •   R01 = x₀                                      R06 to R09: temp             • Rn+10 to R(n²+5.n+16)/2 = the ( n² + 3.n - 2 ) / 2 coefficients a_i,j ; b_i ; c_i
                            •   R02 = y₀                                      R10 = k₁                                                                            which are to be stored as shown below:
                            •   R03 = h = stepsize                        R11 = k₂
                            •   R04 = N = number of steps          ............
                            •   R05 = n = number of stages        Rn+9 = k_n

• Rn+10 = a_2,1 • Rn+11 = b₂
• Rn+12 = a_3,1 • Rn+13 = a_3,2 • Rn+14 = b₃

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

• R(n²+n+18)/2 = a_n,1 ................ • R(n²+3n+14)/2 = a_n,n-1 • R(n²+3n+16)/2 = b_n

• R(n²+3n+18)/2 = c₁ ............................................................ • R(n²+5n+16)/2 = c_n

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 ( n²+5.n+18 )/2 )

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

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

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

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

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

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

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

-Initialize registers R14 thru R26.

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

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

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

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

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

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

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

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

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

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

•   R00: f name   •   R01 = x₀                                      R07 to R11: temp             • R2n+12 to R(n²+7n+20)/2 = the ( n² + 3.n - 2 ) / 2 coefficients a_i,j ; b_i ; c_i
                            •   R02 = y₀                                      R12 to Rn+11 = k_i^y                                                            which are to be stored as shown below:
                            •   R03 = z₀                                      Rn+12 to R2n+11 = k_i^z
                            •   R04 = h   = stepsize
                            •   R05 = N = number of steps
                            •   R06 = n   = number of stages

• R2n+12 = a_2,1 • R2n+13 = b₂
• R2n+14 = a_3,1 • R2n+15 = a_3,2 • R2n+16 = b₃

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

• R(n²+3n+22)/2 = a_n,1 ................ • R(n²+5n+18)/2 = a_n,n-1 • R(n²+5n+20)/2 = b_n

• R(n²+5n+22)/2 = c₁ ............................................................ • R(n²+7n+20)/2 = c_n

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 ( n²+ 7.n + 22 )/2 )

STACK	INPUTS	OUTPUTS
Z	/	z(x₀+N.h)
Y	/	y(x₀+N.h)
X	/	x₀+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 = x₀ • R03 = h = stepsize R05 = error estimate
• R02 = y₀ • R04 = N = number of steps R06 thru R10: temp

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

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

01 LBL "RKF"
02 CLX
03 STO 05
04 LBL 10
05 RCL 04
06 STO 06
07 LBL 01
08 RCL 02
09 RCL 01
10 XEQ IND 00
11 RCL 03
12 *
13 STO 07
14 ST+ X
15 9
16 /
17 RCL 02
18 +
19 RCL 03
20 ST+ X
21 9
22 /
23 RCL 01
24 +
25 XEQ IND 00
26 RCL 03
27 *
28 4
29 /
30 STO 08
31 RCL 07

32 12
33 /
34 +
35 RCL 02
36 +
37 RCL 03
38 3
39 /
40 RCL 01
41 +
42 XEQ IND 00
43 RCL 03
44 *
45 STO 09
46 270
47 *
48 RCL 08
49 972
50 *
51 -
52 RCL 07
53 69
54 *
55 +
56 128
57 /
58 RCL 02
59 +
60 RCL 03
61 .75
62 *

63 RCL 01
64 +
65 XEQ IND 00
66 RCL 03
67 ST+ 01
68 *
69 64
70 *
71 STO 10
72 RCL 07
73 85
74 *
75 -
76 60
77 /
78 RCL 08
79 RCL 09
80 5
81 /
82 -
83 27
84 *
85 +
86 RCL 02
87 +
88 RCL 01
89 XEQ IND 00
90 RCL 03
91 *
92 15
93 *

94 RCL 10
95 +
96 STO 10
97 RCL 09
98 351
99 *
100 +
101 RCL 08
102 540
103 *
104 -
105 RCL 07
106 65
107 *
108 +
109 432
110 /
111 RCL 02
112 +
113 RCL 01
114 RCL 03
115 6
116 /
117 -
118 XEQ IND 00
119 RCL 03
120 *
121 72
122 *
123 RCL 10
124 -

125 RCL 09
126 9
127 *
128 +
129 RCL 07
130 ST+ X
131 -
132 3
133 1/X
134 %
135 ST- 05
136 RCL 10
137 RCL 09
138 81
139 *
140 +
141 PI
142 R-D
143 /
144 RCL 07
145 9
146 /
147 +
148 ST+ 02
149 DSE 06
150 GTO 01
151 RCL 02
152 RCL 01
153 RTN
154 GTO 10
155 END

( 204 bytes / SIZE 011 )

STACK	INPUTS	OUTPUTS
Y	/	y(x₀+N.h)
X	/	x₀+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 = x₀      •   R04 = h = stepsize                         R06 = y-error estimate           R08 to R16: temp
                            •   R02 = y₀      •   R05 = N = number of steps           R07 = z-error estimate
                            •   R03 = z₀
Flags: /
Subroutine: a program calculating f(x;y;z) in Y-register and g(x;y;z) in X-register
                       assuming x is in X-register , y is in Y-register and z is in Z-register upon entry.

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

01 LBL "RKFB"
02 CLX
03 STO 06
04 STO 07
05 LBL 10
06 RCL 05
07 STO 16
08 LBL 01
09 RCL 03
10 RCL 02
11 RCL 01
12 XEQ 00
13 STO 12
14 ST+ X
15 9
16 /
17 RCL 03
18 +
19 X<>Y
20 STO 08
21 ST+ X
22 9
23 /
24 RCL 02
25 +
26 RCL 04
27 ST+ X
28 9
29 /
30 RCL 01
31 +
32 XEQ 00
33 STO 13
34 4
35 /
36 RCL 12
37 12
38 /
39 +
40 RCL 03
41 +
42 X<>Y
43 STO 09

44 4
45 /
46 RCL 08
47 12
48 /
49 +
50 RCL 02
51 +
52 RCL 04
53 3
54 /
55 RCL 01
56 +
57 XEQ 00
58 STO 14
59 270
60 *
61 RCL 13
62 243
63 *
64 -
65 RCL 12
66 69
67 *
68 +
69 128
70 /
71 RCL 03
72 +
73 X<>Y
74 STO 10
75 270
76 *
77 RCL 09
78 243
79 *
80 -
81 RCL 08
82 69
83 *
84 +
85 128
86 /

87 RCL 02
88 +
89 RCL 04
90 .75
91 *
92 RCL 01
93 +
94 XEQ 00
95 64
96 ST* Z
97 *
98 STO 15
99 RCL 14
100 324
101 *
102 -
103 RCL 13
104 405
105 *
106 +
107 RCL 12
108 85
109 *
110 -
111 60
112 /
113 RCL 03
114 +
115 X<>Y
116 STO 11
117 RCL 10
118 324
119 *
120 -
121 RCL 09
122 405
123 *
124 +
125 RCL 08
126 85
127 *
128 -
129 60

130 /
131 RCL 02
132 +
133 RCL 01
134 RCL 04
135 +
136 STO 01
137 XEQ 00
138 15
139 ST* Z
140 *
141 RCL 15
142 +
143 STO 15
144 RCL 14
145 351
146 *
147 +
148 RCL 13
149 135
150 ST* 09
151 *
152 -
153 RCL 12
154 65
155 *
156 +
157 432
158 /
159 RCL 03
160 +
161 X<>Y
162 RCL 11
163 +
164 STO 11
165 RCL 10
166 351
167 *
168 +
169 RCL 09
170 -
171 RCL 08
172 65

173 *
174 +
175 432
176 /
177 RCL 02
178 +
179 RCL 04
180 6
181 /
182 RCL 01
183 X<>Y
184 -
185 XEQ 00
186 72
187 ST* Z
188 *
189 RCL 15
190 -
191 RCL 14
192 9
193 *
194 +
195 RCL 12
196 ST+ X
197 -
198 3
199 1/X
200 %
201 ST- 07
202 R^
203 RCL 11
204 -
205 RCL 10
206 9
207 *
208 +
209 RCL 08
210 ST+ X
211 -
212 3
213 1/X
214 %
215 ST- 06

216 RCL 15
217 RCL 14
218 81
219 ST* 10
220 *
221 +
222 PI
223 R-D
224 /
225 RCL 12
226 9
227 ST/ 08
228 /
229 +
230 ST+ 03
231 RCL 11
232 RCL 10
233 +
234 PI
235 R-D
236 /
237 RCL 08
238 +
239 ST+ 02
240 DSE 16
241 GTO 01
242 RCL 03
243 RCL 02
244 RCL 01
245 RTN
246 GTO 10
247 LBL 00
248 XEQ IND 00
249 RCL 04
250 ST* Z
251 *
252 RTN
253 END

( 343 bytes / SIZE 017)

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

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

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

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

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

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

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

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

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

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

-All n stage explicit Runge-Kutta formulae with error estimate are determined by ( n² + 5.n - 2 ) / 2 coefficients a_i,j ; b_i ; c_i ; d_i

k_n = h.f(x+b_n.h;y+a_n,1.k₁+a_n,2.k₂+ .......... +a_n,n-1.k_n-1)

then, y(x+h) = y(x) + c₁.k₁+ ................ + c_n.k_n and the difference between the higher-order and the lower-order formulae gives an error estimate:

err-estim = d₁.k₁+ ................ + d_n.k_n

-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(n²+7.n+18)/2 are to be initialized before executing "RKE" )

•   R00: f name   •   R01 = x₀                                     R06 = y-error             • Rn+11 to R(n²+7.n+18)/2 = the ( n² + 5.n - 2 ) / 2 coefficients a_i,j ; b_i ; c_i ; d_i
                            •   R02 = y₀                                      R07 to R10: temp                                                         which are to be stored as shown below:
                            •   R03 = h = stepsize                        R11 to Rn+10 = k_i
                            •   R04 =N = number of steps
                            •   R05 = n = number of stages

• Rn+11 = a_2,1 • Rn+12 = b₂
• Rn+13 = a_3,1 • Rn+14 = a_3,2 • Rn+15 = b₃

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

• R(n²+n+20)/2 = a_n,1 ................ • R(n²+3n+16)/2 = a_n,n-1 • R(n²+3n+18)/2 = b_n

• R(n²+3n+20)/2 = c₁ ............................................................ • R(n²+5n+18)/2 = c_n

• R(n²+5n+20)/2 = d₁ ............................................................ • R(n²+7n+18)/2 = d_n

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

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

01 LBL "RKE"
02 CLX
03 STO 06
04 LBL 10
05 RCL 04
06 STO 07
07 LBL 01
08 RCL 05
09 11.9
10 STO 08
11 +
12 STO 09
13 RCL 05
14 10
15 +
16 .1
17 %
18 12

19 +
20 STO 10
21 RCL 02
22 RCL 01
23 XEQ IND 00
24 RCL 03
25 *
26 STO 11
27 LBL 02
28 RCL 08
29 INT
30 E3
31 /
32 11
33 +
34 STO 08
35 CLX
36 LBL 03

37 RCL IND 08
38 RCL IND 09
39 *
40 +
41 ISG 09
42 ISG 08
43 GTO 03
44 RCL 02
45 +
46 RCL 03
47 RCL IND 09
48 *
49 RCL 01
50 +
51 XEQ IND 00
52 RCL 03
53 *
54 STO IND 10

55 ISG 09
56 ISG 10
57 GTO 02
58 RCL 05
59 ST- 10
60 RCL 09
61 +
62 0
63 LBL 04
64 RCL IND Y
65 RCL IND 10
66 *
67 ST- 06
68 CLX
69 RCL IND 09
70 LASTX
71 *
72 +

73 ISG Y
74 ISG 09
75 ISG 10
76 GTO 04
77 ST+ 02
78 RCL 03
79 ST+ 01
80 DSE 07
81 GTO 01
82 RCL 02
83 RCL 01
84 RTN
85 GTO 10
86 END

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

STACK	INPUTS	OUTPUTS
Y	/	y(x₀+N.h)
X	/	x₀+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(n²+9.n+24)/2 are to be initialized before executing "RKEB" )

•   R00: f name   •   R01 = x₀                                   R07 = y-error             • R2n+14 to R(n²+9.n+24)/2 = the ( n² + 5.n - 2 ) / 2 coefficients a_i,j ; b_i ; c_i ; d_i
                            •   R02 = y₀                                   R08 = z-error                                                                which are to be stored as shown below:
                            •   R03 = y₀
                            •   R04 = h = stepsize                        R09 to R13: temp
                            •   R05 =N = number of steps           R14    to Rn+13 = k_i(y)
                            •   R06 = n = number of stages         R14+n to R2n+13 = k_i(z)

• R2n+14 = a_2,1 • R2n+15 = b₂
• R2n+16 = a_3,1 • R2n+17 = a_3,2 • R2n+18 = b₃

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

• R(n²+3n+26)/2 = a_n,1 ................ • R(n²+5n+22)/2 = a_n,n-1 • R(n²+5n+24)/2 = b_n

• R(n²+5n+26)/2 = c₁ ............................................................ • R(n²+7n+24)/2 = c_n

• R(n²+7n+26)/2 = d₁ ............................................................ • R(n²+9n+24)/2 = d_n

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

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

01 LBL "RKEB"
02 CLX
03 STO 07
04 STO 08
05 LBL 10
06 RCL 05
07 STO 09
08 LBL 01
09 14.9
10 STO 10
11 RCL 06
12 +
13 STO 11
14 LASTX
15 STO 13
16 +
17 STO 12
18 RCL 03
19 RCL 02
20 RCL 01
21 XEQ IND 00
22 RCL 04

23 ST* Z
24 *
25 STO IND 11
26 X<>Y
27 STO 14
28 DSE 13
29 LBL 02
30 RCL 10
31 INT
32 .1
33 %
34 14
35 +
36 STO 10
37 RCL 06
38 +
39 STO 11
40 CLST
41 XEQ 03
42 RCL 03
43 +
44 X<>Y

45 RCL 02
46 +
47 RCL 04
48 RCL IND 12
49 *
50 RCL 01
51 +
52 XEQ IND 00
53 RCL 04
54 ST* Z
55 *
56 STO IND 11
57 X<>Y
58 STO IND 10
59 ISG 12
60 DSE 13
61 GTO 02
62 1.001
63 RCL 06
64 -
65 ST+ 10
66 ST+ 11

67 CLST
68 XEQ 03
69 ST+ 03
70 X<>Y
71 ST+ 02
72 RCL 06
73 ST- 10
74 ST- 11
75 CLST
76 XEQ 03
77 ST- 08
78 X<>Y
79 ST- 07
80 RCL 04
81 ST+ 01
82 DSE 09
83 GTO 01
84 RCL 03
85 RCL 02
86 RCL 01
87 RTN
88 GTO 10

89 LBL 03
90 X<>Y
91 RCL IND 10
92 RCL IND 12
93 *
94 +
95 X<>Y
96 RCL IND 11
97 RCL IND 12
98 *
99 +
100 ISG 12
101 ISG 11
102 CLX
103 ISG 10
104 GTO 03
105 RTN
106 END

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

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

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

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

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

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

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

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

7°) Second-order Equations without dy/dx

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

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

y(x+h) = y(x) + h ( y'(x) + k₁/6 + 2k₂/3 )
y'(x+h) = y'(x) + k₁/6 + 2k₂/3 + k₃/6 where k₁ = h.f (x,y) ; k₂ = h.f(x+h/2,y+h.y'/2+h.k₁/8) ; k₃ = h.f(x+h,y+h.y'+h.k₂/2)

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

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

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

01 LBL "RKS1"
02 RCL 05
03 STO 08
04 LBL 01
05 RCL 02
06 RCL 01
07 XEQ IND 00
08 RCL 04
09 ST+ 01
10 *
11 STO 06
12 2

13 /
14 RCL 03
15 +
16 RCL 04
17 *
18 RCL 02
19 +
20 RCL 01
21 XEQ IND 00
22 RCL 04
23 ST+ 01
24 *

25 STO 07
26 RCL 03
27 +
28 RCL 04
29 ST+ X
30 *
31 RCL 02
32 +
33 RCL 01
34 XEQ IND 00
35 RCL 04
36 *

37 RCL 06
38 +
39 RCL 07
40 ST+ X
41 ST+ 06
42 ST+ X
43 +
44 3
45 ST/ 06
46 /
47 X<> 03
48 ST+ 03

49 RCL 06
50 +
51 RCL 04
52 ST+ X
53 *
54 ST+ 02
55 DSE 08
56 GTO 01
57 RCL 03
58 RCL 02
59 RCL 01
60 END

( 85 bytes / SIZE 009 )

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

Example: y(0) = 1 , y'(0) = 1 , y" = - y ( x² + y² )^1/2

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

-Load this short routine:

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

"T" ASTO 00

0 STO 01 STO 03 1 STO 02

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

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

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

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

Notes:

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

Reference:

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