Simulating solar system with Newton's law

Question

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;
        
    
}

Answer 1

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

Answer 2

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.

Answer 3

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