Simulation of a dynamic system of N particles under the gravitational force in 3 dimensions.
The system of N particles is descrived by a 3N coupled ODEs of second order and 3N initial conditions of position and velocity. The ODEs are determined using 2nd Newton's Law and Newton's law of universal gravitation. This system can be converted into 6N ODEs of first order. Assuming that the mass of each particle is constant, the 2 ODEs for each dimension of each particle are
where G is the gravitational constant (6.67408 × 10^-11 m^3 kg^-1 s^-2), v_i and r_i are the velocity and position vectors of the particle i, m_j is the mass of the particle j, |a| is the modulus of vector a, and t the time coordinate.
In order to remove the units and work with values close to 1, the following substitution as been used
where alpha is the cartesian coordinate of a vector (alpha = {x, y, z}, p.e. v_alpha = {v_x, v_y, v_z}), and r_A and v_A are the characteristical position and speed of the system (in this case the position and velocity in the aphelion). Then, the 6N ODEs system is transformed into the following expression:
The method used to calculate the evolution of the system is the Runge-Kutta 4 (RK4), which is an iterative and explicit method. The method applied to this system results in
with initial conditions
and
where the superindex (n) refers to the value of the variable in the n-iteration, and h is the time separation between two consequent iterations (h = (total time)/(number of iterations)).
Click here to see the RK4 code.
The estimation of errors has been done using the conservation of both energy and angular momentum. The error is estimated by comparing the energy and angular momentum's modulus of the system in the initial state and each iteration.
Click here to see the code for the estimation of errors.
Apart from the errors due to computation, the mathematical model of ODEs has also its approximations that should be taken into account:
-
TWO PARTICLES COLLIDING. The formulas descrived above do not descrive the effects of two particles colliding as they are considered as points. Once they get very close, the code cannot calculate accurately the next iteration and thus a jump in the error can be observed (function PLOT_ERROR).
-
RELATIVISTIC EFFECTS. If a particle is very close to another one of great mass, relativistic effects should be considered. This is the example of Mercury in the solar system (because it is very close to the Sun).
The main application of this code is centered on the solar system. Some examples, from simple to complex, are:
- Determine the orbit's variables of the planets (aphelion, perihelion, period...)
- Check the 2nd Law of Kepler
- Determine the mean anomaly of the planets
- Determine if an asteroid is NEO and/or PHO
In this webpage, it is possible to find data for the initial conditions. To get the cartesian coordinates of the objects, follow:
- Click on change next to Ephemeris Type and select Vector table. Then, click on Use Selection Above.
- Click on change next to Target Body and select the body from which you want the coordinates.
- Click on change next to Coordinate Origin and search "@0" to select Solar System Barycenter (SSB).
- Click on change next to Time span and select the date from which you want to start the simulation (all data from all bodies should be taken from the same date).
- Click on Generate Ephemeris
- Go to the line where your chosen date and X, Y, Z appear. For example, if 2019-Apr-18 is selected:
2458591.500000000 = A.D. 2019-Apr-18 00:00:00.0000 TDB
X =-2.089197670392227E-01 Y = 1.581041484991311E+00 Z = 3.802086917814641E-02
VX=-1.335337752542632E-02 VY=-6.379473197539670E-04 VZ= 3.142747436519678E-04
LT= 9.213328042903418E-03 RG= 1.595238299695790E+00 RR= 1.124040429122816E-03
- The cartesian coordinates corresponding to the initial conditions are the ones next to: X, Y, Z, VX, VY, VZ (in AU and AU/day).