Full linear multistep methods as root-finders

TL;DR

This paper analyzes the use of full linear multistep methods as root-finders, demonstrating their convergence properties, limitations, and practical implementation with guaranteed convergence.

Contribution

It provides a general theoretical analysis of LMM-based root-finders, establishes a fundamental convergence barrier, and offers a robust implementation with guaranteed convergence.

Findings

Full LMM-based methods perform well in numerical tests.

A fundamental barrier limits the convergence rate of LMM root-finders.

A robust Brent-based implementation guarantees convergence.

Abstract

Root-finders based on full linear multistep methods (LMMs) use previous function values, derivatives and root estimates to iteratively find a root of a nonlinear function. As ODE solvers, full LMMs are typically not zero-stable. However, used as root-finders, the interpolation points are convergent so that such stability issues are circumvented. A general analysis is provided based on inverse polynomial interpolation, which is used to prove a fundamental barrier on the convergence rate of any LMM-based method. We show, using numerical examples, that full LMM-based methods perform excellently. Finally, we also provide a robust implementation based on Brent's method that is guaranteed to converge.

Tables

Table 1. Table 1: The first few convergence rates for s 𝑠 s points using only function values.

$s$	$p$
$2$	$1.62$
$3$	$1.84$
$4$	$1.92$
$5$	$1.97$

Table 2. Table 2: The first few convergence rates for Adams-Bashforth root-finder method using s 𝑠 s points.

$s$	$p$
$1$	$2$
$2$	$2.41$
$3$	$2.55$
$4$	$2.59$
$5$	$2.61$

Table 3. Table 3: The first few convergence rates for s 𝑠 s points (vertical) using function values and the first d 𝑑 d derivatives (horizontal).

$s\backslash d$	$1$	$2$	$3$	$4$
$1$	$2$	$3$	$4$	$5$
$2$	$2.73$	$3.79$	$4.82$	$5.85$
$3$	$2.91$	$3.95$	$4.97$	$5.98$
$4$	$2.97$	$3.99$	$4.99$	$5.996$

Table 4. Table 4: Test cases with iterations taken for Newton’s method (subscript N 𝑁 N ) and the LMM-based method with s = 2 𝑠 2 s=2 (subscript 2 2 2 ) and s = 3 𝑠 3 s=3 method (subscript 3 3 3 ).

function	root	$x_0$	$\#$itsN	$\#$its2	$\#$its3	$p_2$	$p_3$
$x+e^x$	$-0.57$	$1.50$	11	8	8	$2.73$	$2.93$
$\sqrt{x}-\cos \left(x\right)$	$0.64$	$0.50$	9	7	8	$2.74$	$2.91$
$e^x-x^2+3x-2$	$0.26$	$0.00$	10	8	7	$2.72$	$2.94$
$x^4-3x^2-3$	$1.95$	$1.30$	17	14	14	$2.73$	$2.92$
$x^3-x-1$	$1.32$	$1.00$	12	9	9	$2.73$	$2.64$
$e^{-x}-x^3$	$0.77$	$2.00$	13	10	10	$2.73$	$2.92$
$5\left(\sin \left(x\right)+\cos \left(x\right)\right)-x$	$2.06$	$1.50$	11	9	9	$2.73$	$2.92$
$x-\cos \left(x\right)$	$0.74$	$1.00$	9	7	7	$2.72$	$2.93$
$\log \left(x-1\right)+\cos \left(x-1\right)$	$1.40$	$1.60$	12	9	9	$2.73$	$2.92$
$\sqrt{1+x}-x$	$1.62$	$1.00$	9	7	7	$2.73$	$2.92$
$\sqrt{e^x-x}-2x$	$0.54$	$1.00$	11	8	7	$2.73$	$2.92$
Total number of iterations			$124$	$96$	$95$

Table 5. Table 5: Convergence history for tanh ⁡ ( x ) 𝑥 \tanh(x) of the three methods: Newton (subscript N 𝑁 N ) and the LMM-based methods with s = 2 𝑠 2 s=2 and s = 3 𝑠 3 s=3 . Note that the root of tanh ⁡ ( x ) 𝑥 \tanh(x) is at x = 0 𝑥 0 x=0 .

$x_N$	$x_{\left(s=2\right)}$	$x_{\left(s=3\right)}$
$1.239$	$1.239$	$1.239$
$-1.719$	$-1.719$	$-1.719$
$6.059$	$0.8045$	$0.8045$
$-4.583\cdot {10}^4$	$0.7925$	$-0.6806$
$\mathrm{Inf}$	$-0.7386$	$1.377$
	-$6.783\cdot {10}^{-3}$	$-0.7730$
	$9.323\cdot {10}^{-6}$	$3.466\cdot {10}^{-2}$
		$-3.032\cdot {10}^{-4}$
		$1.831\cdot {10}^{-11}$

Table 6. Table 6: Convergence history for h ​ ( x ) ℎ 𝑥 h(x) from ( 31 ) of the three methods: Newton (subscript N 𝑁 N ) and the LMM-based methods with s = 2 𝑠 2 s=2 and s = 3 𝑠 3 s=3 . Note that the root is at x = 0 𝑥 0 x=0 .

