The Museum of
HP CalculatorsThis 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
-The purpose of the following programs is to solve differential equations with various Runge-Kutta methods:
1°) RK4 is the classical 4th order Runge-Kutta
formula
2°) RKF is a Runge-Kutta-Fehlberg method of
order 4 ( embedded within 5th order )
3°) RK6 uses a 6th order Runge-Kutta -----
4°) RK8 is an 8th order method
5°) ERK is a general Runge-Kutta program suitable
to all explicit formulae
6°) IRK8 uses an implicit Runge-Kutta
method of order 8.
-All these methods are applied to solve:
a) A first order differential equation:
dy/dx = f(x,y)
with the initial value: y(x0) = y0
b) A system of 2 differential equations:
dy/dx = f(x,y,z) ; dz/dx = g(x,y,z)
------------------- y(x0) = y0
; z(x0) = z0
1°) The classical 4th order Runge-Kutta method
A) First order differential
equations: dy/dx = f (x;y)
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 serie through 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 preceding formulae are now generalized to a system of 2 differential equations ( see also "systems of differential equations for the HP-41" ):
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' + 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
-In the following program, a fourth order Runge-Kutta
method is imbedded within a fifth order Runge-Kutta method:
the difference of the results provides an error estimate for
the 4th order method.
2°) A 4th order Runge-Kutta-Fehlberg method
A) First order differential
equations: dy/dx = f (x;y)
Data Registers: Registers R00 to R05 are to be initialized before executing "RKF"
• R00: f name • R01
= x0 • R03 = h = stepsize
• R05 = the error estimate ( clear this register before
the first execution )
• R02 = y0 •
R04 = N = number of steps
R06 thru R11: temp
Flags: /
Subroutine: a program calculating f(x;y)
assuming x is in X-register and y is in Y-register upon entry.
001 LBL "RKF"
002 RCL 04
003 STO 06
004 LBL 01
005 RCL 02
006 RCL 01
007 XEQ IND 00
008 RCL 03
009 *
010 STO 07
011 ST+ X
012 9
013 /
014 RCL 02
015 +
016 RCL 03
017 ST+ X
018 9
019 /
020 RCL 01
021 +
022 XEQ IND 00
023 RCL 03
024 *
025 4
026 /
027 STO 08
028 RCL 07
029 12
030 /
031 +
032 RCL 02
033 +
034 RCL 03
035 3
036 /
037 RCL 01
038 +
039 XEQ IND 00
040 RCL 03
041 *
042 STO 09
043 270
044 *
045 RCL 08
046 972
047 *
048 -
049 RCL 07
050 69
051 *
052 +
053 128
054 /
055 RCL 02
056 +
057 RCL 03
058 .75
059 *
060 RCL 01
061 +
062 XEQ IND 00
063 RCL 03
064 ST+ 01
065 *
066 64
067 *
068 STO 10
069 RCL 07
070 85
071 *
072 -
073 60
074 /
075 RCL 08
076 RCL 09
077 5
078 /
079 -
080 27
081 *
082 +
083 RCL 02
084 +
085 RCL 01
086 XEQ IND 00
087 RCL 03
088 *
089 15
090 *
091 RCL 10
092 +
093 STO 11
094 RCL 09
095 351
096 *
097 +
098 RCL 08
099 540
100 *
101 -
102 RCL 07
103 65
104 *
105 +
106 432
107 /
108 RCL 02
109 +
110 RCL 01
111 RCL 03
112 6
113 /
114 -
115 XEQ IND 00
116 RCL 03
117 *
118 72
119 *
120 RCL 11
121 -
122 RCL 09
123 9
124 *
125 +
126 RCL 07
127 ST+ X
128 -
129 3
130 1/X
131 %
132 ST- 05
if you replace this line by ABS ST+ 05
you'll get a less optimistic error estimate ( in absolute value )
133 RCL 11
134 RCL 09
135 81
136 *
137 +
138 PI
139 R-D
140 /
141 RCL 07
142 9
143 /
144 +
145 ST+ 02
146 DSE 06
147 GTO 01
a three-byte GTO
148 RCL 02
149 RCL 01
150 END
( 197 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
>>> Don't forget to clear register R05: CLX STO 05
"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 132 )
-The error is actually -1.8 10-7
-Press R/S to continue with the same h and N.
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 R07 are to be initialized before executing RKFB ( clear registers R06 & R07 before the first execution )
• R00: f name • R01
= x0 • R04 = h = stepsize
• R06 = y-error estimate
R08 to R18: 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.
001 LBL "RKFB"
002 RCL 05
003 STO 18
004 GTO 01
005 LBL 00
006 XEQ IND 00
007 RCL 04
008 ST* Z
009 *
010 RTN
011 LBL 01
012 RCL 03
013 RCL 02
014 RCL 01
015 XEQ 00
016 STO 13
017 ST+ X
018 9
019 /
020 RCL 03
021 +
022 X<>Y
023 STO 08
024 ST+ X
025 9
026 /
027 RCL 02
028 +
029 RCL 04
030 ST+ X
031 9
032 /
033 RCL 01
034 +
035 XEQ 00
036 STO 14
037 4
038 /
039 RCL 13
040 12
041 /
042 +
043 RCL 03
044 +
045 X<>Y
046 STO 09
047 4
048 /
049 RCL 08
050 12
051 /
052 +
053 RCL 02
054 +
055 RCL 04
056 3
057 /
058 RCL 01
059 +
060 XEQ 00
061 STO 15
062 270
063 *
064 RCL 14
065 243
066 *
067 -
068 RCL 13
069 69
070 *
071 +
072 128
073 /
074 RCL 03
075 +
076 X<>Y
077 STO 10
078 270
079 *
080 RCL 09
081 243
082 *
083 -
084 RCL 08
085 69
086 *
087 +
088 128
089 /
090 RCL 02
091 +
092 RCL 04
093 .75
094 *
095 RCL 01
096 +
097 XEQ 00
098 64
099 ST* Z
100 *
101 STO 16
102 RCL 15
103 324
104 *
105 -
106 RCL 14
107 405
108 *
109 +
110 RCL 13
111 85
112 *
113 -
114 60
115 /
116 RCL 03
117 +
118 X<>Y
119 STO 11
120 RCL 10
121 324
122 *
123 -
124 RCL 09
125 405
126 *
127 +
128 RCL 08
129 85
130 *
131 -
132 60
133 /
134 RCL 02
135 +
136 RCL 04
137 RCL 01
138 +
139 STO 01
140 XEQ 00
141 15
142 ST* Z
143 *
144 RCL 16
145 +
146 STO 17
147 RCL 15
148 351
149 *
150 +
151 RCL 14
152 135
153 ST* 09
154 *
155 -
156 RCL 13
157 65
158 *
159 +
160 432
161 /
162 RCL 03
163 +
164 X<>Y
165 RCL 11
166 +
167 STO 12
168 RCL 10
169 351
170 *
171 +
172 RCL 09
173 -
174 RCL 08
175 65
176 *
177 +
178 432
179 /
180 RCL 02
181 +
182 RCL 04
183 6
184 /
185 RCL 01
186 X<>Y
187 -
188 XEQ 00
189 72
190 ST* Z
191 *
192 RCL 17
193 -
194 RCL 15
195 9
196 *
197 +
198 RCL 13
199 ST+ X
200 -
201 3
202 1/X
203 %
204 ST- 07
if you replace this line by ABS ST+ 07
you'll get a less optimistic error estimate ( in absolute value )
205 R^
206 RCL 12
207 -
208 RCL 10
209 9
210 *
211 +
212 RCL 08
213 ST+ X
214 -
215 3
216 1/X
217 %
218 ST- 06
if you replace this line by ABS ST+ 06
you'll get a less optimistic error estimate ( in absolute value )
219 RCL 17
220 RCL 15
221 81
222 ST* 10
223 *
224 +
225 PI
226 R-D
227 /
228 RCL 13
229 9
230 ST/ 08
231 /
232 +
233 ST+ 03
234 RCL 12
235 RCL 10
236 +
237 PI
238 R-D
239 /
240 RCL 08
241 +
242 ST+ 02
243 DSE 18
244 GTO 01
a three-byte GTO
245 RCL 03
246 RCL 02
247 RCL 01
248 END
( 342 bytes / SIZE 019)
| 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
>>> Don't forget to clear registers R06 & R07: CLX STO 06 STO 07
"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.
3°) A 6th order method
A) First order differential
equations: dy/dx = f (x;y)
-The formula duplicates the Taylor serie through 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.
001 LBL "RK6"
002 RCL 04
003 STO 05
004 LBL 01
005 RCL 02
006 RCL 01
007 XEQ IND 00
008 RCL 03
009 *
010 STO 06
011 3
012 /
013 RCL 02
014 +
015 RCL 03
016 3
017 /
018 RCL 01
019 +
020 XEQ IND 00
021 RCL 03
022 *
023 STO 07
024 1.5
025 /
026 RCL 02
027 +
028 RCL 03
029 1.5
030 /
031 RCL 01
032 +
033 XEQ IND 00
034 RCL 03
035 *
036 STO 08
037 CHS
038 RCL 07
039 4
040 *
041 +
042 RCL 06
043 +
044 12
045 /
046 RCL 02
047 +
048 RCL 03
049 3
050 /
051 RCL 01
052 +
053 XEQ IND 00
054 RCL 03
055 *
056 STO 09
057 90
058 *
059 RCL 08
060 35
061 *
062 +
063 RCL 07
064 110
065 *
066 -
067 RCL 06
068 25
069 *
070 +
071 48
072 /
073 RCL 02
074 +
075 RCL 03
076 1.2
077 /
078 RCL 01
079 +
080 XEQ IND 00
081 RCL 03
082 *
083 STO 10
084 10
085 /
086 RCL 09
087 2
088 /
089 +
090 RCL 08
091 8
092 /
093 -
094 RCL 07
095 55
096 *
097 RCL 06
098 18
099 *
100 -
101 120
102 /
103 -
104 RCL 02
105 +
106 RCL 03
107 6
108 /
109 RCL 01
110 +
111 XEQ IND 00
112 RCL 03
113 ST+ 01
114 *
115 STO 11
116 80
117 *
118 RCL 09
119 118
120 *
121 -
122 RCL 07
123 99
124 *
125 +
126 39
127 /
128 RCL 10
129 128
130 *
131 RCL 08
132 215
133 *
134 +
135 RCL 06
136 783
137 *
138 -
139 780
140 /
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 RCL 10
154 +
155 32
156 *
157 +
158 RCL 09
159 RCL 08
160 +
161 55
162 *
163 +
164 .5
165 %
166 ST+ 02
167 DSE 05
168 GTO 01
a three-byte GTO
169 RCL 02
170 RCL 01
171 END
( 226 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 89s )
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 R18: 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.
001 LBL "RK6B"
002 RCL 05
003 STO 18
004 GTO 01
005 LBL 00
006 XEQ IND 00
007 RCL 04
008 ST* Z
009 *
010 RTN
011 LBL 01
012 RCL 03
013 RCL 02
014 RCL 01
015 XEQ 00
016 STO 12
017 3
018 /
019 RCL 03
020 +
021 X<>Y
022 STO 06
023 3
024 /
025 RCL 02
026 +
027 RCL 04
028 3
029 /
030 RCL 01
031 +
032 XEQ 00
033 STO 13
034 1.5
035 /
036 RCL 03
037 +
038 X<>Y
039 STO 07
040 1.5
041 /
042 RCL 02
043 +
044 RCL 04
045 1.5
046 /
047 RCL 01
048 +
049 XEQ 00
050 STO 14
051 CHS
052 RCL 13
053 4
054 *
055 +
056 RCL 12
057 +
058 12
059 /
060 RCL 03
061 +
062 X<>Y
063 STO 08
064 CHS
065 RCL 07
066 4
067 *
068 +
069 RCL 06
070 +
071 12
072 /
073 RCL 02
074 +
075 RCL 04
076 3
077 /
078 RCL 01
079 +
080 XEQ 00
081 STO 15
082 90
083 *
084 RCL 14
085 35
086 *
087 +
088 RCL 13
089 110
090 *
091 -
092 RCL 12
093 25
094 *
095 +
096 48
097 /
098 RCL 03
099 +
100 X<>Y
101 STO 09
102 90
103 *
104 RCL 08
105 35
106 *
107 +
108 RCL 07
109 110
110 *
111 -
112 RCL 06
113 25
114 *
115 +
116 48
117 /
118 RCL 02
119 +
120 RCL 04
121 1.2
122 /
123 RCL 01
124 +
125 XEQ 00
126 STO 16
127 10
128 /
129 RCL 15
130 2
131 /
132 +
133 RCL 14
134 8
135 /
136 -
137 RCL 13
138 11
139 *
140 24
141 /
142 -
143 RCL 12
144 .15
145 *
146 +
147 RCL 03
148 +
149 X<>Y
150 STO 10
151 10
152 /
153 RCL 09
154 2
155 /
156 +
157 RCL 08
158 8
159 /
160 -
161 RCL 07
162 11
163 *
164 24
165 /
166 -
167 RCL 06
168 .15
169 *
170 +
171 RCL 02
172 +
173 RCL 04
174 6
175 /
176 RCL 01
177 +
178 XEQ 00
179 STO 17
180 80
181 *
182 RCL 15
183 118
184 *
185 -
186 RCL 13
187 99
188 *
189 +
190 20
191 *
192 RCL 16
193 128
194 *
195 +
196 RCL 14
197 215
198 *
199 +
200 RCL 12
201 783
202 *
203 -
204 780
205 /
206 RCL 03
207 +
208 X<>Y
209 STO 11
210 80
211 *
212 RCL 09
213 118
214 *
215 -
216 RCL 07
217 99
218 *
219 +
220 20
221 *
222 RCL 10
223 128
224 *
225 +
226 RCL 08
227 215
228 *
229 +
230 RCL 06
231 783
232 *
233 -
234 780
235 /
236 RCL 02
237 +
238 RCL 04
239 RCL 01
240 +
241 STO 01
242 XEQ 00
243 RCL 12
244 +
245 13
246 *
247 RCL 17
248 RCL 16
249 +
250 32
251 *
252 +
253 RCL 15
254 RCL 14
255 +
256 55
257 *
258 +
259 .5
260 %
261 ST+ 03
262 R^
263 RCL 06
264 +
265 13
266 *
267 RCL 11
268 RCL 10
269 +
270 32
271 *
272 +
273 RCL 09
274 RCL 08
275 +
276 55
277 *
278 +
279 .5
280 %
281 ST+ 02
282 DSE 18
283 GTO 01 a
three-byte GTO
284 RCL 03
285 RCL 02
286 RCL 01
287 END
( 390 bytes / SIZE 019 )
| 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 2mn44s )
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 )
4°) A 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 serie through 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.
001 LBL "RK8"
002 RCL 04
003 STO 52
004 GTO 01
005 LBL 00
006 XEQ IND 00
007 RCL 03
008 *
009 RTN
010 LBL 01
011 RCL 02
012 RCL 01
013 XEQ 00
014 STO 05
015 RCL 12
016 *
017 RCL 02
018 +
019 RCL 03
020 RCL 12
021 *
022 RCL 01
023 +
024 XEQ 00
025 STO 06
026 RCL 05
027 +
028 RCL 16
029 *
030 RCL 02
031 +
032 RCL 03
033 RCL 12
034 *
035 RCL 01
036 +
037 XEQ 00
038 STO 07
039 RCL 17
040 *
041 RCL 06
042 RCL 18
043 *
044 +
045 RCL 05
046 RCL 19
047 /
048 +
049 RCL 02
050 +
051 RCL 03
052 RCL 11
053 *
054 RCL 01
055 +
056 XEQ 00
057 STO 08
058 RCL 20
059 *
060 RCL 07
061 RCL 21
062 *
063 +
064 RCL 05
065 RCL 22
066 *
067 +
068 RCL 02
069 +
070 RCL 03
071 RCL 11
072 *
073 RCL 01
074 +
075 XEQ 00
076 STO 09
077 RCL 23
078 *
079 RCL 08
080 RCL 24
081 *
082 +
083 RCL 07
084 RCL 25
085 *
086 +
087 RCL 05
088 RCL 26
089 *
090 +
091 RCL 02
092 +
093 RCL 03
094 RCL 12
095 *
096 RCL 01
097 +
098 XEQ 00
099 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.
001 LBL "RK8B"
002 RCL 05
003 STO 59
004 GTO 01
005 LBL 00
006 XEQ IND 00
007 RCL 04
008 ST* Z
009 *
010 RTN
011 LBL 01
012 RCL 03
013 RCL 02
014 RCL 01
015 XEQ 00
016 STO 12
017 RCL 19
018 *
019 RCL 03
020 +
021 X<>Y
022 STO 06
023 RCL 19
024 *
025 RCL 02
026 +
027 RCL 04
028 RCL 19
029 *
030 RCL 01
031 +
032 XEQ 00
033 STO 13
034 RCL 12
035 +
036 RCL 23
037 *
038 RCL 03
039 +
040 X<>Y
041 STO 07
042 RCL 06
043 +
044 RCL 23
045 *
046 RCL 02
047 +
048 RCL 04
049 RCL 19
050 *
051 RCL 01
052 +
053 XEQ 00
054 STO 14
055 RCL 24
056 *
057 RCL 13
058 RCL 25
059 *
060 +
061 RCL 12
062 RCL 26
063 /
064 +
065 RCL 03
066 +
067 X<>Y
068 STO 08
069 RCL 24
070 *
071 RCL 07
072 RCL 25
073 *
074 +
075 RCL 06
076 RCL 26
077 /
078 +
079 RCL 02
080 +
081 RCL 04
082 RCL 18
083 *
084 RCL 01
085 +
086 XEQ 00
087 STO 15
088 RCL 27
089 *
090 RCL 14
091 RCL 28
092 *
093 +
094 RCL 12
095 RCL 29
096 *
097 +
098 RCL 03
099 +
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 )
5°) A general Explicit-Runge-Kutta program
A) First order differential equations: dy/dx = f (x;y)
-All n stage explicit Runge-Kutta formulae 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 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
-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+13 to R(n2+7.n+22)/2 are to be initialized before executing "ERKB"
• R00: f name • R01
= x0
R07 to R12: temp
• R2n+13 to R(n2+7n+22)/2 = the ( n2
+ 3.n - 2 ) / 2 coefficients ai,j ; bi ; ci
• R02 = y0
R13 to Rn+12 = kiy
which are to be stored as shown below:
• R03 = z0
Rn+13 to R2n+12 = kiz
• R04 = h = stepsize
• R05 = N = number of steps
• R06 = n = number of stages
• R2n+13 = a2,1
• R2n+14 = b2
• R2n+15 = a3,1
• R2n+16 = a3,2
• R2n+17 = b3
......................................................................... ...................
• R(n2+3n+24)/2 = an,1 ................ • R(n2+5n+20)/2 = an,n-1 • R(n2+5n+22)/2 = bn
• R(n2+5n+24)/2 = c1 ............................................................
• R(n2+7n+22)/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.
001 LBL "ERKB"
002 RCL 05
003 STO 07
004 LBL 01
005 13.9
006 STO 08
007 RCL 06
008 +
009 STO 09
010 LASTX
011 +
012 STO 10
013 LASTX
014 12
015 +
016 .1
017 %
018 14
019 +
020 STO 11
021 RCL 06
022 +
023 STO 12
024 RCL 03
025 RCL 02
026 RCL 01
027 XEQ IND 00
028 RCL 04
029 ST* Z
030 *
031 STO IND 09
032 X<>Y
033 STO 13
034 LBL 02
035 RCL 08
036 INT
037 .1
038 %
039 13
040 +
041 STO 08
042 RCL 06
043 +
044 STO 09
045 CLST
046 LBL 03
047 X<>Y
048 RCL IND 08
049 RCL IND 10
050 *
051 +
052 X<>Y
053 RCL IND 09
054 RCL IND 10
055 *
056 +
057 ISG 10
058 ISG 09
059 CLX
060 ISG 08
061 GTO 03
062 RCL 03
063 +
064 X<>Y
065 RCL 02
066 +
067 RCL 04
068 RCL IND 10
069 *
070 RCL 01
071 +
072 XEQ IND 00
073 RCL 04
074 ST* Z
075 *
076 STO IND 12
077 X<>Y
078 STO IND 11
079 ISG 12
080 CLX
081 ISG 10
082 ISG 11
083 GTO 02
084 RCL 06
085 ST- 11
086 ST- 12
087 CLST
088 LBL 04
089 X<>Y
090 RCL IND 10
091 RCL IND 11
092 *
093 +
094 X<>Y
095 RCL IND 10
096 RCL IND 12
097 *
098 +
099 ISG 12
100 CLX
101 ISG 10
102 ISG 11
103 GTO 04
104 ST+ 03
105 X<>Y
106 ST+ 02
107 RCL 04
108 ST+ 01
109 DSE 07
110 GTO 01
a three-byte GTO
111 RCL 03
112 RCL 02
113 RCL 01
114 END
( 167 bytes / SIZE ( n2+ 7.n + 24 )/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
R21
R22
0 1/2
1/2 are to be stored into
R23 R24 R25
respectively
0 0 1
1
R26 R27 R28 R29
1/6 1/3 1/3 1/6
R30 R31 R32 R33
-With the 7 stage method used in "RK6B": 7 STO 06 and the 34 coefficients:
1/3
1/3
R27
R28
0 2/3
2/3
R39 R30
R31
1/12 1/3 -1/12
1/3
R32 R33 R34
R35
25/48 -55/24 35/48
15/8
5/6 are to be stored into
R36 R37 R38 R39
R40
3/20 -11/24 -1/8
1/2 1/10
1/6
R41 R42 R43 R44 R45
R46
-261/260 33/13 43/156 -118/39 32/195
80/39 1
R47 R48 R49 R50 R51 R52 R53
13/200 0 11/40
11/40 4/25 4/25
13/200
R54 R55 R56 R57 R58 R59 R60
................ 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 )
6°) An implicit Runge-Kutta method of order 8
A) First order differential
equations: dy/dx = f (x;y)
-The formulae used by "RK4" to "RK8" are explicit Runge-Kutta
methods in that the successive k1 , k2 , .... , kn
may be directly computed.
-In the following program , we have k1 = h.f ( x +
b1h ; y + a1,1 k1 + a1,2
k2 + ..... + a1,5 k5 )
.....................................................................................
where ai,j # 0 for j > i-1
k5 = h.f ( x + b5h ; y + a5,1 k1
+ a5,2 k2 + ..... + a5,5 k5
)
and ci = a5,i
-The advantage is that only 5 stages are required to obtain an 8th order
method: the program is shorter and less data registers are needed.
-But we have to solve a 5x5 non-linear system within each step! Speed
is not a characteristic of implicit methods, especially with an HP-41...
-"IRK8" uses an iterative algorithm based on the "fixed point theorem"
-Some implicit Runge-Kutta formulae are also satisfactory for stiff
problems, but the elementary iterative method used by "IRK8" will seldom
lead to convergence
in such a case, unless h is very small: another algorithm
( like the Newton method ) would be better to solve the 5x5 non-linear
system, if you're in the mood...
Data Registers: Registers R00 to R04 and R16 to R45 are to be initialized before executing "IRK8"
• R00: f name • R01
= x0 • R03 = h = stepsize
R05 to R10: temp ; R11 = k1 ......
R15 = k5
• R02 = y0 •
R04 = N = number of steps •
R16 thru R45 = the 30 following numbers:
R16 = 1
R26 = 29/180
R36 = 29/180
R17 = 1/20
R27 = -3/140
R37 = -0.06901154103
the decimal numbers are actually
R18 = 49/180
R28 = 1/2
R38 = 0.05200216599
rational function of sqrt(21)
R19 = 16/45
R29 = 1/20
R39 = -3/140
( all decimals are correct )
R20 = 49/180
R30 = 0.2813091833
R40 = 0
R21 = 1/20
R31 = 73/360
R41 = 1/20
R22 = 0.8273268354
R32 = -0.05283696110
R42 = -7/60
R23 = 1/20
R33 = 3/160
R43 = 2/15
R24 = 0.2702200562
R34 = 0.1726731646
R44 = -7/60
R25 = 0.3674242394
R35 = 1/20
R45 = 1/20
Flags: /
Subroutine: a program calculating f(x;y)
assuming x is in X-register and y is in Y-register upon entry.
01 LBL "IRK8"
02 RCL 04
03 STO 09
04 CLX
05 STO 11
06 STO 12
07 STO 13
08 STO 14
09 STO 15
10 LBL 01
11 11.015
12 STO 05
13 45
14 STO 07
15 CLX
16 STO 08
17 LBL 02
18 15.01
19 STO 06
20 CLX
21 LBL 03
22 RCL IND 06
23 RCL IND 07
24 *
25 +
26 DSE 07
27 DSE 06
28 GTO 03
29 RCL 02
30 +
31 STO 10
32 RCL IND 07
33 RCL 03
34 *
35 RCL 01
36 +
37 XEQ IND 00
38 RCL 03
39 *
40 ENTER^
41 X<> IND 05
42 -
43 ABS
44 ST+ 08
45 DSE 07
46 ISG 05
47 GTO 02
48 RCL 08
49 VIEW X
50 E-9
or another "small" number:
51 X<Y?
52 GTO 01
53 RCL 03
54 ST+ 01
55 RCL 10
56 STO 02
57 DSE 09
58 GTO 01
59 RCL 01
60 CLD
61 END
( 99 bytes / SIZE 046 )
| STACK | INPUTS | OUTPUTS |
| Y | / | y(x0+N.h) |
| X | / | x0+N.h |
Example: y(0) = 1 and dy/dx = -2 x.y Evaluate y(0.5)
-Initialize registers R16 thru R45.
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
5 STO 04 XEQ "IRK8"
yields x = 0.5
( in 5mn57 )
X<>Y
y(0.5) = 0.7788007830 ( exact result
is 0.7788007831 )
-If you replace lines 50-51 with X#0? , you'll get the exact value 0.7788007831 but this severe test might sometimes lead to an infinite loop!
-The HP-41 displays | k1 - k'1 | + ............
+ | k5 - k'5 | where ki
and k'i are 2 successive k-approximations:
-At each step, this sum must tend to zero. Otherwise, reduce the stepsize
h and start over!
-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)
-Now, we have to solve a 10x10 non-linear system within each step.
Data Registers: Registers R00 to R05 and R25 to R54 are to be initialized before executing "IRK8B"
• R00: f name • R01
= x0 • R04 = h = stepsize
R06 to R14: temp ; R15 , ..... , R24 = ki
• R02 = y0 •
R05 = N = number of steps •
R25 thru R54 = the same 30 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 "IRK8B"
02 RCL 05
03 STO 12
04 15.024
05 CLRGX
if you don't have an HP-41CX , replace lines 04-05 with CLX
STO 15 STO 16 STO 17 STO 18 STO 19
06 LBL 01
STO 20 STO 21 STO 22 STO 23 STO 24
07 15.019
08 STO 06
09 20
10 STO 07
11 54
12 STO 10
13 CLX
14 STO 11
15 LBL 02
16 19.014
17 STO 08
18 24
19 STO 09
20 CLST
21 LBL 03
22 X<>Y
23 RCL IND 08
24 RCL IND 10
25 *
26 +
27 X<>Y
28 RCL IND 09
29 RCL IND 10
30 *
31 +
32 DSE 10
33 DSE 09
34 DSE 08
35 GTO 03
36 RCL 03
37 +
38 STO 14
39 X<>Y
40 RCL 02
41 +
42 STO 13
43 RCL IND 10
44 RCL 04
45 *
46 RCL 01
47 +
48 XEQ IND 00
49 RCL 04
50 ST* Z
51 *
52 ENTER^
53 X<> IND 07
54 -
55 ABS
56 X<>Y
57 ENTER^
58 X<> IND 06
59 -
60 ABS
61 +
62 ST+ 11
63 SIGN
64 ST+ 07
65 ST- 10
66 ISG 06
67 GTO 02
68 RCL 11
69 VIEW X
70 E-9
or another "small" number
71 X<Y?
72 GTO 01
73 RCL 04
74 ST+ 01
75 RCL 14
76 STO 03
77 RCL 13
78 STO 02
79 DSE 12
80 GTO 01 a three-byte
GTO
81 RCL 01
82 CLD
83 END
( 138 bytes / SIZE 055 )
| 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(0.5) and z(0.5)
-Initialize registers R25 thru R54.
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
5 STO 05 XEQ "IRK8B"
yields x = 0.5
execution time = 12mn
RDN
y(0.5) = 0.7788007830
( exact result is 0.7788007831 )
RDN
z(0.5) = -0.7788007830 ( exact
result is -0.7788007831 )
-The HP-41 displays | k1 - k'1 | + ............
+ | k10 - k'10 | where ki
and k'i are 2 successive k-approximations:
-At each step, this sum must tend to zero. Otherwise, reduce the stepsize
h and start over!
-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 )
Reference: J.C. Butcher - "The Numerical Analysis of Ordinary Differential Equations" - John Wiley & Sons - ISBN 0-471-91046-5
Go back to the HP-41 software library
Go back to the general software library
Go
back to the main exhibit hall