我正在尝试使用 matlab 编写一个 3 体问题。我得到了月球在太空旋转框架中的轨迹公式。它基本上是 ydotdot, xdotdot=GM/(x^2+y^2)^3/2 公式。公式是什么，并不重要。
我面临的问题是，我应该编写一个程序，以数值方式求解月球的轨迹方程。我使用 ODE45 进行比较，因为我的目标是获得与 ODE45 相同的结果。我的最终问题是，我想以天为单位迭代时间，所以 tspan= [0 365]。最重要的是当我将引力常数转换为秒然后执行 tspace= [0 365] 我得到完全不同的结果然后如果我要执行 [0 365*3600*34] 表示一年中的秒数并且 G= 6.67e -11。看来我的单位很奇怪。
我想知道是否有人可以解释为什么在我使用 ODE 45 时会发生这种情况。为什么我不能使用 ODE45 清楚地将秒转换为天？我需要做一个额外的步骤吗？我的问题中唯一的其他变量是半径、距离和 3 个物体的质量。
太感谢了。我已经为此工作了很长时间。任何帮助将非常感激。
沿每个轴的重力加速度公式是不正确的。
将质量 M_e 的地球置于原点，将月球（质量 M_m）置于 (x,y)。那么地月距离由下式给出：
R_em = sqrt(x^2 + y^2)
地月总力由下式给出：
F_em = G*M_e*M_m/R_em^2
地球重力引起的总加速度由下式给出：
a_em = F_em/M_m = G*M_e/R_em^2
并且指向原点。然后沿每个轴的加速度为：
xdotdot = -F_em*cos(theta) = -F_em*x/R_em = -G*M_e*x/R_em^(3/2)
ydotdot = -F_em*sin(theta) = -F_em*y/R_em = -G*M_e*y/R_em^(3/2)
请注意您所说的公式中缺少的 x 和 y 因子。
我不确定“将引力常数转换为秒”是什么意思。您用于 G 的值的单位为牛顿米^2/kg^2。所以它已经用 MKS（米-千克-秒）系统表示，使用这个值计算的加速度的单位是米/秒^2。
使用 (x_s, y_s) 处的第三个物体（例如，太阳），您计算一个新的 R_s 表示日月距离，并使用太阳的质量 M_s 计算新的加速度矢量（除了加速度现在在(x_s, y_s) 的方向，而不是 (0,0))。一旦将一切都放入一个共同的坐标系（这里是地心坐标——尽管日心坐标可能是一个更方便的选择，如果你在模拟日-地-月系统）。加上初始位置和速度，应该是您在下一个时间步计算位置和速度所需的全部内容。