Search code examples
pythonsympyderivative

Why can't sympy converge on a solution in Python?


I have an equation I'm trying to find where the derivative of this equation is equal to 1.

I use sympy's solve and subtract 1 from the equation to solve for 0, but it never converges on a solution which in this case should be x=260.806 for f'(x)=1.

Why can't it find this solution and what can I change without rewriting for something besides sympy? It throws no errors, just tries indefinitely.

import sympy as sym
from sympy import symbol, symbols, solve, init_printing, diff, lambdify, exp
import matplotlib.pyplot as plt
import numpy as np

x_sym = symbols('x')
init_printing(use_unicode=True)
y_sym_orig = x_sym*exp(0.4*(1-x_sym/550))
y_sym_deriv = diff(y_sym_orig, x_sym, 1)
print('Orig=', y_sym_orig)
print('Deriv=', y_sym_deriv)

y_orig = sym.lambdify(x_sym, y_sym_orig)
y_deriv = sym.lambdify(x_sym, y_sym_deriv)

x_sol=solve(y_sym_deriv-1, x_sym)
print('The derivative has y=1 at: ',x_sol)

plt.figure()
x1 = np.arange(0, 300, .1)
y_graph = y_orig(x1)
y_deriv = y_deriv(x1)
plt.ylabel('y')
plt.xlabel('x')
plt.grid(True, which='both')
plt.ylim(0,1)
plt.plot(x1, y_graph, 'r', label='Original')
plt.plot(x1, y_deriv, 'b', label='Derivitive')
plt.title('Original and Derivitives')
plt.legend()
plt.show()

Solution

  • SymPy strongly prefers rational numbers to floats. Add from sympy import Rational and replace 0.4 by Rational('0.4') in the formula.

    y_sym_orig = x_sym*exp(Rational('0.4')*(1-x_sym/550))
    y_sym_deriv = diff(y_sym_orig, x_sym, 1)
    x_sol = solve(y_sym_deriv-1, x_sym)
    print('The derivative has y=1 at: ', x_sol)
    

    prints

    The derivative has y=1 at:  [-1375*LambertW(exp(3/5)) + 1375]
    

    Notes:

    • Rational(2, 5) is another way to represent 2/5 in SymPy. Rational(0.4) or Rational(2/5) would not help: both versions create a Python float first, whose rational form is not 2/5, but rather 3602879701896397/9007199254740992.

    • Printing out the equation before solving it (as you did) is a good idea. Here the formula for the derivative is -0.00108496341646638*x*exp(-0.000727272727272727*x) + 1.49182469764127*exp(-0.000727272727272727*x) in the original form and -x*exp(-x/1375 + 2/5)/1375 + exp(-x/1375 + 2/5) in the rational form. Solving equations symbolically when they are full of floats is a tall or impossible task, considering how different the floating point arithmetics is from the rules of symbolic mathematics.

    • See Python numbers vs. SymPy Numbers for more.