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
,但它需要一个种子值来开始数字运算。根据种子,获得 0
、40.50
或 87.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()
我正在尝试使用高度方程找出发射的火箭的最大高度。我已经设置了导数,所以现在我需要使用导数求解零。我试图求解 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
,但它需要一个种子值来开始数字运算。根据种子,获得0
、40.50
或87.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()