$x_N$	$x_{\left(s=2\right)}$	$x_{\left(s=3\right)}$
$0.1147$	$0.1147$	$0.1147$
$0.2589$	$-0.2589$	$-0.2589$
$1.0402$	$0.1016$	$0.1016$
$1.6084$	$9.993\cdot {10}^{-2}$	$-5.648\cdot {10}^{-2}$
$1.9407$	$-0.2581$	$0.1959$
$2.2102$	$9.840\cdot {10}^{-2}$	$-0.1611$
$2.4445$	$9.810\cdot {10}^{-2}$	$5.021\cdot {10}^{-2}$
$2.6549$	$-0.2344$	$-7.190\cdot {10}^{-2}$
$2.8478$	$6.602\cdot {10}^{-2}$	$4.947\cdot {10}^{-2}$
$3.0270$	$6.021\cdot {10}^{-2}$	$-3.777\cdot {10}^{-3}$
$3.1953$	$-4.939\cdot {10}^{-2}$	$3.027\cdot {10}^{-4}$
$3.3543$	$-4.019\cdot {10}^{-4}$	$-6.875\cdot {10}^{-6}$
$3.5056$	$1.288\cdot {10}^{-4}$	$1.216\cdot {10}^{-9}$
$3.6502$	$2.028\cdot {10}^{-10}$	$-4.652\cdot {10}^{-15}$
$3.7889$	$-5.308\cdot {10}^{-15}$
$3.9225$

Table 7. Table 7: Test cases with iterations taken for Brent’s method and the LMM-based method. Subscript B 𝐵 B represents Brent while subscript L ​ M ​ M 𝐿 𝑀 𝑀 LMM represents the LMM-based method.

function	root	$[a,b]$	$\#$itsB	$\#$itsLMM
$x+e^x$	$-0.57$	$[-1,1]$	$6$	$4$
$\sqrt{x}-\cos \left(x\right)$	$0.64$	$[0,2]$	$8$	$4$
$e^x-x^2+3x-2$	$0.26$	$[-1,1]$	$5$	$3$
$x^4-3x^2-3$	$1.95$	$[1,3]$	$10$	$8$
$x^3-x-1$	$1.32$	$[0,2]$	$29$	$6$
$e^{-x}-x^3$	$0.77$	$[0,2]$	$9$	$4$
$5\left(\sin \left(x\right)+\cos \left(x\right)\right)-x$	$2.06$	$[0,4]$	$45$	$6$
$x-\cos \left(x\right)$	$0.74$	$[0,1]$	$7$	$3$
$\log \left(x-1\right)+\cos \left(x-1\right)$	$1.39$	$[1.2,1.6]$	$31$	$4$
$\sqrt{1+x}-x$	$1.62$	$[0,2]$	$5$	$3$
$\sqrt{e^x-x}-2x$	$0.54$	$[-1,2]$	$9$	$4$
Total number of iterations			$164$	$49$
Total number of function evaluations			$164$	$98$

Equations

$$f(x)=0.$$

$$\frac{\operatorname{\mathrm{d}\!}x}{\operatorname{\mathrm{d}\!}y}=[f^{-1}]^{\prime}(y)=\frac{1}{f^{\prime}\big{(}x\big{)}}=F(x),$$

$$x(0)=x_{0}+\int_{y_{0}}^{0}F\big{(}x(y)\big{)}\operatorname{\mathrm{d}\!}y,$$

$$x_{n+s}+\sum _{k=0}^{s-1}{a}_k^{(n)}{x}_{n+k}=h_{n+s}\sum _{k=0}^{s-1}{b}_k^{(n)}F(x_{n+k}),$$

$$h_{n+k}=f(x_{n+k})-f(x_{n+k-1})\text{for }k=1,\dots ,s-1.$$

$$x_{n+s}+\sum _{k=0}^{s-1}{a}_k^{(n)}{x}_{n+k}=h_{n+s}\sum _{k=0}^{s-1}{b}_k^{(n)}F(x_{n+k})+h_{n+s}^2\sum _{k=0}^{s-1}{c}_k^{(n)}{F}^{\prime }(x_{n+k})+\dots$$

$$p^s=\sum _{k=0}^{s-1}{p}^k(d+{\sigma }_k),$$

$$f^{-1}(y)-H(y)=\frac{[f^{-1}{]}^{(N+1)}(\upsilon )}{(N+1)!}\prod _{k=0}^{s-1}(y-y_{n+k}{)}^{d+{\sigma }_k},$$

$$|x_{n+s}-\alpha |=\frac{|[f^{-1}{]}^{(N+1)}(\upsilon )|}{(N+1)!}\prod _{k=0}^{s-1}|y_{n+k}{|}^{d+{\sigma }_k}.$$

$$|{\epsilon }_{n+s}|\approx A_0|{\epsilon }_{n+s-1}{|}^{d+{\sigma }_{s-1}}\cdots |{\epsilon }_n{|}^{d+{\sigma }_0},$$

$$|{\epsilon }_{l+1}|=C|{\epsilon }_l{|}^p,$$

$$|{\epsilon }_n{|}^{p^s}\approx A_1|{\epsilon }_n{|}^{{\sum }_{k=0}^{s-1}{p}^k(d+{\sigma }_k)},$$

$$p^s\le (d+1)\sum _{k=0}^{s-1}{p}^k=(d+1)\frac{p^s-1}{p-1}.$$

$$p^{s+1}-(d+2)p^s+d+1\le 0.$$

$$p-(d+2)\le -\frac{d+1}{p^s}<0,$$

$$d+2-p=\frac{d+1}{p^s}.$$

$$\frac{d+1}{p^{s+1}}<\frac{d+1}{p^s},$$

$$p_{s+1}>p_s$$

$$p^{s+1}-2p^s+1=0,$$

$$p^s=\sum _{k=0}^{s-1}{p}^k=\frac{p^s-1}{p-1}.$$

$$p^{s+1}-3p^s+p^{s-1}+1=0,$$

$$p^s=p^{s-1}+\sum _{k=0}^{s-1}{p}^k=p^{s-1}+\frac{p^s-1}{p-1}.$$

$$p^2-3p+1=-\frac{1}{p^{s-1}}<0,$$

