Моделирование солнечной системы с помощью закона Ньютона

Я сделал симуляцию на С++ с законом Ньютона и протестировал ее, сравнивая положения планет с положением из Калькулятора Солнечной системы Дона Кросса (который я преобразовал из JavaScript в С++)
http://cosinekitty.com/solar_system.html

То, что я делаю, это каждый временной шаг (обычно 1 секунда, но шаг 0,2 секунды очень похож на шаг 10 секунд):

1. Рассчитать ускорение ($=$Ньютон Форзес$\times$дельта-время)
   
2. Скорость обновления и позиции
   
3. Сравните позиции с результатами алгоритмов солнечного калькулятора Don Cross.
   Но после 10 дней моделирования я получаю это отклонение расстояния (до калькулятора от Don Cross):
   

$\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}$

Как вы видите, у некоторых планет отклонения небольшие, а у некоторых большие, поэтому мой вопрос: могут ли Ньютоновы быть точными? или Don Cross Солнечной системы калькулятор не? Или в этом районе есть черная материя? Или что еще?

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