# Numerical methods challenge

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