Python - Solving equation equal to 0 using sympy

Question

I'm trying to find the max height of a rocket launched using the equation for the height. I have the derivative already set, so now I need to solve for zero using the derivative. The equation I'm trying to solve for 0 is

-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())

I tried inserting the direct equation into the Eq, like so...

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

thinking the problem might be with the return type of velocity, but that didn't work. I also tried playing around with making them class functions, but that didn't seem to help either.

The result I'm getting is that it runs indefinitely until I stop it. When I do stop it, I get the following errors

    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)

Any help would be greatly appreciated.

Answer 1

There are several problems:

• You can't use ^ for exponentiation. In Python you need **. See sympy - gotchas.
• Another problem is that you have 3 different variables t. This should be just one variable. If there is only one variable in an equation, equation.diff() automatically uses that one. If there are multiple, you also need to pass the correct variable to your velocity function.
• As your equations uses floats, sympy gets very confused as it tries to find exact symbolic solutions, which doesn't work well in floats which are by definition only approximations. Especially the float in the exponent is hard for sympy. To cope, sympy uses a numerical solver nsolve, but which needs a seed value to start number crunching. Depending on the seed, either 0, 40.50 or 87.55 is obtained.

Here is how the updated code could look like:

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())

It could help to draw a plot (using lambdify() to make the functions accessible in 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()