How to extract the function from sympy.diff?

I'm trying to calculate a derivative function using sympy.diff, and then solve the 'x' in this function letting the roots of derivative function func(x) = 0. However, solve this derivative function is very slow, because it returns five solutions, but I only need the solution most close to a fixed value x0.

import sympy

def diff_dist_func(a, b, c):
    x = sympy.Symbol('x')
    x0 = sympy.Symbol('x0')
    y0 = sympy.Symbol('y0')

    dist = ((x - x0) ** 2 + (a * x ** 3 + b * x ** 2 + c * x - y0) ** 2) ** (1 / 2)
    return sympy.diff(dist, x)

a = -0.00020129919480721813
b=0.10107634020780536
c=-12.305150031126267
shortest_dist = diff_dist_func(a, b, c)
    
x0=252.3007982720215
y0=96.55526056735049
solve_shortest_dist = shortest_dist.evalf(subs={'x0': x0, 'y0': y0})  # build the derivative function
solve_x = sympy.solve(solve_shortest_dist, sympy.Symbol('x'), simplify=False, rational=False)  # Here solve the derivative function is very slow.

To speed up the solution, I try to use scipy.optimize.fsolve, which is able to give a starting estimate x0 for the roots of func(x) = 0. So I replace sympy.solve with fsolve(solve_shortest_dist, np.array([x0])), but an error occured TypeError: 'Mul' object is not callable. How to extract the derivative function from the output of sympy.diff, enable it to be solved by scipy.optimize.fsolve? Or are there any way to speed up the solving process?

I'm using the SymPy 1.12, and solve is really fast, no need to move onto numerical libraries for now:

def diff_dist_func(a, b, c):
    x = sympy.Symbol('x')
    x0 = sympy.Symbol('x0')
    y0 = sympy.Symbol('y0')

    dist = sympy.sqrt((x - x0) ** 2 + (a * x ** 3 + b * x ** 2 + c * x - y0) ** 2)
    return sympy.diff(dist, x)

a = -0.00020129919480721813
b = 0.10107634020780536
c = -12.305150031126267
shortest_dist = diff_dist_func(a, b, c)

x0 = 252.3007982720215
y0 = 96.55526056735049
sols = sympy.solve(shortest_dist.subs({'x0': x0, 'y0': y0}), x)
sols
# out: [-5.99737132198727,
#  76.9370574126698,
#  252.303674956197,
#  256.81160524887 - 13.1689464528248*I,
#  256.81160524887 + 13.1689464528248*I]

Note the complex solutions. Let's compute the error wrt x0:

sols = np.array(sols, dtype=complex)
error = np.abs(sols - x0)
error
# out: array([2.58298170e+02, 1.75363741e+02, 2.87668418e-03, 1.39200765e+01,
       1.39200765e+01])

Finally, extract the closest solution to x0:

idx = np.argmin(error)
sols[idx]
# out: 252.303674956197

If you do need the speed and want a numerical solution, scipy.optimize.root (fsolve is a legacy function, so I prefer using the newer root function) can work on the result of the lambdify function. lambdify converts a sympy function into a numpy or scipy function (there are options to use other packages). Once converted, the scipy root-finder can work just fine.

from scipy.optimize import root

# your other code here

solve_shortest_dist_f = sympy.lambdify(sympy.Symbol("x"), solve_shortest_dist)
res = root(solve_shortest_dist_f, x0)
print(res.x[0]) # 252.30367495620214

The res variable contains some information about the solve, but res.x[0] is your root.

Timing this, it takes 6.01ms vs sympy solve which takes 22.2ms. Though, sympy is finding all the roots while scipy is just finding one.