# 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/6 )
y'(t+h) = y'(t) + k1/6 + 2k2/3 + k3/6                 where  k1 = h.f (y)  ;  k2 = h.f(y+h.y'/2+h.k1/8)  ;  k3 = h.f(y+h.y'+h.k2/2)

Alternatively, we can also use multisteps methods if we know 2 or more initial values ( at t , t-h , ... )  without knowing the speeds, for instance:

-Numerov's method ( of order 5 ):     ym+1 = 2.ym - ym-1 + (h2/12).( fm+1 + 10.fm + fm-1 )
-A multistep method of order 7:         ym+1 = ym + ym-2 - ym-3  + (h2/240).( 17. fm+1 + 232.fm + 222.fm-1 + 232.fm-2 + 17.fm-3 )

1°) A Fourth-Order Runge-Kutta Method

a)  Inertial Frame of Reference

Let   M1( x1, y1,z1) ; ....... ; Mn( xn, yn,zn)  be n celestial bodies  with initial velocities   V1( x'1, y'1,z'1) ; ....... ; Vn( x'n, y'n,z'n)  at an instant t
and    m1 , ...... , mn   the n masses of these bodies.
We want to know their future positions and velocities at an instant  t + N.h    ( h = step size ; N = number of steps )
So, we have to solve the system:

d2xi /dt2 = SUMj#i   Gmj ( xj - xi ) / [ ( xi - xj )2+ ( yi - yj )2 + ( zi - zj )2 ]3/2
d2yi /dt2 = SUMj#i   Gmj ( yj - yi ) / [ ( xi - xj )2 + ( yi - yj )2 + ( zi - zj )2 ]3/2          i , j = 1 , ..... , n
d2zi /dt2 = SUMj#i   Gmj ( zj - zi ) / [ ( xi - xj )2 + ( yi - yj )2 + ( zi - zj )2 ]3/2

where G is the gravitational constant.
It may seem very large but it's not too large for the HP-41:
If "GM" is executed directly from extended memory, it can solve the 15-body problem.

Data Registers:    The following registers are to be initialized before executing "GM":

R00 = n = number of bodies                                   R27+n = x1                    R27+4n = x'1
R24 = G = constant of gravitation                           R28+n = y1                    R28+4n = y'1
R25 = h = step size                                                R29+n = z1                     R29+4n = z'1
R26 = N = number of steps
................                       ...................
R27 = m1
R28 = m2                                                               R24+4n = xn                  R24+7n = x'n
...............                                                               R25+4n = yn                  R25+7n = y'n
R26+n =  mn                                                          R26+4n = zn                  R26+7n = z'n

>>> Thus, the n masses, the 3n position coordinates and the 3n velocity coordinates are to be stored into contiguous registers, starting at R27
>>> Registers R27+16n thru R26+19n must be cleared before the first execution  ( you can XEQ CLRG before storing numbers )

-Other registers are used for temporary data storage ( control numbers , partial sums , k-numbers of the Runge-Kutta scheme ... etc ... )

N.B.

G = 6.67259 10-11 m3 kg-1 s-2
but if the units are Astronomical Unit , Solar mass and day , G = k2  where  k = 0.01720209895 is Gaussian gravitational constant.

Program Listing

Note:    If you don't have an X-Function Module,

replace lines 191 to 204 by          and replace lines 18 to 34 by:

RCL 22                                              RCL 22
RCL 21                                              RCL 21
3                                                         4
*                                                         *
+                                                         +
RCL 22                                               RCL 22
RCL 21                                               RCL 21
ST+ X                                                .1
.1                                                         %
%                                                        +
+                                                         ST+ Y
ST+ Y                                                 LBL 02
LBL 11                                               CLX
CLX                                                    RCL IND Y
RCL IND Z                                         STO IND Z
STO IND Y                                        DSE Z
DSE Z                                                 DSE Y
DSE Y                                                GTO 02
GTO 11

