The Precession of Mercury’s Perihelion (simulation)

One of the classical test of general relativity is the perihelion precession of Mercury (also gravitational lensing & gravitational redshift).

I1_1
Image source: http://www.relativity.li/en/epstein2/read/i0_en/i1_en/

Under Newtonian physics, a two-body system consisting of a lone object orbiting a spherical mass would trace out an ellipse with the spherical mass at a focus. The point of closest approach, called the periapsis (or, as the central body in our Solar System is the sun, perihelion), is fixed. A number of effects in our solar system cause the perihelia of planets to precess (rotate) around the sun. The principal cause is the presence of other planets which perturb each other’s orbit. Another (much less significant) effect is solar oblateness.

In 1859, the French mathematician and astronomer Urbain Le Verrier reported that the slow precession of Mercury’s orbit around the Sun could not be completely explained by “Newtonian mechanics” and “perturbations by the known planets“.

Observation: The perihelion precession of Mercury is 574.10±0.65 arcseconds per century relative to the inertial ICFR (International Celestial Reference Frame: barycenter of solar system).

Classic calculation: Newtonian mechanics, taking into account all the effects from the other planets, predicts a precession 531.63±0.69 arcseconds per century.

In the early 20th century, Albert Einstein’s General Theory of Relativity provided the explanation for the observed precession. The effect is very small: the Mercurian relativistic perihelion advance excess is just 42.98 arcseconds per century. Similar, but much smaller, effects operate for other planets: 8.62 arcseconds per century for Venus, 3.84 for Earth, 1.35 for Mars, and 10.05 for 1566 Icarus.

source: wikipedia hahahaha

Sooo…. Lets try this out! :) with numerical simulation..

First, of course the theoretical background, you can find the derivation of the General relativity correction in literature. :D ahaha
Ok, we can calculate that in gravitational system consist of many particles, the Post-Newtonian equation of motion is [ujug-ujug.. mak bedunduk]:
PN
then, become [klothaa..ak]:
(the derivation is a lot more complicated, LOL)
Post-newtonian
Source: Benitez and Gallardo, Celest Mech Dyn Astr 2008
I really want to calculate the gravity using above equation, but if we just want to know Mercury’s precession, that equation is overkill.. [halah]. (maybe next time, e.g. for binary Neutron star, globular cluster, or S-star around supermassive black hole :) )

Assume that only the central body give effect of general relativity, the general relativity correction given by this simple equation:
PN-correction

WOKE, lets do our simulation.. using Rebound integrator. Oh wait, there is no GR-correction in this integrator yet, so we can add it,
simply give “additional_force” like this:

void relativity_correction_force(){
    for (int i=1; i<N; i++){ // for all planet 
        struct particle* com = &(particles[i]);
    
        // relative to the Sun
        struct particle sun = particles[0];
        double dx = com->x - sun.x;
        double dy = com->y - sun.y;
        double dz = com->z - sun.z;
    
        double dvx = com->vx - sun.vx;
        double dvy = com->vy - sun.vy;
        double dvz = com->vz - sun.vz;

        double r = sqrt(dx*dx + dy*dy + dz*dz);
        double _r = 1./r; // reciprocal
    
        // Benitez and Gallardo 2008
        // alpha = k^2/(r^3*c^2); beta  = 4k^2/r - v.v; gamma = 4*(v.r)
        // delta_a = alpha*(beta*r + gamma*v)
        double alpha = G*_r*_r*_r*_c2;
        double beta  = 4.*G*_r - (dvx*dvx + dvy*dvy + dvz*dvz);
        double gamma = 4.*(dvx*dx + dvy*dy + dvz*dz);
    
        // add acceleration to comet
        com->ax += alpha * (beta*dx + gamma*dvx);
        com->ay += alpha * (beta*dy + gamma*dvy);
        com->az += alpha * (beta*dz + gamma*dvz);
    }
}

or similar with this code. :)

using this program and some initial data of planet from HORIZON JPL NASA we can make simulation of “Solar System” for example in 10000 years. We want to know how much the different of the Mercury’s orbit with and without GR correction, hopefully I can get the value given above. ahahaha

Result
Calculation:
Package Rebound
Integrator: IAS15 (similar with RADAU15)
time step: 0.1 days
integration time: 10000 years
including all planets, moon, and 3 massive asteroid.. LOL (why these asteroids here? because I’m too lazy to delete it)

computing time (OpenMP 8 threads)
with GR-correction 1011.007421 s
w/o (Classic) 948.855602 s

First of all lets remind you (if u already know) about Keplerian coordinate, hehe, its easier to look how the orbit change using this coordinate. :D
400px-Orbit1.svg

Mercury’s orbit change after 10000 years:
mercury_a

mercury_e

mercury_i

mercury_omega

mercury_Omega

