Я пытаюсь написать ie некоторый код для симуляции системы из n-тел гравитационно связанных объектов. Я использую схему численного интегрирования Эйлера-Кромера. Когда алгоритм запускается, я знаю, что ожидается колебание скорости (кинетическая c энергия) относительно среднего значения, , однако, когда я выполняю код, я получаю net тенденцию к снижению скорости со временем , Как это можно исправить. Есть ли у кого-то опыт с такой проблемой?
Я считаю, что моя реализация верна, и у меня сложилось впечатление, что я пренебрегаю пренебрежимым эффектом в том, как python обрабатывает числа во время итераций (около 1 миллиона итераций).
График зависимости скорости от времени находится под кодом. Любые предложения относительно того, почему происходит это отклонение, приветствуются. Вот мой код:
import scipy as sp
import numpy as np
import math
from scipy import constants
import matplotlib.pyplot as plt
G = sp.constants.G
m_mars = 6.4185*(10**23)
m_phobos = 1.06*(10**16)
init_r_phobos = 9.3773*(10**6)
init_v_phobos = math.sqrt((G*m_mars)/init_r_phobos)
delta_t = 0.0005
r1 = []
r2 = []
v = []
r_phobos = np.array([init_r_phobos,0])
v_phobos = np.array([0,init_v_phobos])
t = np.arange(0,500,delta_t)
#----------------------------------------------------------------------------
for delta_t in t: ##The problem is here##
r_norm = np.linalg.norm(r_phobos)
a_phobos = -(G*m_mars*r_phobos/(r_norm**3))
v_phobos = v_phobos + (a_phobos*delta_t)
#----------------------------------------------------------------------------
r_phobos = r_phobos + (v_phobos*delta_t)
r1.append(r_phobos[0])
r2.append(r_phobos[1])
v.append(np.linalg.norm(v_phobos))
f, ax = plt.subplots()
points = ax.scatter(t,v)
Скорость в зависимости от графика времени: где большое число итераций и небольшой интервал времени делают колебания похожими на одну толстую линию. 
Здесь более длительное отклонение:
