太阳系模拟中行星偏离轨道、日距增大问题排查求助

我正在尝试模拟太阳系，但行星并未沿轨道运行，反而与太阳的距离不断增大。我无法找出代码中的问题，以下是最小可复现示例（MRE）：

import numpy as np
import matplotlib.pyplot as plt


GRAVITATIONAL_CONSTANT = 6.674e-11
EARTH_MASS = 5.972 * (10**24)
SUN_MASS = 332954.355179 * EARTH_MASS
MERCURY_MASS = 0.06 * EARTH_MASS
AU2M = 1.495979e11
AUD2MS = AU2M / 86400


class SolarSystemBody:
    def __init__(self, name, mass, position, velocity):
        self.name = name
        self.mass = mass
        self.position = np.array(position, dtype=float) * AU2M
        self.velocity = np.array(velocity, dtype=float) * AUD2MS
        self.acceleration = np.zeros(3, dtype=float)

    def update_position(self, dt):
        self.position += self.velocity * dt + 0.5 * self.acceleration * dt**2

    def update_velocity(self, dt):
        self.velocity += self.acceleration * dt

    def gravitational_force(self, sun):
        r = sun.position - self.position
        distance = np.linalg.norm(r)
        direction = r / distance
        force_magnitude = GRAVITATIONAL_CONSTANT * self.mass * sun.mass / distance**2
        return force_magnitude * direction

    def calculate_acceleration(self, sun):
        force = self.gravitational_force(sun)
        self.acceleration = force / self.mass


mercury_position = [0.1269730114807624, 0.281031132701101, 0.01131924496141172]
mercury_velocity = [-0.03126433724097832, 0.01267637703164289, 0.00390363008183905]

sun = SolarSystemBody("Sun", SUN_MASS, [0, 0, 0], [0, 0, 0])
mercury = SolarSystemBody("Mercury", MERCURY_MASS, mercury_position, mercury_velocity)

dt = 3600 * 24
total_time = 365 * dt

pos = np.zeros((365, 3), dtype=float)
i = 0
for t in np.arange(0, total_time, dt):
    print(np.linalg.norm(sun.position - mercury.position))
    pos[i, :] = mercury.position
    mercury.calculate_acceleration(sun)
    mercury.update_velocity(dt)
    mercury.update_position(dt)
    i += 1


fig = plt.figure()
ax = fig.add_subplot(111, projection="3d")
ax.plot(pos[:, 0], pos[:, 1], pos[:, 2])
plt.show()

初始数据通过astroquery的HorizonsClass获取：

from astropy.time import Time
from astroquery.jplhorizons import Horizons
import numpy as np

sim_start_date = "2023-01-01"

data = []
planet_id = 199  # Mercury

obj = Horizons(id=planet_id, location="@sun", epochs=Time(sim_start_date).jd, id_type='id').vectors()
name = obj["targetname"].data[0].split('(')[0].strip()
r = [np.double(obj[xi]) for xi in ['x', 'y', 'z']]
v = [np.double(obj[vxi]) for vxi in ['vx', 'vy', 'vz']]

更新：
显然水星出现了利用太阳的**引力弹弓（slingshotting）**效应，导致其与太阳的距离突然增大。问题根源是dt步长过大，加速度无法及时更新。减小dt即可解决该问题：

dt = 360
total_time = 365 * 240 * dt

内容的提问来源于stack exchange，提问作者O_o