Overview
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
b) Heliocentric Coordinates
3°) A Multistep Method of order 7 ( X-Functions Module required )
a) Inertial Frame of Reference
b) Heliocentric Coordinates
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/6
)
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 )
So, we have to solve the system:
d2xi /dt2 = SUMj#i
Gmj ( xj - xi ) / [ ( xi -
xj )2+ ( yi - yj )2
+ ( zi - zj )2 ]3/2
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
-Lines 83-127-190 are three-byte GTOs
| 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
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 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 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
-Lines 83-163-220 are three-byte GTOs
| 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 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 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 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: /
-Lines 156-166 are three-byte GTOs
| 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 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 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 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 2 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
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
-Execution time = 4mn58s
-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: /
-Lines 155-165-308 are three-byte GTOs
| 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 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 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 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 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
x2 = 0.195600654 = R25
RDN y1 = 0.299249451
= R23
y2 = 3.994996468 = R26
RDN z1 = -0.993676113
= R24
z2 = -0.896746353 = R27
-Execution time = 3mn21s
-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: /
-Lines 202-221-231 are three-byte GTOs
| 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 10 30 XEQ 10 31 XEQ 10 32 XEQ 06 33 LBL 01 34 RCL 16 35 STO 11 36 LBL 02 37 1.004012 38 RCL 00 39 * 40 17.017 41 + 42 REGMOVE 43 3 44 RCL 00 45 * 46 STO Y 47 LASTX 48 16 49 + 50 .1 51 % 52 + 53 + |
54 STO 01
55 + 56 STO 02 57 + 58 STO 03 59 + 60 STO 04 61 + 62 STO 05 63 + 64 + 65 STO 06 66 + 67 STO 07 68 + 69 STO 08 70 LBL 03 71 RCL IND 06 72 RCL IND 08 73 + 74 8 75 * 76 RCL IND 07 77 44 78 * 79 + 80 RCL 15 81 X^2 82 * 83 3 84 / 85 RCL IND 05 86 - 87 RCL IND 02 88 RCL IND 04 89 + 90 16 91 * 92 - 93 RCL IND 03 94 34 95 * 96 + 97 STO IND 01 98 CLX 99 SIGN 100 ST- 02 101 ST- 03 102 ST- 04 103 ST- 05 104 ST- 06 105 ST- 07 106 ST- 08 |
107 DSE 01
108 GTO 03 109 LBL 04 110 RCL 00 111 4 112 * 113 16 114 + 115 STO 05 116 LASTX 117 RCL 00 118 + 119 E3 120 / 121 + 122 STO 04 123 RCL 00 124 19 125 * 126 16 127 + 128 STO 06 129 XEQ 06 130 STO Y 131 RCL 00 132 16 133 + 134 .1 135 % 136 + 137 + 138 STO 01 139 + 140 STO 02 141 + 142 + 143 STO 03 144 + 145 STO 04 146 + 147 STO 05 148 + 149 STO 06 150 + 151 STO 07 152 + 153 STO 08 154 + 155 STO 09 156 CLX 157 LBL 05 158 RCL IND 05 159 RCL IND 09 |
160 +
161 17 162 * 163 RCL IND 06 164 RCL IND 08 165 + 166 232 167 * 168 + 169 RCL IND 07 170 222 171 * 172 + 173 RCL 15 174 X^2 175 * 176 240 177 / 178 RCL IND 04 179 - 180 RCL IND 03 181 + 182 RCL IND 02 183 + 184 ENTER^ 185 X<> IND 01 186 - 187 ABS 188 ST+ Y 189 SIGN 190 ST- 02 191 ST- 03 192 ST- 04 193 ST- 05 194 ST- 06 195 ST- 07 196 ST- 08 197 ST- 09 198 X<>Y 199 DSE 01 200 GTO 05 201 X#0? 202 GTO 04 203 RCL 05 204 INT 205 RCL 06 206 INT 207 1 208 ST+ Z 209 + 210 E3 211 / 212 + |
213 RCL 00
214 12 215 * 216 E6 217 / 218 + 219 REGMOVE 220 DSE 11 221 GTO 02 222 RCL 01 223 3 224 + 225 RCL IND X 226 DSE Y 227 RCL IND Y 228 DSE Z 229 RCL IND Z 230 RTN 231 GTO 01 232 LBL 06 233 CLX 234 STO 01 235 STO 02 236 STO 03 237 LBL 07 238 RCL 04 239 INT 240 RCL 05 241 X=Y? 242 GTO 09 243 RCL IND 10 244 RCL IND 05 245 RCL IND 04 246 - 247 STO 08 248 X^2 249 DSE 04 250 DSE 05 251 RCL IND 05 252 RCL IND 04 253 - 254 STO 09 255 X^2 256 + 257 DSE 04 258 DSE 05 259 RCL IND 05 260 RCL IND 04 261 - 262 STO T 263 X^2 264 + 265 ENTER^ |
266 SQRT
267 * 268 / 269 ST* 08 270 ST* 09 271 * 272 ST+ 01 273 RCL 09 274 ST+ 02 275 RCL 08 276 ST+ 03 277 2 278 ST+ 04 279 LBL 08 280 DSE 05 281 DSE 10 282 GTO 07 283 RCL 14 284 ST* 01 285 ST* 02 286 RCL 03 287 * 288 STO IND 06 289 DSE 06 290 RCL 02 291 STO IND 06 292 DSE 06 293 RCL 01 294 STO IND 06 295 DSE 06 296 2 297 ST- 04 298 RCL 00 299 ST+ 10 300 3 301 * 302 ST+ 05 303 DSE 04 304 GTO 06 305 RTN 306 LBL 09 307 2 308 ST- 05 309 GTO 08 310 LBL 10 311 XEQ 06 312 ST- 05 313 E3 314 / 315 ST- 04 316 END |
( 473 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
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
-Execution time = 5mn07s
-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 201 )
-To avoid unuseful computations of the complicated LBL 06, you can
fix the number of iterations:
Store this value in R12
Replace line 201
by DSE 13
Delete line 198
Replace lines 184 to 188 by STO IND 01 CLX
Delete line 156
Add RCL 12 STO 13 after line 108
-In the example above, 1 iteration is enough to produce the same accuracy ( to at least 9 decimals ) and the execution time is reduced to 2mn47s
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 R13n'+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: /
-Lines 201-220-230-373 are three-byte GTOs
| 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 12 30 XEQ 12 31 XEQ 12 32 XEQ 06 33 LBL 01 34 RCL 18 35 STO 13 36 LBL 02 37 1.004012 38 RCL 00 39 * 40 20.02 41 + 42 REGMOVE 43 3 44 RCL 00 45 * 46 STO Y 47 LASTX 48 19 49 + 50 .1 51 % 52 + 53 + 54 STO 01 55 + 56 STO 02 57 + 58 STO 03 59 + 60 STO 04 61 + 62 STO 05 63 + 64 + 65 STO 06 |
66 +
67 STO 07 68 + 69 STO 08 70 LBL 03 71 RCL IND 06 72 RCL IND 08 73 + 74 8 75 * 76 RCL IND 07 77 44 78 * 79 + 80 RCL 17 81 X^2 82 * 83 3 84 / 85 RCL IND 05 86 - 87 RCL IND 02 88 RCL IND 04 89 + 90 16 91 * 92 - 93 RCL IND 03 94 34 95 * 96 + 97 STO IND 01 98 CLX 99 SIGN 100 ST- 02 101 ST- 03 102 ST- 04 103 ST- 05 104 ST- 06 105 ST- 07 106 ST- 08 107 DSE 01 108 GTO 03 109 LBL 04 110 RCL 00 111 4 112 * 113 19 114 + 115 STO 08 116 LASTX 117 RCL 00 118 + 119 E3 120 / 121 + 122 STO 07 123 RCL 00 124 19 125 ST* Y 126 + 127 STO 09 128 XEQ 06 129 STO Y 130 RCL 00 |
131 19
132 + 133 .1 134 % 135 + 136 + 137 STO 01 138 + 139 STO 02 140 + 141 + 142 STO 03 143 + 144 STO 04 145 + 146 STO 05 147 + 148 STO 06 149 + 150 STO 07 151 + 152 STO 08 153 + 154 STO 09 155 CLX 156 LBL 05 157 RCL IND 05 158 RCL IND 09 159 + 160 17 161 * 162 RCL IND 06 163 RCL IND 08 164 + 165 232 166 * 167 + 168 RCL IND 07 169 222 170 * 171 + 172 RCL 17 173 X^2 174 * 175 240 176 / 177 RCL IND 04 178 - 179 RCL IND 03 180 + 181 RCL IND 02 182 + 183 ENTER^ 184 X<> IND 01 185 - 186 ABS 187 ST+ Y 188 SIGN 189 ST- 02 190 ST- 03 191 ST- 04 192 ST- 05 193 ST- 06 194 ST- 07 195 ST- 08 |
196 ST- 09
197 X<>Y 198 DSE 01 199 GTO 05 200 X#0? 201 GTO 04 202 RCL 05 203 INT 204 RCL 06 205 INT 206 1 207 ST+ Z 208 + 209 E3 210 / 211 + 212 RCL 00 213 12 214 * 215 E6 216 / 217 + 218 REGMOVE 219 DSE 13 220 GTO 02 221 RCL 01 222 3 223 + 224 RCL IND X 225 DSE Y 226 RCL IND Y 227 DSE Z 228 RCL IND Z 229 RTN 230 GTO 01 231 LBL 06 232 CLX 233 STO 04 234 STO 05 235 STO 06 236 LBL 07 237 RCL IND 10 238 RCL IND 08 239 STO 11 240 X^2 241 DSE 08 242 RCL IND 08 243 STO 12 244 X^2 245 + 246 DSE 08 247 RCL IND 08 248 STO T 249 X^2 250 + 251 ENTER^ 252 SQRT 253 * 254 / 255 ST* 11 256 ST* 12 257 * 258 ST- 04 259 RCL 12 260 ST- 05 |
261 RCL 11
262 ST- 06 263 DSE 08 264 DSE 10 265 GTO 07 266 RCL 00 267 ST+ 10 268 3 269 * 270 ST+ 08 271 LBL 08 272 RCL 04 273 STO 01 274 RCL 05 275 STO 02 276 RCL 06 277 STO 03 278 RCL 19 279 RCL IND 07 280 STO 11 281 X^2 282 DSE 07 283 RCL IND 07 284 STO 12 285 X^2 286 + 287 DSE 07 288 RCL IND 07 289 STO T 290 X^2 291 + 292 ENTER^ 293 SQRT 294 * 295 / 296 ST* 11 297 ST* 12 298 * 299 ST- 01 300 RCL 12 301 ST- 02 302 RCL 11 303 ST- 03 304 2 305 ST+ 07 306 LBL 09 307 RCL 07 308 INT 309 RCL 08 310 X=Y? 311 GTO 11 312 RCL IND 10 313 RCL IND 08 314 RCL IND 07 315 - 316 STO 11 317 X^2 318 DSE 07 319 DSE 08 320 RCL IND 08 321 RCL IND 07 322 - 323 STO 12 324 X^2 325 + |
326 DSE 07
327 DSE 08 328 RCL IND 08 329 RCL IND 07 330 - 331 STO T 332 X^2 333 + 334 ENTER^ 335 SQRT 336 * 337 / 338 ST* 11 339 ST* 12 340 * 341 ST+ 01 342 RCL 12 343 ST+ 02 344 RCL 11 345 ST+ 03 346 2 347 ST+ 07 348 LBL 10 349 DSE 08 350 DSE 10 351 GTO 09 352 RCL 16 353 ST* 01 354 ST* 02 355 RCL 03 356 * 357 STO IND 09 358 DSE 09 359 RCL 02 360 STO IND 09 361 DSE 09 362 RCL 01 363 STO IND 09 364 DSE 09 365 2 366 ST- 07 367 RCL 00 368 ST+ 10 369 3 370 * 371 ST+ 08 372 DSE 07 373 GTO 08 374 RTN 375 LBL 11 376 2 377 ST- 08 378 GTO 10 379 LBL 12 380 XEQ 06 381 ST- 08 382 E3 383 / 384 ST- 07 385 END |
( 576 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
x2 = 0.195600616 = R25
RDN y1 = 0.299249432
= R23
y2 = 3.994996461 = R26
RDN z1 = -0.993676071
= R24
z2 = -0.896746334 = R27
-Execution time = 2mn57s
-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 06.
-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 200 )
-Alternatively, you can fix the number of iterations: store this
value in R14
Replace line 200
by DSE 15
Delete line 197
Replace lines 183 to 188 by STO IND 01
CLX SIGN
Delete line 155
Add RCL 14 STO 15 after line 108
-In the example above, 1 iteration is enough to produce the same accuracy
and the execution time is reduced to 2mn07s
-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
with h = 0.5 day, the error is of the order
of 5.8 10 -9 AU after 88
days on an HP-48
but 1.6 10 -8 AU after
88 days 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!
Execution time:
-With n' = 9 ( i-e n = 10 - for instance the Sun and 9 planets )
each XEQ 06 requires 2m57s
-With one iteration ( in R14 ) , n' = 9 and N = 1, the first
execution lasts 16m16s whereas the next one lasts 4mn34s
-With 2 iterations ( in R14 ) , n' = 9 and N = 1,
the first execution lasts 20mn whereas the next
one lasts 8mn18s
Multistep Methods of order 7 for y" = f ( x , y ):
-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 2 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.
-The formula given by E = 0 - which is explicit but highly unstable - is used as a predictor to get the first ym+1 approximation:
ym+1 = -16.ym + 34.ym-1 - 16.ym-2 - ym-3 + (h2/3).( 8.fm + 44.fm-1 + 8.fm-2 )
-Formula (F) is then employed as a corrector.
Reference:
[1] "Numerical Analysis" - F. Scheid - Mc Graw Hill
ISBN 07-055197-9