Function solve_itp

pub fn solve_itp(
    f: impl FnMut(f64) -> f64,
    a: f64,
    b: f64,
    epsilon: f64,
    n0: usize,
    k1: f64,
    ya: f64,
    yb: f64,
) -> f64

Solve an arbitrary function for a zero-crossing.

This uses the ITP method, as described in the paper An Enhancement of the Bisection Method Average Performance Preserving Minmax Optimality.

The values of ya and yb are given as arguments rather than computed from f, as the values may already be known, or they may be less expensive to compute as special cases.

It is assumed that ya < 0.0 and yb > 0.0, otherwise unexpected results may occur.

The value of epsilon must be larger than 2^-63 times b - a, otherwise integer overflow may occur. The a and b parameters represent the lower and upper bounds of the bracket searched for a solution.

The ITP method has tuning parameters. This implementation hardwires k2 to 2, both because it avoids an expensive floating point exponentiation, and because this value has been tested to work well with curve fitting problems.

The n0 parameter controls the relative impact of the bisection and secant components. When it is 0, the number of iterations is guaranteed to be no more than the number required by bisection (thus, this method is strictly superior to bisection). However, when the function is smooth, a value of 1 gives the secant method more of a chance to engage, so the average number of iterations is likely lower, though there can be one more iteration than bisection in the worst case.

The k1 parameter is harder to characterize, and interested users are referred to the paper, as well as encouraged to do empirical testing. To match the paper, a value of 0.2 / (b - a) is suggested, and this is confirmed to give good results.

When the function is monotonic, the returned result is guaranteed to be within epsilon of the zero crossing. For more detailed analysis, again see the paper.