$$f(x) = \sin(x^2)$$
We can do it simply, with python and scipy, using, for example, the brentq function, that solve such problem with the Brent method, wich is considered the best of rootingfinding routines of scipy.
This method finds the first zero of the given function, in an interval $[a,b]$, where the function has different sign at the bounds.
from scipy.optimize import brentq
import numpy as np
from matplotlib import pyplot as plt
x = np.linspace(0,np.pi,1000)
def y(x):
return np.sin(x**2)
x0_old = x[0]
roots = []
for x0 in x:
if np.sign(y(x0))!=np.sign(y(x0_old)):
roots.append(brentq(y,x0_old,x0))
x0_old = x0
print roots
plt.plot(x,y(x))
--
[0.0, 1.7724538509055159, 2.506628274631493, 3.0699801238394655]
But what we can do, if we doesn't know the exact definition of the function at all? Let's suppose now, we know just the x and y arrays. Well, we can use the interpolation method built in scipy, to define a function object, then all the proceeding is the same as before:
from scipy.optimize import brentq from scipy.interpolate import interp1d import numpy as np from matplotlib import pyplot as plt x = np.linspace(0,np.pi,1000) y = np.sin(x**2) y = interp1d(x,y) x0_old = x[0] roots = [] for x0 in x: if np.sign(y(x0))!=np.sign(y(x0_old)): roots.append(brentq(y,x0_old,x0)) x0_old = x0 print roots plt.plot(x,y(x)) -- [0.0, 1.7724531933700884, 2.506628127081023, 3.0699798592999334]
Obviously, whether we know the function has only one root, it doesn't make any sense to use the for statement, but we can apply directly the Brent method to the nonlinear function to increase the code speed.
Although we can find all the roots, the code might be really slow, and inefficient, if the x array is large, even if the function has only one root, because the complexity is related through the x array length.
No comments:
Post a Comment