I have been unable to find any real code implementation of ITP method - Wikipedia. However, it claims to be the bee’s knees - just read the first paragraph of the Wikipedia article:
“In numerical analysis, the ITP method , short for Interpolate Truncate and Project , is the first root-finding algorithm that achieves the superlinear convergence of the secant method while retaining the optimal worst-case performance of the bisection method. It is also the first method with guaranteed average performance strictly better than the bisection method under any continuous distribution. In practice it performs better than traditional interpolation and hybrid based strategies (Brent’s Method, Ridders, Illinois), since it not only converges super-linearly over well behaved functions but also guarantees fast performance under ill-behaved functions where interpolations fail.”
I implemented the pseudocode in the Wikipedia article, and thought I should share it. If anyone wants to put this in a package like
NonlinearSolve.jl, feel free to - I’d be thrilled.
EDIT: provided defaults, and changed argument order to make having defaults more sensible.
EDIT 2: Put the function argument first, diverging further from the Wikipedia pseudocode, but converging towards the conventional (and sensible) ordering. Also added reference to the source of hyperparameter defaults.
""" ITP(f, a, b, ϵ=eps(), κ₁=0.2/(b-a), κ₂=2, n₀=1) Find the root of the funtion `f` on the interval [a, b], given a tolerance ϵ and hyperparameters κ₁, κ₂, n₀. The hyperparameters are subject to the following constraints: 0 < κ₁ < Inf 1 ≤ κ₂ < 1 + (1 + √5)/2 0 < n₀ < Inf A Julia implementation of the Interpolate, Truncate Project algorithm, as detailed in `https://en.wikipedia.org/wiki/ITP_method`. Hyperparameter defaults are set on the basis of the reasoning provided in `https://docs.rs/kurbo/0.8.1/kurbo/common/fn.solve_itp.html` """ function ITP(f, a, b, ϵ=eps(), κ₁=0.2/(b-a), κ₂=2, n₀=1) 0 < κ₁ < Inf || error("κ₁ must be between 0 and ∞") 1 ≤ κ₂ < 1 + (1 + √5)/2 || error("κ₂ must be between 1 and 1+ (1 + √5)/2 (1 + the golden ratio)") 0 < n₀ < Inf || error("n₀ must be between 0 and ∞") n_1div2 = ceil(Int, log2((b-a)/2ϵ)) nₘₐₓ = n_1div2 + n₀ y_a = f(a) y_b = f(b) sign(y_a) == sign(y_b) && error("sign(f(a)) = sign(f(b)). There is no guaranteed root in the given interval.") j = 0 while b-a > 2ϵ # Calculating parameters: x_1div2 = (a+b)/2 r = ϵ*2^(nₘₐₓ - j) - (b-a)/2 δ = κ₁*(b-a)^κ₂ # Interpolation: x_f = /(y_b*a - y_a*b, y_b-y_a) # Truncation: σ = sign(x_1div2 - x_f) δ ≤ abs(x_1div2 - x_f) ? (x_t=x_f+σ*δ) : (x_t = x_1div2) # Projection: abs(x_t - x_1div2) ≤ r ? (x_ITP = x_t) : (x_ITP = x_1div2 - σ*r) # Updating Interval: y_ITP = f(x_ITP) if y_ITP > 0 b = x_ITP y_b = y_ITP elseif y_ITP < 0 a = x_ITP y_a = y_ITP else a = b = x_ITP end j += 1 end return (a+b)/2 end
One could concider adding the following constraint:
ϵ ≤ big(1e-78) && error("A tolerance this small has lead to indefinite iterations in the past. Aborting.")