当两个物体都在运动时，如何求解牛顿引力方程

Question

我正在努力使用scipy.integrate.solve_ivp数值求解牛顿重力方程。当其中一个物体被假定为静止时，我可以让编码工作。然而，正如图中所示，当我使用下面的代码时，我试图解释太阳对木星的作用力和木星对太阳的作用力，太阳的轨道是不正确的。它似乎没有摇晃，但只是离开了COM

G = 6.67408e-11
msun = 1988500e24
m_jupiter = 1898.13e24

x_pos_j = 2.939022030646856E+00*au
y_pos_j = -4.169544536356319E+00*au
z_pos_j = -4.843813063988785E-02*au

x_vel_j = 6.083727195546033E-03*v_factor
y_vel_j = 4.708328571996785E-03*v_factor
z_vel_j = -1.556609498459863E-04*v_factor

def f_grav(t, y):
    x1, x2, x3, v1, v2, v3 = y
    # x1_star, x2_star, x3_star = x_other
    dydt = [v1,
            v2,
            v3,
            -(x1-x1_star)*G*m/((x1-x1_star)**2+(x2-x2_star)**2+(x3-x3_star)**2)**(3/2),
            -(x2-x2_star)*G*m/((x1-x1_star)**2+(x2-x2_star)**2+(x3-x3_star)**2)**(3/2),
            -(x3-x3_star)*G*m/((x1-x1_star)**2+(x2-x2_star)**2+(x3-x3_star)**2)**(3/2)]
    return dydt

x_jupiter = np.array([x_pos_j, y_pos_j, z_pos_j])
x_sun = np.array([0, 0, 0])
v_sun = np.array([0, 0, 0])

tStart = 0e0
t_End = 20*year
dt = t_End/3000
# minstep = dt/100
domain = (tStart, dt)
temp_end = dt

x1_star, x2_star, x3_star = [0, 0, 0]
init_j = [x_pos_j, y_pos_j, z_pos_j, x_vel_j, y_vel_j, z_vel_j]
init_sun = [20, -20, 20, 15, -15, 0]

while tStart < t_End:
    
    m=msun
    ans_j = solve_ivp(fun=f_grav, t_span=domain, y0=init_j)
    x1_star, x2_star, x3_star = ans_j['y'][0:3, -1]
    v1_star, v2_star, v3_star = ans_j['y'][3:6, -1]
    x_jupiter = np.vstack((x_jupiter, (ans_j['y'][0:3, -1])))
    init_j = [x1_star, x2_star, x3_star, v1_star, v2_star, v3_star]
    # print(init_j[0:3])
    
    m=m_jupiter
    ans_sun = solve_ivp(fun=f_grav, t_span=domain, y0=init_sun)
    x1_star, x2_star, x3_star = ans_sun['y'][0:3, -1]
    v1_star, v2_star, v3_star = ans_sun['y'][3:6, -1]
    
    v_sun = np.vstack((v_sun, (ans_sun['y'][3:6, -1])))
    x_sun = np.vstack((x_sun, (ans_sun['y'][0:3, -1])))
    init_sun = [x1_star, x2_star, x3_star, v1_star, v2_star, v3_star]
    
    tStart += dt
    temp_end = tStart + dt
    domain = (tStart,temp_end)

plt.plot(x_jupiter[:,0], x_jupiter[:,1])

plt.plot(x_sun[:,0], x_sun[:,1])

plt.show()

在任何给定时刻求解方程的代码中，我假设另一颗恒星是静止的，但只适用于dt，我认为这不会影响结果。这就是为什么我的图是错误的，如果是的话，当两个物体都在移动时，我怎样才能更好地解这个方程。木星的轨道（右图为蓝色）看起来不太正确，但与太阳的轨道不同（两张图上都是橙色）