11

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.

Charles Brunet
  • 211
  • 2
  • 5

2 Answers2

7

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.

Pedro
  • 9,573
  • 1
  • 36
  • 45
0

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.

Brian Borchers
  • 18,719
  • 1
  • 39
  • 70
  • 1
    These kinds of zeros are clearly impossible to find, but @Charles seemed to be interested in, at worst, black box functions with jump discontinuities, but not so-called removable discontinuities. – Bill Barth Mar 23 '13 at 20:07
  • 1
    Even if you limit yourself to jump discontinuities and even if you limit yourself to continuous functions, if the function isn't Lipschitz continuous on known intervals, then finding all of the zeroes from evaluations at a finite number of points won't ensure that you get all of the roots. – Brian Borchers Mar 24 '13 at 03:46
  • In particular, consider the function $\sin(1/x)$ as an example where finding all zeros on the interval $[0,1]$ will be difficult. – Wolfgang Bangerth Mar 24 '13 at 15:50
  • The OP was willing to specify an $\epsilon$. If the function is pathological, it will find many zeros, but it seems like that's life. It's possible that he may also have to set a maximum number of intervals for searching to avoid such pathologies. – Bill Barth Mar 25 '13 at 12:26