物理月球-地球-太阳系统

我没有仔细审查您的代码或测试其行为，但众所周知，欧拉方法在此类模拟中会获得能量，因此通常使用"梯形"
https://en.wikipedia.org/wiki/Trapezoidal_rule

如果我理解正确的话，你的方法是位置和速度之间的二阶关系（我假设是为了避免能量增益），但它似乎是加速度和速度之间关系的一阶。

你可能可以保持越级设计（而不是梯形），但在尝试这样做的过程中，你应该更好地使额外的状态明确，以消除二阶导数的使用。

编辑：我重读了维基百科上关于蛙跳的页面，重读了你的代码，很明显你没有使用蛙跳（正如我之前所说，因为你只将它应用于位置\速度，而不是速度\加速度）。请注意，该页面还向您展示了如何在加速度取决于位置而非速度的情况下简化蛙跳方程：https://en.wikipedia.org/wiki/Leapfrog_integration

如果你想继续使用显式集成，也许你可以尝试使用Runge-Kutta方法：https://en.wikipedia.org/wiki/Runge%E2%80%93Kutta_methodsRK4通常就足够了。