$$3-p=\frac{1}{p}+\frac{1}{p^s}.$$

$$\frac{1}{p}+\frac{1}{p^{s+1}}<\frac{1}{p}+\frac{1}{p^s},$$

$$q_k^{(n)}=\frac{f(x_{n+k})}{f(x_{n+s-1})},k=0,\dots ,s-2,$$

$$x(y_{n+k})=x(y_{n+s})+\sum _{m=1}^{\infty }\frac{1}{m!}(-h_{n+s}{q}_k{)}^m{x}^{(m)}(y_{n+s}),$$

$$x^{\prime }(y_{n+k})=x^{\prime }(y_{n+s})+\sum _{m=1}^{\infty }\frac{1}{m!}(-h_{n+s}{q}_k{)}^m{x}^{(m+1)}(y_{n+s}),$$

$$x(y_{n+s})\left[1+\sum _{k=0}^{s-1}{a}_k\right]-h_{n+s}{x}^{\prime }(y_{n+s})\left[\sum _{k=0}^{s-1}{a}_k{q}_k+b_k\right]$$

$$+\sum _{m=2}^{\infty }\frac{1}{(m-1)!}(-h_{n+s}{)}^m{x}^{(m)}(y_{n+s})\sum _{k=0}^{s-1}\left[\frac{1}{m}{q}_k^m{a}_k+q_k^{m-1}{b}_k\right]=0.$$

$$\sum _{k=0}^{s-1}{a}_k=-1,$$

Full text

1 Introduction

Suppose we are given a sufficiently smooth nonlinear function $f:\mathbb{R}\to\mathbb{R}$ and we are asked to solve the equation

[TABLE]

This archetypical problem is ubiquitous in all fields of mathematics, science and engineering. For example, ray tracing techniques in optics and computer graphics need to accurately calculate intersection points between straight lines, rays, and objects of varying shapes and sizes [1, 2]. Implicit ODE solvers are often formulated like (1), after which a root-finder of some kind is applied [3].

Depending on the properties of the function $f$, there are several methods that present themselves. Sometimes the derivative is not available for various reasons, in which case the secant method will prove useful. If higher-order convergence is desired, inverse quadratic interpolation may be used [4]. If the derivative of $f$ exists and is available, Newton’s method is a solid choice, especially if $f$ is also convex.

Recently, a new interpretation of root-finding methods in terms of ODEs has been introduce by Grau-Sánchez et al. [5, 6, 7]. Their idea is to consider the inverse function derivative rule as an ODE, so that any explicit ODE solver may be converted to a root-finding method. Indeed, Grau-Sánchez et al. have successfully introduced root-finders based on Adams-type multistep and Runge-Kutta integrators. It goes without saying that only explicit ODE solvers can be usefully converted to root-finding methods. However, predictor-corrector pairs are possible, as those methods are indeed explicit.

We argue that the ODE approach can be interpreted as inverse interpolation with (higher) derivatives. Indeed, any linear integration method is based on polynomial interpolation. Thus, the ODE approach can be seen as a generalisation of inverse interpolation methods such as the secant method or inverse quadratic interpolation. The analysis can thus be combined into a single approach based on inverse polynomial interpolation.

Our main theoretical result is a theorem on the convergence rate of root-finders based on explicit linear multistep methods. We furthermore prove a barrier on the convergence rate of LMM-based root-finders. It turns out that adding a few history points quickly boosts the convergence rate close to the theoretical bound. However, adding many history points ultimately proves an exercise in futility due to diminishing returns in the convergence rate. Two LMM-based methods are constructed explicitly, one using two history points with a convergence rate of $1+\sqrt{3}\approx 2.73$ and another with three history points that converges with rate $2.91$.

In terms of the efficiency measure, defined as $p^{\frac{1}{w}}$ with $p$ the order and $w$ the number of evaluations [4], our method holds up even compared to several optimal memoryless root-finders. The Kung-Traub conjecture famously states the optimal order of memoryless root-finders using $w$ evaluations is $2^{w-1}$, so that the efficiency is $2^{\frac{w-1}{w}}$. Our method requires two evaluations per iteration, provided function and derivative values are stored, leading to an efficiency of $1.74$ for the three-point method. As a consequence, anything up to an eighth-order method has a lower efficiency measure. For instance, optimal fourth-order methods with $\sqrt[3]{4}\approx 1.59$ [8, 9, 10], or eighth-order methods with $\sqrt[4]{8}\approx 1.68$ [11, 12, 13]. Only 16th-order methods, e.g. Geum and Kim’s [14], with $\sqrt[5]{16}\approx 1.74$ would in theory be more efficient.

Using several numerical examples, we show that the LMM-based methods indeed achieve this higher convergence rate. Furthermore, pathological examples where Newton’s method fails to converge are used to show increased stability. We also construct a robust LMM-based method combined with bisection to produce a method that can been seen as an extension of Brent’s [15]. Similar to Brent’s method, whenever an enclosing starting bracket is provided, an interval $[a,b]$ with $f(a)f(b)<0$, our method is guaranteed to converge.

This article is organised in the following way. First, we find the convergence rate of a wide class of root-finders in Section 2 and prove a barrier on the convergence rates. Next, in Section 3 we derive new root-finders based on full linear multistep methods and show that such methods are stable when the initial guess is sufficiently close to the root. After this, some results are presented in Section 4 that verify our earlier theoretical treatment. Finally, we present our robust implementation in Section 5, after which we give our conclusions in Section 6.

2 Barriers on LMM root-finders

Root-finding methods based on the ODE approach of Grau-Sánchez et al. can be derived by assuming that the function $f$ is sufficiently smooth and invertible in the vicinity of the root. Under these assumptions, the chain rule gives

