求解系统微分方程(太阳和木星轨迹)

Question

我在试着解一个微分方程组，找出太阳和木星的轨迹。但我没有一个很好的轨迹，只有一些点。你能帮忙吗？(“太阳”指的是太阳)

​

​

这是我的代码

import numpy as np
import matplotlib.pyplot as plt
from scipy.integrate import odeint   
from mc_deriv import deriv

start = 0
end   = 14*365
nbpas = end/10
t = np.linspace(start,end,nbpas)
M = M_Soleil + M_Jupiter

x0    = x_Jupiter  - x_Soleil
y0    = y_Jupiter  - y_Soleil
vx0   = vx_Jupiter - vx_Soleil
vy0   = vy_Jupiter - vy_Soleil

syst_CI = [x0,y0,vx0,vy0]               
Sols=odeint(deriv,syst_CI,t,args=(M,))  

x  = Sols[:, 0]           
y  = Sols[:, 1]
vx = Sols[:, 2]
vy = Sols[:, 3]

初始化

x_Soleil    = -7.139143380212696e-03       # (UA)
y_Soleil    = -2.792019770161695e-03       # (UA)
x_Jupiter   = +3.996321311604079e+00       # (UA)
y_Jupiter   = +2.932561211517850e+00       # (UA)

vx_Soleil   = -7.139143380212696e-03       # (UA*j^-1)
vy_Soleil   = -2.792019770161695e-03       # (UA*j^-1)
vx_Jupiter  = +3.996321311604079e+00       # (UA*j^-1)
vy_Jupiter  = +2.932561211517850e+00       # (UA*j^-1)

M_Soleil    = 2e30                         # masse Soleil  (kg)
M_Jupiter   = 1.9e27                       # masse Jupiter (kg)
r_Soleil    = 696e6                        # rayon Soleil  (m)  

和外部函数

def deriv(syst,t,M):
    G     = 6.67e-11 
    x     = syst[0]
    y     = syst[1]
    vx    = syst[2]
    vy    = syst[3]
    dxdt  = vx
    dydt  = vy
    dvxdt = -(G*M*x)/((x**2+y**2)**(3/2))
    dvydt = -(G*M*y)/((x**2+y**2)**(3/2))
    return dxdt,dydt,dvxdt,dvydt

剧情

plt.figure(figsize=(7, 5))
plt.title("Trajectoires Soleil-Jupiter")
#plt.xlabel("UA)")
#plt.ylabel("UA)")
plt.plot(x,  y,  '-', color="red")
plt.show() 

​

​

​

​

绘图的结果是：

​

灵光乍现！

​

​

Answer 1

目前，我在您的代码中发现了以下问题，这些问题使您的观察结果无法重现(除了缺少参考值之外)：

• 在初始数据中，长度单位是天文单位，时间单位是一天。引力常数的单位是m^3 s^-1 kg^-2，所以要将它们组合在一个公式中，你需要将AU转换为m，该因子约为150e+09，其中的立方体将被分割。一天是24*3600秒，必须乘以，积分时间间隔也应该以年计算，在你认为天的那一刻，第三个单位没有适当的换算因子。解算后，un jour=one
• 从除法的时间节点中看起来就像是用了Python2。然后指数3/2在整数除法中求值为1，你可以直接在指数中使用1.5，它是一个二进制浮点格式的精确值。其实python 3，那么第一个除法应该是显式整数
• 在复制初始数据时犯了复制粘贴错误，初始位置和速度都有相同的数字，而实数应该形成垂直矢量。没有在代码中解决，图像有正确的速度寻找适合你位置的在线数据，美国宇航局的地平线系统给我2011-Nov-11 04:00的木星位置和速度为

pos: 3.996662712108880E+00，2.938301820497121E+00，-1.017177623308866E-01，vel：-4.560191659347578E-03，6.440946682361135E-03，7.529386668190383E-05

• The归一化到重心框架需要应用动量守恒，木星的质量足够大，仅仅减去速度可能会给出物理上错误的结果。如果不解决，初始数据应该已经是重心的，没有校正应该是necessary

• The的，改变物理常数的精度也会引入误差，从而导致偏离参考位置。目前可见的最“脏”的常数是引力常数和质量，之后是年份类型的不确定性。对于任何(正确的)计算结果，您只能获得可靠的前两位数。