The Interpolate, Truncate, Project (ITP) root-finding algorithm was developed in Oliveira and Takahashi (2020). It’s performance compares favourably with existing methods on both well-behaved functions and ill-behaved functions while retaining the worst-case reliability of the bisection method. For details see the authors’ Kudos summary and the Wikipedia article ITP method.
The itp
function implements the ITP method to find a root \(x^*\) of the function \(f: \mathbb{R} \rightarrow \mathbb{R}\) in the interval \([a, b]\), where \(f(a)f(b) < 0\). If \(f\) is continuous over \([a, b]\) then \(f(x^*) = 0\). If \(f\) is discontinuous over \([a, b]\) then \(x^*\) may be an estimate of a point of discontinuity at which the sign of \(f\) changes, that is, \(f(x^* - \delta)f(x^* + \delta) \leq 0\), where \(0 \leq \delta \leq \epsilon\) for some tolerance value \(\epsilon\).
We use some of the examples presented in Table 1 of Oliveira and Takahashi (2020) to illustrate the use of this function and run the examples using the uniroot
function in the stats
package as a means of comparison, using a convergence tolerance of \(10^{-10}\) in both cases. The itp
function uses the following default values of the tuning parameters: \(\kappa_1 = 0.2 / (b - a)\), \(\kappa_2 = 2\) and \(n_0 = 1\), but these may be changed by the user. See Sections 2 and 3 of Oliveira and Takahashi (2020) for information. The following function prints output from uniroot
in the same style as the output from itp
.
# Method to print part of uniroot output
<- function(x, digits = max(3L, getOption("digits") - 3L)) {
print.list names(x)[1:3] <- c("root", "f(root)", "iterations")
print.default(format(x[1:3], digits = digits), print.gap = 2L, quote = FALSE)
}
library(itp)
These functions are infinitely differentiable and contain only one simple root over \([−1, 1]\).
# Lambert
<- function(x) x * exp(x) - 1
lambert itp(lambert, c(-1, 1))
#> root f(root) iterations
#> 0.5671 2.048e-12 8
uniroot(lambert, c(-1, 1), tol = 1e-10)
#> root f(root) iterations
#> 0.5671 -8.623e-12 8
# Trigonometric 1
<- function(x) tan(x - 1 /10)
trig1 itp(trig1, c(-1, 1))
#> root f(root) iterations
#> 0.1 -8.238e-14 8
uniroot(trig1, c(-1, 1), tol = 1e-10)
#> root f(root) iterations
#> 0.1 6.212e-14 8
This function has a non-simple root at \(10^{-6}\), with a multiplicity of 3.
# Polynomial 3
<- function(x) (x * 1e6 - 1) ^ 3
poly3 itp(poly3, c(-1, 1))
#> root f(root) iterations
#> 1e-06 -1.25e-13 36
# Using n0 = 0 leads to (slightly) fewer iterations, in this example
<- function(x) (x * 1e6 - 1) ^ 3
poly3 itp(poly3, c(-1, 1), n0 = 0)
#> root f(root) iterations
#> 1e-06 -6.739e-14 35
uniroot(poly3, c(-1, 1), tol = 1e-10)
#> root f(root) iterations
#> 1e-06 1.771e-17 65
This function has discontinuities, including one at the location of the root.
# Staircase
<- function(x) ceiling(10 * x - 1) + 1 / 2
staircase itp(staircase, c(-1, 1))
#> root f(root) iterations
#> 7.404e-11 0.5 31
uniroot(staircase, c(-1, 1), tol = 1e-10)
#> root f(root) iterations
#> -3.739e-11 -0.5 31
This function has multiple roots, but we only seek one of them.
# Warsaw
<- function(x) ifelse(x > -1, sin(1 / (x + 1)), -1)
warsaw itp(warsaw, c(-1, 1))
#> root f(root) iterations
#> -0.6817 -5.472e-11 11
uniroot(warsaw, c(-1, 1), tol = 1e-10)
#> root f(root) iterations
#> -0.6817 -3.216e-16 9
In terms of a naive comparison based on the number of iterations itp
and uniroot
perform similarly, except in the repeated-root “Polynomial 3” example, where itp
requires fewer iterations.