The Museum of HP Calculators

# The Gravitational n-body problem for the HP-41

This program is Copyright © 2002-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.

## 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 ( 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  o  s. 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