# Best method to find the zero of a decreasing function numerically

I need to find the zero of a function $$f(\lambda)$$ which is of the form $$\sum \frac{c_i^2}{(1+\lambda d_i)^2} -1$$. I tried using Newton's method, and it works sometimes, but it is highly dependent on the initial choice, and the chances of divergence are very high since the function is "almost" flat.

What would be the best method to find the zero of such a function?

For now, I replaced Newton with the secant method, and it seems to work better.

Off the top of my head, there are a number of things that may be going wrong. You might want to verify the following

• Are the function values or derivatives you are computing numerically stable? Is the rounding error in their evaluation smooth? You can verify this visually by plotting an extremely small interval around a known root such that the function values are $\pm \varepsilon_\mathsf{mach}$. If your function or derivative changes sign near the root, Newton's method won't work.

• What termination criteria are you using? If you're looking for a root $f(\lambda)=0$ with $\lambda \in [a,b]$, do you stop when $|b-a|<\tau$ or when $|f(\lambda)|<\tau$ where $\tau$ is your tolerance? In either case, what is your rationale for your choice of $\tau$ and how does it fit with the error/interval you observe in the previous point?

• Newton's method is a good hammer, but there are many functions that will look a lot like your thumb. The Secant method is a good choice for functions that are almost linear in the search interval, but will perform poorly as soon as you depart from that assumption. A good alternative is the so-called Illinois algorithm. A major advantage of interval methods is that they don't rely on the function or any of its derivatives being numerically smooth as you can run them until $fl(f(\lambda))=0$ or until the search interval $[a,b]$ is numerically zero, i.e. $fl(a)=fl(b)$.

Can you bound your zero? If you can come up with a good heuristic for the approximate location of the zero and then compute bounds for it over an interval that is not too large, the performance of many methods is improved. Here's some Matlab code I use to perform basic root bracketing for monotonic functions.

In terms of interval methods, @Pedro's suggestion to use the Illinois algorithm is a good simple possibility. If you want a big hammer, there is of course Brent's method, which is more complex to implement. It's the basis of Matlab's fzero.

• The Golden section search is an excellent minimum/maximum finder, but not a root finder. Oct 17, 2013 at 19:37
• @Pedro: You're correct -I switched over to general optimization in my head. Oct 17, 2013 at 19:48

Using a line search or trust region method along with Newton's method typically does the trick of finding the zero reliably. Check any advanced numerical methods or optimization book.

Assume $$d_i$$ is an increasing strictly-positive sequence so that there is always a local minima for your function between $$-1/d_i$$ and $$-1/d_{i+1}$$. Since $$f^{\prime\prime}(\lambda)\ge0$$, these local minima are unique. Finding this local minima can always be handled using bracketing methods applied to $$f^\prime(\lambda)$$, which in your case is

$$f^\prime(\lambda)=\sum_i\frac{-2d_ic_i^2}{(1+\lambda d_i)^3}.$$

Let $$\lambda_i\in(-1/d_i,-1/d_{i+1})$$ be these minimas. If $$f(\lambda_i)=0$$, then that is the only root in $$(-1/d_i,-1/d_{i+1})$$. Otherwise if $$f(\lambda_i)<0$$, then only roots exist and are between $$-1/d_i$$ and $$\lambda_i$$, and $$\lambda_i$$ and $$-1/d_{i+1}$$, which can again be found using bracketing methods.

This covers all roots between the $$-1/d_i$$, leaving only $$\lambda\to\pm\infty$$ cases to be explored. Since $$f(\lambda)\to-1$$ as $$\lambda\to\pm\infty$$, we know 2 more roots are guaranteed to exist in $$(-\infty,-1/d_0)\cup(-1/d_n,\infty)$$, which can again be found using bracketing methods. Furthermore, $$f^\prime(\lambda)\ne0$$ in that region, which guarantees no other roots exist in those regions.

Similar cases may be constructed if $$d_i$$ have both positive and negative signs but is a bit messier.

Some example code implementing the above algorithm.