Python - 使用 sympy 求解等于 0 的方程

Python - Solving equation equal to 0 using sympy

我正在尝试使用高度方程找出发射的火箭的最大高度。我已经设置了导数,所以现在我需要使用导数求解零。我试图求解 0 的方程式是

-0.0052t^3 + 4.26t + 0.000161534t^3.751 The related code is as follows

    def velocity(equation):
        time = Symbol('t')
        derivative = equation.diff(time)
        return derivative

    def max_height():
        time = Symbol('t')
        equ = 2.13 * (time ** 2) - 0.0013 * (time ** 4) + 0.000034 * (time ** 4.751)
        return solve(Eq(velocity(equ), 0))

    if __name__ == '__main__':
        t = Symbol('t')
        print(max_height())

我试过将直接方程式插入方程式,就像这样...

return solve(Eq(-0.0052t^3 + 4.26t + 0.000161534t^3.751, 0))

认为问题可能出在 return 类型的速度上,但这没有用。我也试过让它们 class 发挥作用,但这似乎也无济于事。

我得到的结果是它无限期地运行,直到我停止它。当我停止它时,出现以下错误

    Traceback (most recent call last):
  File "C:\Users\...\main.py", line 40, in <module>
    print(max_height())
  File "C:\Users\...\main.py", line 29, in max_height
    return solve(Eq(velocity(equ), 0))
  File "C:\Users\...\venv\lib\site-packages\sympy\solvers\solvers.py", line 1095, in solve
    solution = _solve(f[0], *symbols, **flags)
  File "C:\Users\...\venv\lib\site-packages\sympy\solvers\solvers.py", line 1675, in _solve
    u = unrad(f_num, symbol)
  File "C:\Users\...\venv\lib\site-packages\sympy\solvers\solvers.py", line 3517, in unrad
    neq = unrad(eq, *syms, **flags)
  File "C:\Users\...\venv\lib\site-packages\sympy\solvers\solvers.py", line 3300, in unrad
    eq = _mexpand(eq, recursive=True)
  File "C:\Users\...\venv\lib\site-packages\sympy\core\function.py", line 2837, in _mexpand
    was, expr = expr, expand_mul(expand_multinomial(expr))
  File "C:\Users\...\venv\lib\site-packages\sympy\core\function.py", line 2860, in expand_mul
    return sympify(expr).expand(deep=deep, mul=True, power_exp=False,
  File "C:\Users\...\venv\lib\site-packages\sympy\core\cache.py", line 72, in wrapper
    retval = cfunc(*args, **kwargs)
  File "C:\Users\...\venv\lib\site-packages\sympy\core\expr.py", line 3630, in expand
    expr, _ = Expr._expand_hint(
  File "C:\Users\...\venv\lib\site-packages\sympy\core\expr.py", line 3555, in _expand_hint
    arg, arghit = Expr._expand_hint(arg, hint, **hints)
  File "C:\Users\...\venv\lib\site-packages\sympy\core\expr.py", line 3563, in _expand_hint
    newexpr = getattr(expr, hint)(**hints)
  File "C:\...\venv\lib\site-packages\sympy\core\mul.py", line 936, in _eval_expand_mul
    n, d = fraction(expr)
  File "C:\...\venv\lib\site-packages\sympy\simplify\radsimp.py", line 1113, in fraction
    return Mul(*numer, evaluate=not exact), Mul(*denom, evaluate=not exact)
  File "C:\...\venv\lib\site-packages\sympy\core\cache.py", line 72, in wrapper
    retval = cfunc(*args, **kwargs)
  File "C:\...\venv\lib\site-packages\sympy\core\operations.py", line 85, in __new__
    c_part, nc_part, order_symbols = cls.flatten(args)
  File "C:\...\venv\lib\site-packages\sympy\core\mul.py", line 523, in flatten
    c_part.append(p)

如有任何帮助,我们将不胜感激。

有几个问题:

  • 您不能使用 ^ 求幂。 Python 你需要 **。参见 sympy - gotchas
  • 另一个问题是你有 3 个不同的变量 t。这应该只是一个变量。如果方程式中只有一个变量,equation.diff() 会自动使用该变量。如果有多个,您还需要将正确的变量传递给您的 velocity 函数。
  • 由于您的方程式使用浮点数,因此 sympy 会变得非常困惑,因为它试图找到精确的符号解,这在根据定义仅为近似值的浮点数中效果不佳。特别是指数中的浮动很难得到同情。为了应对,sympy 使用数值求解器 nsolve,但它需要一个种子值来开始数字运算。根据种子,获得 040.5087.55

更新后的代码如下所示:

from sympy import Symbol, Eq, nsolve

def velocity(equation):
    derivative = equation.diff()
    return derivative

def max_height():
    time = Symbol('t')
    equ = 2.13 * (time ** 2) - 0.0013 * (time ** 4) + 0.000034 * (time ** 4.751)
    return nsolve(Eq(velocity(equ), 0), 30)

print(max_height())

它可以帮助绘制绘图(使用 lambdify() 使函数在 matplotlib 中可访问):

from sympy import Symbol, Eq, nsolve, lambdify

def velocity(equation, time):
    derivative = equation.diff(time)
    return derivative

def get_equation(time):
    return 2.13 * (time ** 2) - 0.0013 * (time ** 4) + 0.000034 * (time ** 4.751)

def max_height(equ, time):
    return [nsolve(Eq(velocity(equ, time), 0), initial_guess) for initial_guess in [0, 30, 500]]

time = Symbol('t')
equ = get_equation(time)
max_heigths = max_height(equ, time)

equ_np = lambdify(time, equ)
vel_np = lambdify(time, velocity(equ, time))

import matplotlib.pyplot as plt
import numpy as np

xs = np.linspace(0, 105, 2000)
ys = equ_np(xs)
max_heigths = np.array(max_heigths)
plt.plot(xs, equ_np(xs), label='equation')
plt.plot(xs, vel_np(xs), label='velocity')
plt.axhline(0, ls='--', color='black')
plt.scatter(max_heigths, equ_np(max_heigths), s=100, color='red', alpha=0.5, label='extremes')
plt.legend()
plt.show()