-Lines 83-127-190 are three-byte GTOs

 01  LBL "GM"   02  LBL 00   03  RCL 00   04  26   05  +   06  .1   07  %   08  +   09  RCL 00   10  3   11  *   12  STO 21   13  +   14  STO 22   15  RCL 26           16  STO 23   17  LBL 01   18  RCL 22   19  INT   20  1   21  +   22  STO Y   23  RCL 00    24  9   25  *   26  +   27  RCL 21   28   E3   29  /   30  +   31   E3   32  /   33  +   34  REGMOVE   35  CLX   36  STO 07   37  STO 08   38  RCL 25   39  6   40  STO 10   41  / 42  STO 09   43  XEQ 04   44  RCL 25   45  2   46  ST* 09   47  /   48  STO 07   49  4   50  ST/ 10   51  /   52  RCL 25   53  *   54  STO 08   55  XEQ 04   56  RCL 25   57  STO 07   58  CLX   59  STO 09           60  4   61  ST* 08   62  ST* 10   63  XEQ 04   64  4   65  RCL 21    66  ST* Y   67  RCL 22   68  +   69  ST+ Y   70  ENTER^   71  LBL 03   72  CLX   73  X<> IND Z   74  RCL 25   75  *   76  ST+ IND Y   77  DSE Z   78  DSE Y   79  GTO 03   80  DSE 23   81  GTO 01   82  RTN 83  GTO 00   84  LBL 04   85  26.026   86  STO 11   87  3 E-5   88  RCL 22   89  +   90  STO 12   91  RCL 21   92  +   93  STO 14   94  LASTX   95  +   96  STO 16   97  LASTX   98  +   99  STO 18 100  LASTX 101  + 102  STO 19        103  LASTX 104  + 105  STO 20  106  LBL 05 107  RCL 00 108  ST+ 11 109  RCL 22 110  INT 111  STO 13 112  RCL 21 113  + 114  STO 15 115  LASTX  116  + 117  STO 17 118  CLX 119  STO 04 120  STO 05 121  STO 06 122  LBL 06 123  RCL 12 124  INT 125  RCL 13 126  X=Y? 127  GTO 08 128  XEQ 09 129  STO 03 130  1 131  ST- 12 132  ST- 13 133  ST- 14 134  ST- 15 135  ST- 16 136  ST- 17 137  XEQ 09 138  STO 02 139  1 140  ST- 12 141  ST- 13 142  ST- 14 143  ST- 15 144  ST- 16 145  ST- 17 146  RCL IND 11 147  XEQ 09 148  STO 01 149  X^2 150  RCL 02 151  X^2 152  RCL 03 153  X^2 154  + 155  + 156  ENTER^ 157  SQRT 158  * 159  / 160  ST* 01 161  ST* 02 162  ST* 03 163  RCL 01 164  ST+ 04 165  RCL 02 166  ST+ 05 167  RCL 03 168  ST+ 06 169  2 170  ST+ 12 171  ST+ 14 172  ST+16 173  SIGN 174  LBL 07 175  ST- 13 176  ST- 15 177  ST- 17 178  DSE 11 179  GTO 06 180  RCL 06 181  XEQ 10 182  RCL 05 183  XEQ 10 184  RCL 04 185  XEQ 10 186  3 187  ST- 14 188  ST- 16 189  DSE 12 190  GTO 05 191  RCL 22        192  RCL 21 193  ST+ X 194  1 195  + 196  .1 197  % 198  + 199  + 200  RCL 21 201   E6 202  / 203  + 204  REGMOVE 205  RTN 206  LBL 08 207  3 208  GTO 07 209  LBL 09 210  RCL IND 13 211  RCL IND 12 212  - 213  RCL IND 15 214  RCL IND 14 215  - 216  RCL 07 217  * 218  + 219  RCL IND 17 220  RCL IND 16 221  - 222  RCL 08 223  * 224  + 225  RTN 226  LBL 10 227  RCL 24 228  * 229  ENTER^ 230  STO IND 18 231  RCL 09 232  * 233  ST+ IND 19 234  CLX 235  RCL 10 236  / 237  ST+ IND 20 238  1 239  ST- 18 240  ST- 19 241  ST- 20 242  RTN 243  END

( 375 bytes / SIZE 27+19n )

Example:

-Let's imagine a system of three stars      M1 ( 2 ; 0 ;0 )       M2 ( 0 ; 4 ; 0 )      M3 ( 0 ; 0 ; 1 )                     unit = 1 AU
with  initial velocities                 V1 ( 0 ; 0.03 ; 0 )  V2 ( 0 ; 0 ; 0.01 )  V3 ( -0.02 ; 0 ; 0 )               unit = 1 AU per day
and masses                        m1 = 2  ;  m2 = 1  ;  m3 = 3                                                           unit = 1 Solar mass

-Calculate the positions and velocities 10 days later.

- SIZE 084 ( or greater )
XEQ CLRG  ( or  75.083 XEQ CLRGX if you have an HP-41 CX )

3     bodies                    >>>                                        3    STO 00
constant of gravitation    >>>    17.20209895 E-3    X^2    STO 24
step size   h = 10 days   >>>                                       10    STO 25
number of steps:    1      >>>                                         1    STO 26

then, we store the 3 masses and the 18 position and velocity coordinates as shown below  and  XEQ "GM"

>>>  the new positions and velocities have replaced the previous ones in registers R30 thru R47  ( next to last column )       Execution time = 1mn52s

 Register input h = 10 ; N=1 h = 5 ; N=2 m R27=m1 R28=m2 R29=m3 2    1    3 undisturbed undisturbed p  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