[SciPy-user] All the roots of a function in an interval

[Gael Varoquaux]

Hi Anne,

Thanks a lot for your reply, enlightening as usual.

On Sat, Sep 29, 2007 at 03:51:30PM -0400, Anne Archibald wrote:
> The easiest and most reliable method here is to simply sample the
> function at equispaced points the nearest distance apart they can be;
> then every sign change gives you an interval bracketing a root.

I did exactly as you suggested. It is indeed the best solution.
Unfortunately it is slowish, and cannot be well implemented with arrays.
Here is the code I got to, just in case it can be of some use for
somebody later on:

def find_x(m_K, m_Rb, delta_x):
    """ Returns the 2D array of x for which:
        m_K * sin(k_Rb * x) - m_Rb * sin(k_K *(x + delta_x)) == 0
        with delta_x a 1D array.
    """
    delta_x = c_[delta_x]
    f = lambda x, delta_x: m_K*sin(k_Rb*x) - m_Rb*sin(k_K*(x + delta_x))
    interval = linspace(-5*Delta_x, 5*Delta_x, 150)
    scan = f(interval, delta_x)
    sign_changes_mask = (scan[:, 1:]*scan[:, :-1])<0
    min_num_x = (sign_changes_mask.sum(axis=-1)).min()
    x = empty((len(delta_x), min_num_x))
    indexes = arange(len(interval))
    for index, this_delta_x in enumerate(delta_x.flat):
        #this_x = _find_x(m_K, m_Rb, this_delta_x)
        # Keep only the min_num_x inner most solutions:
        sign_changes_index = indexes[sign_changes_mask[index, :]]
        sign_changes_index = sign_changes_index[
                                int((len(sign_changes_index)-min_num_x)/2):
                                int((len(sign_changes_index)+min_num_x)/2)]
        this_x = []
        for column_index in sign_changes_index:
            a = interval[column_index]
            b = interval[column_index + 1]
            this_x.append(optimize.brentq(f, a, b, args=(this_delta_x,)))
        x[index, :] = this_x
    return x

One of the difficulties was that, depending on my parameter "delta_x",
the number of root changed, so I couldn't easily fit this in an array.
The good news is that the further away these roots where from 0, the less
important they were in the calculation (I have a gaussian prefactor), so
I just dump the extra ones, and thus build an array out of a fixed number
of roots, which speeds up a lot the calculations I do on these roots.

This code has served me well over the week end, and produced good
results. I think I just figured out a way to address my problem with a
totally different approach, so no more root finding! For those
interested, it is an inverse problem, for which I have an analytical
model, and it seems that as the root finding is expansive, and the model
get harder and harder to derive analytically as I refine it, the brute
force approach of calculating the forward problem for all the values of
my only unknown parameter, and matching it to the data scales better than
using the analytical model to do the inverse problem. I don't know if I
am clear, but if not just forget it.

Any way, thanks a lot for your help, this mailing list is a gold mine !

Gaël