Skip to main content

Nicolás' blog

Numerical methods challenge

Nicolás Guarín-Zapata

Comments

Also available in:

Español

For the next month 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.

Day 1: Bisection method

The first method to be considered is the bisection method. This method is used to solve the equation \(f(x) = 0\) for \(x\) real, and \(f\) continuous. It starts with an interval \([a,b]\), where \(f(a)\) and \(f(b)\) should have opposite signs. The method then proceeds by halving the interval and selecting the one where the solution appears, i.e., the sign of the function changes.

We will use the function \(f(x) = \cos(x) - x^2\) to test the codes, and the initial interval is [0, 1].

The following are the codes:

Python

from __future__ import division, print_function
from numpy import log2, ceil, abs, cos

def bisection(fun, a, b, xtol=1e-6, ftol=1e-12):
    if fun(a) * fun(b) > 0:
        c = None
        msg = "The function should have a sign change in the interval."
    else:
        nmax = int(ceil(log2((b - a)/xtol)))
        for cont in range(nmax):
            c = 0.5*(a + b)
            if abs(fun(c)) < ftol:
                msg = "Root found with desired accuracy."
                break
            elif fun(a) * fun(c) < 0:
                b = c
            elif fun(b) * fun(c) < 0:
                a = c
            msg = "Maximum number of iterations reached."
    return c, msg

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

print(bisection(fun, 0.0, 1.0))

With result

(0.824131965637207, 'Maximum number of iterations reached.')

Julia

function bisection(fun, a, b, xtol=1e-6, ftol=1e-12)
    if fun(a) * fun(b) > 0
        c = nothing
        msg = "The function should have a sign change in the interval."
    else
        nmax = ceil(log2((b - a)/xtol))
        for cont = 1:nmax
            c = 0.5*(a + b)
            if abs(fun(c)) < ftol
                msg = "Root found with desired accuracy."
                break
            elseif fun(a) * fun(c) < 0
                b = c
            elseif fun(b) * fun(c) < 0
                a = c
            end
            msg = "Maximum number of iterations reached."
        end
    end
    return c, msg
end

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

println(bisection(fun, 0.0, 1.0))

With result

(0.824131965637207, "Maximum number of iterations reached.")

In this case, both codes are really close to each other. The Python code has 25 lines, while the Julia one has 27. As expected, the results given by them are the same.

Edit (2017-10-02)

As suggested by Edward Villegas, I decided to compare execution times. I used %timeit for IPython and @benchmark (from BenchmarkTools) for Julia.

In IPython, we have

%timeit bisection(fun, 0.0, 2.0)

with result

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

And in Julia, we have

@benchmark bisection(fun, 0.0, 2.0)

with result

BenchmarkTools.Trial:
  memory estimate:  48 bytes
  allocs estimate:  2
  --------------
  minimum time:     1.505 μs (0.00% GC)
  median time:      1.548 μs (0.00% GC)
  mean time:        1.604 μs (0.00% GC)
  maximum time:     38.425 μs (0.00% GC)
  --------------
  samples:          10000
  evals/sample:     10

And it seems that the Julia version is around 100 times faster than the Python counterpart.

Comments

Comments powered by Disqus