[TABLE]

which we may interpret as an autonomous ODE for the inverse. Integrating (2) from an initial guess $y_{0}=f(x_{0})$ to $y=0$ yields

[TABLE]

where $x(0)$ is the location of the root. Immediately, we see that applying the forward Euler method to (2) gives Newton’s method. From (3), we see that the step size of the integrator should be taken as $0-y_{0}=-f(x_{0})$. However, Newton’s method may also be interpreted as an inverse linear Taylor method, i.e., a method where the inverse function is approximated by a first-order Taylor polynomial. Indeed, any linear numerical integration method applied to (2) can be interpreted as an inverse polynomial interpolation.

As such, explicit linear multistep methods applied to (2) will also produce a polynomial approximation to the inverse function. Such a method has the form

[TABLE]

where indeed $b_{s}^{(n)}=0$, otherwise we end up with an implicit root-finder, which would not be very useful. The coefficients of the method, $\{a_{k}^{(n)}\}_{k=0}^{s-1}$ and $\{b_{k}^{(n)}\}_{k=0}^{s-1}$, will depend on the previous step sizes and will therefore be different each step. The step sizes are given by $h_{n+k}=y_{n+k}-y_{n+k-1}$, the differences in $y$-values. Since we wish to find the root, we set $y_{n+s}=0$, leading to $h_{n+s}=y_{n+s}-y_{n+s-1}=-y_{n+s-1}$. Furthermore, the $y$-values are of course given by the function values of the root estimates, i.e.,

[TABLE]

Like an ODE solver, we may use an implicit LMM in tandem with an explicit LMM to form a predictor-corrector pair, the whole forming an explicit method. Unlike an ODE solver, we may construct derivative-free root-finders based on the LMM approach by setting all $b_{k}^{(n)}=0$ for $k=0,\ldots s-1$ and for all $n>0$, e.g., the secant method. For an ODE solver this would obviously not make sense. Similar to ODE solvers, we may introduce higher derivatives by using

[TABLE]

The following theorem provides the maximal convergence rate for any method of the form (6). Furthermore, it provides a fundamental barrier on the convergence rates of LMM-based root-finders. Under certain conditions, it also rewards our intuition in the sense that methods using more information converge faster to the root.

Let us introduce some notation first. We denote $d$ the number of derivatives used in the method (6). Higher derivatives of the inverse are found by iteratively applying the inverse function derivative rule. Methods defined by (4) are the special case of (6) with $d=1$. We also introduce coefficients $\sigma_{k}$ that indicate whether the coefficients $a_{k}^{(n)}$ are arbitrarily fixed from the outset or left free to maximise the order of convergence, i.e., $\sigma_{k}=1$ if $a_{k}^{(n)}$ is free and $\sigma_{k}=0$ otherwise.

Theorem 1**.**

For simple roots, the convergence rate $p$ for any method of the form (6), where the coefficients are chosen so as to give the highest order of convergence, is given by the largest real root of

[TABLE]

*for all $s\geq 1$ and $d\geq 1$, or $s\geq 2$ and $d=0$. The convergence rate is bounded by $p<d+2$.

Additionally, if $\sigma_{k}=1$ for all $k=0,\ldots,s-1$, the convergence rates form a monotonically increasing sequence in $s$ and $p\to d+2$ as $s\to\infty$.*

Proof.

1. Any method of the form (6) implicitly uses inverse polynomial (Hermite) interpolation applied to the inverse function $f^{-1}$, let us call the resulting interpolation $H$. Let $y_{n+k}$, $k=0,\ldots,s-1$ be the interpolation points. At each point $y_{n+k}$, there are $d+\sigma_{k}$ values are interpolated, the inverse function value $x_{k}$ if $\sigma_{k}=1$ and $d$ derivative values. Thus, the polynomial interpolation error formula gives

[TABLE]

where $\upsilon$ is in the interval spanned by the interpolation points and $N=sd+\sum_{k=0}^{s-1}\sigma_{k}$. The approximation to the root is then computed as $x_{n+s}=H(0)$. Let us denote the exact value of the root as $\alpha$, then

[TABLE]

Define $\varepsilon_{n+k}=x_{n+k}-\alpha$ and recognise that $f(x_{n+k})=f(\alpha+\varepsilon_{n+k})=f^{\prime}(\alpha)\varepsilon_{n+k}+\mathcal{O}(\varepsilon_{n+k}^{2})$, where $f^{\prime}(\alpha)\neq 0$. Thus, we find

[TABLE]

where $A_{0}>0$ is a constant depending on $[f^{-1}]^{(N+1)}(\upsilon)$, $s$ and $f^{\prime}(\alpha)$. The error behaviour is of the form

[TABLE]

asymptotically as $l\to\infty$. Here, $C>0$ is a constant. Applying $(\ast)$ $s$ times on the left and $s-1$ times on the right-hand side leads to

[TABLE]

where all the constants have been absorbed into $A_{1}$. Thus, (7) is established.

2. For methods that only use a single point, we have $s=1$ and $\sigma_{0}=1$ so that (7) simplifies to $p=d+1$. Hence, also in the case $s=1$, we have $p<d+2$.

3. Finally, by its definition we can bound $\sigma_{k}\leq 1$, so that we obtain

[TABLE]

Simplifying, we obtain

[TABLE]

Note that $p=1$ is always a solution if we impose equality. However, the maximal convergence rate is given by the largest real root, so that we look for solutions $p>1$. Dividing by $p^{s}$ yields,

[TABLE]

which holds for all $s\geq 1$. Hence, we obtain $p<d+2$.

4. Suppose now that $\sigma_{k}=1$ for all $k=0,\ldots,s-1$, so that the convergence rate satisfies $(\star)$ with equality. Rewriting, we obtain

