# Numerical methods challenge: Day 25

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.

## Finite element method

Today we have the Finite element method to solve the equation:

\begin{equation*} \nabla^2 u = -f(x) \end{equation*}

with

\begin{equation*} u(\sqrt{3}, 1) = 0 \end{equation*}

As in the Ritz method we form a functional that is equivalent to the differential equation, propose an approximation as a linear combination of a set of basis functions and find the best set of coefficients for that combination. That best solution is found minimizing the functional.

The functional for this differential equation is

\begin{equation*} \Pi[u] = -\int_\Omega \left(\nabla u\right)^2 d\Omega -\int_\Omega u f(x) d\Omega \end{equation*}

Using linear triangles, the local stiffness matrices read

\begin{equation*} K_\text{local} = \frac{1}{2|J|} \begin{bmatrix} 2 & -1 &1\\ -1 & 1 &0\\ -1 & 0 &1\\ \end{bmatrix} \end{equation*}

and

\begin{equation*} b_\text{local} = -|J| f(x_m) \mathbf{1}_3\, , \end{equation*}

where $|J|$ is the Jacobian determinant of the transformation and $x_m$ is the centroid of the triangle. I am skipping a great deal about assembling, but it would be just too extensive to describe the complete process.

We will solve the problem in a mesh with 6 triangles forming a regular hexagon

Following are the codes.

### Python

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

def assem(coords, elems, source):
ncoords = coords.shape
stiff = np.zeros((ncoords, ncoords))
rhs = np.zeros((ncoords))
for el_cont, elem in enumerate(elems):
stiff_loc, jaco = local_stiff(coords[elem])
rhs[elem] += jaco*np.mean(source[elem])
for row in range(3):
for col in range(3):
row_glob, col_glob = elem[row], elem[col]
stiff[row_glob, col_glob] += stiff_loc[row, col]
return stiff, rhs

def local_stiff(coords):
dHdr = np.array([[-1, 1, 0], [-1, 0, 1]])
jaco = dHdr.dot(coords)
stiff = 0.5*np.array([[2, -1, -1], [-1, 1, 0], [-1, 0, 1]])
return stiff/np.linalg.det(jaco), np.linalg.det(jaco)

sq3 = np.sqrt(3)
coords = np.array([
[sq3, -1],
[0, 0],
[2*sq3, 0],
[0, 2],
[2*sq3, 2],
[sq3, 3],
[sq3, 1]])
elems = np.array([
[1, 0, 6],
[0, 2, 6],
[2, 4, 6],
[4, 5, 6],
[5, 3, 6],
[3, 1, 6]])
source = -np.ones(7)
stiff, rhs = assem(coords, elems, source)
free = range(6)
sol = np.linalg.solve(stiff[np.ix_(free, free)], rhs[free])
sol_c = np.zeros(coords.shape)
sol_c[free] = sol
plt.tricontourf(coords[:, 0], coords[:, 1], sol_c, cmap="hot")
plt.colorbar()
plt.axis("image")
plt.show()


### Julia

using PyPlot

function assem(coords, elems, source)
ncoords = size(coords)
nelems = size(elems)
stiff = zeros(ncoords, ncoords)
rhs = zeros(ncoords)
for el_cont = 1:nelems
elem = elems[el_cont, :]
stiff_loc, jaco = local_stiff(coords[elem, :])
rhs[elem] += jaco*mean(source[elem])
for row = 1:3
for col = 1:3
row_glob = elem[row]
col_glob = elem[col]
stiff[row_glob, col_glob] += stiff_loc[row, col]
end
end
end
return stiff, rhs
end

function local_stiff(coords)
dHdr = [-1 1 0; -1 0 1]
jaco = dHdr * coords
stiff = 0.5*[2 -1 -1; -1 1 0; -1 0 1]
return stiff/det(jaco), det(jaco)
end

sq3 = sqrt(3)
coords =[sq3 -1;
0 0;
2*sq3 0;
0 2;
2*sq3 2;
sq3 3;
sq3 1]
elems =[2 1 7;
1 3 7;
3 5 7;
5 6 7;
6 4 7;
4 2 7]
source = -ones(7)
stiff, rhs = assem(coords, elems, source)
free = 1:6
sol = stiff[free, free] \ rhs[free]
sol_c = zeros(size(coords))
sol_c[free] = sol
tricontourf(coords[:, 1], coords[:, 2], sol_c, cmap="hot")
colorbar()
axis("image")
show()


Both have the same result, as follows ### Comparison Python/Julia

Regarding number of lines we have: 51 in Python and 56 in Julia. The comparison in execution time is done with %timeit magic command in IPython and @benchmark in Julia. For the test we are just comparing the time it takes to generate the matrices.

For Python:

%timeit assem(coords, elems, source)


with result

1000 loops, best of 3: 671 µs per loop


For Julia:

@benchmark assem(coords, elems, source)


with result

BenchmarkTools.Trial:
memory estimate:  13.30 KiB
allocs estimate:  179
--------------
minimum time:     7.777 μs (0.00% GC)
median time:      8.934 μs (0.00% GC)
mean time:        10.810 μs (14.54% GC)
maximum time:     797.432 μs (95.85% GC)
--------------
samples:          10000
evals/sample:     4


In this case, we can say that the Python code is about 70 times slower than Julia code.