Instantly share code, notes, and snippets.

# np-8/LICENCE.md

Last active Nov 24, 2022
Find all zeroes (roots) of a function with python (1d)

# Zero / root finder using scipy.optimize.fsolve (Python)

• For functions that have only one tunable variable (other arguments are fixed)
• It can find any roots from interval (start, stop).
• Use relatively small stepsize step to find all the roots. This is used as stepsize for changing the x0 for the fsolve().
• Can only search for zeroes in one dimension (other dimensions must be fixed). If you have other needs, I would recommend using sympy for calculating the analytical solution.
• Can handle singularities (skips them)

### Usage example:

def f(u, V=90, ell=5):
# Example from https://stackoverflow.com/q/14878110/3015186
w = np.sqrt(V ** 2 - u ** 2)

jl = scipy.special.jn(ell, u)
jl1 = scipy.special.yn(ell - 1, u)
kl = scipy.special.kn(ell, w)
kl1 = scipy.special.kn(ell - 1, w)

return jl / (u * jl1) + kl / (w * kl1)

r = RootFinder(1, 20, 0.01)
roots = r.find(f, 90, 5)

Output: Roots: [ 2.84599497 8.82720551 12.38857782 15.74736542 19.02545276]

Note: The example picture seems to cross zero multiple times, also when there is no marked zero. This is a plotting artifact since the example function f() has singularities.

This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters. Learn more about bidirectional Unicode characters
 """ Simple zero/root finder for functions with one variable. Finds all the zeroes within a range of values, using user-defined searchstep and tolerance. Utilizes the scipy.optimize.fsolve: https://docs.scipy.org/doc/scipy/reference/generated/scipy.optimize.fsolve.html LICENCE: MIT """ import numpy as np import scipy from scipy.optimize import fsolve from matplotlib import pyplot as plt class RootFinder: def __init__(self, start, stop, step=0.01, root_dtype="float64", xtol=1e-9): self.start = start self.stop = stop self.step = step self.xtol = xtol self.roots = np.array([], dtype=root_dtype) def add_to_roots(self, x): if (x < self.start) or (x > self.stop): return # outside range if any(abs(self.roots - x) < self.xtol): return # root already found. self.roots = np.append(self.roots, x) def find(self, f, *args): current = self.start for x0 in np.arange(self.start, self.stop + self.step, self.step): if x0 < current: continue x = self.find_root(f, x0, *args) if x is None: # no root found. continue current = x self.add_to_roots(x) return self.roots def find_root(self, f, x0, *args): x, _, ier, _ = fsolve(f, x0=x0, args=args, full_output=True, xtol=self.xtol) if ier == 1: return x[0] return None def f(u, V=90, ell=5): # Example from https://stackoverflow.com/q/14878110/3015186 w = np.sqrt(V ** 2 - u ** 2) jl = scipy.special.jn(ell, u) jl1 = scipy.special.yn(ell - 1, u) kl = scipy.special.kn(ell, w) kl1 = scipy.special.kn(ell - 1, w) return jl / (u * jl1) + kl / (w * kl1) if __name__ == "__main__": r = RootFinder(1, 20, 0.01) args = (90, 5) roots = r.find(f, *args) print("Roots: ", roots) # plot results u = np.linspace(1, 20, num=600) fig, ax = plt.subplots() ax.plot(u, f(u, *args)) ax.scatter(roots, f(np.array(roots), *args), color="r", s=10) ax.grid(color="grey", ls="--", lw=0.5) plt.show()

Permission is hereby granted, free of charge, to any person obtaining a copy of this software and associated documentation files (the "Software"), to deal in the Software without restriction, including without limitation the rights to use, copy, modify, merge, publish, distribute, sublicense, and/or sell copies of the Software, and to permit persons to whom the Software is furnished to do so, subject to the following conditions:

The above copyright notice and this permission notice shall be included in all copies or substantial portions of the Software.

THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE.

### eperez1700 commented Jul 12, 2022

Thank you very much Niko!