[TABLE]

Thus, the convergence rate is given by the intersection point of a straight line and an inverse power law, both as a function of $p$. First, we observe that there is always an intersection at $p=1$. Furthermore, the slope of the straight line is $-1$, while the slope of the inverse power law at $p=1$ is given by $-s(d+1)$. Therefore, the slope of the right-hand side is smaller than the slope of the left-hand side for all $s\geq 1$ and $d\geq 1$, or $d=0$ and $s\geq 2$. Thus, to the right of $p=1$, the inverse power law is below the straight line. Combined with the fact that $\frac{d+1}{p^{s}}$ is a convex function for $p>0$ and it goes asymptotically to zero, we see that the second intersection point exists and is unique. Fix $s$ and call the second intersection point, i.e., the convergence rate, $p_{s}$. Finally, we note that for $p>1$, we have

[TABLE]

so that $p_{s+1}$ will be moved towards the right compared to $p_{s}$. Thus, we find

[TABLE]

for all $d\geq 1$ and $s\geq 1$, or $d=0$ and $s\geq 2$. The result now follows from the Monotone Convergence Theorem [16]. ∎

From Theorem 1, we find several special cases, such as the derivative-free interpolation root-finders, i.e., $d=0$. Note that derivative-free root-finders with $s=1$ simply do not exist, as then the inverse is approximated with a constant function.

Corollary 1**.**

Inverse polynomial interpolation root-finders, i.e., $d=0$ resulting in all $b_{k}^{(n)}=0$ in (4), can attain at most a convergence rate that is quadratic. Their convergence rates are given by the largest real root of

[TABLE]

for all $s\geq 2$. The convergence rates are bounded by $p<2$ and form a monotonically increasing sequence in $s$, with $p\to 2$ as $s\to\infty$.

Proof.

The coefficients $\{a_{k}^{(n)}\}_{k=0}^{s-1}$ are chosen to maximise the order of convergence, so that $\sigma_{k}=1$ for all $k=0,\ldots,s-1$, while $d=0$, leading to

[TABLE]

Simplifying yields (8). Furthermore, the condition that $\sigma_{k}=1$ for all $k=0,\ldots,s-1$ is satisfied so that the convergence rates form a monotonically increasing sequence that converges to $d+2=2$. ∎

Inverse polynomial root-finders such as the secant method ($s=2$) or inverse quadratic interpolation ($s=3$) are derivative-free, so that their highest convergence rate is $2$ according to Theorem 1. The first few convergence rates for derivative-free inverse polynomial interpolation methods are presented in Table 1. The well known convergence rates for the secant method and the inverse quadratic interpolation method are indeed reproduced. As becomes clear from the table, the rates quickly approach $2$ but never quite get there. The increase in convergence rate becomes smaller and smaller as we increase the number of interpolation points.

Next, we cover the Adams-Bashforth methods also discussed in [5]. As ODE solvers, Adams-Bashforth methods are explicit integration methods that have order of accuracy $s$ [17]. However, as Theorem 1 suggests, as root-finders they will have a convergence rate that is smaller than cubic, since $d=1$. In fact, the convergence rate of Adams-Bashforth root-finders is bounded by $\frac{3+\sqrt{5}}{2}=2.62$ as was proven by Grau-Sánchez et al. [5]. The following corollary is a generalisation of their result.

Corollary 2**.**

The Adams-Bashforth root-finder methods with $s\geq 2$ exhibit convergence rates given by the largest real root of

[TABLE]

for all $s\geq 1$. The convergence rates are bounded by $p<\frac{3+\sqrt{5}}{2}$ and form a monotonically increasing sequence in $s$, with $p\to\frac{3+\sqrt{5}}{2}$ as $s\to\infty$.

Proof.

1. Adams-Bashforth methods have $a_{k}^{(n)}=0$ for $k=0,\ldots,s-2$, resulting in $\sigma_{k}=0$ for $k=0,\ldots,s-2$ and $\sigma_{s-1}=1$. We may write $\sigma_{k}$ for simplicity as a Kronecker delta, i.e., $\sigma_{k}=\delta_{k,s-1}$. Furthermore, the methods use a single derivative of $f^{-1}$ so that $d=1$. The $s=1$ method is equal to Newton’s method, which has a quadratic convergence rate. For $s\geq 2$, we find from Theorem 1 that

[TABLE]

Simplifying yields (9). Again, we assume that $p>1$ and we divide by $p^{s-1}$, so that

[TABLE]

which holds for all $s\geq 2$. This implies that $p<\frac{3+\sqrt{5}}{2}$ for all $s\geq 2$.

2. The proof that the convergence rates make up a monotonically increasing sequence is similar to the one given for Theorem 1. First, we note that (9) always has a root at $p=1$. Next, we rewrite it to read

[TABLE]

The left-hand side is a straight line with slope $-1$ while the right-hand side is a convex function that has slope $-3$ at $p=1$. Thus, to the right of $p=1$, the convex function is below the straight line. Furthermore, the right-hand side goes asymptotically to zero for large $p$. This implies that there is a unique intersection point with $p>1$, which is the convergence rate. For fixed $s$, call the second intersection point $p_{s}$. Finally, we note that

[TABLE]

for $p>1$, from which we see that the intersection point is moved to the right for larger $s$, i.e., $p_{s+1}>p_{s}$. The result again follows from the Monotone Convergence Theorem. ∎

The first few convergence rates for the Adams-Bashforth root-finder methods are given in Table 2 and agree with the rates found by Grau-Sánchez et al. As becomes clear from the table, the convergence rates quickly draw near the bound of $2.62$. Yet again we are met with steeply diminishing returns as we increase the number of history points $s$.

