I am struggling to numerically solve Newton's equation for gravity using scipy.integrate.solve_ivp. I can get the coding working when one of the bodies is assumed stationary. However as shown in the picture when I use the code below where I am trying to account for the force of the sun on Jupiter and that of Jupiter on the sun, the suns orbit is not correct. It doesn't seem to wobble but instead just heads off away from the 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()
In the code at any given moment to solve the equation I assume that the other star is stationary but only for dt, which I don't think should effect the results. Is this why my plots are wrong and if so how can I better solve the equation when both bodies are moving. The orbit of Jupiter (blue on right image) looks more less correct but not that of the sun (on both images and is orange).
