hp41programs

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

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/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 )
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
 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

-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