# Numerical methods challenge: Day 24

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:

with

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

The main difference is that we use a piecewise interpolation for the basis functions,

this leads to the system of equations

where the local stiffness matrices read

and

where \(|J|\) is the Jacobian determinant of the transformation. I am skipping a great deal about assembling, but it would be just too extensive to describe the complete process.

We will test the implementation with the function \(f(x) = x^3\), that leads to the solution

Following are the codes.

### Python

from __future__ import division, print_function import numpy as np import matplotlib.pyplot as plt def FEM1D(coords, source): N = len(coords) stiff_loc = np.array([[2.0, -2.0], [-2.0, 2.0]]) eles = [np.array([cont, cont + 1]) for cont in range(0, N - 1)] stiff = np.zeros((N, N)) rhs = np.zeros(N) for ele in eles: jaco = coords[ele[1]] - coords[ele[0]] rhs[ele] = rhs[ele] + jaco*source(coords[ele]) for cont1, row in enumerate(ele): for cont2, col in enumerate(ele): stiff[row, col] = stiff[row, col] + stiff_loc[cont1, cont2]/jaco return stiff, rhs N = 100 fun = lambda x: x**3 x = np.linspace(0, 1, N) stiff, rhs = FEM1D(x, fun) sol = np.zeros(N) sol[1:-1] = np.linalg.solve(stiff[1:-1, 1:-1], -rhs[1:-1]) #%% Plotting plt.figure(figsize=(4, 3)) plt.plot(x, sol) plt.plot(x, x*(x**4 - 1)/20, linestyle="dashed") plt.xlabel(r"$x$") plt.ylabel(r"$y$") plt.legend(["Ritz solution", "Exact solution"]) plt.tight_layout() plt.show()

### Julia

using PyPlot function FEM1D(coords, source) N = length(coords) stiff_loc = [2.0 -2.0; -2.0 2.0] eles = [[cont, cont + 1] for cont in 1:N-1] stiff = zeros(N, N) rhs = zeros(N) for ele in eles jaco = coords[ele[2]] - coords[ele[1]] rhs[ele] = rhs[ele] + jaco*source(coords[ele]) stiff[ele, ele] = stiff[ele, ele] + stiff_loc/jaco end return stiff, rhs end N = 100 fun(x) = x.^3 x = linspace(0, 1, N) stiff, rhs = FEM1D(x, fun) sol = zeros(N) sol[2:end-1] = -stiff[2:end-1, 2:end-1]\rhs[2:end-1] #%% Plotting figure(figsize=(4, 3)) plot(x, sol) plot(x, x.*(x.^4 - 1)/20, linestyle="dashed") xlabel(L"$x$") ylabel(L"$y$") legend(["FEM solution", "Exact solution"]) tight_layout() show()

Both have the same result, as follows

### Comparison Python/Julia

Regarding number of lines we have: 37 in Python and 35 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 FEM1D(x, fun)

with result

100 loops, best of 3: 2.15 ms per loop

For Julia:

@benchmark FEM1D(x, fun)

with result

BenchmarkTools.Trial: memory estimate: 183.73 KiB allocs estimate: 1392 -------------- minimum time: 60.045 μs (0.00% GC) median time: 70.445 μs (0.00% GC) mean time: 98.276 μs (25.64% GC) maximum time: 4.269 ms (96.70% GC) -------------- samples: 10000 evals/sample: 1

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

## Comments

Comments powered by Disqus