# Inverse quadratic interpolation in Julia

numerical-analysis root-finding 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