Find all the roots of a function in a given interval

Question

I need to find all the roots of a scalar function in a given interval. The function may have discontinuities. The algorithm can have a precision of ε (e.g. it is ok if the algorithm doesn't find two distinct roots that are closer than ε).

Does such algorithm exists? Could you point me papers about that?

---

Actually, I have a function to find a zero in a given interval using Brent's algorithm, and a function to find a minimum in a given interval. Using those two functions, I built my own algorithm, but I was wondering if a better algorithm exists. My algorithm is like that:

I start with an interval [a,b] and a function f. If sign(f(a+ε)) ≠ sign(f(b-ε)), I know there is at least one zero between a and b, and I find z = zero(]a,b[). I test if z really is a zero (it could be a discontinuity), by looking a the value of z-εand z+ε. If it is, I add it to the list of found zeros. If f(a+ε) and f(b-ε) both are positive, I search m = min(]a, b[). If f(m) still is positive, I search m = max(]a,b[) because there could be a discontinuity between a and b. I do the opposite if f(a+ε) and f(b-ε) were negatives.

Now, from the point I found (z or m) I build a stack containing the zeros, discontinuities, and inflection points of my function. After the first iteration, the stack now looks like [a, z, b]. I start again the algorithm from intervals ]a,z[ and ]z,b[. When, between two points a and b, the extrema have the same sign than both interval ends, and there is no discontinuities at both extrema, I remove the interval from the stack. The algorithm ends when there is no more intervals.

Answer 1

If you're using Matlab, you may want to try the Chebfun system (disclaimer: I used to be an active developer of this project). It can find all the roots of a one-dimensional function in a closed or open interval to machine precision.

The main idea behind the Chebfun root-finder is to use a combination of recursive bisection and the Colleague Matrix, an analogue of the Companion Matrix, on the coefficients of an interpolant of the target function.

I have a simplified version of the code here. The function chebroots takes an anonymous function as its first input, the finite interval as a second and third argument, and a degree N as it's fourth and final argument. For reasonable results, you can set N to 100.

Answer 2

In general, this is a hopeless quest- without some information about the continuity and/or differentiability of the function, anything could happen. Consider for example the MATLAB function defined on the interval from 0 to 1:

function y=f(x)

y=1.0;

if (x==0.01)

y=0.0;

end

if (x==0.013)

y=0.0;

end

if (x==0.753124)

y=0.0;

end

Treating this function as a block box, there's no way to see that it has zeroes at these three points and no other points in the interval from 0 to 1 without checking every floating point number between 0 and 1.