The Museum of HP Calculators
The Museum of HP Calculators
This program is Copyright © 2007 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.
1°) Systems of first-order Differential Equations
a) The "Classical" 4th-order Runge-Kutta
Method
b) A Runge-Kutta Method of order 6 (
new )
c) A Runge-Kutta
Method of order 8 ( new )
d) Runge-Kutta-Fehlberg
4-5 Method ( new )
e) Bulirsch-Stoer
Method ( new )
2°) Systems of second-order Conservative Equations
a) Systems of 2
Equations ( new )
b) Systems of 3
Equations ( new )
c) Systems of n
Equations ( new )
-The programs listed in paragraph 1°) solve a system of n first-order differential equations:
dy1/dx = f1(x,y1,
... ,yn)
dy2/dx = f2(x,y1,
... ,yn)
..................
with the initial values: yi(x0)
i = 1 , ... , n
dyn/dx = fn(x,y1, ... ,yn)
-In paragraph 2°) the HP-41 solves systems of second-order equations where the first derivatives y'i(x) do not appear explicitly:
d2y1/dx2 =
f1(x,y1, ... ,yn)
d2y2/dx2
= f2(x,y1, ... ,yn)
..................
with the initial values: yi(x0)
, y'i(x0)
i = 1 , ... , n
d2yn/dx2 = fn(x,y1, ... ,yn)
-All these programs use 2 instructions from the X-functions module:
REGMOVE and REGSWAP - Except in §2°)a)b)
1°) Systems of first-order Differential Equations
a) The "Classical" 4th-order
Runge-Kutta Method
Data Registers: • R00 = subroutine name ( Registers R00 thru R03 and R10 thru R10+n are to be initialized before executing "RK4" )
• R01 = n = number of equations = number of functions
R04 thru R09: temp
• R02 = h = stepsize
• R03 = N = number of steps
• R10 = x0
• R11 = y1(x0)
R11+n thru R10+4n: temp
• R12 = y2(x0)
............... ( initial
values )
• R10+n = yn(x0)
Flags: F06-F07
Subroutine: A program which calculates
and stores f1 ( x, y1 , ... , yn)
, ... , fn ( x, y1 , ... , yn) in
registers R11+n, ... , R10+2n respectively
with x, y1, ... , yn in regiters R10, R11,
... , R10+n
If you don't have an HP-41 CX, line 24 can be replaced by the 5 lines:
0
LBL 00
STO IND Y
ISG Y
GTO 00
01 LBL "RK4"
02 RCL 03
03 STO 09
04 RCL 01
05 RCL 01
06 RCL 01
07 10.01
08 +
09 STO 04
10 INT
11 +
12 STO 05
13 +
14 STO 06
15 +
16 STO 07
17 RCL 05
18 RCL 06
19 1
20 +
21 E3
22 /
23 +
24 CLRGX
25 11
26 LASTX
27 +
28 RCL 01
29 E6
30 /
31 +
32 STO 08
33 GTO 03
34 LBL 01
35 XEQ IND 00
36 LBL 02
37 RCL IND 05
38 RCL 02
39 *
40 ST+ IND 06
41 FS? 06
42 ST+ IND 06
43 FC? 07
44 ST+ X
45 RCL IND 07
46 +
47 STO IND 04
48 DSE 07
49 DSE 06
50 DSE 05
51 DSE 04
52 GTO 02
53 RCL 01
54 ST+ 04
55 ST+ 05
56 ST+ 06
57 ST+ 07
58 RTN
59 LBL 03
60 RCL 08
61 REGMOVE
62 CF 06
63 SF 07
64 XEQ 01
65 SF 06
66 RCL 02
67 ST+ 10
68 XEQ 01
69 CF 07
70 XEQ 01
71 CF 06
72 RCL 02
73 ST+ 10
74 XEQ 01
75 RCL 08
76 REGSWAP
77 RCL 04
78 RCL 06
79 3
80 SIGN
81 LBL 04
82 CLX
83 X<> IND Y
84 LASTX
85 /
86 ST+ IND Z
87 DSE Y
88 DSE Z
89 GTO 04
90 DSE 09
91 GTO 03
92 RCL 13
93 RCL 12
94 RCL 11
95 RCL 10
96 END
( 155 bytes / SIZE 4n+11 )
| STACK | INPUTS | OUTPUTS |
| T | / | y3(x0+N.h) |
| Z | / | y2(x0+N.h) |
| Y | / | y1(x0+N.h) |
| X | / | x0+N.h |
Example: Let's solve the system
dy1/dx = - y1y2y3
dy2/dx = x ( y1+y2-y3
)
dy3/dx = xy1- y2y3
with the initial values: y1(0) = 1 ; y2(0) = 1 ; y3(0) = 2
x ; y1 ; y2 ; y3 will be in R10
; R11 ; R12 ; R13 respectively, and we have to write a program to compute
the 3 functions
and store the results in R14 ; R15 ; R16 ( just after the
"initial-value registers" ).
For instance, the following subroutine will do the job - let's call
it "T"
01 LBL "T"
02 RCL 11
03 RCL 12
04 RCL 13
05 *
06 *
07 CHS
08 STO 14
09 RCL 11
10 RCL 12
11 +
12 RCL 13
13 -
14 RCL 10
15 *
16 STO 15
17 RCL 10
18 RCL 11
19 *
20 RCL 12
21 RCL 13
22 *
23 -
24 STO16
25 RTN
We store the name of the subroutine in R00 "T"
ASTO 00
the numbers of equations in R01
3 STO 01
then the initial values: x
in R10
0 STO 10
y1 in R11
1 STO 11
y2 in R12
1 STO 12
y3 in R13
2 STO 13
Suppose we want to know y1(1) , y2(1) , y3(1) with a step size h = 0.1 and therefore N = 10
h/2 is stored in R02
0.05 STO 02
N is stored in R03
10 STO 03 and XEQ "RK4"
about 2mn24s later, the HP-41 gives:
x = 1
in X and R10
y1(1) = 0.2582093855 in Y and R11
y2(1) = 1.157619553 in Z and
R12
y3(1) = 0.8421786516 in T and R13
-To continue with the same h and N, simply press R/S and you'll have y1(2) y2(2) y3(2)
-An estimation of the accuracy may be obtained by doing the calculation again with a smaller step size like h=0.05
It yields: y1(1) = 0.2582079989 y2(1) = 1.157623732 y3(1) = 0.8421783424
Theoretically, the results tend to the exact solution as h tends to
zero,
but naturally, roundoff errors ( and execution time ) will limit the
accuracy attainable!
b) A Runge-Kutta Method of
Order 6
Data Registers: • R00 = subroutine name ( Registers R00 thru R03 and R20 thru R20+n are to be initialized before executing "RK6" )
• R01 = n = number of equations = number of functions
R04 thru R13: temp R14 thru R19:
unused
• R02 = h = stepsize
• R03 = N = number of steps
• R20 = x0
• R21 = y1(x0)
R21+n thru R20+8n: temp
• R22 = y2(x0)
............... ( initial
values )
• R20+n = yn(x0)
Flags: /
Subroutine: A program which calculates
and stores f1 ( x, y1 , ... , yn)
, ... , fn ( x, y1 , ... , yn) in
registers R21+n, ... , R20+2n respectively
with x, y1, ... , yn in regiters R20, R21,
... , R20+n
01 LBL "RK6"
02 RCL 03
03 STO 04
04 21.021
05 RCL 01
06 E6
07 /
08 +
09 RCL 01
10 7
11 *
12 +
13 STO 05
14 RCL 01
15 RCL 01
16 RCL 01
17 20.02
18 +
19 STO 06
20 +
21 STO 07
22 +
23 STO 08
24 +
25 STO 09
26 +
27 STO 10
28 +
29 STO 11
30 +
31 STO 12
32 LBL 01
33 RCL 20
34 STO 13
35 RCL 05
36 REGSWAP
37 REGMOVE
38 XEQ IND 00
39 LBL 02
40 RCL IND 07
41 RCL 02
42 *
43 STO IND 08
44 3
45 /
46 ST+ IND 06
47 DSE 08
48 DSE 07
49 DSE 06
50 GTO 02
51 RCL 01
52 ST+ 06
53 ST+ 07
54 ST+ 08
55 RCL 02
56 LASTX
57 /
58 ST+ 20
59 XEQ IND 00
60 RCL 05
61 REGMOVE
62 LBL 03
63 RCL IND 07
64 RCL 02
65 *
66 STO IND 09
67 1.5
68 /
69 ST+ IND 06
70 DSE 09
71 DSE 07
72 DSE 06
73 GTO 03
74 RCL 01
75 ST+ 06
76 ST+ 07
77 ST+ 09
78 RCL 02
79 3
80 /
81 ST+ 20
82 XEQ IND 00
83 RCL 05
84 REGMOVE
85 LBL 04
86 RCL IND 07
87 RCL 02
88 *
89 STO IND 10
90 RCL IND 09
91 4
92 *
93 -
94 RCL IND 08
95 -
96 12
97 /
98 ST- IND 06
99 DSE 10
100 DSE 09
101 DSE 08
102 DSE 07
103 DSE 06
104 GTO 04
105 RCL 01
106 ST+ 06
107 ST+ 07
108 ST+ 08
109 ST+ 09
110 ST+ 10
111 RCL 02
112 3
113 /
114 ST- 20
115 XEQ IND 00
116 RCL 05
117 REGMOVE
118 LBL 05
119 RCL IND 07
120 RCL 02
121 *
122 STO IND 11
123 18
124 *
125 RCL IND 10
126 7
127 *
128 +
129 RCL IND 09
130 22
131 *
132 -
133 RCL IND 08
134 5
135 *
136 +
137 9.6
138 /
139 ST+ IND 06
140 DSE 11
141 DSE 10
142 DSE 09
143 DSE 08
144 DSE 07
145 DSE 06
146 GTO 05
147 RCL 01
148 ST+ 06
149 ST+ 07
150 ST+ 08
151 ST+ 09
152 ST+ 10
153 ST+ 11
154 RCL 02
155 2
156 /
157 ST+ 20
158 XEQ IND 00
159 RCL 05
160 REGMOVE
161 LBL 06
162 RCL IND 07
163 RCL 02
164 *
165 STO IND 12
166 5
167 /
168 RCL IND 11
169 +
170 RCL IND 10
171 4
172 /
173 -
174 RCL IND 08
175 18
176 *
177 RCL IND 09
178 55
179 *
180 -
181 60
182 /
183 +
184 2
185 /
186 ST+ IND 06
187 DSE 12
188 DSE 11
189 DSE 10
190 DSE 09
191 DSE 08
192 DSE 07
193 DSE 06
194 GTO 06
195 RCL 01
196 ST+ 06
197 ST+ 07
198 ST+ 08
199 ST+ 09
200 ST+ 10
201 ST+ 11
202 ST+ 12
203 RCL 02
204 6
205 /
206 RCL 13
207 +
208 STO 20
209 XEQ IND 00
210 RCL 05
211 REGMOVE
212 LBL 07
213 RCL IND 07
214 RCL 02
215 *
216 80
217 X<>Y
218 ST* Y
219 X<> IND 09
220 99
221 *
222 +
223 RCL IND 11
224 118
225 *
226 -
227 20
228 *
229 RCL IND 12
230 ST+ IND 09
231 128
232 *
233 +
234 RCL IND 10
235 ST+ IND 11
236 215
237 *
238 +
239 RCL IND 08
240 783
241 *
242 -
243 780
244 /
245 ST+ IND 06
246 DSE 12
247 DSE 11
248 DSE 10
249 DSE 09
250 DSE 08
251 DSE 07
252 DSE 06
253 GTO 07
254 RCL 01
255 ST+ 06
256 ST+ 07
257 ST+ 08
258 ST+ 09
259 ST+ 10
260 ST+ 11
261 ST+ 12
262 RCL 02
263 RCL 13
264 +
265 STO 20
266 XEQ IND 00
267 RCL 05
268 REGMOVE
269 LBL 08
270 RCL IND 07
271 RCL 02
272 *
273 RCL IND 08
274 +
275 13
276 *
277 RCL IND 09
278 32
279 *
280 +
281 RCL IND 11
282 55
283 *
284 +
285 .5
286 %
287 ST+ IND 06
288 DSE 11
289 DSE 09
290 DSE 08
291 DSE 07
292 DSE 06
293 GTO 08
294 RCL 01
295 ST+ 06
296 ST+ 07
297 ST+ 08
298 ST+ 09
299 ST+ 11
300 DSE 04
301 GTO 01
a three-byte GTO ( or replace this line by GTO 15 and
line 32 by LBL 15 )
302 RCL 23
303 RCL 22
304 RCL 21
305 RCL 20
306 END
( 498 bytes / SIZE 21+8n )
| STACK | INPUTS | OUTPUTS |
| T | / | y3(x0+N.h) |
| Z | / | y2(x0+N.h) |
| Y | / | y1(x0+N.h) |
| X | / | x0+N.h |
Example:
dy/dx = -y.z.u
y(0) = 1 Evaluate
y(1) , z(1) , u(1) with h = 0.1 ( whence
N = 10 )
dz/dx = x.( y + z - u )
z(0) = 1
du/dx = x.y - z.u
u(0) = 2
LBL "T" *
RCL 21 -
RCL 20 RCL 23
RTN "T"
ASTO 00
RCL 21 *
RCL 22 RCL 20
RCL 21 *
RCL 22 CHS
+
*
*
-
RCL 23 STO 24
RCL 23 STO 25
RCL 22 STO 26
Initial values: 0 STO
20
1 STO 21 STO 22
2 STO 23
n = 3 STO 01 (
3 functions )
h = 0.1 STO 02
N = 10 STO 03 XEQ "RK6" >>>>
x = 1
= R20
( in 6mn28s )
RDN y(1) = 0.258207889 = R21
y(1) = 0.258207906
RDN z(1) = 1.157623948 = R22
the exact results, rounded to 9D are
z(1) = 1.157623981
RDN u(1) = 0.842178329 = R23
u(1) = 0.842178312
-To continue the calculations, simply press R/S ( after changing
h & N in registers R02 & R03 if needed )
c) A Runge-Kutta Method of
Order 8
Data Registers: • R00 = subroutine name ( Registers R00 to R03 , R14 to R48 and R50 to R50+n are to be initialized before executing "RK8" )
• R01 = n = number of equations = number of functions
R04 thru R13: temp R49: temp
• R02 = h = stepsize
• R03 = N = number of steps
• R14 thru R48 = Constants
• R50 = x0
• R51 = y1(x0)
R51+n thru R50+9n: temp
• R52 = y2(x0)
............... ( initial
values )
• R50+n = yn(x0)
•
R14 = 0.8961811934
• R23 = 0.1996369936
• R32 = 0.9576053440
• R41 = 0.9401094520
•
R15 = -0.2117115009
• R24 = 0.09790045616
• R33 = -0.6325461607
• R42 = -2.229158210
•
R16 = 0.8273268354
• R25 = 0.3285172131
• R34 = 0.1527777778
• R43 = 7.553840442
•
R17 = 0.06514850915
• R26 = -0.3497052863
• R35 = -0.009086961101
• R44 = -7.164951553
•
R18 = 0.5766714727
• R27 = -0.03302551131
• R36 = 1.062498447
• R45 = 2.451380432
•
R19 = 0.1855068535
• R28 = 0.1289862930
• R37 = -1.810863083
• R46 = 0.05
•
R20 = 0.3865246267
• R29 = 0.1726731646
• R38 = 2.031083139
• R47 = 0.2722222222
•
R21 = -0.4634553896
• R30 = -0.01186868389
• R39 = -0.6379313502
• R48 = 0.3555555556
•
R22 = 0.3772937693
• R31 = 0.002002165993
• R40 = -0.5512205631
Flags: /
Subroutine: A program which calculates
and stores f1 ( x, y1 , ... , yn)
, ... , fn ( x, y1 , ... , yn) in
registers R51+n, ... , R50+2n respectively
with x, y1, ... , yn in regiters R50, R51,
... , R50+n
01 LBL "RK8"
02 51.051
03 RCL 01
04 E6
05 /
06 +
07 RCL 01
08 8
09 *
10 +
11 STO 05
12 RCL 01
13 RCL 01
14 RCL 01
15 50.05
16 +
17 STO 06
18 +
19 STO 07
20 +
21 STO 08
22 +
23 STO 09
24 +
25 STO 10
26 +
27 STO 11
28 +
29 STO 12
30 +
31 STO 13
32 RCL 03
33 STO 49
34 GTO 01
35 LBL 00
36 XEQ IND 00
37 RCL 05
38 REGMOVE
39 RTN
40 LBL 01
41 RCL 50
42 STO 04
43 RCL 05
44 REGSWAP
45 REGMOVE
46 XEQ IND 00
47 LBL 02
48 RCL IND 07
49 RCL 02
50 *
51 STO IND 08
52 2
53 /
54 ST+ IND 06
55 DSE 08
56 DSE 07
57 DSE 06
58 GTO 02
59 RCL 01
60 ST+ 06
61 ST+ 07
62 ST+ 08
63 RCL 02
64 LASTX
65 /
66 ST+ 50
67 XEQ 00
68 LBL 03
69 RCL IND 07
70 RCL 02
71 *
72 STO IND 09
73 RCL IND 08
74 +
75 4
76 /
77 ST+ IND 06
78 DSE 09
79 DSE 08
80 DSE 07
81 DSE 06
82 GTO 03
83 RCL 01
84 ST+ 06
85 ST+ 07
86 ST+ 08
87 ST+ 09
88 XEQ 00
89 LBL 04
90 RCL IND 07
91 RCL 02
92 *
93 STO IND 10
94 RCL 14
95 *
96 RCL IND 09
97 RCL 15
98 *
99 +
100 RCL IND 08
101 7
102 /
103 +
104 ST+ IND 06
105 DSE 10
106 DSE 09
107 DSE 08
108 DSE 07
109 DSE 06
110 GTO 04
111 RCL 01
112 ST+ 06
113 ST+ 07
114 ST+ 08
115 ST+ 09
116 ST+ 10
117 RCL 02
118 RCL 16
119 *
120 RCL 04
121 +
122 STO 50
123 XEQ 00
124 LBL 05
125 RCL IND 07
126 RCL 02
127 *
128 STO IND 11
129 RCL 17
130 *
131 RCL IND 10
132 RCL 18
133 *
134 +
135 RCL IND 08
136 RCL 19
137 *
138 +
139 ST+ IND 06
140 DSE 11
141 DSE 10
142 DSE 08
143 DSE 07
144 DSE 06
145 GTO 05
146 RCL 01
147 ST+ 06
148 ST+ 07
149 ST+ 08
150 ST+ 10
151 ST+ 11
152 XEQ 00
153 LBL 06
154 RCL IND 07
155 RCL 02
156 *
157 STO IND 12
158 RCL 20
159 *
160 RCL IND 11
161 RCL 21
162 *
163 +
164 RCL IND 10
165 RCL 22
166 *
167 +
168 RCL IND 08
169 RCL 23
170 *
171 +
172 ST+ IND 06
173 DSE 12
174 DSE 11
175 DSE 10
176 DSE 08
177 DSE 07
178 DSE 06
179 GTO 06
180 RCL 01
181 ST+ 06
182 ST+ 07
183 ST+ 08
184 ST+ 10
185 ST+ 11
186 ST+ 12
187 RCL 02
188 2
189 /
190 RCL 04
191 +
192 STO 50
193 XEQ 00
194 LBL 07
195 RCL IND 07
196 RCL 02
197 *
198 STO IND 13
199 RCL 24
200 *
201 RCL IND 12
202 RCL 25
203 *
204 +
205 RCL IND 11
206 RCL 26
207 *
208 +
209 RCL IND 10
210 RCL 27
211 *
212 +
213 RCL IND 08
214 RCL 28
215 *
216 +
217 ST+ IND 06
218 DSE 13
219 DSE 12
220 DSE 11
221 DSE 10
222 DSE 08
223 DSE 07
224 DSE 06
225 GTO 07
226 RCL 01
227 ST+ 06
228 ST+ 07
229 ST+ 08
230 ST+ 10
231 ST+ 11
232 ST+ 12
233 ST+ 13
234 RCL 02
235 RCL 29
236 *
237 RCL 04
238 +
239 STO 50
240 XEQ 00
241 LBL 08
242 RCL IND 07
243 RCL 02
244 *
245 STO IND 09
246 9
247 /
248 RCL IND 13
249 RCL 30
250 *
251 +
252 RCL IND 12
253 RCL 31
254 *
255 +
256 RCL IND 08
257 14
258 /
259 +
260 ST+ IND 06
261 DSE 13
262 DSE 12
263 DSE 09
264 DSE 08
265 DSE 07
266 DSE 06
267 GTO 08
268 RCL 01
269 ST+ 06
270 ST+ 07
271 ST+ 08
272 ST+ 09
273 ST+ 12
274 ST+ 13
275 XEQ 00
276 LBL 09
277 RCL IND 07
278 RCL 02
279 *
280 STO IND 10
281 RCL 32
282 *
283 RCL IND 09
284 RCL 33
285 *
286 +
287 RCL IND 13
288 RCL 34
289 *
290 +
291 RCL IND 12
292 RCL 35
293 *
294 +
295 RCL IND 08
296 32
297 /
298 +
299 ST+ IND 06
300 DSE 13
301 DSE 12
302 DSE 10
303 DSE 09
304 DSE 08
305 DSE 07
306 DSE 06
307 GTO 09
308 RCL 01
309 ST+ 06
310 ST+ 07
311 ST+ 08
312 ST+ 09
313 ST+ 10
314 ST+ 12
315 ST+ 13
316 RCL 02
317 2
318 /
319 RCL 04
320 +
321 STO 50
322 XEQ 00
323 LBL 10
324 RCL IND 07
325 RCL 02
326 *
327 STO IND 11
328 RCL 36
329 *
330 RCL IND 10
331 RCL 37
332 *
333 +
334 RCL IND 09
335 RCL 38
336 *
337 +
338 RCL IND 13
339 RCL 39
340 *
341 +
342 RCL IND 12
343 9
344 /
345 +
346 RCL IND 08
347 14
348 /
349 +
350 ST+ IND 06
351 DSE 13
352 DSE 12
353 DSE 11
354 DSE 10
355 DSE 09
356 DSE 08
357 DSE 07
358 DSE 06
359 GTO 10
360 RCL 01
361 ST+ 06
362 ST+ 07
363 ST+ 08
364 ST+ 09
365 ST+ 10
366 ST+ 11
367 ST+ 12
368 ST+ 13
369 RCL 02
370 RCL 16
371 *
372 RCL 04
373 +
374 STO 50
375 XEQ 00
376 LBL 11
377 RCL IND 07
378 RCL 02
379 *
380 X<> IND 12
381 RCL 40
382 *
383 RCL IND 12
384 RCL 41
385 *
386 +
387 RCL IND 11
388 RCL 42
389 *
390 +
391 RCL IND 10
392 RCL 43
393 *
394 +
395 RCL IND 09
396 RCL 44
397 *
398 +
399 RCL IND 13
400 RCL 45
401 *
402 +
403 ST+ IND 06
404 DSE 13
405 DSE 12
406 DSE 11
407 DSE 10
408 DSE 09
409 DSE 07
410 DSE 06
411 GTO 11
412 RCL 01
413 ST+ 06
414 ST+ 07
415 ST+ 09
416 ST+ 10
417 ST+ 11
418 ST+ 12
419 ST+ 13
420 RCL 02
421 RCL 04
422 +
423 STO 50
424 XEQ 00
425 LBL 12
426 RCL IND 07
427 RCL 02
428 *
429 RCL IND 08
430 +
431 RCL 46
432 *
433 RCL IND 10
434 RCL IND 12
435 +
436 RCL 47
437 *
438 +
439 RCL IND 11
440 RCL 48
441 *
442 +
443 ST+ IND 06
444 DSE 12
445 DSE 11
446 DSE 10
447 DSE 08
448 DSE 07
449 DSE 06
450 GTO 12
451 RCL 01
452 ST+ 06
453 ST+ 07
454 ST+ 08
455 ST+ 10
456 ST+ 11
457 ST+ 12
458 DSE 49
459 GTO 01
a three-byte GTO
460 RCL 53
461 RCL 52
462 RCL 51
463 RCL 50
464 END
( 765 bytes / SIZE 51+9n )
| STACK | INPUTS | OUTPUTS |
| T | / | y3(x0+N.h) |
| Z | / | y2(x0+N.h) |
| Y | / | y1(x0+N.h) |
| X | / | x0+N.h |
Example:
dy/dx = -y.z.u
y(0) = 1 Evaluate
y(1) , z(1) , u(1) with h = 0.1 and N =
10
dz/dx = x.( y + z - u )
z(0) = 1
du/dx = x.y - z.u
u(0) = 2
LBL "T" *
RCL 51 -
LASTX RCL 53
RTN
"T"
ASTO 00
RCL 51 *
RCL 52 RCL 50
RCL 51 *
RCL 52 CHS
+
*
*
-
RCL 53 STO 54
RCL 53 STO 55
RCL 52 STO 56
Initial values: 0 STO
50
1 STO 51 STO 52
2 STO 53
n = 3 STO 01 (
3 functions )
h = 0.1 STO 02
N = 10 STO 03 XEQ "RK8" >>>>
x = 1
= R50
( in 9mn31s )
RDN y(1) = 0.258207906 = R51
y(1) = 0.258207906
RDN z(1) = 1.157623979 = R52
the exact results, rounded to 9D are
z(1) = 1.157623981
RDN u(1) = 0.842178313 = R53
u(1) = 0.842178312
-To continue the calculations, simply press R/S ( after changing
h & N in registers R02 & R03 if needed )
-If you replace lines 431 thru 442 by 9 * RCL IND
10 RCL IND 12 + 49 * + RCL IND 11
64 * + PI R-D /
registers R46-R47-R48 will be unused.
d) Runge-Kutta-Fehlberg 4-5
Method
-The following program uses a 4th-order formula to compute the solution,
and the difference between this 4th-order formula and a 5th-order
formula provides an error estimate.
Data Registers: • R00 = subroutine name ( Registers R00 thru R03 and R20 thru R20+n are to be initialized before executing "RKF" )
• R01 = n = number of equations = number of functions
R04 thru R14: temp R15 thru R19:
unused
• R02 = h = stepsize
• R03 = N = number of steps
• R20 = x0
at the end:
• R21 = y1(x0)
R21+n = err(y1)
R21+n thru R20+8n: temp
• R22 = y2(x0)
R22+n = err(y2)
...............
......................
• R20+n = yn(x0) R20+2n = err(yn)
Where err(yi) is an error-estimate of yi ( during the calculations, these error-estimates are moved in registers R21+2n .... R20+3n )
Flags: /
Subroutine: A program which calculates
and stores f1 ( x, y1 , ... , yn)
, ... , fn ( x, y1 , ... , yn) in
registers R21+n, ... , R20+2n respectively
with x, y1, ... , yn in regiters R20, R21,
... , R20+n
01 LBL "RKF"
02 21.021
03 RCL 01
04 E6
05 /
06 +
07 RCL 01
08 7
09 *
10 +
11 STO 05
12 FRC
13 RCL 01
14 .1
15 %
16 +
17 +
18 RCL 01
19 +
20 21
21 +
22 STO 06
23 RCL 01
24 ST+ X
25 E3
26 /
27 RCL 01
28 +
29 21.02
30 +
31 CLRGX
If you don't have an HP-41 CX, replace this line by 0
LBL 00 STO IND Y ISG Y GTO 00
32 RCL 01
33 RCL 01
34 LASTX
35 1
36 -
37 +
38 STO 07
39 +
40 STO 08
41 +
42 STO 13
43 +
44 STO 09
45 +
46 STO 10
47 +
48 STO 11
49 +
50 STO 12
51 LBL 10
52 RCL 03
53 STO 04
54 RCL 06
55 REGSWAP
56 LBL 01
57 RCL 20
58 STO 14
59 RCL 05
60 REGSWAP
61 REGMOVE
62 XEQ IND 00
63 LBL 02
64 RCL IND 08
65 RCL 02
66 *
67 STO IND 09
68 ST+ X
69 9
70 /
71 ST+ IND 07
72 DSE 09
73 DSE 08
74 DSE 07
75 GTO 02
76 RCL 01
77 ST+ 07
78 ST+ 08
79 ST+ 09
80 RCL 02
81 ST+ X
82 LASTX
83 /
84 ST+ 20
85 XEQ IND 00
86 RCL 05
87 REGMOVE
88 LBL 03
89 RCL IND 08
90 RCL 02
91 *
92 4
93 /
94 STO IND 10
95 RCL IND 09
96 12
97 /
98 +
99 ST+ IND 07
100 DSE 10
101 DSE 09
102 DSE 08
103 DSE 07
104 GTO 03
105 RCL 01
106 ST+ 07
107 ST+ 08
108 ST+ 09
109 ST+ 10
110 RCL 02
111 9
112 /
113 ST+ 20
114 XEQ IND 00
115 RCL 05
116 REGMOVE
117 LBL 04
118 RCL IND 08
119 RCL 02
120 *
121 STO IND 11
122 270
123 *
124 RCL IND 10
125 972
126 *
127 -
128 RCL IND 09
129 69
130 *
131 +
132 128
133 /
134 ST+ IND 07
135 DSE 11
136 DSE 10
137 DSE 09
138 DSE 08
139 DSE 07
140 GTO 04
141 RCL 01
142 ST+ 07
143 ST+ 08
144 ST+ 09
145 ST+ 10
146 ST+ 11
147 RCL 02
148 2.4
149 /
150 ST+ 20
151 XEQ IND 00
152 RCL 05
153 REGMOVE
154 LBL 05
155 RCL IND 08
156 RCL 02
157 *
158 64
159 *
160 STO IND 12
161 RCL IND 09
162 85
163 *
164 -
165 60
166 /
167 RCL IND 10
168 RCL IND 11
169 5
170 /
171 -
172 27
173 *
174 +
175 ST+ IND 07
176 DSE 12
177 DSE 11
178 DSE 10
179 DSE 09
180 DSE 08
181 DSE 07
182 GTO 05
183 RCL 01
184 ST+ 07
185 ST+ 08
186 ST+ 09
187 ST+ 10
188 ST+ 11
189 ST+ 12
190 RCL 02
191 ST+ 14
192 RCL 14
193 STO 20
194 XEQ IND 00
195 RCL 05
196 REGMOVE
197 LBL 06
198 RCL IND 08
199 RCL 02
200 *
201 15
202 *
203 ST+ IND 12
204 RCL IND 12
205 RCL IND 11
206 351
207 *
208 +
209 RCL IND 10
210 540
211 *
212 -
213 RCL IND 09
214 65
215 *
216 +
217 432
218 /
219 ST+ IND 07
220 DSE 12
221 DSE 11
222 DSE 10
223 DSE 09
224 DSE 08
225 DSE 07
226 GTO 06
227 RCL 01
228 ST+ 07
229 ST+ 08
230 ST+ 09
231 ST+ 10
232 ST+ 11
233 ST+ 12
234 RCL 02
235 6
236 /
237 ST- 20
238 XEQ IND 00
239 RCL 05
240 REGMOVE
241 LBL 07
242 RCL IND 08
243 RCL 02
244 *
245 72
246 *
247 RCL IND 12
248 -
249 RCL IND 11
250 9
251 *
252 +
253 RCL IND 09
254 ST+ X
255 -
256 3
257 1/X
258 %
259 ST- IND 13
260 RCL IND 12
261 RCL IND 11
262 81
263 *
264 +
265 PI
266 R-D
267 /
268 RCL IND 09
269 9
270 /
271 +
272 ST+ IND 07
273 DSE 13
274 DSE 12
275 DSE 11
276 DSE 09
277 DSE 08
278 DSE 07
279 GTO 07
280 RCL 01
281 ST+ 07
282 ST+ 08
283 ST+ 09
284 ST+ 11
285 ST+ 12
286 ST+ 13
287 RCL 14
288 STO 20
289 DSE 04
290 GTO 01
a three-byte GTO
291 RCL 06
292 REGSWAP
293 RCL 23
294 RCL 22
295 RCL 21
296 RCL 20
297 RTN
298 GTO 10
a three-byte GTO
299 END
( 481 bytes / SIZE 21+8n )
| STACK | INPUTS | OUTPUTS |
| T | / | y3(x0+N.h) |
| Z | / | y2(x0+N.h) |
| Y | / | y1(x0+N.h) |
| X | / | x0+N.h |
Example:
dy/dx = -y.z.u
y(0) = 1 Evaluate
y(1) , z(1) , u(1) with h = 0.1 & N
= 10
dz/dx = x.( y + z - u )
z(0) = 1
du/dx = x.y - z.u
u(0) = 2
LBL "T" *
RCL 21 -
RCL 20 RCL 23
RTN "T"
ASTO 00
RCL 21 *
RCL 22 RCL 20
RCL 21 *
RCL 22 CHS
+
*
*
-
RCL 23 STO 24
RCL 23 STO 25
RCL 22 STO 26
Initial values: 0 STO
20
1 STO 21 STO 22
2 STO 23
n = 3 STO 01 (
3 functions )
h = 0.1 STO 02
N = 10 STO 03 XEQ "RKF" >>>>
x = 1
= R20
( in 5mn40s )
RDN y(1) = 0.258207319 = R21
err(y) = -603 E-9 = R24
RDN z(1) = 1.157624973 = R22
and the error estimates
err(z) = 1062 E-9 = R25
RDN u(1) = 0.842178529 = R23
err(u) = 1768 E-9 = R26
-Actually, the true errors are -587 E-9 for y(1) 992 E-9 for z(1) 217 E-9 for u(1)
-To continue the calculations, simply press R/S ( after changing h & N in registers R02 & R03 if needed )
-If you want to use the 5th-order formula to compute yi(x)
, replace line 259 with ST+ IND 07 ST+ IND 13 but the
errors
will be probably overestimated.
-In the above example, it gives: y(1) = 0.258207897
z(1) = 1.157624056
u(1) = 0.842178340
e) Bulirsch-Stoer Method
-This program uses the formulae given in "The Bulirsch-Stoer Method
for the HP-41" and - first of all - in "Numerical Recipes"
Data Registers: • R00 = subroutine name ( Registers R00-R01-R02 and R20 thru R20+n are to be initialized before executing "BST" )
• R01 = n = number of equations = number of functions
R03 thru R18: temp R19: unused
• R02 = tolerance = a small positive number
that specifies the desired accuracy
• R20 = x0
• R21 = y1(x0)
R21+n thru R20+14n: temp
• R22 = y2(x0)
............... ( initial
values )
• R20+n = yn(x0)
Flags: F09-F10
Subroutine: A program which calculates
and stores f1 ( x, y1 , ... , yn)
, ... , fn ( x, y1 , ... , yn) in
registers R21+n, ... , R20+2n respectively
with x, y1, ... , yn in regiters R20, R21,
... , R20+n
01 LBL "BST"
02 STO 03
03 1
04 STO 04
05 9
06 STO 05
07 21.021
08 RCL 01
09 E6
10 /
11 +
12 RCL 01
13 5
14 *
15 +
16 STO 08
17 RCL 01
18 ST+ X
19 ST- Y
20 .1
21 %
22 -
23 -
24 STO 16
25 20.02
26 RCL 01
27 +
28 STO 10
29 INT
30 .1
31 %
32 RCL 01
33 STO T
34 +
35 +
36 STO 11
37 +
38 STO 12
39 +
40 STO 13
41 +
42 STO 14
43 GTO 01
44 LBL 00
45 2
46 ST/ 04
47 GTO 02
48 LBL 01
49 RCL 20
50 STO 15
51 RCL 08
52 REGSWAP
53 REGMOVE
54 XEQ IND 00
55 RCL 16
56 REGMOVE
57 6
58 RCL 05
59 X>Y?
60 GTO 02
61 RCL 04
62 ST+ 04
63 LBL 02
64 SF 09
65 RCL 04
66 STO 17
67 ABS
68 RCL 03
69 RCL 15
70 -
71 ABS
72 X<=Y?
73 CF 09
74 X>Y?
75 X<>Y
76 LASTX
77 SIGN
78 *
79 STO 04
80 SF 10
81 CLX
82 STO 05
83 RCL 08
84 INT
85 STO 09
86 LBL 03
87 RCL 08
88 REGMOVE
89 2
90 ST+ 05
91 18
92 RCL 05
93 X=Y?
94 GTO 00
95 RCL 01
96 ST+ 09
97 SIGN
98 -
99 E3
100 /
101 ISG X
102 STO 07
103 RCL 04
104 RCL 05
105 /
106 STO 06
107 LBL 04
108 RCL IND 12
109 RCL 06
110 *
111 RCL IND 10
112 STO IND 13
113 +
114 STO IND 14
115 DSE 14
116 DSE 13
117 DSE 12
118 DSE 10
119 GTO 04
120 RCL 01
121 ST+ 10
122 ST+ 12
123 ST+ 13
124 ST+ 14
125 LBL 05
126 RCL 08
127 RCL 01
128 -
129 REGMOVE
130 RCL 06
131 RCL07
132 INT
133 *
134 RCL 15
135 +
136 STO 20
137 XEQ IND 00
138 LBL 06
139 RCL IND 11
140 RCL 06
141 ST+ X
142 *
143 RCL IND 14
144 X<> IND 13
145 +
146 STO IND 14
147 DSE 14
148 DSE 13
149 DSE 11
150 GTO 06
151 RCL 01
152 ST+ 11
153 ST+ 13
154 ST+ 14
155 ISG 07
156 GTO 05
157 RCL 08
158 RCL 01
159 -
160 REGMOVE
161 RCL 04
162 RCL 15
163 +
164 STO 20
165 XEQ IND 00
166 LBL 07
167 RCL IND 11
168 RCL 06
169 *
170 RCL IND 13
171 +
172 ST+ IND 10
173 2
174 ST/ IND 10
175 SIGN
176 ST- 11
177 ST- 13
178 DSE 10
179 GTO 07
180 21
181 RCL 09
182 E3
183 /
184 +
185 RCL 01
186 ST+ 10
187 ST+ 11
188 ST+ 13
189 E6
190 /
191 +
192 REGMOVE
193 FS?C 10
194 GTO 03
a three-byte GTO
195 RCL 09
196 RCL 01
197 +
198 DSE X
199 E3
200 /
201 RCL 09
202 +
203 STO 06
204 CLX
205 STO 18
206 LBL 08
207 RCL 05
208 2
209 ST- Y
210 /
211 STO 07
212 LBL 09
213 RCL 06
214 RCL 07
215 RCL 01
216 *
217 -
218 RCL IND 06
219 ENTER^
220 X<> IND Z
221 STO Z
222 -
223 RCL 05
224 X^2
225 ST* Y
226 RCL 07
227 ST+ X
228 X^2
229 -
230 /
231 +
232 STO IND 06
233 DSE 07
234 GTO 09
235 LASTX
236 ABS
237 RCL 18
238 X<Y?
239 X<>Y
240 STO 18
241 ISG 06
242 GTO 08
243 RCL 02
244 X<Y?
245 GTO 03
a three-byte GTO
246 RCL 09
247 RCL 08
248 FRC
249 +
250 REGMOVE
251 FS? 09
252 GTO 01
a three-byte GTO
253 RCL 23
254 RCL 22
255 RCL 21
256 RCL 20
257 RTN
258 STO 03
259 RCL 04
260 RCL 17
261 X=Y?
262 GTO 01
a three-byte GTO
263 STO 04
264 9
265 STO 05
266 GTO 01
a three-byte GTO
267 END
( 395 bytes / SIZE 21+14n )
| STACK | INPUTS | OUTPUTS |
| T | / | y3(x1) |
| Z | / | y2(x1) |
| Y | / | y1(x1) |
| X | x1 | x1 |
Example:
dy/dx = -y.z.u
y(0) = 1 Evaluate
y(1) , z(1) , u(1) with a tolerance = 10 -7
dz/dx = x.( y + z - u )
z(0) = 1
du/dx = x.y - z.u
u(0) = 2
LBL "T" *
RCL 21 -
LASTX RCL 23
RTN "T"
ASTO 00
RCL 21 *
RCL 22 RCL 20
RCL 21 *
RCL 22 CHS
+
*
*
-
RCL 23 STO 24
RCL 23 STO 25
RCL 22 STO 26
Initial values: 0 STO
20
1 STO 21 STO 22
2 STO 23
n = 3 STO 01
( 3 functions )
tol = E-7 STO 02
1 XEQ "BST" >>>>
x = 1
= R20
( in 12mn06s )
RDN y(1) = 0.258207909 = R21
y(1) = 0.258207906
RDN z(1) = 1.157623986 = R22
the exact results, rounded to 9D are
z(1) = 1.157623981
RDN u(1) = 0.842178304 = R23
u(1) = 0.842178312
-The long execution time is due to the fact that the method is first
used with H = 1
but the desired accuracy is not satisfied. So the stepsize is
divided by 2.
-To compute y(2) z(2) u(2) key in 2 R/S it yields ( in 6m46s )
y(2) = 0.106363294
y(2) = 0.106363329
z(2) = 3.886706181
the exact results, rounded to 9D are
z(2) = 3.886706159
u(2) = 0.196515847
u(2) = 0.196515847
Notes:
-The initial stepsize is always +/-1 ( line 03 )
-If you want to choose your own initial stepsize, place it in Y-register
after replacing line 03 by X<>Y
-The modified midpoint rule ( and Richarson extrapolation ) are used
with h/2 , h/4 , h/6 , ..... , h/16 ( at most ) until the desired
accuracy is satisfied.
-If this is not satisfied, even with h/16, the stepsize is divided
by 2 ( LBL 00 )
-If - for instance - you want to stop the extrapolation with
h/14, replace line 91 by 16 ( instead of 18 )
-Thus, the SIZE may be reduced to 21+13n
-On the other hand, if you want to continue the extrapolation with
h/18, replace line 91 by 20
-In this case, the SIZE must be at least 21+15n
2°) Systems of second-order
Conservative Equations
-The 3 following programs use the 4th-order Runge-Kutta formula:
y(x+h) = y(x) + h ( y'(x) + k1/6 + 2k2/3
)
y'(x+h) = y'(x) + k1/6 + 2k2/3 +
k3/6
where k1 = h.f (x,y) ; k2 = h.f(x+h/2,y+h.y'/2+h.k1/8)
; k3 = h.f(x+h,y+h.y'+h.k2/2)
a)
Systems of 2 Equations
-"RKS2" solves the system:
y" = f(x,y,z)
y(x0) = y0 y'(x0)
= y'0
z" = g(x,y,z)
z(x0) = z0 z'(x0)
= z'0
Data Registers: • R00: subroutine name ( Registers R00 thru R07 are to be initialized before executing "RKS2" )
• R01 = x0 •
R04 = h/2 = half of the stepsize •
R06 = y'0
R08 to R12: temp
• R02 = y0 •
R05 = N = number of steps
• R07 = z'0
• R03 = z0
Flags: /
Subroutine: a program which calculates 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 "RKS2"
02 RCL 05
03 STO 12
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 09
14 2
15 /
16 RCL 07
17 +
18 RCL 04
19 *
20 RCL 03
21 +
22 X<>Y
23 STO 08
24 2
25 /
26 RCL 06
27 +
28 RCL 04
29 *
30 RCL 02
31 +
32 RCL 01
33 XEQ IND 00
34 RCL 04
35 ST+ 01
36 ST* Z
37 *
38 STO 11
39 RCL 07
40 +
41 X<>Y
42 STO 10
43 RCL 06
44 +
45 RCL 04
46 ST+ X
47 ST* Z
48 *
49 RCL 03
50 ST+ Z
51 CLX
52 RCL 02
53 +
54 RCL 01
55 XEQ IND 00
56 RCL 04
57 ST* Z
58 *
59 RCL 09
60 +
61 RCL 11
62 ST+ X
63 ST+ 09
64 ST+ X
65 +
66 3
67 ST/ 09
68 /
69 X<> 07
70 ST+ 07
71 RCL 09
72 +
73 X<>Y
74 RCL 08
75 +
76 RCL 10
77 ST+ X
78 ST+ 08
79 ST+ X
80 +
81 3
82 ST/ 08
83 /
84 X<> 06
85 ST+ 06
86 RCL 08
87 +
88 RCL 04
89 ST+ X
90 ST* Z
91 *
92 ST+ 02
93 X<>Y
94 ST+ 03
95 DSE 12
96 GTO 01
a three-byte GTO
97 RCL 03
98 RCL 02
99 RCL 01
100 END
( 139 bytes / SIZE 013 )
| STACK | INPUTS | OUTPUTS |
| Z | / | z(x0+N.h) |
| Y | / | y(x0+N.h) |
| X | / | x0+N.h |
Example:
d2y/dx2 = -y.z
y(0) = 2 y'(0) = 1
Evaluate y(1) , z(1) , y'(1) , z'(1) with
h = 0.1 & N = 10
d2z/dx2 = x.( y + z )
z(0) = 1 z'(0) = 1
LBL "T" CHS
*
"T" ASTO 00
RDN
ST* Z RTN
STO Z -
X<>Y R^
0 STO 01
1 STO 06
2 STO 02
1 STO 07
1 STO 03
0.05 STO 04
10 STO 05
XEQ "RKS2" >>>> x =
1
( in 41 seconds )
RDN y(1) = 1.531358015 and
R06 = y'(1) = -2.312838895
RDN z(1) = 2.620254480
R07 = z'(1) = 2.941751649
-With h = 0.05 it yields y(1) =
1.531356736 y'(1) = -2.312840085
z(1) = 2.620254295 z'(1) =
2.941748608
-Actually, the exact results - rounded to 9D - are
y(1) = 1.531356646 y'(1) = -2.312840137
z(1) = 2.620254281 z'(1) =
2.941748399
-To continue the calculations, simply press R/S ( after changing
h/2 & N in registers R04 & R05 if needed )
b)
Systems of 3 Equations
-"RKS3" solves the system:
y" = f(x,y,z,u)
y(x0) = y0 y'(x0)
= y'0
z" = g(x,y,z,u)
z(x0) = z0 z'(x0)
= z'0
u" = h(x,y,z,u)
u(x0) = u0 u'(x0)
= u'0
Data Registers: • R00: subroutine name ( Registers R00 thru R09 are to be initialized before executing "RKS3" )
• R01 = x0 •
R05 = h/2 = half of the stepsize •
R07 = y'0
R10 to R16: temp
• R02 = y0 •
R06 = N = number of steps
• R08 = z'0
• R03 = z0
• R09 = u'0
• R04 = u0
Flags: /
Subroutine: a program which calculates 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.
01 LBL "RKS3"
02 RCL 06
03 STO 16
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 12
16 2
17 /
18 RCL 09
19 +
20 RCL 05
21 *
22 RCL 04
23 +
24 X<>Y
25 STO 11
26 2
27 /
28 RCL 08
29 +
30 RCL 05
31 *
32 RCL 03
33 +
34 R^
35 STO 10
36 2
37 /
38 RCL 07
39 +
40 RCL 05
41 *
42 RCL 02
43 +
44 RCL 01
45 XEQ IND 00
46 RCL 05
47 ST+ 01
48 ST* Z
49 ST* T
50 *
51 STO 15
52 RCL 09
53 +
54 X<>Y
55 STO 14
56 RCL 08
57 +
58 R^
59 STO 13
60 RCL 07
61 +
62 RCL 05
63 ST+ X
64 ST* Z
65 ST* T
66 *
67 RCL 04
68 ST+ T
69 CLX
70 RCL 03
71 ST+ Z
72 CLX
73 RCL 02
74 +
75 RCL 01
76 XEQ IND 00
77 RCL 05
78 ST* Z
79 ST* T
80 *
81 RCL 12
82 +
83 RCL 15
84 ST+ X
85 ST+ 12
86 ST+ X
87 +
88 3
89 ST/ 12
90 /
91 X<> 09
92 ST+ 09
93 RCL 12
94 +
95 X<>Y
96 RCL 11
97 +
98 RCL 14
99 ST+ X
100 ST+ 11
101 ST+ X
102 +
103 3
104 ST/ 11
105 /
106 X<> 08
107 ST+ 08
108 RCL 11
109 +
110 R^
111 RCL 10
112 +
113 RCL 13
114 ST+ X
115 ST+ 10
116 ST+ X
117 +
118 3
119 ST/ 10
120 /
121 X<> 07
122 ST+ 07
123 RCL 10
124 +
125 RCL 05
126 ST+ X
127 ST* Z
128 ST* T
129 *
130 ST+ 02
131 X<>Y
132 ST+ 03
133 R^
134 ST+ 04
135 DSE 16
136 GTO 01
a three-byte GTO
137 RCL 04
138 RCL 03
139 RCL 02
140 RCL 01
141 END
( 194 bytes / SIZE 017 )
| STACK | INPUTS | OUTPUTS |
| T | / | u(x0+N.h) |
| Z | / | z(x0+N.h) |
| Y | / | y(x0+N.h) |
| X | / | x0+N.h |
Example:
d2y/dx2 = -y.z.u
y(0) = 1 y'(0) = 1
Evaluate y(1) , z(1) , u(1) , y'(1) , z'(1) , u'(1)
with h = 0.1 & N = 10
d2z/dx2 = x.( y + z - u )
z(0) = 1 z'(0) = 1
d2u/dx2 = x.y - z.u
u(0) = 2 u'(0) = 1
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 STO 07
STO 08 STO 09
2 STO 04
0.05 STO 05
10 STO 06 XEQ "RKS3"
>>>> x = 1
( in 65 seconds )
RDN y(1) = 0.439528419
y'(1) = -2.101120400 = R07
RDN z(1) = 2.070938499
z'(1) = 1.269599239 = R08
RDN u(1) = 1.744522976
u'(1) = -1.704232092 = R09
-With h = 0.05 it yields y(1) =
0.439524393
y'(1) = -2.101122784
z(1) = 2.070940521
z'(1) = 1.269597110
u(1) = 1.744524843
u'(1) = -1.704234567
-Actually, the exact results rounded to 9D are:
y(1) = 0.439524100
y'(1) = -2.101122880
z(1) = 2.070940654
z'(1) = 1.269596950
u(1) = 1.744524964
u'(1) = -1.704234756
-Press R/S to continue the calculations ( after changing h &
N in registers R05 & R06 if needed )
-The subroutine may be difficult to write, but "RKS3" is faster than
"RKS" below, for 3 equations.
c)
Systems of n Equations
Data Registers: • R00 = subroutine name ( Registers R00 thru R03 and R10 thru R10+2n are to be initialized before executing "RKS" )
• R01 = n = number of equations = number of functions
R04 thru R09: temp
• R02 = h = stepsize
• R03 = N = number of steps
• R10 = x0
• R11 = y1(x0)
• R11+n = y'1(x0)
R11+2n thru R10+6n: temp
• R12 = y2(x0)
• R12+n = y'2(x0)
...............
......................
• R10+n = yn(x0) • R10+2n = y'n(x0)
( during the calculations, y'i are moved in registers R11+2n .... R10+3n )
Flags: /
Subroutine: A program which calculates
and stores f1 ( x, y1 , ... , yn)
, ... , fn ( x, y1 , ... , yn) in
registers R11+n, ... , R10+2n respectively
with x, y1, ... , yn in regiters R10, R11,
... , R10+n
01 LBL "RKS"
02 RCL 03
03 STO 04
04 11.011
05 RCL 01
06 E6
07 /
08 +
09 RCL 01
10 5
11 *
12 +
13 STO 05
14 RCL 01
15 RCL 01
16 RCL 01
17 10.01
18 +
19 STO 06
20 +
21 STO 07
22 +
23 STO 08
24 +
25 STO 09
26 LBL 01
27 RCL 05
28 REGSWAP
29 REGMOVE
30 FRC
31 11
32 +
33 RCL 01
34 .1
35 %
36 ST+ X
37 +
38 +
39 REGMOVE
40 XEQ IND 00
41 LBL 02
42 RCL IND 07
43 RCL 02
44 *
45 STO IND 09
46 4
47 /
48 RCL IND 08
49 +
50 RCL 02
51 *
52 2
53 /
54 ST+ IND 06
55 DSE 09
56 DSE 08
57 DSE 07
58 DSE 06
59 GTO 02
60 RCL 01
61 ST+ 06
62 ST+ 07
63 ST+ 08
64 ST+ X
65 ST+ 09
66 RCL 02
67 LASTX
68 /
69 ST+ 10
70 XEQ IND 00
71 RCL 05
72 REGMOVE
73 LBL 03
74 RCL IND 07
75 RCL 02
76 *
77 STO IND 09
78 2
79 /
80 RCL IND 08
81 +
82 RCL 02
83 *
84 ST+ IND 06
85 DSE 09
86 DSE 08
87 DSE 07
88 DSE 06
89 GTO 03
90 RCL 01
91 ST+ 06
92 ST+ 08
93 ST+ 09
94 3
95 *
96 ST+ 07
97 RCL 02
98 2
99 /
100 ST+ 10
101 XEQ IND 00
102 RCL 05
103 REGMOVE
104 LBL 04
105 RCL IND 09
106 ST+ X
107 ST+ IND 07
108 X<> IND 07
109 ST+ IND 07
110 6
111 /
112 RCL IND 08
113 +
114 RCL 02
115 *
116 ST+ IND 06
117 RCL 08
118 RCL 01
119 -
120 RCL IND 07
121 RCL IND Y
122 RCL 02
123 *
124 +
125 6
126 /
127 RCL IND 08
128 +
129 STO IND Y
130 DSE 09
131 DSE 08
132 DSE 07
133 DSE 06
134 GTO 04
135 RCL 01
136 ST+ 06
137 ST- 07
138 ST+ 08
139 DSE 04
140 GTO 01
a three-byte GTO
141 RCL 13
142 RCL 12
143 RCL 11
144 RCL 10
145 END
( 225 bytes / SIZE 6n+11 )
| STACK | INPUTS | OUTPUTS |
| T | / | y3(x0+N.h) |
| Z | / | y2(x0+N.h) |
| Y | / | y1(x0+N.h) |
| X | / | x0+N.h |
Example:
d2y/dx2 = -y.z.u
y(0) = 1 y'(0) = 1
Evaluate y(1) , z(1) , u(1) , y'(1) , z'(1) , u'(1)
with h = 0.1 & N = 10
d2z/dx2 = x.( y + z - u )
z(0) = 1 z'(0) = 1
d2u/dx2 = x.y - z.u
u(0) = 2 u'(0) = 1
LBL "T" *
RCL 11 -
RCL 10 RCL 13
RTN "T"
ASTO 00
RCL 11 *
RCL 12 RCL 10
RCL 11 *
RCL 12 CHS
+
*
*
-
RCL 13 STO 14
RCL 13 STO 15
RCL 12 STO 16
Initial values: 0 STO
10
1 STO 11 STO 12 STO 14
STO 15 STO 16
2 STO 13
n = 3 STO 01 (
3 functions )
h = 0.1 STO 02
N = 10 STO 03 XEQ "RKS" >>>>
x = 1
= R10
( in 2mn19s )
RDN y(1) = 0.439528419 = R11
y'(1) = -2.101120401 = R14
RDN z(1) = 2.070938499 = R12
z'(1) = 1.269599239 = R15
RDN u(1) = 1.744522977 = R13
u'(1) = -1.704232091 = R16
-Press R/S to continue the calculations ( after changing h &
N in registers R02 & R03 if needed )
References:
[1] J.C. Butcher - "The Numerical Analysis of Ordinary Differential
Equations" - John Wiley & Sons - ISBN 0-471-91046-5
[2] Abramowitz and Stegun - "Handbook of Mathematical Functions"
- Dover Publications - ISBN 0-486-61272-4
[3] F. Scheid - "Numerical Analysis" - Mc Graw Hill
ISBN 07-055197-9
[4] Press , Vetterling , Teukolsky , Flannery - "Numerical Recipes
in C++" - Cambridge University Press - ISBN 0-521-75033-4
Go back to the HP-41 software library
Go back to the general software library
Go
back to the main exhibit hall