The Adams-Bashforth root-finder methods cannot attain a convergence rate higher than $2.62$, which is still some way off the cubic bound given by Theorem 1. However, the theorem also tells us that using full linear multistep methods will result in convergence rates closer to cubic, at least in the limit of large $s$. For ODE solvers, trying to obtain a higher convergence rate by increasing the number of points often leads to instabilities. Indeed, the stability regions of many LMMs tend to shrink as the order is increased [3]. In general, polynomial interpolation on equispaced points is unstable, e.g., Runge’s phenomenon [18].

We must note that the stability issues in ODE solvers arise from the fact that polynomial interpolation is applied on an equispaced grid. Root-finders are designed to home in on a root, and when convergent, the step sizes will decrease rapidly. Runge’s phenomenon can be countered by placing more nodes closer to the boundary, for example Gauß or Lebesgue nodes. As we will see, inverse polynomial interpolation is stable on the set of points generated by the root-finder itself, provided the starting guess is sufficiently close.

Let us inspect the convergence rates of different LMM-based root-finders using Theorem 1, see Table 3. These convergence rates are computed under the assumption that all derivatives and point values are used, i.e., $\sigma_{k}=1$ for $k=0,\ldots,s-1$ in Theorem 1. The convergence rate of a $d$-derivative method can be boosted by at most $1$, and the table shows that this mark is attained very quickly indeed. Adding a few history points raises the convergence rate significantly, but finding schemes with $s>3$ is likely to be a waste of time.

Thus, provided that the root-finders are stable, a higher convergence rate can be achieved by adding history points, as well as derivative information.

3 Full LMM-based root-finders

Let us investigate full LMM-based root-finders that use a single derivative, thus methods of the form (4). The current step size is then given by $h_{n+s}=-f(x_{n+s-1})$. Let us define $q_{k}^{(n)}$ as

[TABLE]

so that $h_{n+s}q_{k}^{(n)}=-f(x_{n+k})$ is the total step between $y_{n+k}$ and $y_{n+s}=0$. The Taylor expansions of $x(y_{n+k})$ and $x^{\prime}(y_{n+k})$ about $y_{n+s}$ are then given by

[TABLE]

where we have dropped the superscript $(n)$ for brevity. Substituting these into (4), we obtain

[TABLE]

The consistency conditions then are given by

[TABLE]

This gives us two equations for $2s$ coefficients, so that we can eliminate another $2s-2$ leading order terms, resulting in the conditions

[TABLE]

where $m=2,\ldots,2s-1$.

3.1 The $s=2$ method

The $s=2$ LMM-based method is given by

[TABLE]

where we have again suppressed the superscript $(n)$ on the coefficients. Here, $h_{n+2}=-f(x_{n+1})$ so that we may write $q=q_{0}$, i.e.

[TABLE]

Applying (14) and (15), we find a set of linear equations, i.e.,

[TABLE]

These equations may be solved, provided $q\neq 1$, to yield

[TABLE]

The condition $q\neq 1$ is equivalent to $f(x_{n+1})\neq f(x_{n})$. This condition is not very restrictive, as stronger conditions are needed to ensure convergence.

The above method may also be derived from the inverse polynomial interpolation perspective, using the ansatz

[TABLE]

where $h_{i}$, $i=0,1,2,3$ are undetermined coefficients. The coefficients are fixed by demanding that $H$ interpolates $f^{-1}$ and its derivative at $y=f(x_{n+1})$ and $y=f(x_{n})$, i.e.,

[TABLE]

Solving for $h_{i}$, $i=0,1,2,3$ and setting $y=0$, we find the same update $x_{n+2}$ as (16).

The stability of the $s=2$ LMM method depends on the coefficients of the LMM in much the same way as an ODE solver. Indeed, we can set the sequence $\tilde{x}_{n}=x_{n}+z_{n}$ where $x_{n}$ is the sequence generated by exact arithmetic we wish to find while $z_{n}$ is a parasitic mode. It can be shown that the parasitic mode satisfies

[TABLE]

so that it may grow unbounded if the roots are greater than $1$ in modulus. Using the ansatz $z_{n}=B\lambda^{n}$, we find the characteristic polynomial of the $s=2$ method, i.e.,

[TABLE]

where the roots cans simply be read off. Stability of the root-finder is ensured if the stability polynomial of the method has a single root with $\lambda=1$, while the other roots satisfy $|\lambda|<1$. This property is called zero-stability for linear multistep ODE solvers. Thus, to suppress parasitic modes we need

[TABLE]

This reduces to $q$ being either $q<0$, or $q>3$, so that $|q|>3$ is a sufficient condition. Thus, if the sequence $\{|f(x_{n})|\}_{n=1}^{\infty}$ is decreasing fast enough, any parasitic mode is suppressed. We may estimate $q$ as a ratio of errors, since $f(x_{n})=f^{\prime}(\alpha)\varepsilon_{n}+\mathcal{O}(\varepsilon_{n}^{2})$, so that

[TABLE]

Using $\varepsilon_{n+1}=C\varepsilon_{n}^{p}$ with $p=1+\sqrt{3}$, we find that

[TABLE]

with $C_{1}=\left(\frac{1}{3C}\right)^{\frac{1}{\sqrt{3}}}$. We conclude that the method will be stable if the initial errors are smaller than the constant $C_{1}$, which depends on the details of the function $f$ in the vicinity of the root. This condition translates to having the starting values sufficiently close to the root. This is a rather typical stability condition for root-finders.

3.2 $s=3$ method

We may again apply (14) - (15) to find a method with $s=3$, this time we have $6$ coefficients, given by

[TABLE]

