Visual Basic Program: "VB Mercury numerical"

Introduction and Purpose

The purpose of the program "VB Mercury numerical" is to calculate the numerical solutions of the analytical equations in chapter 15 of the book "Introducing Einsteins Relativity" by Ray D'Inverno.
Specific are studied the paragraph To get a copy select: VB Mercury numerical.zip


Program Operation

The program consists of 2 Forms: Control Form and Display Form
Mercury numerical
Picture 1A
  • Picture 1A shows the control display
    Without selection the picture shows the initial display.
    After selection the picture shows the situation after the Start button is selected.
  • The Control Form uses the following control parameters:
    Test , Sub Test , Eccentricity 1, l , m , k , rev 1 and Phi0
    This are the parameters shown in Yellow. The parameters only can be modified when they are in Yellow.

In order to control the excution of the program following control parameters are used:


Operation Test 1 - Select Start

After selecting Start the background color of 8 parameters will change in color Yellow. These fields represents the parameters that can be changed.
TestEcc 1lm Chvr vukepsSl R
1.1.71002.00714.1421 .09899.00699.994534.06099.99
1.2.51002.00514.1421 .07071.00500.994534.06099.99
1.3.31002.00314.1421 .04242.00300.994534.06099.99
1.4.71502.00717.3205 .09899.00698.994534.040149.99
1.5.72002.00720.0000 .08082.00466.994534.030199.99
1.6.71001.5.00712.2474 .08573.00699.994534.044100.00
1.7.71003.00717.3205 .12124.00699.994534.09099.99
The above table shows the calculated values


Test 1 Classical 15.7

                     ^2
    d^2R        d phi      mu
    ------- - R -----  = - --
     dt^2        dt        R^2 

       Equation 15.7
Forward Backward
Picture 1B
For i = 0 To imax 
  For j = 1 To di
    Select Case test
    Case Is = 1
     
      ' 15.7
      vdPhi = h / (r * r)           ' dphi/dt = 15.6
      ac = r * vdPhi * vdPhi - mu / (r * r)
      vr = vr + ac * ddt
      phi = phi + vdPhi * ddt
      r = r + vr * ddt
     
      TESTMAXR eccr

      If t2 > 0 Then t2 = 0
    End Select
  next j

    ' Display      

    xxx = r * Cos(phi)
    yyy = r * Sin(phi)
    px = px0 + fac * xxx
    py = py0 - fac * yyy
    Form2.PSet (px, py), QBColor(kleur)

next i 


Test 2 Classical 15.10 and 15.11

Test 1 calculates the equations 15.10 and 15.11.
Equation 15.10 is the same as equation 15.25 (See test 3) with the factor "3 * mu^2" being equal to zero.
    d^2u          mu 
    ------- + U = ---
    d phi^2       h^2

       Equation 15.10
Picture 1C


Test 3 Mercury 15.25

    d^2u           m 
    ------- + U = --- + 3 mu^2
    d phi^2       h^2

       Equation 15.25
Equation 15.25

For i = 0 To imax
  For j = 1 To di
    Select Case test
    Case Is = 2, 3
      '  15.25 page 196
      ac = m / h2 + C3 * m * u * u - u
      vu = vu + ac * ddPhi
      u = u + vu * ddPhi
      r = 1 / u
      phi = Phi00 * pi / 180 + i * dPhi + j * ddPhi
    
      TESTMAXR eccr
    
      If t2 > 0 Then t2 = 0
    End Select
  next j

    ' Display      

    xxx = r * Cos(phi)
    yyy = r * Sin(phi)
    px = px0 + fac * xxx
    py = py0 - fac * yyy
    Form2.PSet (px, py), QBColor(kleur)

next i 
Equation 15.25 is very easy to implement.
The function PSet stores one pixel on the display. The next equation 15.24 (and 15.23) are much more difficult.


Test 4 Mercury 15.24

      2        2      
   du      2  k -1   2m              3 
  ----  + U = ---- + --- * u + 2 * mu
  dphi         h^2   h^2

       Equation 15.24
      2    2      
   du     k -1         2m         2 
  ----  = ---- + U * (--- + 2 * mu  -u) 
  dphi     h^2        h^2

       Equation 15.24
