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:
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
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.
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] - coords[ele] 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(["FEM solution", "Exact solution"]) plt.tight_layout() plt.show()
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] - coords[ele] 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
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.
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.