where $q_{0}=\frac{f(x_{n})}{f(x_{n+2})}$ and $q_{1}=\frac{f(x_{n+1})}{f(x_{n+2})}$. Here, we have the conditions $q_{0}\neq 1$ and $q_{1}\neq 1$, reducing to the condition that all $y$-values must be unique. Again, this condition is not very restrictive for reasons detailed above.

Methods with a greater number of history points are possible, however, the gain in convergence rate from $s=3$ to $s=4$ is rather slim, as indicated by Table 3. If such methods are desirable, they can be derived by selecting coefficients that satisfy (14) - (15).

4 Results

Like the secant method, the $s=2$ full LMM root-finding method needs two starting points for the iteration. However, as the analytical derivative is available, we choose to simply start Newton’s method with one point, say $x_{0}$, giving $x_{1}$. The $s=3$ method needs three starting values, therefore the next value $x_{2}$ is obtained from the $s=2$ LMM method. The LMM-based methods can be efficiently implemented by storing the function value and derivative value of the previous step, thus resulting in a need for only one function and one derivative evaluation per iteration.

The efficiency measure is defined as $p^{\frac{1}{w}}$ with $p$ the order of convergence and $w$ the number of evaluations per iteration [4]. Assuming the function itself and the derivative cost the same to evaluate, the $s=3$ LMM-based method has an efficiency measure of $\sqrt{2.91}\approx 1.71$. Compared to Newton’s method, with an efficiency measure of $\sqrt{2}\approx 1.41$, this is certainly an improvement.

4.1 Numerical examples

Here, we provide a number of test cases and show how many iterations LMM-based root-finders take versus the number needed by Newton’s method, see Table 4. We have used a selection of different test cases with polynomials, exponentials, trigonometric functions, square roots and combinations thereof. For each of the test problems shown, the methods converged within a few iterations. Some problems were deliberately started near a maximum or minimum to see the behaviour when the derivatives are small.

The test computations were performed using the variable-precision arithmetic of MATLAB’s Symbolic Math Toolbox. The number of digits was set to 300 while the convergence criterion used was

[TABLE]

with $\eta=250$. The numerical convergence rates were computed with the error behaviour

[TABLE]

asymptotically as $l\to\infty$. The limiting value of the estimates for $p$ is displayed in Table 4.

Overall, Newton’s method consistently displayed a quadratic convergence rate, which is the reason we did not display it in the table. The LMM-based methods, on the other hand, generally have a higher convergence rate that may vary somewhat from problem to problem. This is due to the fact that the step sizes may vary slightly while the convergence rates of the LMM-based methods only holds asymptotically, even with so many digits.

4.2 Pathological functions

A classical example of a pathological function for Newton’s method is the hyperbolic tangent $\tanh(x)$. Students often believe that Newton’s method converges for any monotone function until they are asked to find the root of $\tanh(x)$. Hence, we have used this function as a test case using standard double precision floating point arithmetic and a convergence criterion reading

[TABLE]

with $\epsilon$ the machine precision. Newton’s method fails to converge for starting values with approximately $|x_{0}|\geq 1.089$, see Table 5. Our $s=2$ LMM-based method extends this range somewhat more and converges for any starting value with roughly $|x_{0}|\leq 1.239$. The behaviour of the $s=3$ LMM-based method is similar, though it does take two more iterations to converge.

Starting at $x_{0}=1.239$, Newton’s method diverges quickly, returned as $-\mathrm{inf}$ after only 4 iterations. The LMM-based method, on the other hand, bounces around positive and negative values for about 5 iterations until it is close enough to the root. After that, the asymptotic convergence rate sets in and the root is quickly found, reaching the root within machine precision at 7 iterations.

Donovan et al.[19] developed another test function for which Newton’s method fails, it gives a false convergence result to be precise. The test function is given by

[TABLE]

which is, in fact, infinitely steep near the root $x=0$, yet smooth, see Figure 1. Again, we used double precision arithmetic and (30) as a stopping criterion. Newton’s method diverges for any starting value except the exact root itself. However, Newton’s method eventually gives a false convergence result as the increment $|x_{l+1}-x_{l}|$ falls below the tolerance. The $s=2$ and $s=3$ LMM-based methods converge when starting with $|x_{0}|\leq 0.1147$ for this problem, see Table 6.

Starting at the maximal $x_{0}=0.1147$ for instance, the LMM-based methods bounce several times between positive and negative $x$-values without making much headway. After that, the root is close enough and the asymptotic convergence rate sets in, reaching the root to within machine precision in a few steps.

We believe the reason that the LMM-based method has increased stability is due to the fact that it uses two points to evaluate the function and its derivative. In both cases, the iterations jump between positive and negative values, enclosing the root. In this fashion, the LMM-based method acts much like the regula falsi method. Once the iterates are close enough to the root, the asymptotic convergence rate sets in and the iterates converge in but a few steps.

5 A robust implementation

As with most open root-finding algorithms, the conditions under which the method is guaranteed to converge are rather restrictive. Therefore, we have designed a bracketing version of the LMM-based method that is guaranteed to converge. The algorithm is based on Brent’s method, using similar conditions to catch either slow convergence or runaway divergence. This version of the LMM-based method does, however, require an enclosing bracket $[a,b]$ on which the function changes sign, i.e., $f(a)f(b)<0$. Alternatively, such a method can start out as an open method, switching to the bracketing method when a sign change is detected.

The algorithm consists of a cascade of methods increasing in accuracy but decreasing in robustness, similar to Brent’s method. At the lowest level stands the most robust method, bisection, guarding against steps outside the current search bracket. Bisection is guaranteed to converge, but does so rather slowly. On the highest level we use the full $s=3$ LMM-based method discussed in the previous section. Thus, in the best possible case, the method will converge with a rate of $2.91$. The method is, by virtue of the bisection method, guaranteed to converge to a root.

