Regula falsi method in Julia

Regula falsi or false position method is a root-finding algorithm very similar to the secant method. In the regula falsi method, the range \([x_0,x_1]\) where the root is found is redefined in each iteration, depending on the sign of the function evaluation in the new \(x\), this will be set as the minimum or maximum of the new range. Unlike the secant method, regular position method always converges and usually faster than the bisection method.

In each iteration, after calculate the new \(x\) as

\[ x = x_1 - f(x_1) {x_1-x_0 \over f(x_1)-f(x_0)} \ , \]

this will be set as the new \(x_0\) value in the case of \(f(x)f(x_0)\) is positive and as the new \(x_1\) value otherwise.

In Julia:

function falsi(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
        x = x1 - y1 * (x1-x0)/(y1-y0)
        # 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(y0*y)==1
            x0 = x
            y0 = y
        else
            x1 = x
            y1 = y
        end
    end
    error("Max iteration exceeded")
end

The main differences in the structure of the regula falsi method implementation with the secant method implementation are that we reassign the variables according to the sign of \(f(x_0)f(x)\) and that the first function evaluations are outside the loop (since we need to check the sign of the function evaluated in \(x\) there's no point in evaluating the function again, just save the result in a variable and use it again at the beginning of the loop).

Also, along with the accuracy in the \(x\) axis (xtol, the minimum difference between two evaluated \(x\)) we added accuracy in the \(y\) axis (ytol, the minimum to be assumed as zero), which by default is set to the machine precission times two (the same default used in root finding algorithms implementations in Scipy).

As an example, for the equation \(f(x)=x^3-23\), the solution is \(x=\sqrt[3]{23}\):

julia> (f(x)=x^3-23; x0=1; x1=5; falsi(f,x0,x1))
2.843859313381865
julia> cbrt(23)
2.8438669798515654

Like in the other root finding algorithms implementations in this blog, extra arguments and complex numbers can be used.