The Museum of HP Calculators
The Museum of HP Calculators
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°) A 4th-order Runge-Kutta method
a) Inertial Frame of Reference
b) Heliocentric Coordinates
2°) Numerov's Method ( X-Functions Module required )
a) Inertial Frame of Reference (
new )
b) Heliocentric Coordinates (
new )
3°) A Multistep Method of order 7 ( X-Functions Module required )
a) Inertial Frame of Reference (
new )
b) Heliocentric Coordinates (
new )
Two kinds of programs are presented to solve the gravitational n-body
problem.
-In the first one ( "GM" ) , the rectangular coordinates x,y,z are
reckoned to an inertial frame of reference.
-In the second one ( "PLN" ) , one of the celestial bodies ( for instance
the Sun ) is chosen as the origin, which is more natural for planetary
motions.
The problem is treated in the Newtonian theory of gravitation and the
relativistic effects are not taken into account.
We have to solve a system of 3n second order differential equations
of the type: d2y/dt2 = f (y)
The time t and the first derivative y' = dy/dt do not appear
explicitly in this trajectory problem.
So, we can use a special 4th order Runge-Kutta method, namely:
y(t+h) = y(t) + h ( y'(t) + k1/6 + 2k2/3
)
y'(t+h) = y'(t) + k1/6 + 2k2/3 + k3/6
where k1 = h.f (y) ; k2 = h.f(y+h.y'/2+h.k1/8)
; k3 = h.f(y+h.y'+h.k2/2)
Alternatively, we can also use multisteps methods if we know 2 or more initial values ( at t , t-h , ... ) without knowing the speeds, for instance:
-Numerov's method ( of order 5 ): ym+1
= 2.ym - ym-1 + (h2/12).( fm+1
+ 10.fm + fm-1 )
-A multistep method of order 7:
ym+1 = ym + ym-2 - ym-3
+ (h2/240).( 17. fm+1 + 232.fm + 222.fm-1
+ 232.fm-2 + 17.fm-3 )
1°) A Fourth-Order Runge-Kutta Method
a) Inertial Frame of Reference
Let M1( x1, y1,z1)
; ....... ; Mn( xn, yn,zn)
be n celestial bodies with initial velocities V1(
x'1, y'1,z'1) ; ....... ; Vn(
x'n, y'n,z'n) at an instant t
and m1 , ...... , mn
the n masses of these bodies.
We want to know their future positions and velocities at an instant
t + N.h ( h = step size ; N = number of steps )
d2xi /dt2 = SUMj#i
Gmj ( xj - xi ) / [ ( xi -
xj )2+ ( yi - yj )2
+ ( zi - zj )2 ]3/2
So, we have to solve the system:
d2yi /dt2 = SUMj#i
Gmj ( yj - yi ) / [ ( xi -
xj )2 + ( yi - yj )2
+ ( zi - zj )2 ]3/2
i , j = 1 , ..... , n
d2zi /dt2 = SUMj#i
Gmj ( zj - zi ) / [ ( xi -
xj )2 + ( yi - yj )2
+ ( zi - zj )2 ]3/2
where G is the gravitational constant.
It may seem very large but it's not too large for the HP-41:
If "GM" is executed directly from extended memory, it can solve
the 15-body problem.
Data Registers:
The following registers are to be initialized before executing "GM":
R00 = n = number of bodies
R27+n = x1
R27+4n = x'1
R24 = G = constant of gravitation
R28+n = y1
R28+4n = y'1
R25 = h = step size
R29+n = z1
R29+4n = z'1
R26 = N = number of steps
................
...................
R27 = m1
R28 = m2
R24+4n = xn
R24+7n = x'n
...............
R25+4n = yn
R25+7n = y'n
R26+n = mn
R26+4n = zn
R26+7n = z'n
>>> Thus, the n masses, the 3n position coordinates and the 3n velocity
coordinates are to be stored into contiguous registers, starting at R27
>>> Registers R27+16n thru R26+19n must be cleared before the first
execution ( you can XEQ CLRG before storing numbers )
-Other registers are used for temporary data storage ( control
numbers , partial sums , k-numbers of the Runge-Kutta scheme ... etc ...
)
N.B. G = 6.67259 10-11 m3
kg-1 s-2
but if
the units are Astronomical Unit , Solar mass and day , G = k2
where k = 0.01720209895 is Gaussian gravitational constant.
Program Listing
Note: If you don't have an X-Function Module,
replace lines 191 to 204 by and replace lines 18 to 34 by:
RCL 22
RCL 22
RCL 21
RCL 21
3
4
*
*
+
+
RCL 22
RCL 22
RCL 21
RCL 21
ST+ X
.1
.1
%
%
+
+
ST+ Y
ST+ Y
LBL 02
LBL 11
CLX
CLX
RCL IND Y
RCL IND Z
STO IND Z
STO IND Y
DSE Z
DSE Z
DSE Y
DSE Y
GTO 02
GTO 11
01 LBL "GM"
02 LBL 00
03 RCL 00
04 26
05 +
06 .1
07 %
08 +
09 RCL 00
10 3
11 *
12 STO 21
13 +
14 STO 22
15 RCL 26
16 STO 23
17 LBL 01
18 RCL 22
19 INT
20 1
21 +
22 STO Y
23 RCL 00
24 9
25 *
26 +
27 RCL 21
28 E3
29 /
30 +
31 E3
32 /
33 +
34 REGMOVE
35 CLX
36 STO 07
37 STO 08
38 RCL 25
39 6
40 STO 10
41 /
42 STO 09
43 XEQ 04
44 RCL 25
45 2
46 ST* 09
47 /
48 STO 07
49 4
50 ST/ 10
51 /
52 RCL 25
53 *
54 STO 08
55 XEQ 04
56 RCL 25
57 STO 07
58 CLX
59 STO 09
60 4
61 ST* 08
62 ST* 10
63 XEQ 04
64 4
65 RCL 21
66 ST* Y
67 RCL 22
68 +
69 ST+ Y
70 ENTER^
71 LBL 03
72 CLX
73 X<> IND Z
74 RCL 25
75 *
76 ST+ IND Y
77 DSE Z
78 DSE Y
79 GTO 03
80 DSE 23
81 GTO 01
82 RTN
83 GTO 00
( a synthetic three-byte GTO )
84 LBL 04
85 26.026
86 STO 11
87 3 E-5
88 RCL 22
89 +
90 STO 12
91 RCL 21
92 +
93 STO 14
94 LASTX
95 +
96 STO 16
97 LASTX
98 +
99 STO 18
100 LASTX
101 +
102 STO 19
103 LASTX
104 +
105 STO 20
106 LBL 05
107 RCL 00
108 ST+ 11
109 RCL 22
110 INT
111 STO 13
112 RCL 21
113 +
114 STO 15
115 LASTX
116 +
117 STO 17
118 CLX
119 STO 04
120 STO 05
121 STO 06
122 LBL 06
123 RCL 12
124 INT
125 RCL 13
126 X=Y?
127 GTO 08
( a synthetic three-byte GTO )
128 XEQ 09
129 STO 03
130 1
131 ST- 12
132 ST- 13
133 ST- 14
134 ST- 15
135 ST- 16
136 ST- 17
137 XEQ 09
138 STO 02
139 1
140 ST- 12
141 ST- 13
142 ST- 14
143 ST- 15
144 ST- 16
145 ST- 17
146 RCL IND 11
147 XEQ 09
148 STO 01
149 X^2
150 RCL 02
151 X^2
152 RCL 03
153 X^2
154 +
155 +
156 ENTER^
157 SQRT
158 *
159 /
160 ST* 01
161 ST* 02
162 ST* 03
163 RCL 01
164 ST+ 04
165 RCL 02
166 ST+ 05
167 RCL 03
168 ST+ 06
169 2
170 ST+ 12
171 ST+ 14
172 ST+16
173 SIGN
174 LBL 07
175 ST- 13
176 ST- 15
177 ST- 17
178 DSE 11
179 GTO 06
180 RCL 06
181 XEQ 10
182 RCL 05
183 XEQ 10
184 RCL 04
185 XEQ 10
186 3
187 ST- 14
188 ST- 16
189 DSE 12
190 GTO 05
( a synthetic three-byte GTO )
191 RCL 22
192 RCL 21
193 ST+ X
194 1
195 +
196 .1
197 %
198 +
199 +
200 RCL 21
201 E6
202 /
203 +
204 REGMOVE
205 RTN
206 LBL 08
207 3
208 GTO 07
209 LBL 09
210 RCL IND 13
211 RCL IND 12
212 -
213 RCL IND 15
214 RCL IND 14
215 -
216 RCL 07
217 *
218 +
219 RCL IND 17
220 RCL IND 16
221 -
222 RCL 08
223 *
224 +
225 RTN
226 LBL 10
227 RCL 24
228 *
229 ENTER^
230 STO IND 18
231 RCL 09
232 *
233 ST+ IND 19
234 CLX
235 RCL 10
236 /
237 ST+ IND 20
238 1
239 ST- 18
240 ST- 19
241 ST- 20
242 RTN
243 END
( 375 bytes / SIZE 27+19n )
Example: Let's imagine a system of three stars
M1 ( 2 ; 0 ;0 ) M2
( 0 ; 4 ; 0 ) M3 ( 0 ; 0 ; 1 )
unit = 1 AU
with initial velocities V1
( 0 ; 0.03 ; 0 ) V2 ( 0 ; 0 ; 0.01 )
V3 ( -0.02 ; 0 ; 0 )
unit = 1 AU per day
and masses m1 = 2
; m2 = 1 ; m3 = 3
unit = 1 Solar mass
-Calculate the positions and velocities 10 days later.
- SIZE 084 ( or greater )
XEQ CLRG ( or 75.083 XEQ CLRGX if you have
an HP-41 CX )
3 bodies
>>>
3 STO 00
constant of gravitation >>>
17.20209895 E-3 X^2 STO 24
step size h = 10 days >>>
10 STO 25
number of steps: 1
>>>
1 STO 26
then, we store the 3 masses and the 18 position and velocity coordinates as shown below and XEQ "GM"
>>> the new positions and velocities have replaced the previous
ones in registers R30 thru R47 ( next to last column )
Execution time = 1mn52s
| Register | input | h = 10 ; N=1 | h = 5 ; N=2 | |
| m | R27=m1
R28=m2 R29=m3 |
2
1 3 |
undisturbed | undisturbed |
|
p
|
R30=x1
R31=y1 R32=z1 -------- R33=x2 R34=y2 R36=z2 -------- R36=x3 R37=y3 R38=z3 |
2
0 0 ----- 0 4 0 ----- 0 0 1 |
1.992077590
0.300333861 0.003673761 -------------- 0.000661665 3.996080594 0.100603408 -------------- -0.194938948 0.001083895 0.997349690 |
1.992077585
0.300333570 0.003673682 -------------- 0.000661669 3.996080575 0.100603412 -------------- -0.194938946 0.001084095 0.997349741 |
| v
e l . |
R39=x'1
R40=y'1 R41=z'1 -------- R42=x'2 R43=y'2 R44=z'2 -------- R45=x'3 R46=y'3 R47=z'3 |
0
0.03 0 ----- 0 0 0.01 ----- -0.02 0 0 |
-0.001550090
0.030038159 0.000706688 -------------- 0.000132598 -0.000790384 0.010117548 -------------- -0.019010806 0.000238022 -0.000510308 |
-0.001550083
0.030038158 0.000706684 -------------- 0.000132598 -0.000790385 0.010117549 -------------- -0.019010811 0.000238023 -0.000510306 |
-To estimate the accuracy of the results, we can perform the calculation
with a half step size ( h = 5 days ) and N = 2 instead of 1:
We obtain the results in the last column: errors are smaller
than 0.000001 AU
- To continue with the same h and N , simply press R/S
b) Heliocentric Coordinates
When computing planetary trajectories for instance, it's often better
to know directly the heliocentric coordinates by choosing the Sun as the
origin.
In this case however, the reference frame is no more inertial , the
equations become more complex and therefore, the program is longer.
On the other hand, it executes faster and requires less data registers
( the 16-body problem can be solved ). The equations are now:
d2xi /dt2 = - Gm0 xi/ri3
- SUMj Gmj xj/rj3
+ SUMj#i Gmj ( xj -
xi ) / rij3
d2yi /dt2 = - Gm0 yi/ri3
- SUMj Gmj yj/rj3
+ SUMj#i Gmj ( yj -
yi ) / rij3
i , j = 1 , ..... , n-1
d2zi /dt2 = - Gm0 zi/ri3
- SUMj Gmj zj/rj3
+ SUMj#i Gmj ( zj -
zi ) / rij3
where ri = ( xi2 + yi2
+ zi2 ) 1/2 and
rij = [ ( xi - xj )2 + ( yi
- yj )2 + ( zi - zj )2
] 1/2
Data Registers: The following registers are to be initialized before executing "PLN":
R00 = n-1 = number of bodies minus 1
R30+n = x1
R27+4n = x'1
R27 = G = constant of gravitation
R31+n = y1
R28+4n = y'1
R28 = h = step size
R32+n = z1
R29+4n = z'1
R29 = N = number of steps
................
...................
R30 = m0 = mass of the origin
R31 = m1
R24+4n = xn-1
R21+7n = x'n-1
...............
R25+4n = yn-1
R22+7n = y'n-1
R26+4n = zn-1
R23+7n = z'n-1
R29+n = mn-1
>>> Thus, the n masses, the 3n-3 position coordinates and the 3n-3 velocity
coordinates are to be stored into contiguous registers, starting at R30
>>> Registers R15+16n thru R11+19n must be cleared before
the first execution ( you can XEQ CLRG before storing numbers
)
-Other registers are used for temporary data storage ( control numbers , partial sums , k-numbers of the Runge-Kutta scheme ... etc ... )
N.B. G = 6.67259 10-11 m3
kg-1 s-2
but if
units are Astronomical Unit , Solar mass and day , G = k2
where k = 0.01720209895 is Gaussian gravitational constant.
Program Listing
Note: If you don't have an X-Function Module,
replace lines 221 to 234 by and replace lines 18 to 34 by:
RCL 25
RCL 25
RCL 24
RCL 24
3
4
*
*
+
+
RCL 25
RCL 25
RCL 24
RCL 24
ST+ X
.1
.1
%
%
+
+
ST+ Y
ST+ Y
LBL 02
LBL 13
CLX
CLX
RCL IND Y
RCL IND Z
STO IND Z
STO IND Y
DSE Z
DSE Z
DSE Y
DSE Y
GTO 02
GTO 13
01 LBL "PLN"
02 LBL 00
03 RCL 00
04 30
05 +
06 .1
07 %
08 +
09 RCL 00
10 3
11 *
12 STO 24
13 +
14 STO 25
15 RCL 29
16 STO 26
17 LBL 01
18 RCL 25
19 INT
20 1
21 +
22 STO Y
23 RCL 00
24 9
25 *
26 +
27 RCL 24
28 E3
29 /
30 +
31 E3
32 /
33 +
34 REGMOVE
35 CLX
36 STO 07
37 STO 08
38 RCL 28
39 6
40 STO 10
41 /
42 STO 09
43 XEQ 04
44 RCL 28
45 2
46 ST* 09
47 /
48 STO 07
49 4
50 ST/ 10
51 /
52 RCL 28
53 *
54 STO 08
55 XEQ 04
56 RCL 28
57 STO 07
58 CLX
59 STO 09
60 4
61 ST* 08
62 ST* 10
63 XEQ 04
64 4
65 RCL 24
66 ST* Y
67 RCL 25
68 +
69 ST+ Y
70 ENTER^
71 LBL 03
72 CLX
73 X<> IND Z
74 RCL 28
75 *
76 ST+ IND Y
77 DSE Z
78 DSE Y
79 GTO 03
80 DSE 26
81 GTO 01
82 RTN
83 GTO 00
( a synthetic three-byte GTO )
84 LBL 04
85 30.03
86 RCL 00
87 +
88 STO 11
89 RCL 25
90 STO 12
91 RCL 24
92 +
93 STO 14
94 LASTX
95 +
96 STO 16
97 LASTX
98 +
99 STO 18
100 LASTX
101 +
102 STO 19
103 LASTX
104 +
105 STO 20
106 CLX
107 STO 04
108 STO 05
109 STO 06
110 LBL 05
111 XEQ 11
112 STO 03
113 1
114 ST- 12
115 ST- 14
116 ST- 16
117 XEQ 11
118 STO 02
119 1
120 ST- 12
121 ST- 14
122 ST- 16
123 XEQ 11
124 STO 01
125 1
126 ST- 12
127 ST- 14
128 ST- 16
129 XEQ 12
130 DSE 11
131 GTO 05
132 RCL 04
133 STO 21
134 RCL 05
135 STO 22
136 RCL 06
137 STO 23
138 RCL 24
139 ST+ 12
140 ST+ 14
141 ST+ 16
142 LBL 06
143 RCL 00
144 ST+ 11
145 RCL 25
146 STO 13
147 RCL 24
148 +
149 STO 15
150 LASTX
151 +
152 STO 17
153 RCL 21
154 STO 04
155 RCL 22
156 STO 05
157 RCL 23
158 STO 06
159 LBL 14
160 RCL 12
161 RCL 13
162 X=Y?
163 GTO 08 (
a synthetic three-byte GTO )
164 XEQ 09
165 STO 03
166 1
167 ST- 12
168 ST- 13
169 ST- 14
170 ST- 15
171 ST- 16
172 ST- 17
173 XEQ 09
174 STO 02
175 1
176 ST- 12
177 ST- 13
178 ST- 14
179 ST- 15
180 ST- 16
181 ST- 17
182 XEQ 09
183 STO 01
184 XEQ 12
185 2
186 ST+ 12
187 ST+ 14
188 ST+ 16
189 SIGN
190 LBL 07
191 ST- 13
192 ST- 15
193 ST- 17
194 DSE 11
195 GTO 14
196 XEQ 11
197 STO 03
198 1
199 ST- 12
200 ST- 14
201 ST- 16
202 XEQ 11
203 STO 02
204 1
205 ST- 12
206 ST- 14
207 ST- 16
208 XEQ 11
209 STO 01
210 XEQ 12
211 RCL 06
212 XEQ 10
213 RCL 05
214 XEQ 10
215 RCL 04
216 XEQ 10
217 ST- 14
218 ST- 16
219 DSE 12
220 GTO 06
( a synthetic three-byte GTO )
221 RCL 25
222 RCL 24
223 ST+ X
224 1
225 +
226 .1
227 %
228 +
229 +
230 RCL 24
231 E6
232 /
233 +
234 REGMOVE
235 RTN
236 LBL 08
237 3
238 GTO 07
239 LBL 09
240 RCL IND 12
241 RCL IND 13
242 -
243 RCL IND 14
244 RCL IND 15
245 -
246 RCL 07
247 *
248 +
249 RCL IND 16
250 RCL IND 17
251 -
252 RCL 08
253 *
254 +
255 RTN
256 LBL 10
257 RCL 27
258 *
259 CHS
260 ENTER^
261 STO IND 18
262 RCL 09
263 *
264 ST+ IND 19
265 CLX
266 RCL 10
267 /
268 ST+ IND 20
269 1
270 ST- 18
271 ST- 19
272 ST- 20
273 RTN
274 LBL 11
275 RCL IND 12
276 RCL IND 14
277 RCL 07
278 *
279 +
280 RCL IND 16
281 RCL 08
282 *
283 +
284 RTN
285 LBL 12
286 RCL IND 11
287 RCL 01
288 X^2
289 RCL 02
290 X^2
291 RCL 03
292 X^2
293 +
294 +
295 ENTER^
296 SQRT
297 *
298 /
299 ST* 01
300 ST* 02
301 ST* 03
302 RCL 01
303 ST+ 04
304 RCL 02
305 ST+ 05
306 RCL 03
307 ST+ 06
308 RTN
309 END
( 486 bytes / SIZE 12+19n )
Example: let's take the same 3 stars and let's choose the third star as the origin. The coordinates become:
M1 ( 2 ; 0 ; -1 ) M2 ( 0 ; 4 ; -1 ) and V1 ( 0.02 ; 0.03 ; 0 ) V2 ( 0.02 ; 0 ; 0.01 )
-Calculate the positions and velocities 10 days later.
- SIZE 069 ( or greater )
XEQ CLRG ( or 63.068 XEQ CLRGX if you
have an HP-41 CX )
3 bodies
>>>
3-1 = 2 STO 00
constant of gravitation >>>
17.20209895 E-3 X^2 STO 27
step size h = 10 days >>>
10 STO 28
number of steps: 1
>>>
1 STO 29
then, we store the 3 masses and the 12 rectangular coordinates ( position and velocity ) as shown below and XEQ "PLN"
>>> the new postions and velocities have replaced the previous
ones in registers R33 thru R44 ( last column )
Execution time = 1mn28s
| Register | input | h = 10 ; N=1 | |
| m | R30=m0
R31=m1 R32=m2 |
3
2 1 |
undisturbed |
| p
o s. |
R33=x1
R34=y1 R35=z1 -------- R36=x2 R37=y2 R38=z2 |
2
0 -1 ----- 0 4 -1 |
2.187016538
0.299249966 -0.993675929 -------------- 0.195600614 3.994996700 -0.896746283 |
| v
e l. |
R39=x'1
R40=y'1 R41=z'1 -------- R42=x'2 R43=y'2 R44=z'2 |
0.02
0.03 0 ----- 0.02 0 0.01 |
0.017460717
0.029800137 0.001216996 -------------- 0.019143404 -0.001028406 0.010627856 |
- To continue with the same h and N , simply press R/S
Accuracy:
-This 4th order Runge-Kutta method doesn't have a built-in error estimate.
-Therefore, you'll have to do the calculation a second time with a
smaller step size to estimate accuracy.
(When h is divided by 2 , errors are roughly divided by 16).
-If you apply "GM" or "PLN" to the whole Solar system, h = 1day
gives an error of about 7 10-6 AU in the position of Mercury
after 88 days.
Execution time:
Here are the results of a few tests ( with N = 1 ):
| n | GM | PLN |
| 3 | 1mn52 | 1mn28 |
| 5 | 5mn04 | 4mn27 |
| 10 | 19mn50 | 18mn51 |
-Execution time can be saved by using CLX SIGN instead of 1 because
digit entries are very slow.
Note: The "classical" 4th order Runge-Kutta
method could also be used to solve this problem and the programs would
be shorter,
but they would run a little slower and much more data registers would be
needed.
2°) Numerov's Method
a) Inertial Frame of Reference
Data Registers: • R00 = n = number of bodies ( Registers R00 and R13 thru R15+7n are to be initialized before executing "GM2" )
R01 to R10: temp R11-R12: unused
• R13 = G = k2
• R16+n = x1(0)
• R16+4n = x1(-1)
R16+7n thru R15+19n: temp
• R14 = h = stepsize
• R17+n = y1(0)
• R17+4n = y1(-1)
• R15 = N = number of steps
• R18+ n = z1(0)
• R18+4n = z1(-1)
• R16 = m1
................
.................
• R17 = m2
................
.................
.............
................
..................
• R15+n = mn
• R13+4n = xn(0)
• R13+7n = xn(-1)
• R14+4n = yn(0)
• R14+7n = yn(-1)
• R15+4n = zn(0)
• R15+7n = zn(-1)
where mi are the n masses,
xi(0) yi(0) zi(0)
are the positions at an instant t
and xi(-1) yi(-1)
zi(-1) are the positions at the instant
t-h
Flags: /
Subroutines: /
01 LBL "GM2"
02 RCL 00
03 7
04 *
05 15
06 +
07 STO 05
08 LASTX
09 RCL 00
10 4
11 *
12 +
13 E3
14 /
15 +
16 STO 04
17 RCL 00
18 19
19 *
20 15
21 +
22 STO 06
23 LASTX
24 RCL 00
25 +
26 .015
27 +
28 STO 07
29 XEQ 05
30 ST- 05
31 E3
32 /
33 ST- 04
34 XEQ 05
35 RCL 06
36 1
37 +
38 RCL Y
39 E6
40 /
41 +
42 ENTER^
43 INT
44 R^
45 -
46 E3
47 /
48 +
49 REGMOVE
lines 01 thru 49 are executed the first time only
50 LBL 01
51 RCL 15
52 STO 10
53 LBL 02
54 1.004006
55 RCL 00
56 *
57 16.016
58 +
59 REGMOVE
60 RCL 00
61 3
62 *
63 GTO 00
64 LBL 03
65 RCL 00
66 4
67 *
68 15
69 +
70 STO 05
71 LASTX
72 RCL 00
73 +
74 E3
75 /
76 +
77 STO 04
78 RCL 00
79 13
80 *
81 15
82 +
83 STO 06
84 XEQ 05
85 LBL 00
86 STO Y
87 RCL 00
88 15
89 +
90 .1
91 %
92 +
93 +
94 STO 01
95 +
96 STO 02
97 +
98 STO 03
99 +
100 STO 04
101 +
102 STO 05
103 +
104 STO 06
105 CLX
106 LBL 04
107 RCL IND 06
108 RCL IND 05
109 10
110 *
111 +
112 RCL IND 04
113 +
114 RCL 14
115 X^2
116 *
117 12
118 /
119 RCL IND 03
120 -
121 RCL IND 02
122 ST+ X
123 +
124 ENTER^
125 X<> IND 01
126 -
127 ABS
128 +
129 DSE 06
130 DSE 05
131 DSE 04
132 DSE 03
133 DSE 02
134 DSE 01
135 GTO 04
136 X#0?
137 GTO 03
138 RCL 04
139 INT
140 RCL 05
141 INT
142 1
143 ST+ Z
144 +
145 E3
146 /
147 +
148 RCL 00
149 6
150 *
151 E6
152 /
153 +
154 REGMOVE
155 DSE 10
156 GTO 02
a three-byte GTO
157 RCL 01
158 3
159 +
160 RCL IND X
161 DSE Y
162 RCL IND Y
163 DSE Z
164 RCL IND Z
165 RTN
166 GTO 01
a three-byte GTO
167 LBL 05
168 CLX
169 STO 01
170 STO 02
171 STO 03
172 LBL 06
173 RCL 04
174 INT
175 RCL 05
176 X=Y?
177 GTO 08
178 RCL IND 07
179 RCL IND 05
180 RCL IND 04
181 -
182 STO 08
183 X^2
184 DSE 04
185 DSE 05
186 RCL IND 05
187 RCL IND 04
188 -
189 STO 09
190 X^2
191 +
192 DSE 04
193 DSE 05
194 RCL IND 05
195 RCL IND 04
196 -
197 STO T
198 X^2
199 +
200 ENTER^
201 SQRT
202 *
203 /
204 ST* 08
205 ST* 09
206 *
207 ST+ 01
208 RCL 09
209 ST+ 02
210 RCL 08
211 ST+ 03
212
213 ST+ 04
214 LBL 07
215 DSE 05
216 DSE 07
217 GTO 06
218 RCL 13
219 ST* 01
220 ST* 02
221 RCL 03
222 *
223 STO IND 06
224 DSE 06
225 RCL 02
226 STO IND 06
227 DSE 06
228 RCL 01
229 STO IND 06
230 DSE 06
231 2
232 ST- 04
233 RCL 00
234 ST+ 07
235 3
236 *
237 ST+ 05
238 DSE 04
239 GTO 05
240 RTN
241 LBL 08
242 2
243 ST- 05
244 GTO 07
245 END
( 362 bytes / SIZE 19n+16 )
| STACK | INPUTS | OUTPUTS |
| Z | / | z1 |
| Y | / | y1 |
| X | / | x1 |
Example: n = 3 stars / m1 = 2 ; m2 = 1 ; m3 = 3 ( unit = 1 Solar mass ) We have the following coordinates:
t = 0
t = -5 days
masses
positions(0)
positions(-1)
( unit = 1 AU )
2 STO 16
2 STO 19
1.997888568 STO 28
x1
1 STO 17
0 STO 20
-0.149784693 STO 29
y1
3 STO 18
0 STO 21
0.001032468 STO 30
z1
0 STO 22
0.000165879 STO 31
x2
4 STO 23
3.999043454 STO 32
y2
0 STO 24
-0.049838219 STO 33
z2
0 STO 25
0.101352328 STO 34
x3
0 STO 26
0.000175311 STO 35
y3
1 STO 27
0.999257761 STO 36
z3
>Evaluate the coordinates of these 3 stars when t = 10 days
n = 3 STO 00 k2
= 17.20209895 E-3 X^2 STO 13
h = 5 days STO 14 N = 2 STO 15
XEQ "GM2" >>>> x1
= 1.992077642 = R19
( Execution time = 4mn58s ) and
x2 = 0.000661670 = R22 x3
= -0.194938984 = R25
RDN y1 = 0.300333555
= R20
y2 = 3.996080573 = R23 y3
= 0.001084105 = R26
RDN z1 = 0.003673650
= R21
z2 = 0.100603410 = R24 z3
= 0.997349763 = R27
-The errors are of the order of 6 10 -8 AU
-The positions at t = 5 days are in registers R28 thru R36
-Press R/S or XEQ 01 to continue ( after changing N in register R15
if needed )
-You can also press XEQ "GM2" but it would needlessly re-calculate
values which are already stored in the required registers.
Variant:
-The Numerov's formula is applied recursively until 2 successive approximations
are equal ( line 136 )
-Therefore, the complicated function ( LBL 05 ) is uselessly calculated
at least once.
-Alternatively, you can fix the number of iterations: store this
value in R11
Replace line 136
by DSE 12
Replace lines 124 to 128 by STO IND 01
Delete line 105
Add RCL 11 STO 12 after line 59
Note:
-When h is divided by 2, the errors are approximately divided by 16
= 24
b) Heliocentric Coordinates
Data Registers: • R00 = n-1 = number of bodies minus 1 = n' ( Registers R00 and R16 thru R19+7n' are to be initialized before executing "PLN2" )
R01 to R13: temp R14-R15: unused
• R16 = G = k2
• R20+n' = x1(0)
• R20+4n' = x1(-1)
R20+7n' thru R19+19n': temp
• R17 = h = stepsize
• R21+n' = y1(0)
• R21+4n' = y1(-1)
• R18 = N = number of steps
• R22+ n' = z1(0)
• R22+4n' = z1(-1)
• R19 = m0 = mass of the origin
................
.................
• R20 = m1
................
.................
( n' = n-1 )
.............
................
..................
• R19+n' = mn'
• R17+4n' = xn'(0)
• R17+7n' = xn'(-1)
• R18+4n' = yn'(0)
• R18+7n' = yn'(-1)
• R19+4n' = zn'(0)
• R19+7n' = zn'(-1)
where mi are the n masses,
xi(0) yi(0) zi(0)
are the positions at an instant t
and xi(-1) yi(-1)
zi(-1) are the positions at the instant
t-h
Flags: /
Subroutines: /
01 LBL "PLN2"
02 RCL 00
03 7
04 *
05 19
06 +
07 STO 08
08 LASTX
09 RCL 00
10 4
11 *
12 +
13 E3
14 /
15 +
16 STO 07
17 RCL 00
18 19
19 ST* Y
20 +
21 STO 09
22 LASTX
23 RCL 00
24 +
25 .019
26 +
27 STO 10
28 XEQ 05
29 ST- 08
30 E3
31 /
32 ST- 07
33 XEQ 05
34 RCL 09
35 1
36 +
37 RCL Y
38 E6
39 /
40 +
41 ENTER^
42 INT
43 R^
44 -
45 E3
46 /
47 +
48 REGMOVE
lines 01 thru 48 are executed the first time only
49 LBL 01
50 RCL 18
51 STO 13
52 LBL 02
53 1.004006
54 RCL 00
55 *
56 20.02
57 +
58 REGMOVE
59 RCL 00
60 3
61 *
62 GTO 00
63 LBL 03
64 RCL 00
65 4
66 *
67 19
68 +
69 STO 08
70 LASTX
71 RCL 00
72 +
73 E3
74 /
75 +
76 STO 07
77 RCL 00
78 13
79 *
80 19
81 +
82 STO 09
83 XEQ 05
84 LBL 00
85 STO Y
86 RCL 00
87 19
88 +
89 .1
90 %
91 +
92 +
93 STO 01
94 +
95 STO 02
96 +
97 STO 03
98 +
99 STO 04
100 +
101 STO 05
102 +
103 STO 06
104 CLX
105 LBL 04
106 RCL IND 06
107 RCL IND 05
108 10
109 *
110 +
111 RCL IND 04
112 +
113 RCL 17
114 X^2
115 *
116 12
117 /
118 RCL IND 03
119 -
120 RCL IND 02
121 ST+ X
122 +
123 ENTER^
124 X<> IND 01
125 -
126 ABS
127 +
128 DSE 06
129 DSE 05
130 DSE 04
131 DSE 03
132 DSE 02
133 DSE 01
134 GTO 04
135 X#0?
136 GTO 03
137 RCL 04
138 INT
139 RCL 05
140 INT
141 1
142 ST+ Z
143 +
144 E3
145 /
146 +
147 RCL 00
148 6
149 *
150 E6
151 /
152 +
153 REGMOVE
154 DSE 13
155 GTO 02
a three-byte GTO
156 RCL 01
157 3
158 +
159 RCL IND X
160 DSE Y
161 RCL IND Y
162 DSE Z
163 RCL IND Z
164 RTN
165 GTO 01
a three-byte GTO
166 LBL 05
167 CLX
168 STO 04
169 STO 05
170 STO 06
171 LBL 06
172 RCL IND 10
173 RCL IND 08
174 STO 11
175 X^2
176 DSE 08
177 RCL IND 08
178 STO 12
179 X^2
180 +
181 DSE 08
182 RCL IND 08
183 STO T
184 X^2
185 +
186 ENTER^
187 SQRT
188 *
189 /
190 ST* 11
191 ST* 12
192 *
193 ST- 04
194 RCL 12
195 ST- 05
196 RCL 11
197 ST- 06
198 DSE 08
199 DSE 10
200 GTO 06
201 RCL 00
202 ST+ 10
203 3
204 *
205 ST+ 08
206 LBL 07
207 RCL 04
208 STO 01
209 RCL 05
210 STO 02
211 RCL 06
212 STO 03
213 RCL 19
214 RCL IND 07
215 STO 11
216 X^2
217 DSE 07
218 RCL IND 07
219 STO 12
220 X^2
221 +
222 DSE 07
223 RCL IND 07
224 STO T
225 X^2
226 +
227 ENTER^
228 SQRT
229 *
230 /
231 ST* 11
232 ST* 12
233 *
234 ST- 01
235 RCL 12
236 ST- 02
237 RCL 11
238 ST- 03
239 2
240 ST+ 07
241 LBL 08
242 RCL 07
243 INT
244 RCL 08
245 X=Y?
246 GTO 10
247 RCL IND 10
248 RCL IND 08
249 RCL IND 07
250 -
251 STO 11
252 X^2
253 DSE 07
254 DSE 08
255 RCL IND 08
256 RCL IND 07
257 -
258 STO 12
259 X^2
260 +
261 DSE 07
262 DSE 08
263 RCL IND 08
264 RCL IND 07
265 -
266 STO T
267 X^2
268 +
269 ENTER^
270 SQRT
271 *
272 /
273 ST* 11
274 ST* 12
275 *
276 ST+ 01
277 RCL 12
278 ST+ 02
279 RCL 11
280 ST+ 03
288 2
282 ST+ 07
283 LBL 09
284 DSE 08
285 DSE 10
286 GTO 08
287 RCL 16
288 ST* 01
289 ST* 02
290 RCL 03
291 *
292 STO IND 09
293 DSE 09
294 RCL 02
295 STO IND 09
296 DSE 09
297 RCL 01
298 STO IND 09
299 DSE 09
300 2
301 ST- 07
302 RCL 00
303 ST+ 10
304 3
305 *
306 ST+ 08
307 DSE 07
308 GTO 07
a three-byte GTO
309 RTN
310 LBL 10
311 2
312 ST- 08
313 GTO 09
314 END
( 465 bytes / SIZE 19n' + 20 )
( n' = n-1 )
| STACK | INPUTS | OUTPUTS |
| Z | / | z1 |
| Y | / | y1 |
| X | / | x1 |
Example: n = 3 stars / m0 = 3 ; m1 = 2 ; m2 = 1 ( unit = 1 Solar mass ) We have the following coordinates:
t = 0
t = -5 days
masses
positions(0)
positions(-1)
( unit = 1 AU )
3 STO 19
2 STO 22
1.896536240 STO 28
x1
2 STO 20
0 STO 23
-0.149960004 STO 29
y1
1 STO 21
-1 STO 24
-0.998225293 STO 30
z1
0 STO 25
-0.101186449 STO 31
x2
4 STO 26
3.998868143 STO 32
y2
-1 STO 27
-1.049095980 STO 33
z2
>Evaluate the coordinates of the 2 stars when t = 10 days
n' = 2 STO 00 k2
= 17.20209895 E-3 X^2 STO 16
h = 5 days STO 17 N = 2 STO 18
XEQ "PLN2" >>>> x1
= 2.187016625 = R22
( Execution time = 3mn21s ) and
x2 = 0.195600654 = R25
RDN y1 = 0.299249451
= R23
y2 = 3.994996468 = R26
RDN z1 = -0.993676113
= R24
z2 = -0.896746353 = R27
-The accuracy is of the order of 10 -7 AU
-The positions at t = 5 days are in registers R28 thru R33
-Press R/S or XEQ 01 to continue ( after changing N in register R18
if needed )
-You can also press XEQ "PLN2" but it would needlessly re-calculate
values which are already stored in the required registers.
Variant:
-The Numerov's formula is applied recursively until 2 successive approximations
are equal ( line 135 )
-Alternatively, you can fix the number of iterations: store this
value in R14
Replace line 135
by DSE 15
Replace lines 123 to 127 by STO IND 01
Delete line 104
Add RCL 14 STO 15 after line 58
-In the example above, 3 iterations are enough and the execution time
becomes 2mn12s
Note:
-For a planet like Mercury with h = 1 day, the error
is of the order of 2.7 10 -5 AU
after 88 days
---- h = 0.5 day, -------------------------- 1.6
10 -6 AU -------------
-3 iterations ( in register R13 ) produce a good accuracy
for this planet with h = 0.5 or 1 day.
-When h is divided by 2, the errors are approximately divided by 16
= 24
3°) A Multistep Method of order 7
a) Inertial Frame of Reference
Data Registers: • R00 = n = number of bodies ( Registers R00 and R14 thru R16+13n are to be initialized before executing "GM3" )
R01 to R11: temp R12-R13: unused R16+7n thru R15+19n: temp
• R14 = G = k2
• R17+n = x1(0)
• R17+4n = x1(-1)
• R17+7n = x1(-2)
• R17+10n = x1(-3)
• R15 = h = stepsize
• R18+n = y1(0)
• R18+4n = y1(-1)
• R18+7n = y1(-2)
• R18+10n = y1(-3)
• R16 = N = number of steps
• R19+ n = z1(0)
• R19+4n = z1(-1)
• R19+7n = z1(-2)
• R19+10n = z1(-3)
• R17 = m1
................
.................
......................
......................
• R18 = m2
................
.................
......................
......................
.............
................
..................
.......................
.......................
• R16+n = mn
• R14+4n = xn(0)
• R14+7n = xn(-1)
• R14+10n = xn(-2)
• R14+13n = xn(-3)
• R15+4n = yn(0)
• R15+7n = yn(-1)
• R15+10n = yn(-2)
• R15+13n = yn(-3)
• R16+4n = zn(0)
• R16+7n = zn(-1)
• R16+10n = zn(-2)
• R16+13n = zn(-3)
where mi are the n masses,
xi(0) yi(0) zi(0)
are the positions at an instant t
xi(-1) yi(-1) zi(-1)
are the positions at the instant t-h
xi(-2) yi(-2) zi(-2)
are the positions at the instant t-2h
xi(-3) yi(-3) zi(-3)
are the positions at the instant t-3h
Flags: /
Subroutines: /
01 LBL "GM3"
02 RCL 00
03 13
04 *
05 16
06 +
07 STO 05
08 LASTX
09 RCL 00
10 10
11 *
12 +
13 E3
14 /
15 +
16 STO 04
17 RCL 00
18 31
19 *
20 16
21 +
22 STO 06
23 LASTX
24 RCL 00
25 +
26 .016
27 +
28 STO 10
29 XEQ 09
30 XEQ 09
31 XEQ 09
32 XEQ 05
33 RCL 06
34 1
35 +
36 RCL Y
37 E6
38 /
39 +
40 ENTER^
41 INT
42 R^
43 -
44 E3
45 /
46 +
47 REGMOVE
lines 01 thru 47 are executed the first time only
48 LBL 01
49 RCL 16
50 STO 11
51 LBL 02
52 1.004012
53 RCL 00
54 *
55 17.017
56 +
57 REGMOVE
58 RCL 00
59 3
60 *
61 GTO 00
62 LBL 03
63 RCL 00
64 4
65 *
66 16
67 +
68 STO 05
69 LASTX
70 RCL 00
71 +
72 E3
73 /
74 +
75 STO 04
76 RCL 00
77 19
78 *
79 16
80 +
81 STO 06
82 XEQ 05
83 LBL 00
84 STO Y
85 RCL 00
86 16
87 +
88 .1
89 %
90 +
91 +
92 STO 01
93 +
94 STO 02
95 +
96 +
97 STO 03
98 +
99 STO 04
100 +
101 STO 05
102 +
103 STO 06
104 +
105 STO 07
106 +
107 STO 08
108 +
109 STO 09
110 CLX
111 LBL 04
112 RCL IND 05
113 RCL IND 09
114 +
115 17
116 *
117 RCL IND 06
118 RCL IND 08
119 +
120 232
121 *
122 +
123 RCL IND 07
124 222
125 *
126 +
127 RCL 15
128 X^2
129 *
130 240
131 /
132 RCL IND 04
133 -
134 RCL IND 03
135 +
136 RCL IND 02
137 +
138 ENTER^
139 X<> IND 01
140 -
141 ABS
142 ST+ Y
143 SIGN
144 ST- 02
145 ST- 03
146 ST- 04
147 ST- 05
148 ST- 06
149 ST- 07
150 ST- 08
151 ST- 09
152 X<>Y
153 DSE 01
154 GTO 04
155 X#0?
156 GTO 03
a three-byte GTO
157 RCL 05
158 INT
159 RCL 06
160 INT
161 1
162 ST+ Z
163 +
164 E3
165 /
166 +
167 RCL 00
168 12
169 *
170 E6
171 /
172 +
173 REGMOVE
174 DSE 11
175 GTO 02
a three-byte GTO
176 RCL 01
177 3
178 +
179 RCL IND X
180 DSE Y
181 RCL IND Y
182 DSE Z
183 RCL IND Z
184 RTN
185 GTO 01
a three-byte GTO
186 LBL 05
187 CLX
188 STO 01
189 STO 02
190 STO 03
191 LBL 06
192 RCL 04
193 INT
194 RCL 05
195 X=Y?
196 GTO 08
197 RCL IND 10
198 RCL IND 05
199 RCL IND 04
200 -
201 STO 08
202 X^2
203 DSE 04
204 DSE 05
205 RCL IND 05
206 RCL IND 04
207 -
208 STO 09
209 X^2
210 +
211 DSE 04
212 DSE 05
213 RCL IND 05
214 RCL IND 04
215 -
216 STO T
217 X^2
218 +
219 ENTER^
200 SQRT
221 *
222 /
223 ST* 08
224 ST* 09
225 *
226 ST+ 01
227 RCL 09
228 ST+ 02
229 RCL 08
230 ST+ 03
231 2
232 ST+ 04
233 LBL 07
234 DSE 05
235 DSE 10
236 GTO 06
237 RCL 14
238 ST* 01
239 ST* 02
240 RCL 03
241 *
242 STO IND 06
243 DSE 06
244 RCL 02
245 STO IND 06
246 DSE 06
247 RCL 01
248 STO IND 06
249 DSE 06
250 2
251 ST- 04
252 RCL 00
253 ST+ 10
254 3
255 *
256 ST+ 05
257 DSE 04
258 GTO 05
259 RTN
260 LBL 08
261 2
262 ST- 05
263 GTO 07
264 LBL 09
265 XEQ 05
266 ST- 05
267 E3
268 /
269 ST- 04
270 END
( 409 bytes / SIZE 17+31n )
| STACK | INPUTS | OUTPUTS |
| Z | / | z1 |
| Y | / | y1 |
| X | / | x1 |
Example: n = 3 stars / m1 = 2 ; m2 = 1 ; m3 = 3 ( unit = 1 Solar mass ) We have the following coordinates:
t = 0
t = -5 days
t = -10 days
t = -15 days
masses
positions(0)
positions(-1)
positions(-2)
positions(-3)
( unit = 1 AU )
2 STO 17
2 STO 20
1.997888568 STO 29
1.991382737 STO 38
1.980240265 STO 47
x1
1 STO 18
0 STO 21
-0.149784693 STO 30
-0.298912394 STO 39
-0.446978169 STO 48
y1
3 STO 19
0 STO 22
0.001032468 STO 31
0.004296703 STO 40
0.010056489 STO 49
z1
0 STO 23
0.000165879 STO 32
0.000666440 STO 41
0.001508330 STO 50
x2
4 STO 24
3.999043454 STO 33
3.996203288 STO 42
3.991522280 STO 51
y2
0 STO 25
-0.049838219 STO 34
-0.099339682 STO 43
-0.148486062 STO 52
z2
0 STO 26
0.101352328 STO 35
0.205522696 STO 44
0.312670380 STO 53
x3
0 STO 27
0.000175311 STO 36
0.000540500 STO 45
0.000811352 STO 54
y3
1 STO 28
0.999257761 STO 37
0.996915425 STO 46
0.992791028 STO 55
z3
>Evaluate the coordinates of these 3 stars when t = 10 days
n = 3 STO 00 k2
= 17.20209895 E-3 X^2 STO 14
h = 5 days STO 15 N = 2 STO 16
XEQ "GM3" >>>> x1
= 1.992077585 = R20
( Execution time = 6mn21s ) and x2
= 0.000661670 = R23 x3 = -0.194938946
= R26
RDN y1 = 0.300333545
= R21
y2 = 3.996080575 = R24 y3
= 0.001084113 = R27
RDN z1 = 0.003673675
= R22
z2 = 0.100603412 = R25 z3
= 0.997349746 = R28
-The errors are of the order of 6 10 -9 AU
-The coordinates corresponding to t = 5 days , 0 , -5 days
are now in registers R29 to R37 , R38 to R46 , R47 to R55 respectively
-Press R/S or XEQ 01 to continue ( after changing N in register R16
if needed )
-You can also press XEQ "GM3" but it would needlessly re-calculate
values which are already stored in the required registers.
Variant:
-The formula is applied recursively until 2 successive approximations
are equal ( line 155 )
-To avoid unuseful computations of the complicated LBL 05, you can
fix the number of iterations:
Store this value in R12
Replace line 155
by DSE 13
Delete line 152
Replace lines 138 to 142 by STO IND 01 CLX
Delete line 110
Add RCL 12 STO 13 after line 57
Note:
-When h is divided by 2, the errors are approximately divided by 64
= 26
-However, roundoff errors will limit the attainable accuracy.
b) Heliocentric Coordinates
Data Registers: • R00 = n-1 = number of bodies minus 1 = n' ( Registers R00 and R16 thru R31n'+19 are to be initialized before executing "PLN3" )
R01 to R13: temp R14-R15: unused
• R16 = G = k2
• R20+n' = x1(0)
• R20+4n' = x1(-1)
• R20+7n' = x1(-2)
• R20+10n' = x1(-3)
• R17 = h = stepsize
• R21+n' = y1(0)
• R21+4n' = y1(-1)
• R21+7n' = y1(-2)
• R21+10n' = y1(-3)
• R18 = N = number of steps
• R22+ n' = z1(0)
• R22+4n' = z1(-1)
• R22+7n' = z1(-2)
• R22+10n' = z1(-3)
• R19 = m0 = mass of the origin
................
.................
....................
........................
• R20 = m1
................
.................
....................
.......................
.............
................
..................
.....................
........................
• R19+n' = mn'
• R17+4n' = xn'(0)
• R17+7n' = xn'(-1)
• R17+10n' = xn'(-2)
• R17+13n' = xn'(-3)
• R18+4n' = yn'(0)
• R18+7n' = yn'(-1)
• R18+10n' = yn'(-2)
• R18+13n' = yn'(-3)
• R19+4n' = zn'(0)
• R19+7n' = zn'(-1)
• R19+10n' = zn'(-2)
• R19+13n' = zn'(-3)
where mi are the n masses,
xi(0) yi(0) zi(0)
are the positions at an instant t
xi(-1) yi(-1) zi(-1)
are the positions at the instant t-h
( n' = n-1 )
xi(-2) yi(-2) zi(-2)
are the positions at the instant t-2h
xi(-3) yi(-3) zi(-3)
are the positions at the instant t-3h
R20+13n' thru R19+31n': temp
Flags: /
Subroutines: /
01 LBL "PLN3"
02 RCL 00
03 13
04 *
05 19
06 +
07 STO 08
08 LASTX
09 RCL 00
10 10
11 *
12 +
13 E3
14 /
15 +
16 STO 07
17 RCL 00
18 31
19 *
20 19
21 +
22 STO 09
23 LASTX
24 RCL 00
25 +
26 .019
27 +
28 STO 10
29 XEQ 11
30 XEQ 11
31 XEQ 11
32 XEQ 05
33 RCL 09
34 1
35 +
36 RCL Y
37 E6
38 /
39 +
40 ENTER^
41 INT
42 R^
43 -
44 E3
45 /
46 +
47 REGMOVE
lines 01 thru 47 are executed the first time only
48 LBL 01
49 RCL 18
50 STO 13
51 LBL 02
52 1.004012
53 RCL 00
54 *
55 20.02
56 +
57 REGMOVE
58 RCL 00
59 3
60 *
61 GTO 00
62 LBL 03
63 RCL 00
64 4
65 *
66 19
67 +
68 STO 08
69 LASTX
70 RCL 00
71 +
72 E3
73 /
74 +
75 STO 07
76 RCL 00
77 19
78 ST* Y
79 +
80 STO 09
81 XEQ 05
82 LBL 00
83 STO Y
84 RCL 00
85 19
86 +
87 .1
88 %
89 +
90 +
91 STO 01
92 +
93 STO 02
94 +
95 +
96 STO 03
97 +
98 STO 04
99 +
100 STO 05
101 +
102 STO 06
103 +
104 STO 07
105 +
106 STO 08
107 +
108 STO 09
109 CLX
110 LBL 04
111 RCL IND 05
112 RCL IND 09
113 +
114 17
115 *
116 RCL IND 06
117 RCL IND 08
118 +
119 232
120 *
121 +
122 RCL IND 07
123 222
124 *
125 +
126 RCL 17
127 X^2
128 *
129 240
130 /
131 RCL IND 04
132 -
133 RCL IND 03
134 +
135 RCL IND 02
136 +
137 ENTER^
138 X<> IND 01
139 -
140 ABS
141 ST+ Y
142 SIGN
143 ST- 02
144 ST- 03
145 ST- 04
146 ST- 05
147 ST- 06
148 ST- 07
149 ST- 08
150 ST- 09
151 X<>Y
152 DSE 01
153 GTO 04
154 X#0?
155 GTO 03
a three-byte GTO
156 RCL 05
157 INT
158 RCL 06
159 INT
160 1
161 ST+ Z
162 +
163 E3
164 /
165 +
166 RCL 00
167 12
168 *
169 E6
170 /
171 +
172 REGMOVE
173 DSE 13
174 GTO 02
a three-byte GTO
175 RCL 01
176 3
177 +
178 RCL IND X
179 DSE Y
180 RCL IND Y
181 DSE Z
182 RCL IND Z
183 RTN
184 GTO 01
a three-byte GTO
185 LBL 05
186 CLX
187 STO 04
188 STO 05
189 STO 06
190 LBL 06
191 RCL IND 10
192 RCL IND 08
193 STO 11
194 X^2
195 DSE 08
196 RCL IND 08
197 STO 12
198 X^2
199 +
200 DSE 08
201 RCL IND 08
202 STO T
203 X^2
204 +
205 ENTER^
206 SQRT
207 *
208 /
209 ST* 11
210 ST* 12
211 *
212 ST- 04
213 RCL 12
214 ST- 05
215 RCL 11
216 ST- 06
217 DSE 08
218 DSE 10
219 GTO 06
220 RCL 00
221 ST+ 10
222 3
223 *
224 ST+ 08
225 LBL 07
226 RCL 04
227 STO 01
228 RCL 05
229 STO 02
230 RCL 06
231 STO 03
232 RCL 19
233 RCL IND 07
234 STO 11
235 X^2
236 DSE 07
237 RCL IND 07
238 STO 12
239 X^2
240 +
241 DSE 07
242 RCL IND 07
243 STO T
244 X^2
245 +
246 ENTER^
247 SQRT
248 *
249 /
250 ST* 11
251 ST* 12
252 *
253 ST- 01
254 RCL 12
255 ST- 02
256 RCL 11
257 ST- 03
258 2
259 ST+ 07
260 LBL 08
261 RCL 07
262 INT
263 RCL 08
264 X=Y?
265 GTO 10
266 RCL IND 10
267 RCL IND 08
268 RCL IND 07
269 -
270 STO 11
271 X^2
272 DSE 07
273 DSE 08
274 RCL IND 08
275 RCL IND 07
276 -
277 STO 12
278 X^2
279 +
280 DSE 07
281 DSE 08
282 RCL IND 08
283 RCL IND 07
284 -
285 STO T
286 X^2
287 +
288 ENTER^
289 SQRT
290 *
291 /
292 ST* 11
293 ST* 12
294 *
295 ST+ 01
296 RCL 12
297 ST+ 02
298 RCL 11
299 ST+ 03
300 2
301 ST+ 07
302 LBL 09
303 DSE 08
304 DSE 10
305 GTO 08
306 RCL 16
307 ST* 01
308 ST* 02
309 RCL 03
310 *
311 STO IND 09
312 DSE 09
313 RCL 02
314 STO IND 09
315 DSE 09
316 RCL 01
317 STO IND 09
318 DSE 09
319 2
320 ST- 07
321 RCL 00
322 ST+ 10
323 3
324 *
325 ST+ 08
326 DSE 07
327 GTO 07
a three-byte GTO
328 RTN
329 LBL 10
330 2
331 ST- 08
332 GTO 09
333 LBL 11
334 XEQ 05
335 ST- 08
336 E3
337 /
338 ST- 07
339 END
( 511 bytes / SIZE 31n' + 20 )
( n' = n-1 )
| STACK | INPUTS | OUTPUTS |
| Z | / | z1 |
| Y | / | y1 |
| X | / | x1 |
Example: n = 3 stars / m0 = 3 ; m1 = 2 ; m2 = 1 ( unit = 1 Solar mass ) We have the following coordinates:
t = 0
t = -5 days
t = -10 days
t = -15 days
masses
positions(0)
positions(-1)
positions(-2)
positions(-3)
( unit = 1 AU )
3 STO 19
2 STO 22
1.896536240 STO 28
1.785860041 STO 34 1.667569885
STO 40
x1
2 STO 20
0 STO 23
-0.149960004 STO 29
-0.299452894 STO 35 -0.447789521
STO 41
y1
1 STO 21
-1 STO 24
-0.998225293 STO 30
-0.992618722 STO 36 -0.982734539
STO 42
z1
0 STO 25
-0.101186449 STO 31
-0.204856256 STO 37 -0.311162050
STO 43
x2
4 STO 26
3.998868143 STO 32
3.995662788 STO 38 3.990710928
STO 44
y2
-1 STO 27
-1.049095980 STO 33
-1.096255107 STO 39 -1.141277090
STO 45
z2
>Evaluate the coordinates of the 2 stars when t = 10 days
n' = 2 STO 00 k2
= 17.20209895 E-3 X^2 STO 16
h = 5 days STO 17 N = 2 STO 18
XEQ "PLN3" >>>> x1
= 2.187016531 = R22
( Execution time = 3mn52s ) and
x2 = 0.195600616 = R25
RDN y1 = 0.299249432
= R23
y2 = 3.994996461 = R26
RDN z1 = -0.993676071
= R24
z2 = -0.896746334 = R27
-The accuracy is of the order of 8 10 -9 AU
-The positions at t = 5 days are in registers R28 thru R33
------------------- 0 --- -------------- R34
thru R39
----------------- -5 ------------------
R40 thru R45
-Press R/S or XEQ 01 to continue ( after changing N in register R18
if needed )
-Execution time is then smaller since the lines before LBL 01 are no
more executed:
these lines contain 4 evaluations of the subroutine LBL 05.
-You could also press XEQ "PLN3" but it would needlessly re-calculate
values which are already stored in the proper registers.
Variant:
-The implicit formula is applied recursively until 2 successive approximations
are equal ( line 154 )
-Alternatively, you can fix the number of iterations: store this
value in R14
Replace line 154
by DSE 15
Delete line 151
Replace lines 137 to 142 by STO IND 01 CLX
SIGN
Delete line 109
Add RCL 14 STO 15 after line 57
-In the example above, 3 iterations are enough to produce the same accuracy
and the execution time is reduced to 2mn58s
-Unfortunately, we can't always know the number of required iterations
in advance!
Note:
-For a planet like Mercury, with h = 1 day, the error
is of the order of 3.6 10 -7 AU
after 88 days
---- h = 0.5 day, -------------------------- 5.8
10 -9 AU ------------- on an HP-48,
but 1.6 10 -8 AU ------------- on an HP-41 because
of roundoff errors.
-When h is divided by 2, the errors are approximately divided by 64
= 26
... if we don't take the roundoff errors into account!
Further Improvements?
-Other multistep methods of order 7 can be used: the formula
ym+1 = a.ym + b.ym-1 + c.ym-2 + d.ym-3 + h2.( A. fm+1 + B.fm +C.fm-1 + D.fm-2 + E.fm-3 )
will duplicate the Taylor polynomial up to the terms in h7 iff the 9 coefficients satisfy:
a = -16 + 240 E
A = E
E is undetermined
b = 34 - 480 E
B = 8/3 - 24 E
c = -16 + 240 E
C = 44/3 - 194 E
d = -1
D = 8/3 - 24 E
-The 4 programs above use E = 17/240 so that b = 0 , it yields:
ym+1 = ym + ym-2 - ym-3 + (h2/240).( 17. fm+1 + 232.fm + 222.fm-1 + 232.fm-2 + 17.fm-3 ) (F)
-Another possible choice is E = 5/72 which leads to
ym+1 = (2.ym + 2.ym-1 + 2.ym-2 )/3 - ym-3 + (h2/72).( 5. fm+1 + 72.fm +86.fm-1 + 72.fm-2 + 5.fm-3 ) (F')
-The accuracy is almost the same, perhaps slightly better in the above
examples but I don't think it's really significant.
-Other values are of course possible, but in order to obtain a stable
formula, the roots r of the polynomial x4
- a.x3 - b.x2 - c.x - d should verify
| r | <= 1 if r is a single root and | r | < 1 if r is a multiple zero.
-Unfortunately, this is impossible! 1 is always a double root
and the method of Numerov has the same defect.
-We can however satisfy these conditions for the other roots.
-Formula (F) - which is implicit - is computed with fm+1
= fm to get the first ym+1 approximation.
-A better predictor could be the formula given by E = 0
- which is explicit - and then formula (F) or (F') as
a corrector.
-It would probably reduce the number of iterations, thus saving execution
time, especially for large n-values.
Reference:
[1] "Numerical Analysis" - F. Scheid - Mc Graw Hill
ISBN 07-055197-9
Go back to the HP-41 software library
Go back to the general software library
Go
back to the main exhibit hall