sympy solve() 无法找到解

Question

我现在有一段比较基础的代码，目的是解决下面这个方程，但它一直没有找到答案，只是在那儿不停地闪烁一个“_”。

from sympy import *

gamma = 1.4
M_a = 1.0

y = Symbol('y', real=True) 
eqn = Eq((1.0/((gamma/2.0)*(y**2.0))) * ((((1.0 + ((gamma - 1.0)/2.0)*(y**2.0))/(1.0 + ((gamma - 1.0)/2.0)*(M_a**2.0)))**(gamma/(gamma - 1.0))) - 1.0) ,(-0.5704/((1.0 - (y**2.0))**(1.0/2.0))))
print solve(eqn, y)

简单地打印这个方程的结果是：

1.42857142857143*y**(-2.0)*((0.166666666666667*y**2.0 + 0.833333333333333)**3.5- 1.0) == -0.5704*(-y**2.0 + 1.0)**(-0.5)

把这个方程放到wolfram或者maple里计算，能得到正确的答案，大约是0.696256。

所以我在想，为什么sympy这个库不能解决这个方程。

这个方程应该长这样：图片

如果sympy不能用来解决这个方程，那我可以用什么呢？

谢谢！

Phil

Answer 1

如果你只关心数字解，那就别用 solve 了。它会努力去找一个符号解，但多项式的符号解通常很复杂，而且一般情况下根本不存在。

如果你只是想要一个数字解，可以用 nsolve：

In [60]: nsolve(eqn, 0)
Out[60]: mpc(real='0.69625557901731519', imag='-3.4211388289180104e-49')

（这里的0是对解的一个猜测）。结果中的虚部很小，可以忽略不计。你可以用 N(solution, chop=True) 来去掉它。

如果你确实需要符号解，另一个建议是避免使用浮点数的指数。要用有理数或整数的指数。为了找到多项式的符号解，首先必须把它转换成一个真正的多项式，也就是要有整数系数。这意味着浮点系数必须转换为有理数。但这样往往会导致很大的数字，从而产生非常高的次数（我怀疑这就是你遇到卡住的原因，考虑一下：

In [61]: gamma
Out[61]: 1.4

In [62]: Rational(gamma)
Out[62]:
3152519739159347
────────────────
2251799813685248