我正在尝试求解微分方程组并找到太阳和木星的轨迹。但我没有一个好的轨迹,只有一些点。 你能帮忙吗? (“Soleil”指的是太阳)
这是我的代码
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()
Eureka ,它有效!!!
最佳答案
目前,我在您的代码中发现以下问题,这些问题导致您的观察结果无法重现(除了缺少引用值之外):
初始数据中,长度单位为天文单位,时间单位为一天。万有引力常数的单位是
m^3 s^-1 kg^-2
,因此要将它们合并到一个公式中,您需要将 AU 转换为 m,因子约为150e+09
,其中的立方体要被分割出来。有一天是24*3600
秒,必须相乘。积分时间间隔也应该以年为单位计算,目前您似乎想到的是天,这是没有适当转换因子的第三个单位。 [已解决,un jour=一天]
从时间节点构造的划分来看,好像您使用的是 python 2。然后指数
3/2
评估为1
在整数除法中,可以直接使用1.5
在指数中,它是二进制浮点格式的精确值。 [实际上是 python 3,那么第一个除法应该显式为整数]在复制初始数据时,您犯了复制粘贴错误,初始位置和速度具有相同的数字,而实数应形成垂直向量。 [未在代码中解决,图像具有正确的速度]寻找适合您位置的在线数据,NASA HORIZON 系统为我提供了
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
对重心框架的标准化需要应用动量守恒,木星的质量足够大,仅减去速度可能会给出物理上错误的结果。 [未解决,初始数据应该已经是重心的,无需修正]
物理常数的不同精度也会引入误差,从而导致偏离引用位置。目前可见的最“脏”常数是引力常数和质量,其次是年份类型的不确定性。您只能获得任何(正确)计算结果的可靠前两位数字。
关于python - 求解系统微分方程(太阳和木星轨迹),我们在Stack Overflow上找到一个类似的问题: https://stackoverflow.com/questions/59518547/