Newton-Raphson method in Julia

Newton-Raphson method (or Newton's method) is a method to find the root of a real function. Knowing the function and its derivative, it will calculate successive approximations to the root from an initial guess, calculating the x-intercept of the tangent line of this guess and using x-intercept value as the new guess for the next iteration. The process will continue till the computed value is accurate enough, according to a tolerance parameter.

On each iteration, an approximation to the root, \(x_{n+1}\), will be calculated. The x-intercept of the tangent line of \(f(0)\) is

\[ 0 = f'(x_n)(x_{n+1}-x_n)+f(x_n). \]

Solving for \(x_{n+1}\) gives the new approximation

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

The convergence may fail due to a bad initial guess or if the derivative function is not continuous at the root or close to the root.

In a stationary point of the function (or close to it) the derivative is zero (or close), then the next approximation would be far from the previous one and the method would be very inefficient and not as quick as expected. However, in this case we will stop the iterations due to a division by zero.

The Newton-Raphson method to find a root of a function for one variable might be implemented in Julia as follows:

function newton(f::Function, x0::Number, fprime::Function, args::Tuple=();
                tol::AbstractFloat=1e-8, maxiter::Integer=50, eps::AbstractFloat=1e-10)
    for _ in 1:maxiter
        yprime = fprime(x0, args...)
        if abs(yprime) < eps
            warn("First derivative is zero")
            return x0
        end
        y = f(x0, args...)
        x1 = x0 - y/yprime
        if abs(x1-x0) < tol
            return x1
        end
        x0 = x1
    end
    error("Max iteration exceeded")
end

As an example, \(f(x)=x^2-2\) has the real roots \(x=\sqrt 2\) and with initial guesses \(x_0=1\) and \(x_0=-1\):

julia> (f(x)=x^2-2; fprime(x)=2x; x0=1; newton(f,x0,fprime))
1.4142135623730951
julia> (f(x)=x^2-2; fprime(x)=2x; x0=-1; newton(f,x0,fprime))
-1.4142135623730951

In this implementation, extra arguments args can be passed to the function f and fprime. An optional keyword argument can be used as the smallest number \(\epsilon\) to use as zero derivative and stop the iterations. Also, x0 can be a complex number, for example, the function \(f(x)=x^2+2\) has the roots \(x=i\sqrt 2\):

julia> (f(x)=x^2+2; fprime(x)=2x; x0=im; newton(f,x0,fprime))
0.0 + 1.4142135623730951im
julia> (f(x)=x^2+2; fprime(x)=2x; x0=-im; newton(f,x0,fprime))
0.0 - 1.4142135623730951im