Like Brent’s method and Dekker’s method, the LMM-based method keeps track of three points $a$, $b$ and $c$. Here, $b$ is the best estimate of the root so far, $c$ is the previous value for $b$ while $a$ is the contrapoint so that $\mathrm{int}[a,b]$ encloses the root. Ideally, all three values are used to compute the next value for $b$. However, extra conditions are added to ensure the inverse actually makes sense on the interval $\mathrm{int}[a,c]$.

Consider the case where the sign of $f^{\prime}(c)$ is not equal to the sign of $\frac{f(b)-f(a)}{b-a}$, but the sign of $f^{\prime}(b)$ is. It follows that there is an extremum between $b$ and $c$, and the inverse function does not exist in the interval $\mathrm{int}[a,c]$, leading to an error if we were to compute the inverse interpolation.

Thus, the following condition should be applied to each derivative value: the sign of the derivative needs to be the same as the sign of the secant slope on $\mathrm{int}[a,b]$, i.e.,

[TABLE]

where $z=a,b,c$. Only when the derivative at a point satisfies (32) can the derivative sensibly contribute to the inverse. Otherwise the derivative information should be discarded, leading to lower-order interpolation. If all derivatives are discarded, the resulting interpolation is inverse quadratic or the secant method.

Ultimately, the method provides an interval on which the function $f$ changes sign with a relative size of some given tolerance $\delta$, i.e.,

[TABLE]

We shall use $\delta=2\epsilon$ in all our examples, with $\epsilon$ the machine precision. As an input, the algorithm has $f$, $f^{\prime}$, $a$ and $b$ such that $f(a)f(b)<0$. The algorithm can be described in the following way:

If all three function values are different, use $s=3$, otherwise use $s=2$. 2. 2.

Check the sign of the derivatives at points $a$, $b$ and $c$, include the derivatives that have the proper sign. 3. 3.

If the interpolation step is worse than a bisection step, or outside the interval $[a,b]$, use bisection. 4. 4.

If the step is smaller than the tolerance, use the tolerance as step size. 5. 5.

If the convergence criterion is met, exit, otherwise go to 1.

The first step determines the number of history points that can be used. The second step determines which derivative values should be taken into account. In effect, only the second step is essentially different from Brent’s method, with all the following steps exactly the same [15]. The extra conditions on the derivatives gives rise to a selection of 12 possible root-finders, including inverse quadratic interpolation and the secant method. Our method can therefore be seen as building another 10 options on top of Brent’s method. Naturally, sufficiently close to the root, the derivative conditions will be satisfied at all three points and the method will use the full LMM method with $s=3$.

5.1 Comparison with Brent’s method

Here we give a few examples of the robust LMM-based root-finding algorithm discussed above compared to Brent’s method. As a performance measure, we use the number of iterations. Standard double precision arithmetic is employed, as that provides sufficient material for comparison. For both methods, the stopping criterion is given by (33), the relative size of the interval must be sufficiently small.

Table 7 shows that for most functions, both Brent’s method and the LMM-based method take a comparable number of iterations. However, in some cases, the difference is considerable. In the worst case considered here, Brent’s method takes $7.5$ times as many iterations to converge. In terms of efficiency index, Brent’s method should be superior with an efficiency index of $1.84$ against $1.71$ of the LMM-based method. Taken over the whole set of test functions, however, Brent’s method takes more than three times as many iterations in total, leading to a significant increase in function evaluations. We conclude therefore that practically, the LMM-based root-finder is a better choice.

6 Conclusions

We have discussed root-finders based on full linear multistep methods. Such LMM-based methods may be interpreted as inverse polynomial (Hermite) interpolation methods, resulting in a simple and general convergence analysis. Furthermore, we have proven a fundamental barrier for LMM-based root-finders: their convergence rate cannot exceed $d+2$, where $d$ is the number of derivatives used in the method.

The results indicate that compared to the Adams-Bashforth root-finder methods of Grau-Sánchez et al. [5], any full LMM-based method with $s\geq 2$ has a higher convergence rate. As ODE solvers, full LMMs are typically not zero-stable and special choices of the coefficients have to be made. Employed as root-finders on the other hand, it turns out that LMMs are stable, due to the rapid decrease of the step size. This allows the usage of full LMMs that are otherwise not zero-stable.

Contrary to the Adams-type methods, the full LMM-based root-finders can achieve the convergence rate of $3$ in the limit that all history points are used. The $s=2$ and $s=3$ methods, $s$ being the number of history points, were explicitly constructed and provide a convergence rate of $2.73$ and $2.91$, respectively. Numerical experiments confirm these predicted convergence rates. Furthermore, application to pathological functions where Newton’s method diverges show that the LMM-based methods also have enhanced stability properties.

Finally, we have implemented a robust LMM-based method that is guaranteed to converge when provided with an enclosing bracket of the root. The resulting robust LMM root-finder algorithm is a cascade of twelve root-finders increasing in convergence rate but decreasing in reliability. At the base sits bisection, so that the method is indeed guaranteed to converge to the root. At the top resides the $s=3$ LMM-based root-finder, providing a maximal convergence rate of $2.91$.

In terms of efficiency index, Brent’s method is theoretically the preferred choice with $1.84$ compared to $1.71$ for the LMM-based method. However, numerical examples show that the increased convergence rate leads to a significant decrease in the total number of function evaluations over a range of test functions. Therefore, in practical situations, provided the derivative is available, the LMM-based method performs better.

Acknowledgements

This work was generously supported by Philips Lighting and the Intelligent Lighting Institute.