hp41programs

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:   d²y/dt² = 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) + k₁/6 + 2k₂/6 )
 y'(t+h) = y'(t) + k₁/6 + 2k₂/3 + k₃/6                 where  k₁ = h.f (y)  ;  k₂ = h.f(y+h.y'/2+h.k₁/8)  ;  k₃ = h.f(y+h.y'+h.k₂/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 ):     y_m+1 = 2.y_m - y_m-1 + (h²/12).( f_m+1 + 10.f_m + f_m-1 )
-A multistep method of order 7:         y_m+1 = y_m + y_m-2 - y_m-3  + (h²/240).( 17. f_m+1 + 232.f_m + 222.f_m-1 + 232.f_m-2 + 17.f_m-3 )
 

1°) A Fourth-Order Runge-Kutta Method
 

   a)  Inertial Frame of Reference
 
 

Let   M₁( x₁, y₁,z₁) ; ....... ; M_n( x_n, y_n,z_n)  be n celestial bodies  with initial velocities   V₁( x'₁, y'₁,z'₁) ; ....... ; V_n( x'_n, y'_n,z'_n)  at an instant t
 and    m₁ , ...... , m_n   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:

              d²x_i /dt² = SUM_j#i   Gm_j ( x_j - x_i ) / [ ( x_i - x_j )²+ ( y_i - y_j )² + ( z_i - z_j )² ]^3/2
              d²y_i /dt² = SUM_j#i   Gm_j ( y_j - y_i ) / [ ( x_i - x_j )² + ( y_i - y_j )² + ( z_i - z_j )² ]^3/2          i , j = 1 , ..... , n
              d²z_i /dt² = SUM_j#i   Gm_j ( z_j - z_i ) / [ ( x_i - x_j )² + ( y_i - y_j )² + ( z_i - z_j )² ]^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 = x₁                    R27+4n = x'₁
      R24 = G = constant of gravitation                           R28+n = y₁                    R28+4n = y'₁
      R25 = h = step size                                                R29+n = z₁                     R29+4n = z'₁
      R26 = N = number of steps
                                                                                    ................                       ...................
      R27 = m₁
      R28 = m₂                                                               R24+4n = x_n                  R24+7n = x'_n
      ...............                                                               R25+4n = y_n                  R25+7n = y'_n
      R26+n =  m_n                                                          R26+4n = z_n                  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 m³ kg^-1 s^-2
 but if the units are Astronomical Unit , Solar mass and day , G = k²  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
 
 

 
    ( 375 bytes / SIZE 27+19n )
 

Example:

