# 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.