Secant method in Julia

The secant method is a root-finding algorithm that recursively calculates the roots of secant lines of points defined in \(f\).

Starting with initial values \(x_0\) and \(x_1\), the equation to find the root of the line between \((x_0, f(x_0))\) and \((x_1, f(x_1))\) is

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

Solving the equation for \(x\) gives

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

The new \(x\) obtained will be the next \(x_1\) in the recurrence process, which can be expressed as

\[ x_n = x_{n-1} - f(x_{n-1}){x_{n-1}-x_{n-2} \over f(x_{n-1})-f(x_{n-2})} \ . \]

The process of solving \(x\) will continue until it reaches a level of tolerance is reached, that is a defined minimum value of the difference between \(x_1\) and \(x\).

In Julia:

function secant(f::Function, x0::Number, x1::Number, args::Tuple=();
                tol::AbstractFloat=1e-5, maxiter::Integer=50)
    for _ in 1:maxiter
        y1 = f(x1, args...)
        y0 = f(x0, args...)
        x = x1 - y1 * (x1-x0)/(y1-y0)
        if abs(x-x1) < tol
            return x
        end
        x0 = x1
        x1 = x
    end
    error("Max iteration exceeded")
end

As an example, the function \(f(x)=x^2-10\) has the roots \(\pm \sqrt{10}\) and with initial guesses \(x_0=1\) and \(x_1=\pm 2\):

julia> (f(x)=x^2-10; x0=1; x1=2; secant(f,x0,x1))
3.162277660040216
julia> (f(x)=x^2-10; x0=1; x1=-2; secant(f,x0,x1))
-3.1622776609633
julia> sqrt(10)
3.1622776601683795

In the same way of the implementation of the Newton-Raphson method, extra arguments args can be passed to the function f and result may be complex numbers if complex numbers are passed as initial guesses. The function \(f(x,c)=x^2+10+c\) with \(c=5\) has the roots \(\pm i \sqrt{15}\), with initial guesses \(x_0=1\) and \(x_1=2i\)

julia> (f(x,c)=x^2+10+c; x0=1; x1=2im; secant(f,x0,x1,(5,)))
-8.268421911988619e-11 + 3.8729833464880765im
julia> sqrt(-15+0im)
0.0 + 3.872983346207417im