Root Finding for Functions with many maxima and minima

Question

I wondering if anyone can provide advice on the best combination of algorithms to find the roots (or any one root) of a function which is "dense" in that it has many local maxima and minima for example consider $\sin\left (x^{2} \right ) + \sin^{2}\left (x \right )$. I have been trying to distinguish between the best methods for this type of function but I can't seem to find any references to functions that behave like that. For the purposes of the answer, I do however have a range in which at least one root is known to lie. Many thanks.

Answer 1

Let's take a look at a graph of the function $f(x)=\sin(x^2)+\sin^2(x)$, which has the salient characteristics of lots of roots and "many local maxima and minima". Here it is over the interval $[0,10]$:

I produced this image with Mathematica, after finding the roots with NSolve, which is quite good at this sort of thing:

roots = x /. NSolve[Sin[x^2] + Sin[x]^2 == 0 && 0 < x < 10, x];
Short[roots]
(* Out: {2.02145, 2.41336, 3.07079, <<25>>, 9.54572, 9.70421, 9.87839} *)

Looks like there are 31 roots lying strictly between 1 and 10. In fact, the technique used by Mathematica is outlined very roughly in this Wolfram blog post, which makes the point that determination of the number of the roots in a region is a key step in the process. As outlined in the post, this can be done quite efficiently using numerical integration in the complex plane. Specifically, the argument principle states that for a simple, closed, contractible curve $\gamma$ that doesn't pass through a pole or zero of a function $f$, we have

$$\frac{1}{2\pi i} \int_{\gamma} \frac{f'(z)}{f(z)} \, dz = Z-P,$$

where $Z$ and $P$ are the number of zeros and poles of $f$ inside $\gamma$ counted according to multiplicity.

For example, if we want to determine the number of roots of our example function $f(z)=\sin(z^2)+\sin(z)^2$ in the interval $(0,10)$, we might let $$\gamma(t) = 5+4.99\cos(t) + i\sin(t)/5.$$ (Note that $f$ has no poles.) If we plot this in the complex plane together with a contour plot of $|f(z)|$, we get something like so:

The integral along $\gamma$ can be evaluated numerically and result rounded, since we're expecting an integer. In Mathematica, you'd do

f[x_] = Sin[x^2] + Sin[x]^2;
g[t_] = 5 + 4.99 Cos[t] + I*Sin[t]/5;
NIntegrate[f'[g[t]]/f[g[t]]*g'[t], {t, 0, 2 Pi}]/(2 Pi*I)
(* Out: 31. - 7.50964*10^-17 I *)

Thus, the interval has 31 roots, as expected. In practice, this process is performed recursively so that we can isolate the roots and then use an iterative process like Newton's method on the isolated portions.

It should be mentioned that the integral for a function like this can be highly oscillatory. Mathematica does very well on this because it implements the Levin technique for oscillatory integrals. The same integral in scipy doesn't go so well.

from numpy import arange, pi, sin, cos
from scipy import integrate
from scipy.misc import derivative
def f(z): return sin(z**2)+sin(z)**2
def g(t): return 5+4.99*cos(t)+1j*sin(t)/5.0
def integrand(t):
    return (derivative(f,g(t),0.0001)/f(g(t)))*derivative(g,t,0.0001)/(2.0*pi*1j)
integrate.quad(integrand, 0, 2*pi)

#Out: (30.99998942349135, 5.6807397269042424e-05)

There's also a long warning issued by scipy. The convergence can be improved by breaking the integral in two parts, but that needs to be done by the user.

---

A completely different approach is taken by the Matlab routine FindRealRoots. It's quite simple - we first approximate the function with a Chebyshev polynomial and then we find the roots of the polynomial. We can apply it to our sample problem like so:

roots = FindRealRoots(@(x) sin(x.^2)+sin(x).^2,0,10,100);
length(roots)
%Out ans=33

The procedure finds all the roots (including the double root at zero). They are not as accurate as the Mathematica version, though, and you've got to experiment with the third parameter (namely, the order of the Chebyshev polynomial) to see if you've got all the roots or not. It's very simple and quite quick, however.