Skip to main content

Nicolás' blog

Numerical methods challenge: Day 3

Nicolás Guarín-Zapata

Comments

Also available in:

Español

During October (2017) I will write a program per day for some well-known numerical methods in both Python and Julia. It is intended to be an exercise then don't expect the code to be good enough for real use. Also, I should mention that I have almost no experience with Julia, so it probably won't be idiomatic Julia but more Python-like Julia.

Newton's method

Today's method is Newton's method. This method is used to solve the equation \(f(x) = 0\) for \(x\) real, and \(f\) a differentiable function. It starts with an initial guess \(x_0\) and it succesively refine it by finding the intercept of the tangent line to the function with zero. The new approximation is computed from the old one using

\begin{equation*} x_k = x_{k-1} - \frac{f(x)}{f'(x)} \end{equation*}

Convergence of this method is generally faster than bisection method. Nevertheless, the convergence is not guaranteed. Another drawback is the need of the derivative of the function.

We will use the function \(f(x) = \cos(x) - x\) to test the codes, and the initial guess is 1.0.

Following are the codes.

Python

from __future__ import division, print_function
from numpy import abs, cos, sin

def newton(fun, grad, x, niter=50, ftol=1e-12, verbose=False):
    msg = "Maximum number of iterations reached."
    for cont in range(niter):
        if abs(grad(x)) < ftol:
            x = None
            msg = "Derivative near to zero."
            break
        if verbose:
            print("n: {}, x: {}".format(cont, x))
        x = x - fun(x)/grad(x)
        if abs(fun(x)) < ftol:
            msg = "Root found with desired accuracy."
            break
    return x, msg


def fun(x):
    return cos(x) - x


def grad(x):
    return -sin(x) - 1.0


print(newton(fun, grad, 1.0))

Julia

function newton(fun, grad, x, niter=50, ftol=1e-12, verbose=false)
    msg = "Maximum number of iterations reached."
    for cont = 1:niter
        if abs(grad(x)) < ftol
            x = nothing
            msg = "Derivative near to zero."
            break
        end
        if verbose
            println("n: $(cont), x: $(x)")
        end
        x = x - fun(x)/grad(x)
        if abs(fun(x)) < ftol
            msg = "Root found with desired accuracy."
            break
        end
    end
    return x, msg
end


function fun(x)
    return cos(x) - x
end


function grad(x)
    return -sin(x) - 1.0
end


println(newton(fun, grad, 1.0))

Comparison

Regarding number of lines we have: 28 in Python and 32 in Julia. The comparison in execution time is done with %timeit magic command in IPython and @benchmark in Julia.

For Python:

%timeit newton(fun, grad, 1.0)

with result

10000 loops, best of 3: 27.3 µs per loop

For Julia:

@benchmark newton(fun, grad, 1.0)

with result

BenchmarkTools.Trial:
  memory estimate:  48 bytes
  allocs estimate:  2
  --------------
  minimum time:     327.925 ns (0.00% GC)
  median time:      337.956 ns (0.00% GC)
  mean time:        351.064 ns (0.80% GC)
  maximum time:     8.118 μs (92.60% GC)
  --------------
  samples:          10000
  evals/sample:     226

Comments

Comments powered by Disqus