python - How to find all zeros of a function using numpy (and scipy)?

Question

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:

1. I evaluate my function on a given number of points
2. I detect whether there is a change of sign
3. I find the zero between the points that are changing sign

4. 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:

1. 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.
2. 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)

Answer 1

Why are you limited to numpy? Scipy has a package that does exactly what you want:

http://docs.scipy.org/doc/scipy/reference/optimize.nonlin.html

One lesson I've learned: numerical programming is hard, so don't do it :)

---

Anyway, if you're dead set on building the algorithm yourself, the doc page on scipy I linked (takes forever to load, btw) gives you a list of algorithms to start with. One method that I've used before is to discretize the function to the degree that is necessary for your problem. (That is, tune delta x so that it is much smaller than the characteristic size in your problem.) This lets you look for features of the function (like changes in sign). AND, you can compute the derivative of a line segment (probably since kindergarten) pretty easily, so your discretized function has a well-defined first derivative. Because you've tuned the dx to be smaller than the characteristic size, you're guaranteed not to miss any features of the function that are important for your problem.

If you want to know what "characteristic size" means, look for some parameter of your function with units of length or 1/length. That is, for some function f(x), assume x has units of length and f has no units. Then look for the things that multiply x. For example, if you want to discretize cos(pi x), the parameter that multiplies x (if x has units of length) must have units of 1/length. So the characteristic size of cos(pi x) is 1/pi. If you make your discretization much smaller than this, you won't have any issues. To be sure, this trick won't always work, so you may need to do some tinkering.