使用SymPy定义两体问题常微分方程时Eq返回False问题咨询

问题原因 #

你代码的问题出在导数定义的参数写错了：Derivative(x, t, 2)中第一个参数是未绑定自变量的函数对象x，而右侧表达式里用的是绑定了自变量的x(t)，二者是完全不同的SymPy对象，因此SymPy会直接判定这个等式不成立，返回False。

修正方案 #

修正时需要把导数作用对象改为x(t)，同时建议做两个调整避免后续运算出问题：

1. 引力为指向原点的向心力，表达式需要加负号，原来的正号对应的是斥力，无法得到椭圆轨道解
2. 不要用1.5这类浮点数作为幂次，替换为SymPy的有理数Rational(3,2)，避免符号运算时浮点数带来的精度问题或者逻辑判断错误

正确的方程定义代码如下：

from sympy import symbols, Function, Derivative, Eq, Rational
t, G, M = symbols('t G M', real = True)
x = Function('x')
y = Function('y')

# x方向运动微分方程
eq_x = Eq(Derivative(x(t), t, 2), - G * M * x(t) / (x(t)**2 + y(t)**2)**Rational(3,2))
# y方向运动微分方程
eq_y = Eq(Derivative(y(t), t, 2), - G * M * y(t) / (x(t)**2 + y(t)**2)**Rational(3,2))

后续求解建议 #

如果后续要解这套方程组，笛卡尔坐标下SymPy的dsolve很难直接给出闭式符号解，建议将坐标系转换为极坐标（定义径向距离r(t)和极角θ(t)），利用角动量守恒简化方程后再求解，会更容易得到符合开普勒定律的椭圆轨道解。如果需要数值模拟，可以将符号方程转换为数值函数后用scipy的常微分方程求解器计算。

内容的提问来源于stack exchange，提问作者Jörg Neulist