# Simulating solar system with Newton's law

I made a simulation in C++ with Newtons law and test it comparing the planets positions with the position from Solar system Calculator Don Cross (which I converted from JavaScript to C++)
http://cosinekitty.com/solar_system.html

What I do is every time step(usually 1 second but step 0.2 second is very similar to 10 seconds step) :

1. Calculate accelaration ( $$=$$ newton forzes $$\times$$ deltatime)
2. Update speed and positions
3. Compare postions with results from Don Cross solar calculator alghorithms
But after 10 days of simulation I get this distance deviation (to the calculator from Don Cross) results :

$$\mathrm{Mercury} \ 4498.7 \ \mathrm{km}$$

$$\mathrm{Venus \ X} \ 1939.8 \ \mathrm{km}$$

$$\mathrm{Earth \ X} \ 10614.6 \ \mathrm{km}$$

$$\mathrm{Moon \ X} \ 7800.2 \ \mathrm{km}$$

$$\mathrm{Mars \ X} \ 445.2 \ \mathrm{km}$$

$$\mathrm{Ceres \ X} \ 129.5 \ \mathrm{km}$$

$$\mathrm{Pallas \ X} \ 432.4 \ \mathrm{km}$$

$$\mathrm{Juno \ X} \ 151.4 \ \mathrm{km}$$

$$\mathrm{Vesta \ X} \ 157.6 \ \mathrm{km}$$

$$\mathrm{Ida \ X} \ 73.6 \ \mathrm{km}$$

$$\mathrm{Gaspra} \ 455.3 \ \mathrm{km}$$

$$\mathrm{9P/T1} \ 241.5 \ \mathrm{km}$$

$$\mathrm{19P/B} \ 402.7 \ \mathrm{km}$$

$$\mathrm{67P/C-G} \ 533.2 \ \mathrm{km}$$

$$\mathrm{81P/W2} \ 110.7 \ \mathrm{km}$$

$$\mathrm{Jupiter} \ 172.3 \ \mathrm{km}$$

$$\mathrm{Saturn} \ 261.2 \ \mathrm{km}$$

$$\mathrm{Uranus} \ 71.4 \ \mathrm{km}$$

$$\mathrm{Neptune} \ 31.3 \ \mathrm{km}$$

$$\mathrm{Pluto \ X \ } \ 45.7 \ \mathrm{km}$$

As you see some planets have little desviations and some bigger, so my question is: Can Newton's be accurate? or Don Cross solar system calculator is not? Or there is black matter in that region? Or what else?

