Ridders' method in Julia

Ridders' method is a root-finding method based on the regula falsi method that uses an exponential function to fit a given function bracketed between \(x_0\) and \(x_1\).

In the algorithm, in each iteration first the function is evaluated at a third point \(x_2 = (x_0+x_1)/2\), then applies the unique exponential function

\[ f(x_0) -2f(x_2)e^Q +f(x_1)e^{2Q} = 0 \]

to transform the function into a straight line.

Solving the quadratic equation for \(e^Q\) gives

\[ e^Q = {f(x_2) + \mathrm{sign}(f(x_1)) \sqrt{f(x_2)^2-f(x_0)f(x_1)} \over f(x_1)} \ , \]

and applying the false position method to \(f(x_0), f(x_2)e^Q\) and \(f(x_1)e^{2Q}\) gives the new guess

\[ x = x_2 + (x_2-x_0) {\mathrm{sign}(f(x_0)-x(x_1)) f(x_2) \over \sqrt{f(x_2)^2-f(x_0)f(x_1)}} \ . \]

This equation converges quadratically, therefore the number of significant digits approximately doubles with each iteration (two function evaluations).

After each iteration, \(x_0\) and \(x_1\) will be redefined, being the new guess one of them and the other one defined according to the sign of the evaluations of \(f\) at \(x_0\), \(x_1\) and \(x_2\).

In Julia:

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

The code follows the same structure as the false position method implementation, supports complex numbers and defined accuracy in \(x\) and \(y\) dimensions.

As an example, for the equation \(f(x)=x^2/12+x-4 \quad (1 \leq x \leq 5)\), the solution is \(x = \sqrt{84}-6\):

julia> (f(x)=x^2/12+x-4; x0=1; x1=5; ridders(f,x0,x1))
3.1651513899117836
julia> sqrt(84)-6
3.1651513899116797