itp {itp}R Documentation

The Interpolate, Truncate, Project (ITP) root-finding algorithm

Description

Performs one-dimensional root-finding using the ITP algorithm of Oliveira and Takahashi (2021). The function itp searches an interval [a, b] for a root (i.e. a zero) of the function f with respect to its first argument. Each iteration results in a bracketing interval for the root that is narrower than the previous interval. If the function is discontinuous then a point of discontinuity at which the function changes sign may be found.

Usage

itp(
  f,
  interval,
  ...,
  a = min(interval),
  b = max(interval),
  f.a = f(a, ...),
  f.b = f(b, ...),
  epsilon = 1e-10,
  k1 = 0.2/(b - a),
  k2 = 2,
  n0 = 1
)

Arguments

f

The function for which the root is sought.

interval

A numeric vector c(a, b) of length 2 containing the end points of the interval to be searched for the root. The function values at the end points must be of opposite signs.

...

additional named or unnamed arguments to be pass to f.

a, b

An alternative way to set the lower and upper end points of the interval to be searched. The function values at these end points must be of opposite signs.

f.a, f.b

The values of f(a) and f(b), respectively.

epsilon

A positive numeric scalar. The desired accuracy of the root. The algorithm continues until the width of the bracketing interval for the root is less than or equal to 2 * epsilon. The value of epsilon should be greater than 2-63(b-a) to avoid integer overflow.

k1, k2, n0

the values of the tuning parameters κ1, κ2, n0. See Details.

Details

Page 8 of Oliveira and Takahashi (2021) describes the ITP algorithm and the roles of the tuning parameters κ1, κ2 and n0. The algorithm is described using pseudocode. The Wikipedia entry for the ITP method provides a summary. If the input function f is continuous over the interval [a, b] then the value of f evaluated at the estimated root is (approximately) equal to 0. If f is discontinuous over the interval [a, b] then the bracketing interval returned after convergence has the property that the signs of the function f at the end points of this interval are different and therefore the estimated root may be a point of discontinuity at which the sign of f changes.

The ITP method requires at most nmax = n1/2 + n0 iterations, where n1/2 is the smallest integer not less than log2((b-a) / 2ε). If n0 = 0 then the ITP method will require no more iterations than the bisection method. Depending on the function f, setting a larger value for n0, e.g. the default setting n0=1 used by the itp function, may result in a smaller number of iterations.

The default values of the other tuning parameters (epsilon = 1e-10, k1 = 0.1, k2 = 2 / (b - a)) are set based on arguments made in Oliveira and Takahashi (2021).

Value

An object (a list) of class "itp" containing the following components:

root

the location of the root, calculated as (a+b)/2, where [a, b] is the bracketing interval after convergence.

f.root

the value of the function evaluated at root.

iter

the number of iterations performed.

a,b

the end points of the bracketing interval [a, b] after convergence.

f.a,f.b

the values of function at a and b after convergence.

estim.prec

an approximate estimated precision for root, equal to the half the width of the final bracket for the root.

If the root occurs at one of the input endpoints a or b then iter = 0 and estim.prec = NA.

References

Oliveira, I. F. D. and Takahashi, R. H. C. (2021). An Enhancement of the Bisection Method Average Performance Preserving Minmax Optimality, ACM Transactions on Mathematical Software, 47(1), 1-24. doi: 10.1145/3423597

See Also

print.itp to print objects of class "itp" returned from itp.

Examples

#### ----- The example used in the Wikipedia entry for the ITP method

wiki <- function(x) x ^ 3 - x - 2
itp(wiki, c(1, 2), epsilon = 0.0005, k1 = 0.1, n0 = 1)
# The default setting (with k1 = 0.2) wins by 1 iteration
itp(wiki, c(1, 2), epsilon = 0.0005, n0 = 1)

#### ----- Some examples from Table 1 of Oliveira and Takahashi (2021)

### Well-behaved functions

# Lambert
lambert <- function(x) x * exp(x) - 1
itp(lambert, c(-1, 1))

# Trigonometric 1
trig1 <- function(x) tan(x - 1 /10)
itp(trig1, c(-1, 1))

# Logarithmic
logarithmic <- function(x) log(abs(x - 10 / 9))
itp(logarithmic, c(-1, 1))

# Linear
linear <- function(x) x
# Solution in one iteration
itp(linear, c(-1, 1))
# Solution at an input endpoint
itp(linear, c(-1, 0))

### Ill-behaved functions

## Non-simple zero

# Polynomial 3
poly3 <- function(x) (x * 1e6 - 1) ^ 3
itp(poly3, c(-1, 1))
# Using n0 = 0 leads to fewer iterations, in this example
poly3 <- function(x) (x * 1e6 - 1) ^ 3
itp(poly3, c(-1, 1), n0 = 0)

## Discontinuous

# Staircase
staircase <- function(x) ceiling(10 * x - 1) + 1 / 2
itp(staircase, c(-1, 1))

## Multiple roots

# Warsaw
warsaw <- function(x) ifelse(x > -1, sin(1 / (x + 1)), -1)
itp(warsaw, c(-1, 1))

[Package itp version 1.0.1 Index]