Finding Roots Numerically

Often times you might need to find when a certain equation is true, or solve for the roots of an ugly function. Sometimes these equations are straight up impossible to solve analytically. Enter numeric solutions!

%matplotlib inline
import matplotlib.pyplot as plt
import seaborn as sbn
import sympy as sym
import numpy as np
sym.init_printing()

Say we had a function like below. There exists a cube root equation, but nobody likes to use it because it is awful:

x = sym.symbols('x',real=True)
expr = sym.pi*x**2 -x**3 - 3*sym.exp(1)

Plotting it looks like below:

xs = np.arange(-2,3,0.01)
plt.plot(xs, sym.lambdify(x,expr)(xs))
plt.show()

Say we want to knew when the function crosses the x-axis, or when $y=0$. For some functions, I could make Sympy solve the system:

sym.solve(expr,x)

$\displaystyle \left[ \frac{\pi}{3} - \frac{\left(- \frac{1}{2} + \frac{\sqrt{3} i}{2}\right) \sqrt[3]{- \pi^{3} + \frac{\sqrt{- 4 \pi^{6} + \left(- 2 \pi^{3} + 81 e\right)^{2}}}{2} + \frac{81 e}{2}}}{3} - \frac{\pi^{2}}{3 \left(- \frac{1}{2} + \frac{\sqrt{3} i}{2}\right) \sqrt[3]{- \pi^{3} + \frac{\sqrt{- 4 \pi^{6} + \left(- 2 \pi^{3} + 81 e\right)^{2}}}{2} + \frac{81 e}{2}}}, \ \frac{\pi}{3} - \frac{\pi^{2}}{3 \left(- \frac{1}{2} - \frac{\sqrt{3} i}{2}\right) \sqrt[3]{- \pi^{3} + \frac{\sqrt{- 4 \pi^{6} + \left(- 2 \pi^{3} + 81 e\right)^{2}}}{2} + \frac{81 e}{2}}} - \frac{\left(- \frac{1}{2} - \frac{\sqrt{3} i}{2}\right) \sqrt[3]{- \pi^{3} + \frac{\sqrt{- 4 \pi^{6} + \left(- 2 \pi^{3} + 81 e\right)^{2}}}{2} + \frac{81 e}{2}}}{3}, \ - \frac{\sqrt[3]{- \pi^{3} + \frac{\sqrt{- 4 \pi^{6} + \left(- 2 \pi^{3} + 81 e\right)^{2}}}{2} + \frac{81 e}{2}}}{3} - \frac{\pi^{2}}{3 \sqrt[3]{- \pi^{3} + \frac{\sqrt{- 4 \pi^{6} + \left(- 2 \pi^{3} + 81 e\right)^{2}}}{2} + \frac{81 e}{2}}} + \frac{\pi}{3}\right]$

Sometimes you can indeed look at these numerically, though you need to pull the individual items out of the list. So if I wanted to look at the first solution:

sym.solve(expr,x)[2].evalf()

$\displaystyle -1.34776929418083$

In general, Sympy is focused on symbolic manipulations in this way. If all you care about is just the final number, you can actually ask Sympy to do so numerically using sym.nsolve:

sym.nsolve(expr,x,1)

$\displaystyle -1.34776929418083$

The third parameter here always has to be an initial guess for the value of the root. Numeric solutions like this will always return only a single answer! So in situations where you may have multiple roots, guessing near the root you want is important!

Scipy's Fsolve

While Sympy does have this numeric solver, in general it is better suited for symbolic mathematics. If Sympy is our symbolic Python library, Scipy and Numpy are our numeric Python libraries. So lets see how Scipy would approach the problem:

from scipy.optimize import fsolve

Two interesting (and different) things about most numeric solvers:

• They assume things are in the form $f(x) = 0$
• You have to give them an initial "guess" as to what the solution might be near. This is why making a quick plot beforehand can be useful!

In this case, I'm already trying to find when the equation equals 0, so nothing needs to be changed there. Based on the plot it looks like the solution is sorta near -1, so I'll guess that:

fsolve(sym.lambdify(x,expr),-1)

array([-1.34776929])

Note that the first input into fsolve is the name of a function! If you have defined a function above using def, you'd just provide the name of the function there. Here, I'd already defined my expression as a Sympy expression above, so I just used lambdify to convert it to a start python function.

As so we arrived at the same place! There are a plethora of other numeric solving functions out there, but I've found fsolve to be a nice combination of ease of use and accuracy.