Different between integration w/o and with GR-correction is clearly visible in argument of perihelion parameter.
mercury_omega2
assume that the other parameters are exactly the same (not change due to GR-correction) we can simply calculate the precession of Mercury’s perihelion using argument of perihelion. After 10000 years the different is 1.22397536º or 44.063 arcseconds/century. Uh.. not exactly as calculated before (by someone).. hahaha. I think 1 arcseconds is tolerable for this simulation :P
This different perhaps due to:
- that last calculation only using argument of perihelion
- need a smaller “integrator_epsilon” (error)
- the precession is not linear! (look in paper above)
- didn’t include correction from Oblateness of the Sun (quadrupole moment of gravity)

:P

Neville’s algorithm (polynomial interpolation and extrapolation)

Salah satu jenis interpolasi (dan ekstrapolasi) yang penting adalah interpolasi polinom Lagrange. Penurunannya dapat dilihat disini dan disitu (halah).

Interpolasi polinom tidak sering digunakan terutama pada kumpulan data yang besar, orang lebih sering menggunakan Cubic Spline untuk “memuluskan” data.. wkwkwk. Namun interpolasi polinom ini banyak digunakan dalam teknik integrasi juga penyelesaian ODE (ya iyalah).. salah satunya metode terkenal algoritma Bulirsch-Stoer-Gragg (ekstrapolasi ke h = 0) :D

Contoh program sederhana dalam C++
neville.cpp

/* Neville's algorithm
 * Copyleft (2014). Ridlo W. Wibowo
 *
 */

#include <iostream>
#include <stdio.h>
#include <stdlib.h>
#include <fstream>
using namespace std;


double x[100], y[100];
int N = -1;

double neville(double input);
void plot();

int main(){
    ifstream read("data.txt");
    while (!read.eof()){
        N += 1;
        read >> x[N] >> y[N];
        //cout << x[N] << " " << y[N] << endl; 
    }
    read.close();
    cout << "Jumlah data : " << N << endl;
    
    /* interpolate input "x" */
    double input[10] = {1.1, 1.6, 2.0, 2.3, 2.8, 3.5, 4.2, 4.5, 4.8, 5.2}; // last data is outside region
    
    /* y-interpolate and extrapolate */        
    ofstream outNeville("data.out");
    for (int i=0; i<10; i++){
        outNeville << input[i] << " " << neville(input[i]) << endl;
    }   
    
    plot();
    return 0;
}

double neville(double input){
    double Q[100][100];    
    // initiate first column    
    for (int i=0; i<N; i++){
        Q[i][0] = y[i];
    }
    
    // iteration
    for (int i=1; i<N; i++){
        for (int j=1; j<i; j++){
            Q[i][j] = ((input-x[i-j])*Q[i][j-1] - (input-x[i])*Q[i-1][j-1])/(x[i] - x[i-j]);
            //cout << Q[i][j] << " ";
        }
        //cout << endl;
    }

    return Q[N-1][N-2]; // output last iteration
}

void plot(){
    ofstream ploter("plot.in");
    ploter << "# gnuplot input file\n";
    ploter << "set xlabel \"x\"\n";
    ploter << "set ylabel \"f(x)\"\n";
    ploter << "plot \"data.txt\" t \"data input\", \"data.out\" t \"interpolation + last data as extrapolation\" \n";
    ploter.close();
    
    system("gnuplot -persist < plot.in");
    system("rm -f plot.in");
}

data.txt

0.9 1.3
1.3 1.5
1.9 1.85
2.1 2.1
2.6 2.6
3.0 2.7
3.9 2.4
4.4 2.15
4.7 2.05
5.0 2.1

neville_algorithm

:)

Tangential of radial vector in velocity direction

How to compute the tangential of radial vector in velocity direction in 3-dimension space?

Problem:
assume there is a particle move in 3D space relative to the origin, if radial vector from origin and velocity vector of particle are known, how to compute tangential of radial vector in velocity direction?

\vec{t} = (\vec{r} \times \vec{v}) \times \vec{r}

is it correct?
LOL

Belajar bahasa pemrograman

Belajar “bahasa” pemrograman baru semakin mudah dengan adanya beberapa website yang menawarkan pembelajaran yang interaktif dan menarik, rekomendasi dari saya antara lain:
http://www.codecademy.com/
https://www.khanacademy.org/computing/cs
http://code.org/

terkadang saya menjumpai anak-anak yang kaget dengan bahasa pemrograman dan akhirnya tidak menyukai dunia programming, padahal kedepannya untuk mempercepat pekerjaan kita membutuhkan bantuan komputer.
Apapun pekerjaan kita, tidak ada yang salah dengan belajar bahasa pemrograman, mengerti komputer dan teknologi saat ini bekerja dan yang paling penting berlatih berpikir secara runut dan efisien (layaknya belajar “bahasa matematika” ☞ kasusnya sama dengan anak-anak yang kemudian tidak menyukai matematika… ^^)