# Numerical methods challenge: Day 21

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.

## Jacobi iteration

Today we have a Finite difference method combined with Jacobi method We are solving the Laplace equation.

\begin{equation*} \nabla^ 2 u = 0 \end{equation*}

with

\begin{align*} u(0, y) = 1 -y,\quad u(x, 0) = 1 - x,\\ u(1, y) = u(x, 1) = 0 \end{align*}

Following are the codes.

### Python

from __future__ import division, print_function
import numpy as np
import matplotlib.pyplot as plt

def jacobi(u, update, tol=1e-6, niter=500):
for n in range(niter):
u_new = update(u)
if np.linalg.norm(u_new - u) < tol:
break
else:
u = u_new.copy()
return u, n

def heat_FDM(u):
u_new = u.copy()
u_new[1:-1, 1:-1] = 0.25*(u_new[0:-2, 1:-1] + u_new[2:, 1:-1] +\
u_new[1:-1, 0:-2] + u_new[1:-1, 2:])
return u_new

nx = 50
ny = 50
x_vec = np.linspace(-0.5, 0.5, nx)
y_vec = np.linspace(-0.5, 0.5, ny)
u0 = np.zeros((nx, ny))
u0[:, 0] = 1 - x_vec
u0[0, :] = 1 - y_vec
nvec = [100, 1000, 10000, 100000]
for num, niter in enumerate(nvec):
u, n = jacobi(u0, heat_FDM, tol=1e-12, niter=niter)
plt.subplot(2, 2, num + 1)
plt.contourf(u, cmap='hot')
plt.xlabel(r"$x$")
plt.ylabel(r"$y$")
plt.title('{} iterations'.format(n + 1))
plt.axis('image')
plt.tight_layout()
plt.show()


### Julia

using PyPlot

function jacobi(u, update; tol=1e-6, niter=500)
num = niter
for n = 1:niter
u_new = update(u)
if norm(u_new - u) < tol
num = n
break
else
u = copy(u_new)
end
end
return u, num
end

function heat_FDM(u)
u_new = copy(u)
u_new[2:end-1, 2:end-1] = 0.25*(u_new[1:end-2, 2:end-1] +
u_new[3:end, 2:end-1] + u_new[2:end-1, 1:end-2] + u_new[2:end-1, 3:end])
return u_new
end

nx = 50
ny = 50
x_vec = linspace(-0.5, 0.5, nx)
y_vec = linspace(-0.5, 0.5, ny)
u0 = zeros(nx, ny)
u0[:, 1] = 1 - x_vec
u0[1, :] = 1 - y_vec
nvec = [100, 1000, 10000, 100000]
for (num, niter) = enumerate(nvec)
u, n = jacobi(u0, heat_FDM, tol=1e-12, niter=niter)
subplot(2, 2, num)
contourf(u, cmap="hot")
xlabel(L"$x$")
ylabel(L"$y$")
title("\$(n) iterations")
axis("image")
end
tight_layout()
show() ### Comparison Python/Julia

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

For Python:

%timeit jacobi(u0, heat_FDM, tol=1e-12, niter=1000)


with result

10 loops, best of 3: 29.6 ms per loop


For Julia:

@benchmark jacobi(u0, heat_FDM, tol=1e-12, niter=1000)


with result

BenchmarkTools.Trial:
memory estimate:  247.89 MiB
allocs estimate:  43002
--------------
minimum time:     196.943 ms (5.66% GC)
median time:      203.230 ms (5.74% GC)
mean time:        203.060 ms (6.01% GC)
maximum time:     208.017 ms (5.51% GC)
--------------
samples:          25
evals/sample:     1


In this case, we can say that the Python code is roughly 10 times faster than Julia.