Inverse quadratic interpolation in Julia

Inverse quadratic interpolation is rarely used directly as a root-finding algorithm, but is is part of the popular Brent's method.

In the inverse quadratic interpolation, a function \(f(x)\) is evaluated in three points (\(x_0\), \(x_1\) and \(x_2\)) within an interval where the root of the function is known to be found. Assuming that the function \(f(x)\) has an inverse quadratic function and applying the Lagrange polynomial for the three points gives the inverse quadratic function

\[ g(y) = {(y-f(x_1))(y-f(x_2)) \over (f(x_0)-f(x_1))(f(x_0)-f(x_2))} \ x_0 + {(y-f(x_0))(y-f(x_2)) \over (f(x_1)-f(x_0))(f(x_1)-f(x_2))} \ x_1 + {(y-f(x_0))(y-f(x_1)) \over (f(x_2)-f(x_0))(f(x_2)-f(x_1))} \ x_2 \ . \]

Then, we replace \(y=0\) to obtain the new guess

\[ x = {f(x_1)f(x_2) \over (f(x_0)-f(x_1))(f(x_0)-f(x_2))} \ x_0 + {f(x_0)f(x_2) \over (f(x_1)-f(x_0))(f(x_1)-f(x_2))} \ x_1 + {f(x_0)f(x_1) \over (f(x_2)-f(x_0))(f(x_2)-f(x_1))} \ x_2 \ . \]

After each iteration, \(x_1\) will be set as \(x_0\), \(x_2\) as \(x_1\) and \(x\) as \(x_2\) until the algorithm converges.

In Julia:

function invquadinterp(f::Function, x0::Number, x1::Number, x2::Number,
                       args::Tuple=(); xtol::AbstractFloat=1e-5,
                       ytol=2eps(Float64), maxiter::Integer=50)
    y0 = f(x0,args...)
    y1 = f(x1,args...)
    y2 = f(x2,args...)
    for _ in 1:maxiter
        x = x0*y1*y2/((y0-y1)*(y0-y2)) +
            x1*y0*y2/((y1-y0)*(y1-y2)) +
            x2*y0*y1/((y2-y0)*(y2-y1))
        # x-tolerance.
        if min(abs(x-x0),abs(x-x1),abs(x-x2)) < xtol
            return x
        end
        y = f(x,args...)
        # y-tolerance.
        if abs(y) < ytol
            return x
        end
        x0 = x1
        y0 = y1
        x1 = x2
        y1 = y2
        x2 = x
        y2 = y
    end
    error("Max iteration exceeded")
end

This implementation has support for extra function arguments and the accuracy can be defined in the \(x\) and \(y\) dimensions.

As an example, for the function \(f(x)=x^4-2x^2+1/4 \quad (0 \leq x \leq 1)\), the solution is \(\sqrt{1-\sqrt{3}/2}\):

julia> (f(x)=x^4-2x^2+1/4; invquadinterp(f,0,.5,1))
0.3660254037449329
julia> sqrt(1-sqrt(3)/2)
0.3660254037844387