Equation 15.24
For i = 0 To imax 
  For j = 1 To di
    Select Case test
    Case Is = 4
      '  15.24 page 196
      C2 = 2
      k1 = (k2 - 1) / h2           ' Eq 15.24
      u2 = u * u
      fac3 = k1 + u * (2 * m / h2 - u + C2 * m * u2)
      
      If fac3 < 0 And Count = 0 Then
         sign1 = sign1 * -1
         Count = 1
      End If
      If sign1 = 1 Then kleur = LGreen Else kleur = Yellow

      If Count > 0 Then
        Count = Count + 1:
        If Count >= maxcount Then Count = 0:
        If sign1 = 1 Then kleur = LRed Else kleur = Black
      End If
      vu = sign1 * Sqr(Abs(fac3))
      u = u + vu * ddPhi
      r = 1 / u
      phi = Phi00 * pi / 180 + i * dPhi + j * ddPhi
           
      TESTMAXR eccr
     
      If t2 > 0 Then t2 = 0
For more technical information about this program select: Description Test 4


Test 5 Mercury 15.23

     2
    r   dphi              
      * ---  = h
         dt
   Equation 15.21 

         -1         -1   2    2      2
k2*(1-2m)     (1-2m)   dr    r  dphi
      --    -    --   *--  -   *----    = 1
       r          r    dt        dt
                 Equation 15.23

      -1   2           -1   2      2
(1-2m)   dr   k2*(1-2m)    r  dphi
   --   *-- =       --   -   *----  - 1
    r    dt          r         dt
                 Equation 15.23
           2              2      2
         dr              r  dphi
 fac1  * -- = k2* fac1 -   *----  - 1
         dt                  dt
              Equation 15.23
                     2      2
                    r  dphi
fac2 = { k2* fac1 -   *----  - 1 } / fac1
                        dt
              Equation 15.23
For i = 0 To imax 
  For j = 1 To di
    Select Case test
    Case Is = 5
      ' 15.23   Page 196
      vdPhi = h / (r * r)              ' 15.21
      
      fac1 = 1 / (1 - 2 * m / r)
      fac2 = (k2 * fac1 - r * r * vdPhi * vdPhi - 1) / fac1
        
        If fac2 < 0 And Count = 0 Then
          sign1 = sign1 * -1
          Count = 1
        End If
        If fac2 < 0 Then
          fac2 = -fac2
        End If
        If Count > 0 Then
            Count = Count + 1:
            If Count = 20 Then Count = 0:
        End If
      v = sign1 * Sqr(fac2)
      phi = phi + vdPhi * ddt
      r = r - v * ddt
     
      TESTMAXR eccr

      If t2 > 0 Then t2 = 0

Equation 15.23
Test 5 is a combination of the equations 15.21, 15.22 and 15.23.
Equation 15.23 is the most difficult.
This equation contains the two derivatives dr/dt^2 and dhi/dt^2.
dhi/dt is calculated in equation 15.21.
The derivative dr/dt^2 is called fac2. That means that dr/dt = sqr(fac2)
The problem is that fac2 can become negative, which has to be taken into account.
Test 5 requires the following parameters:
Ecc (which is used to calculate h), l (which is used to calculate r0), m and k
The purpose of the following table is to show the calculated k value to obtain the requested Eccentricity. This "calculation" is done by try and error.
TestEcclm Chvr vukEpsilonPrecession
5.1.71002.00714.1421 .09899.00699.994534.06025.65
5.2.61002.00614.1421 .08485.00600.99318.06025.65
5.3.51002.00514.1421 .07071.00500.992045.05925.65
5.4.41002 .00414.1421 .05656.00400.991112.05925.65
5.5.31002 .00314.1421 .04242.00300.990412.06025.65
5.6.71003.00717.3205 .12124.00699.99145.09043.80
5.7.71004.00720.0000 .14000.00699.98802.12067.41
5.8.71001.00710.0000 .07000.00699.997301.03011.76


Test 6 Mercury 15.27

              m^2
    Eps = 3 * ---
              h^2  

  2
 d u           m          h^2* U^2
------- + U = --- + Eps * -------- 
d phi^2       h^2             m

            Equation 15.27
  2
 d u       m          h^2* U^2
------- = --- + Eps * -------- - u
d phi^2   h^2             m

            Equation 15.27
    
