How to find all zeros of a function using numpy (and scipy)?
numpy, python, scipy
Solution
The main problem I see with this is if you can actually find all roots --- as have already been mentioned in comments, this is not always possible. If you are sure that your function is not completely pathological (`sin(1/x)` was already mentioned), the next one is what's your tolerance to missing a root or several of them. Put differently, it's about to what length you are prepared to go to make sure you did not miss any --- to the best of my knowledge, there is no general method to isolate all the roots for you, so you'll have to do it yourself. What you show is a reasonable first step already. A couple of comments:
- Brent's method is indeed a good choice here.
- First of all, deal with the divergencies. Since in your function you have Bessels in the denominators, you can first solve for their roots -- better look them up in e.g., Abramovitch and Stegun (Mathworld link). This will be a better than using an ad hoc grid you're using.
- What you can do, once you've found two roots or divergencies, `x_1` and `x_2`, run the search again in the interval `[x_1+epsilon, x_2-epsilon]`. Continue until no more roots are found (Brent's method is guaranteed to converge to a root, provided there is one).
- If you cannot enumerate all the divergencies, you might want to be a little more careful in verifying a candidate is indeed a divergency: given `x` don't just check that `f(x)` is large, check that, e.g. `|f(x-epsilon/2)| > |f(x-epsilon)|` for several values of `epsilon` (1e-8, 1e-9, 1e-10, something like that).
- If you want to make sure you don't have roots which simply touch zero, look for the extrema of the function, and for each extremum, `x_e`, check the value of `f(x_e)`.
Problem
Suppose I have a function `f(x)` defined between `a` and `b`. This function can have many zeros, but also many asymptotes. I need to retrieve all the zeros of this function. What is the best way to do it? Actually, my strategy is the following: - I evaluate my function on a given number of points - I detect whether there is a change of sign - I find the zero between the points that are changing sign I verify if the zero found is really a zero, or if this is an asymptote ``` U = numpy.linspace(a, b, 100) # evaluate function at 100 different points c = f(U) s = numpy.sign(c) for i in range(100-1): if s[i] + s[i+1] == 0: # oposite signs u = scipy.optimize.brentq(f, U[i], U[i+1]) z = f(u) if numpy.isnan(z) or abs(z) > 1e-3: continue print('found zero at {}'.format(u)) ``` This algorithm seems to work, except I see two potential problems: - It will not detect a zero that doesn't cross the x axis (for example, in a function like `f(x) = x**2`) However, I don't think it can occur with the function I'm evaluating. - If the discretization points are too far, there could be more that one zero between them, and the algorithm could fail finding them. Do you have a better strategy (still efficient) to find all the zeros of a function? I don't think it's important for the question, but for those who are curious, I'm dealing with characteristic equations of wave propagation in optical fiber. The function looks like (where `V` and `ell` are previously defined, and `ell` is an positive integer): ``` def f(u): w = numpy.sqrt(V**2 - u**2) jl = scipy.special.jn(ell, u) jl1 = scipy.special.jnjn(ell-1, u) kl = scipy.special.jnkn(ell, w) kl1 = scipy.special.jnkn(ell-1, w) return jl / (u*jl1) + kl / (w*kl1) ```