-Let's imagine a system of three stars      M₁ ( 2 ; 0 ;0 )       M₂ ( 0 ; 4 ; 0 )      M₃ ( 0 ; 0 ; 1 )                     unit = 1 AU
              with  initial velocities                 V₁ ( 0 ; 0.03 ; 0 )  V₂ ( 0 ; 0 ; 0.01 )  V₃ ( -0.02 ; 0 ; 0 )               unit = 1 AU per day
                     and masses                        m₁ = 2  ;  m₂ = 1  ;  m₃ = 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
 
 

 
 -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:
 

               d²x_i /dt² =  - Gm₀ x_i/r_i³ - SUM_j Gm_j x_j/r_j³ +  SUM_j#i   Gm_j ( x_j - x_i ) / r_ij³
               d²y_i /dt² =  - Gm₀ y_i/r_i³ - SUM_j Gm_j y_j/r_j³ +  SUM_j#i   Gm_j ( y_j - y_i ) / r_ij³                  i , j = 1 , ..... , n-1
               d²z_i /dt² =  - Gm₀ z_i/r_i³ - SUM_j Gm_j z_j/r_j³ +  SUM_j#i   Gm_j ( z_j - z_i ) / r_ij³
 

              where  r_i =  ( x_i² + y_i² + z_i² ) ^1/2   and   r_ij = [ ( x_i - x_j )² + ( y_i - y_j )² + ( z_i - z_j )² ] ^1/2
 

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

      R00 = n-1 = number of bodies minus 1                   R30+n = x₁                    R27+4n = x'₁
      R27 = G = constant of gravitation                           R31+n = y₁                    R28+4n = y'₁
      R28 = h = step size                                                R32+n = z₁                     R29+4n = z'₁
      R29 = N = number of steps
                                                                                    ................                       ...................
      R30 = m₀ = mass of the origin
      R31 = m₁                                                               R24+4n = x_n-1               R21+7n = x'_n-1
      ...............                                                               R25+4n = y_n-1               R22+7n = y'_n-1
                                                                                    R26+4n = z_n-1               R23+7n = z'_n-1
     R29+n = m_n-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 m³ kg^-1 s^-2
 but if units are Astronomical Unit , Solar mass and day , G = k²  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
 
 

 
      ( 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:

                   M₁ ( 2 ; 0 ; -1 )    M₂ ( 0 ; 4 ; -1 )   and  V₁ ( 0.02 ; 0.03 ; 0 )    V₂ ( 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
 
 

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

 
-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 = k²                                 •  R16+n = x₁⁽⁰⁾           •  R16+4n = x₁^(-1)           R16+7n thru R15+19n:  temp
                                      •  R14 = h = stepsize                          •  R17+n = y₁⁽⁰⁾           •  R17+4n = y₁^(-1)
                                      •  R15 = N = number of steps            •  R18+ n = z₁⁽⁰⁾           •  R18+4n = z₁^(-1)
                                      •  R16 = m₁                                           ................                     .................
                                      •  R17 = m₂                                           ................                     .................
                                          .............                                            ................                     ..................

                                      •  R15+n = m_n                                    •  R13+4n = x_n⁽⁰⁾           •  R13+7n = x_n^(-1)
                                                                                                 •  R14+4n = y_n⁽⁰⁾           •  R14+7n = y_n^(-1)
                                                                                                 •  R15+4n = z_n⁽⁰⁾           •  R15+7n = z_n^(-1)

   where  m_i  are the n masses,   x_i⁽⁰⁾ y_i⁽⁰⁾ z_i⁽⁰⁾   are the positions at an instant  t
                                    and    x_i^(-1) y_i^(-1) z_i^(-1)  are the positions at the instant  t-h

Flags: /
Subroutines: /

-Lines 156-166 are three-byte GTOs
 
 

 
    ( 362 bytes / SIZE 19n+16 )
 
 

 
Example:       n = 3 stars  /   m₁ = 2  ;  m₂ = 1  ;  m₃ = 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                       x₁
  1  STO 17            0  STO 20                  -0.149784693  STO 29                       y₁
  3  STO 18            0  STO 21                    0.001032468  STO 30                       z₁
                              0  STO 22                    0.000165879  STO 31                       x₂
                              4  STO 23                    3.999043454  STO 32                       y₂
                              0  STO 24                  -0.049838219   STO 33                       z₂
                              0  STO 25                    0.101352328  STO 34                       x₃
                              0  STO 26                    0.000175311  STO 35                       y₃
                              1  STO 27                    0.999257761  STO 36                       z₃

   >Evaluate the coordinates of these 3 stars when  t = 10 days

  n = 3  STO 00      k² = 17.20209895 E-3  X^2   STO 13    h = 5 days  STO 14   N = 2  STO 15
 

  XEQ "GM2"  >>>>      x₁ =  1.992077642  = R19                 x₂ = 0.000661670 = R22     x₃ = -0.194938984 = R25
                       RDN        y₁ =  0.300333555 = R20                  y₂ = 3.996080573 = R23     y₃ =  0.001084105 = R26
                       RDN        z₁ =  0.003673650 = R21                  z₂ = 0.100603410 = R24     z₃ =  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 = 2⁴
 

   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 = k²                                 •  R20+n' = x₁⁽⁰⁾           •  R20+4n' = x₁^(-1)              R20+7n' thru R19+19n':  temp
                                      •  R17 = h = stepsize                          •  R21+n' = y₁⁽⁰⁾           •  R21+4n' = y₁^(-1)
                                      •  R18 = N = number of steps            •  R22+ n' = z₁⁽⁰⁾           •  R22+4n' = z₁^(-1)
                                      •  R19 = m₀ = mass of the origin            ................                     .................
                                      •  R20 = m₁                                           ................                    .................                       ( n' = n-1 )
                                          .............                                            ................                     ..................

                                      •  R19+n' = m_n'                                    •  R17+4n' = x_n'⁽⁰⁾           •  R17+7n' = x_n'^(-1)
                                                                                                 •  R18+4n' = y_n'⁽⁰⁾           •  R18+7n' = y_n'^(-1)
                                                                                                 •  R19+4n' = z_n'⁽⁰⁾           •  R19+7n' = z_n'^(-1)

   where  m_i  are the n masses,   x_i⁽⁰⁾ y_i⁽⁰⁾ z_i⁽⁰⁾   are the positions at an instant  t
                                    and    x_i^(-1) y_i^(-1) z_i^(-1)  are the positions at the instant  t-h

Flags: /
Subroutines: /

-Lines 155-165-308 are three-byte GTOs
 
 

 
    ( 465 bytes / SIZE 19n' + 20 )     ( n' = n-1 )
 
 

 
Example:       n = 3 stars  /    m₀ = 3  ;  m₁ = 2  ;  m₂ = 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                       x₁
  2  STO 20            0  STO 23                  -0.149960004  STO 29                       y₁
  1  STO 21           -1  STO 24                  -0.998225293  STO 30                       z₁
                              0   STO 25                  -0.101186449  STO 31                       x₂
                              4  STO 26                    3.998868143   STO 32                       y₂
                             -1  STO 27                  -1.049095980   STO 33                       z₂
 

   >Evaluate the coordinates of the 2 stars when  t = 10 days

  n' = 2  STO 00      k² = 17.20209895 E-3  X^2   STO 16    h = 5 days  STO 17   N = 2  STO 18
 

  XEQ "PLN2"  >>>>      x₁ =  2.187016625  = R22                      x₂ =  0.195600654 = R25
                        RDN        y₁ =  0.299249451 = R23                       y₂ =  3.994996468 = R26
                        RDN        z₁ = -0.993676113 = R24                       z₂ = -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 = 2⁴
 

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 = k²                                 •  R17+n = x₁⁽⁰⁾           •  R17+4n = x₁^(-1)         •  R17+7n = x₁^(-2)       •  R17+10n = x₁^(-3)
                                      •  R15 = h = stepsize                          •  R18+n = y₁⁽⁰⁾           •  R18+4n = y₁^(-1)         •  R18+7n = y₁^(-2)       •  R18+10n = y₁^(-3)
                                      •  R16 = N = number of steps            •  R19+ n = z₁⁽⁰⁾           •  R19+4n = z₁^(-1)         •  R19+7n = z₁^(-2)       •  R19+10n = z₁^(-3)
                                      •  R17 = m₁                                           ................                     .................                   ......................             ......................
                                      •  R18 = m₂                                           ................                     .................                   ......................             ......................
                                          .............                                            ................                     ..................                  .......................            .......................

                                      •  R16+n = m_n                                    •  R14+4n = x_n⁽⁰⁾           •  R14+7n = x_n^(-1)      •  R14+10n = x_n^(-2)       •  R14+13n = x_n^(-3)
                                                                                                 •  R15+4n = y_n⁽⁰⁾           •  R15+7n = y_n^(-1)      •  R15+10n = y_n^(-2)       •  R15+13n = y_n^(-3)
                                                                                                 •  R16+4n = z_n⁽⁰⁾           •  R16+7n = z_n^(-1)       •  R16+10n = z_n^(-2)       •  R16+13n = z_n^(-3)

   where  m_i  are the n masses,

                x_i⁽⁰⁾ y_i⁽⁰⁾ z_i⁽⁰⁾   are the positions at an instant  t
                x_i^(-1) y_i^(-1) z_i^(-1)  are the positions at the instant  t-h
                x_i^(-2) y_i^(-2) z_i^(-2)  are the positions at the instant  t-2h
                x_i^(-3) y_i^(-3) z_i^(-3)  are the positions at the instant  t-3h

Flags: /
Subroutines: /

-Lines 202-221-231 are three-byte GTOs
 
 

 
    ( 473 bytes / SIZE 17+31n )
 
 

 
Example:     n = 3 stars  /   m₁ = 2  ;  m₂ = 1  ;  m₃ = 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                       x₁
  1  STO 18            0  STO 21                  -0.149784693  STO 30            -0.298912394  STO 39          -0.446978169  STO 48                      y₁
  3  STO 19            0  STO 22                    0.001032468  STO 31             0.004296703  STO 40           0.010056489  STO 49                       z₁
                              0  STO 23                    0.000165879  STO 32             0.000666440  STO 41            0.001508330  STO 50                      x₂
                              4  STO 24                    3.999043454  STO 33             3.996203288  STO 42            3.991522280  STO 51                      y₂
                              0  STO 25                  -0.049838219   STO 34           -0.099339682  STO 43          -0.148486062  STO 52                       z₂
                              0  STO 26                    0.101352328  STO 35             0.205522696  STO 44            0.312670380  STO 53                       x₃
                              0  STO 27                    0.000175311  STO 36             0.000540500  STO 45            0.000811352  STO 54                       y₃
                              1  STO 28                    0.999257761  STO 37             0.996915425  STO 46            0.992791028  STO 55                       z₃

   >Evaluate the coordinates of these 3 stars when  t = 10 days

  n = 3  STO 00      k² = 17.20209895 E-3  X^2   STO 14    h = 5 days  STO 15   N = 2  STO 16
 

  XEQ "GM3"  >>>>      x₁ =  1.992077585  = R20             x₂ = 0.000661670 = R23     x₃ = -0.194938946 = R26
                       RDN        y₁ =  0.300333545 = R21              y₂ = 3.996080575 = R24     y₃ =  0.001084113 = R27
                       RDN        z₁ =  0.003673675 = R22              z₂ = 0.100603412 = R25     z₃ =  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 = 2⁶
-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 = k²                                 •  R20+n' = x₁⁽⁰⁾           •  R20+4n' = x₁^(-1)         •  R20+7n' = x₁^(-2)       •  R20+10n' = x₁^(-3)
                                      •  R17 = h = stepsize                          •  R21+n' = y₁⁽⁰⁾           •  R21+4n' = y₁^(-1)         •  R21+7n' = y₁^(-2)       •  R21+10n' = y₁^(-3)
                                      •  R18 = N = number of steps            •  R22+ n' = z₁⁽⁰⁾           •  R22+4n' = z₁^(-1)         •  R22+7n' = z₁^(-2)       •  R22+10n' = z₁^(-3)
                                      •  R19 = m₀ = mass of the origin            ................                     .................                    ....................               ........................
                                      •  R20 = m₁                                           ................                     .................                   ....................                .......................
                                          .............                                            ................                     ..................                  .....................               ........................

                                      •  R19+n' = m_n'                                   •  R17+4n' = x_n'⁽⁰⁾         •  R17+7n' = x_n'^(-1)        •  R17+10n' = x_n'^(-2)      •  R17+13n' = x_n'^(-3)
                                                                                                 •  R18+4n' = y_n'⁽⁰⁾         •  R18+7n' = y_n'^(-1)        •  R18+10n' = y_n'^(-2)      •  R18+13n' = y_n'^(-3)
                                                                                                 •  R19+4n' = z_n'⁽⁰⁾         •  R19+7n' = z_n'^(-1)         •  R19+10n' = z_n'^(-2)      •  R19+13n' = z_n'^(-3)

   where  m_i  are the n masses,

                x_i⁽⁰⁾ y_i⁽⁰⁾ z_i⁽⁰⁾   are the positions at an instant  t
                x_i^(-1) y_i^(-1) z_i^(-1)  are the positions at the instant  t-h                                                        ( n' = n-1 )
                x_i^(-2) y_i^(-2) z_i^(-2)  are the positions at the instant  t-2h
                x_i^(-3) y_i^(-3) z_i^(-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
 
 

 
    ( 576 bytes / SIZE 31n' + 20 )     ( n' = n-1 )
 
 

 
Example:       n = 3 stars  /    m₀ = 3  ;  m₁ = 2  ;  m₂ = 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                x₁
  2  STO 20            0  STO 23                  -0.149960004  STO 29        -0.299452894  STO 35      -0.447789521  STO 41               y₁
  1  STO 21           -1  STO 24                  -0.998225293  STO 30        -0.992618722  STO 36     -0.982734539  STO 42                z₁
                              0   STO 25                  -0.101186449  STO 31        -0.204856256  STO 37     -0.311162050  STO 43                x₂
                              4  STO 26                    3.998868143   STO 32         3.995662788  STO 38       3.990710928  STO 44                y₂
                             -1  STO 27                  -1.049095980   STO 33        -1.096255107  STO 39     -1.141277090  STO 45                z₂
 

   >Evaluate the coordinates of the 2 stars when  t = 10 days

  n' = 2  STO 00      k² = 17.20209895 E-3  X^2   STO 16    h = 5 days  STO 17   N = 2  STO 18
 

  XEQ "PLN3"  >>>>      x₁ =  2.187016531  = R22               x₂ =  0.195600616 = R25
                        RDN        y₁ =  0.299249432 = R23                y₂ =  3.994996461 = R26
                        RDN        z₁ = -0.993676071 = R24                z₂ = -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 = 2⁶
  ... 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

    y_m+1 = a.y_m + b.y_m-1 + c.y_m-2 + d.y_m-3  + h².( A. f_m+1 + B.f_m +C.f_m-1 + D.f_m-2 + E.f_m-3 )

 will duplicate the Taylor polynomial up to the terms in h⁷ 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:

    y_m+1 = y_m + y_m-2 - y_m-3  + (h²/240).( 17. f_m+1 + 232.f_m + 222.f_m-1 + 232.f_m-2 + 17.f_m-3 )        (F)

-Another possible choice is  E = 5/72  which leads to

  y_m+1 = (2.y_m + 2.y_m-1 + 2.y_m-2 )/3 - y_m-3  + (h²/72).( 5. f_m+1 + 72.f_m +86.f_m-1 + 72.f_m-2 + 5.f_m-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    x⁴ - a.x³ - b.x² - 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  y_m+1 approximation:

   y_m+1 = -16.y_m + 34.y_m-1 - 16.y_m-2 - y_m-3  + (h²/3).(  8.f_m + 44.f_m-1 + 8.f_m-2 )

-Formula (F) is then employed as a corrector.
 

Reference:

[1]   "Numerical Analysis" - F. Scheid - Mc Graw Hill   ISBN 07-055197-9