root finding in python

Question

Edit: A big issue here is that scipy.optimize.brentq requires that the limits of the search interval have opposite sign. If you slice up your search interval into arbitrary sections and run brentq on each section as I do below and as Dan does in the comments, you wind up throwing a lot of useless ValueErrors. Is there a slick way to deal with this in Python?

Original post: I'm repeatedly searching functions for their largest zero in python. Right now I'm using scipy.optimize.brentq to look for a root and then using a brutish search method if my initial bounds don't work:

#function to find the largest root of f
def bigRoot(func, pars):
    try:
        root = brentq(func,0.001,4,pars)
    except ValueError:
        s = 0.1
        while True:
            try:
                root = brentq(func,4-s,4,pars)
                break
            except ValueError:
                s += 0.1
                continue
    return root

There are two big problems with this.

First I assume that if there are multiple roots in an interval, brentq will return the largest. I did some simple testing and I've never seen it return anything except the largest root but I don't know if that's true in all cases.

The second problem is that in the script I'm using this function will always return zero in certain cases even though the function I pass to bigRoot diverges at 0. If I change the step size of the search from 0.1 to 0.01 then it will return a constant nonzero value in those cases. I realize the details of that depend on the function I'm passing to bigRoot but I think the problem might be with the way I'm doing the search.

The question is, what's a smarter way to look for the largest root of a function in python?

---

Thanks Dan; a little more info is below as requested.

The functions I'm searching are well behaved in the regions I'm interested in. An example is plotted below (code at the end of the post).

The only singular point is at 0 (the peak that goes off the top of the plot is finite) and there are either two or three roots. The largest root usually isn't greater than 1 but it never does anything like run away to infinity. The intervals between roots get smaller at the low end of the domain but they'll never become extremely small (I'd say they'll always be larger than 10^-3).

from numpy import exp as e
#this isn't the function I plotted
def V(r):
    return  27.2*(
                23.2*e(-43.8*r) + 
                8.74E-9*e(-32.9*r)/r**6 - 
                5.98E-6*e(-0.116*r)/r**4 + 
                0.0529*( 23*e(-62.5*r) - 6.44*e(-32*r) )/r -
                29.3*e(-59.5*r)
            )

#this is the definition of the function in the plot
def f(r,b,E):
    return 1 - b**2/r**2 - V(r)/E

#the plot is of f(r,0.1,0.06)

Answer 1

Good question, but it's a math problem rather than a Python problem.

In the absence of an analytic formula for the roots of a function, there's no way to guarantee that you've found the largest root of that function, even on a given finite interval. For example, I can construct a function which oscillates between ±1 faster and faster as it approaches 1.

f(x) = sin(1/(1-x))

This would bork any numerical method that tries to find the largest root on the interval [0,1), since for any root there are always larger ones in the interval.

So you'll have to give some background about the characteristics of the functions in question in order to get any more insight into this general problem.

UPDATE: Looks like the functions are well-behaved. The brentq docs suggest there is no guarantee of finding the largest/smallest root in the interval. Try partitioning the intervals and recursively searching for smaller and larger other roots.

from scipy.optimize import brentq

# This function should recursively find ALL the roots in the interval
# and return them ordered from smallest to largest.

from scipy.optimize import brentq
def find_all_roots(f, a, b, pars=(), min_window=0.01):
    try:
        one_root = brentq(f, a, b, pars)
        print "Root at %g in [%g,%g] interval" % (one_root, a, b)
    except ValueError:
        print "No root in [%g,%g] interval" % (a, b)
        return [] # No root in the interval

    if one_root-min_window>a:
        lesser_roots = find_all_roots(f, a, one_root-min_window, pars)
    else:
        lesser_roots = []

    if one_root+min_window<b:
        greater_roots = find_all_roots(f, one_root+min_window, b, pars)
    else:
        greater_roots = []

    return lesser_roots + [one_root] + greater_roots

I tried this on your function and it finds the largest root, at ~0.14.

There's something tricky about brentq, though:

print find_all_roots(sin, 0, 10, ())

Root at 0 in [0,10] interval
Root at 3.14159 in [0.01,10] interval
No root in [0.01,3.13159] interval
No root in [3.15159,10] interval
[0.0, 3.141592653589793]

The sin function should have roots at 0, π, 2π, 3π. But this approach is only finding the first two. I realized that the problem is right there in the docs: f(a) and f(b) must have opposite signs. It appears that all of the scipy.optimize root-finding functions have the same requirement, so partitioning the intervals arbitrarily won't work.