void CGravitator::CalcAceleration(double timeseconds){
unsigned int i,j,iend;
if (sunStatic)iend=m_np-1;
else iend=m_np;
for (i = 0; i < iend; i++) {
m_planetas[i].aceleration.set(0,0,0);
CVector3 totalGravitationalForce;
// Loop through all bodies in the universe to calculate their gravitational pull on the object (TODO: Ignore very far away objects for better performance)

for (j = 0; j < m_np; j++) {
if (i == j) continue; // No need to calculate the gravitational pull of a body on itself as it will always be 0.
double distancia =CVector3::Distancia(m_planetas[i].pos,m_planetas[j].pos);
double force = KGNEWTON * m_planetas[i].masa * m_planetas[j].masa / pow(distancia, 2);
CVector3 forceDirection = CVector3::Normalize(m_planetas[j].pos - m_planetas[i].pos);

totalGravitationalForce += forceDirection * force;
}
CVector3 incspeed = totalGravitationalForce / m_planetas[i].masa ;
m_planetas[i].aceleration += incspeed * timeseconds;

}

• Are you calculating all the forces between the bodies, or just the force between each body and the Sun? Commented Nov 29, 2021 at 4:23
• BTW, it's pretty easy to change your code from Euler integration to a symplectic integrator like Verlet or Leapfrog, and you'll be able to use a much larger time step. Commented Nov 29, 2021 at 4:26
• "Can Newton's be accurate? or Don Cross solar system calculator is not?" When writing new simulation code it's always good to benchmark against something more well established. But when the results differ, the default assumption should be that the fault is with your code.
– jkej
Commented Nov 29, 2021 at 12:29
• This isn't Code Review so I won't make this an answer, but a few tips for numerical stability and performance. #1 You're first multiplying and then dividing by m_planetas[i].masa. Cancel those out. (And change "force" in your variable names to "acceleration".) #2 You're first (presumably) taking a square root inside CVector3::Distancia and then squaring the distance. Cancel those out too. And #3 yes, avoid Euler integration. And #4 please include your actual velocity and position update code in your question — it might have bugs or numerical inaccuracies that affect your results. Commented Nov 30, 2021 at 0:19
• One way that you can cross-check your program is to calculate the total momentum of all objects in your system. This shouldn’t change, so any variation represents some kind of error. You’ll always have some, but it should be very small. Commented Nov 30, 2021 at 1:16

You need to use a better numerical method. Euler’s method is notoriously bad for orbital mechanics because the numerical errors always accumulate. In particular, Euler’s method does not conserve energy, so you get orbits that just magically gain energy and spiral away out of control.

You need to use something like the Verlet method or some other symplectic integrator.

https://en.wikipedia.org/wiki/Verlet_integration

https://en.wikipedia.org/wiki/Symplectic_integrator

• Thanks, but I used The Verlet code from wikipedia and got the same results Commented Nov 29, 2021 at 23:57
• I would recommend using the JPL ephemerides for checking your computations: ssd.jpl.nasa.gov/horizons
– Dale
Commented Nov 30, 2021 at 0:27
• Thanks, I tryed parsing Horizon files but got bigger errors!. Also its minimun time step is 0.5 seconds. Since some planets use the same code(in Don Cross source), and the ones with bigger errors have smaller constats, I will try with a precission library like gmp, but it is tedious. Commented Nov 30, 2021 at 0:56
• I doubt that gmp will help. You will just have the inaccurate positions to high precision.
– Dale
Commented Nov 30, 2021 at 2:01
• You don't need to use a symplectic integrator. JPL does not use one, for example. They use a variable step size, variable order Adams family integrator. Maybe Cowell or Gauss-Jackson? This are in the Adams family of numerical integrator, and both that position and velocity differently. Commented Nov 30, 2021 at 15:12

Besides the deficiencies in the numerical method pointed out in the other answer, simulating Mercury's orbit must take into account Jupiter's gravitation, as well as relativistic effects. This only explains a tiny fraction of the deviation but should be considered when the numerical method gets better.

Apparently the perihelion precession of Mercury due to Jupiter's pull and relativistic effects is about 574 arcseconds in a century, or 1.57E-2 arcseconds/day. With an orbital circumference of about 3.6E8 km and an arcsecond being 1.296E-6 of a turn, that amounts to about 4.3km, or 43 km in 10 days.

While the difference in the position of the perihelion does not directly translate into a location difference it should give an idea of the effect.

• thanks, this give me an idea of the error magnitude Commented Dec 1, 2021 at 19:10

The error came from the data sources. Now I replaced to NAIF cspice.lib and download a general SPK file for only the planets, and the error is really small without other gravity bodies like asteorids and comets. Is really simple to use only 2 functions( furnsh_c and spkezr_c), an include and a lib.
Coded with vc6++ pure win 32. I runned makeall.bat from spice unzipped.(couldn' make it work in MFC)
Also time step can be any fraction of seconds so spice and newton are synchronized. So Newton tell's me that there is no black matter.
Here are the results without asteroids and comets gravity for 10 days and no corrections in cpsice spkezr_c

distance error in KM newton(Verlet) from spice
time = 864000.250
MERCURY BARYCENTER = 5.511496277969209e+002
SATURN BARYCENTER = 8.535731413118873e+001
VENUS BARYCENTER = 2.701394194074592e+002
URANUS BARYCENTER = 8.651056255887706e+001
EARTH = 9.664941717935676e+001
NEPTUNE BARYCENTER = 8.654038466254323e+001
MOON = 1.208560265740111e+002
MARS BARYCENTER = 5.954829440293592e+001
PLUTO BARYCENTER = 8.661640570487361e+001
JUPITER BARYCENTER = 8.256645190275238e+001
TOTAL ERROR DISTANCE = 1.525933904320919e+003
Time Ini GREGORIAN = 2000 JAN 01 12:00:00.000
Time End GREGORIAN = 2000 JAN 11 12:00:00.250