Почему моя симуляция n-тела не производит орбиты (используя Землю и Венеру вокруг Солнца)?

Question

Я пытаюсь сделать симуляцию n-body в python. Я решил протестировать алгоритм в системе из трех тел, состоящей из Земли, Венеры и Солнца, в течение 10 лет с использованием временного шага в один день. Мой алгоритм будет работать до конца, но вывод не кажется правильным;Я ожидаю, что Земля и Венера будут иметь почти круглые пути вокруг Солнца;фактический результат показан ниже:

Я надеюсь, что кто-то может показать мне, где моя реализация неверна или моя логика неверна. Реализация ниже:

import numpy as np
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D # view 3d

class CelestialBody():

    def __init__(self, identifier, facecolor, mass, position, velocity=None, acceleration=None, notes=None):
        """ """
        self.identifier = identifier
        self.facecolor = facecolor
        self.mass = mass
        self.position = np.array(position, dtype=float)
        self.velocity = self.autocorrect_parameter(velocity)
        self.acceleration = self.autocorrect_parameter(acceleration)
        self._previous_positions = [self.position]
        self._previous_velocities = [self.velocity]
        self._previous_accelerations = [self.acceleration]
        self.notes = notes

    @property
    def previous_positions(self):
        return np.array(self._previous_positions)

    @property
    def previous_velocities(self):
        return np.array(self._previous_velocities)

    @property
    def previous_accelerations(self):
        return np.array(self._previous_accelerations)

    def autocorrect_parameter(self, args):
        """ """
        if args is None:
            return np.zeros(self.position.shape, dtype=float)
        else:
            return np.array(args, dtype=float)

Из класса, определенного выше, мы можем создать экземпляры Earth, Venus и Sun. Масса дана в единицах кг, позиции в метрах, скорость в м / с, ускорение в м / с ^ 2. Солнце изначально находится в начале координат;начальное положение Земли полностью в направлении оси x (со скоростью, направленной в направлении оси y, по касательной ...);и начальная позиция Венеры полностью в направлении y (с является скоростью, направленной в направлении x). Используемые значения получены и преобразованы через wolframalpha .

Sun = CelestialBody(identifier='Sun', facecolor='yellow', mass=1.988e30, position=[0, 0, 0])
Earth = CelestialBody(identifier='Earth', facecolor='blue', mass=5.972e24, position=[1.496e11, 0, 0], velocity=[0, 29766.5, 0])
Venus = CelestialBody(identifier='Venus', facecolor='darkorange', mass=4.867e24, position=[0, 1.082e11, 0], velocity=[34995.75, 0, 0])

Ниже приведен метод для запуска моделирования. Во-первых, цикл for используется для постоянного суммирования ускорения, вызванного силой массы j-th, на массу i-th на расстоянии, которое их разделяет (где i =/= j). Из чистого ускорения скорость получается умножением ускорения на временной шаг;положение рассчитывается путем умножения скорости на шаг по времени. Это происходит для каждого временного шага во внешнем цикле;для каждого временного шага вновь рассчитанные положение, скорость и ускорение добавляются в соответствующие атрибуты класса CelestialBody.

class NBodySimulation():

    def __init__(self, bodies, gravitational_constant=6.67e-11, t=0):
        """ """
        self.bodies = bodies
        self.gravitational_constant = gravitational_constant
        self.t = t
        self.collision_times = []

    def update_accelerations(self):
        """ """
        for ith_body, ibody in enumerate(self.bodies):
            ibody.acceleration *= 0
            for jth_body, jbody in enumerate(self.bodies):
                if ith_body != jth_body:
                    # dpos = ibody.position - jbody.position
                    dpos = jbody.position - ibody.position
                    r = np.sum(dpos**2)
                    acc = -1 * self.gravitational_constant * jbody.mass * dpos / np.sqrt(r**3)
                    ibody.acceleration += acc
            ibody._previous_accelerations.append(ibody.acceleration)
            self.bodies[ith_body] = ibody

    def update_velocities_and_positions(self, dt):
        """ """
        for ith_body, ibody in enumerate(self.bodies):
            ibody.velocity += ibody.acceleration * dt
            ibody.position += ibody.velocity * dt
            ibody._previous_velocities.append(ibody.velocity)
            ibody._previous_positions.append(ibody.position)
            self.bodies[ith_body] = ibody

    def run_simulation(self, dt, duration):
        """ """
        nsteps = int(duration / dt)
        for ith_step in range(nsteps):
            self.update_accelerations()
            self.update_velocities_and_positions(dt)
            self.t += dt


    def view_positions(self):
        """ """
        if self.bodies[0].position.size == 2:
            fig, ax = plt.subplots()
            ndim = 2
        elif self.bodies[0].position.size == 3:
            fig = plt.figure()
            ax = fig.add_subplot(1, 1, 1, projection='3d')
            ndim = 3
        npos = len(self.bodies[0].previous_positions)
        for ibody in self.bodies:
            if ndim == 2:
                ax.plot(ibody.previous_positions[0, :], ibody.previous_positions[1, :], color=ibody.facecolor, label=ibody.identifier, alpha=0.5)
            else:
                ax.plot(ibody.previous_positions[0, :], ibody.previous_positions[1, :], ibody.previous_positions[2, :], color=ibody.facecolor, label=ibody.identifier, alpha=0.5)
        # ax.view_init(azim=60, elev=45)
        fig.legend(mode='expand', loc='lower center',     ncol=len(self.bodies))
        plt.show()
        plt.close(fig)

Чтобы запустить код,

bodies = [Sun, Earth, Venus]
dt = 86400 # seconds in one day
duration = 315360000 # seconds in 10 years
NB = NBodySimulation(bodies)
NB.run_simulation(dt, duration)
NB.view_positions()

Мой вывод неверен из-за неправильной реализации или неисправной логики? Любая помощь будет оценена.