For i = 0 To imax 
  For j = 1 To di
    Select Case test
    Case Is = 6
     
      ' 15.27
      ac = m / h2 + Eps * h2 * u * u / m - u
      vu = vu + ac * ddPhi
      u = u + vu * ddPhi
     
      r = 1 / u
      phi = Phi00 * pi / 180 + i * dPhi + j * ddPhi
    
      TESTMAXR eccr
    
      If t2 > 0 Then t2 = 0
    End Select
  next j

    ' Display      

    xxx = r * Cos(phi)
    yyy = r * Sin(phi)
    px = px0 + fac * xxx
    py = py0 - fac * yyy
    Form2.PSet (px, py), QBColor(kleur)

next i 

Test 6 requires the following parameters:
Ecc (which is used to calculate h), l (which is used to calculate u0), m and vu


Test 7 Mercury 15.30

          m
    u0 = --- * (1 + e * cos(phi))
         h^2  

          m      m
    u0 = ---  + ---*e * cos(phi))
         h^2    h^2 


 d^2u          m          h^2* U^2
------- + U = --- + Eps * -------- 
d phi^2       h^2             m

            Equation 15.30

    
Equation 15.30
1 For i = 0 To imax 
2  For j = 1 To di
3   Select test
4   Case Is = 7
     
      ' 15.30
      
5     u0 = (m / h2) * (1 + e * Cos(phi))   ' Error
6     u0 = (m / h2) * (1 + C * Cos(phi))
7     ' u0 = m / h2 + C * Cos(phi)     ' Error
8     ac = (h2 / m) * u0 * u0 - u1
      vu = vu + ac * ddPhi
      u1 = u1 + vu * ddPhi
     
      r = 1 / u1
      phi = Phi00 * pi / 180 + i * dPhi + j * ddPhi

      TESTMAXR eccr
    
      If t2 > 0 Then t2 = 0
    End Select
  next j

    ' Display      

    xxx = r * Cos(phi)
    yyy = r * Sin(phi)
    px = px0 + fac * xxx
    py = py0 - fac * yyy
    Form2.PSet (px, py), QBColor(kleur)

next i 

Test 7 requires the following parameters:
Ecc (which is used to calculate h), l (which is used to calculate u0), m and vu
u0 is the initial value for u1 at t = 0
The numerical calculation causes a problem in line #5, because e (Eccentricity) is used. This value is too large. In stead the parameter C is used which is much smaller.
In Test 7 C is set equal to Ecc/l = 0.7/100 = 0.007
You can try a value for C = 0.001 and C=0.01 and observe the difference.


Sub Test 1 Classical 15.11

       mu 
   U = ---  + C * cos(Phi - Phi0)
       h^2

        Equation 15.11
Equation 15.11


Sub Test 2 Classical 15.12

    l     
    -  = 1 + Ecc * (cos (Phi - Phi0)) 
    R           

             Equation 15.12
Equation 15.12


Sub Test 3 Mercury 15.31

    m                    m*ecc        m*ecc^2
A =---*(1+0.5*ecc^2)  B= -----  C = - -------   
   h^2                    h^2           6*h^2

    u1 =  A +B*phi*sin(phi)+C*cos(2*phi)
           
          Equation 15.31
Equation 15.31

For i = 1 To imax
    phi = Phi0 + i * dPhi + pi / 2
    Select Case subtest
    Case Is = 3
       '  equation 15.31
       u0 = m / h2 * (1 + ecc * Cos(phi))
       u1 = aa + bb * phi * Sin(phi) + CC * Cos(2 * phi)
       u = u0 + Eps * u1
       r = 1 / u
    End Select
    
    TESTMAXR eccr
    
    If t2 > 0 Then t2 = 0
    ' Display    
    
    px = px0 + fac * r * Cos(phi):
    py = py0 - fac * r * Sin(phi)
    Form2.PSet (px, py), QBColor(kleur)

    If maxcount = irev Then GoTo end1
    If stop1 = 1 Then GoTo end2
Next i



SubTest 4 Mercury 15.33

           m                           d^2u       h^2     
    u0 =  --- * (1 + Ecc * cos(phi))   ---- + u = --- * u0^2
          h^2                          dt^2        m

                                           Equation 15.33
Equation 15.33


Sub Test 5 Mercury 15.34A

           m                             
    u0 =  --- * [1 + Ecc * cos(phi) + Eps*Ecc*phi*sin(phi)] 
          h^2                         

            Equation 15.34A
Picture 1B


Sub Test 6 Mercury 15.34B

           m                             
    u0 =  --- * {1 + Ecc * cos[phi(1-Eps)]}  
          h^2                       

              Equation 15.34B
Equation 15.34B m=1
Equation 15.34B m=2


Program Description

The program uses the following subroutines: Main and TESTMAXR. The subroutine Main is called when the Start Command is selected. The subroutine TESTMAXR is used to calculate the Perihelion (Shortest distance) and Aphelion (Farthest distance).
Using that information, after 1 complete revolutions the parameters a, ae, v1 and v2 are calculated and displayed on the Control Form


Reflection 1

The starting point of the the different equations is to consider Mercury as a test particle. This is a very simple approach when you want to calculate the movements of the planets in the Solar System Using General Relativity.
Of course you can claim that in order to perform such a simulation you can also use Newton's Law but than you will miss the precession of the perihelion.
A second issue is, and that is the most important issue, I want to investigate what is involved when you want to use GR. The main lesson is: It is already difficult just to perform a simplified example considering Mercury only. It will be extremely difficult when all planets are considerd, because all the inter actions involved

What makes numerical calculation difficult, specific for example equation 15.23 that the starting point should be to calculate the parameters and initial conditions using observations.
Specific you have to calculate the masses involved and the parameters h and k. The two parameters h and k are called constants, that may be true for each object (mass) involved, but in fact they are very important because they define the trajectory of each individual object. I expect you need as many parameter combinations (h,k) as objects involved In the above I have calculated h directly from eccentricity but in general you cannot follow such an approach.

Observations in GR are specific difficult because proper time is used. GR does not use the concept of simultaneity, which makes the whole issue of observation "handling" difficult.


Reflection 2

The program is based using two concepts: Tests and Subtests. Tests are the numerical implementation of differential (or difference) equation. Subtests are the implementation of analytical solutions. Analytical solutions are functions in phi.
All the implementations of tests (the differential equations 15.23, 15.24, 15.25 and 15.27) show the same result.
For the subtests this is different. That means the analytical solutions are not correct (except equation 15.34B). The same problem arises when you implement the algorithms show in the book "Astronomical Algorithms" by Jean Meeus. Using these algorithms (functions of t) you can calculate the position of the planet Mercury of any day in the future. However the accuracy is not enough to demonstrate the precession of the perihelion. It should be mentioned that nothing is wrong with this book: it is excellent. The same with the book "Astronomy with your personnal computer" by Peter Duffett-Smith.


Reflection 3

Paragaph 15.3 (Advance of the perihelion of Mercury) is an example of a one-body problem in general relativity. This is an extremely simple problem, because only one object (one mass) is considered. This is not even an realistic problem because any realistic problem studied should consists of at least 2 objects or masses.
A much more realistic example consists of at least 3 objects: A BH at the center of our Galaxy, one star and one planet. That means at least 2 moving objects or two moving gravitational fields
Such an example shows much better the issues involved when you calculate a numerical solution using GR.
In order to study the behavior of the planet Mercury also at least three objects should be involved.


Reflection 4

Numerical solutions require three steps. The first step involves observations. The second step consists of a pure mathematical process and the third step again involves observations, that means you have to test the result (predictions) of the mathematical process with actual observations in the future.
The mathematical process consists of a mathematical description of the physical process studied. In our case a collection of astronomical objects.
In the first step and in the third step, because humans are involved to do the observations, light signals are involved. That means the speed and the trajectories of the photons is an issue. In step the most important parameter is the speed of the gravitational radiation.
This distinction at first side seems strange, but in reality the physical laws that govern the movement of astronomical objects and photons have nothing in common. In fact you should study these laws "with your eyes closed". The only exception comes when some of the objects involved (gas clouds) are electrical charged. In that sense you can make a distinction between baryonic and non baryonic matter but not between visible (which emits photons) and non visible or darkmatter.


Feedback

None


Created:1 September 2016

Back to my home page: Contents of This Document