Problem Formulation

This package solves the following 2D ordinary differential equation using the finite elements method:

Strong Formulation

Let:

$\Omega = [0,1] \times [0,1] \subset \mathbb{R}^2$;
$\Gamma$ be the boundary of $\Omega$;
$\hat{\Omega} = \Gamma \cup \Omega$;
$\Delta u(x,y) = u_{xx}(x,y) + u_{yy}(x,y)$.

Given $\alpha > 0$, $\beta \geq 0$, and a function $f : \Omega \to \mathbb{R}$, we want to find $u : \hat{\Omega} \to \mathbb{R}$ such that:

\[\begin{cases} -\alpha \Delta u(x,y) + \beta u(x,y) = f(x,y) & \forall (x,y) \in \Omega \\ \\ u(x,y) = 0 & \forall (x,y) \in \Gamma. \end{cases}\]

Weak Formulation

Let:

$V$ represent the space of test functions;
$H$ represent the space of solution functions that satisfy the given variational formulation;
$v \in V$ be a test function;
$u \in H$ be a solution function;
$\displaystyle\kappa(f,g) = \alpha \int_\Omega \Delta f(x) \Delta g(x) d\Omega + \beta \int_\Omega f(x)g(x) d\Omega$;
$\displaystyle(f,g) = \int_\Omega f(x) g(x) d\Omega$.

Given $\alpha > 0$, $\beta \geq 0$, and a function $f : \Omega \to \mathbb{R}$, we want to find $u \in H, u : \hat{\Omega} \to \mathbb{R}$ such that:

\[\begin{cases} \kappa(u,v) = (f,v) & \forall (x,y) \in \Omega \\ \\ u(x,y) = 0 & \forall (x,y) \in \Gamma. \end{cases}\]

Approximate Problem

Let:

$V_h$ represent the finite-dimensional space of approximate test functions;
$H_h$ represent the finite-dimensional space of approximate solution functions that satisfy the discrete variational formulation;
$v_h \in V_h$ be an approximate test function;
$u_h \in H_h$ be an approximate solution function.

Given $\alpha > 0$, $\beta \geq 0$, and a function $f : \Omega \to \mathbb{R}$, we want to find $u_h \in H_h, u_h : \hat{\Omega} \to \mathbb{R}$ such that:

\[\begin{cases} \kappa(u_h,v_h) = (f,v_h) & \forall (x,y) \in \Omega \\ \\ u_h(x,y) = 0 & \forall (x,y) \in \Gamma. \end{cases}\]

Matrix Formulation

Let

$\varphi_i$ be the basis functions of the finite-dimensional space $H_h = V_h$.

Given a matrix $K$ and a vector $F$, we want to find $C$ such that $KC = F$

\[\begin{align*} \begin{bmatrix} \kappa\big(\varphi_1,\varphi_1\big)& \dots& \kappa\big(\varphi_m,\varphi_1\big) \\ \vdots& \ddots& \vdots \\ \kappa\big(\varphi_1,\varphi_m\big)& \dots& \kappa\big(\varphi_m,\varphi_m\big) \\ \end{bmatrix} \begin{bmatrix} c_1\\ \vdots\\ c_m \end{bmatrix} = \begin{bmatrix} (f,\varphi_1)\\ \vdots\\ (f,\varphi_m)\\ \end{bmatrix}. \end{align*}\]

More details on the calculations can be found at Two Dimensional Analysis and Calculations.

Function Signatures

Functions signatures from the FiniteElements2dDirichlet package.

EQ, LG and m Initializers

FiniteElements2dDirichlet.init_LG_matrix — Method

init_LG_matrix(Nx, Ny)

Return the matrix that correlates the local to the global elements.

The size of the matrix is 4×Nx*Ny

Arguments

Nx::Int: the number of elements in the x axis.
Ny::Int: the number of elements in the y axis.

Examples

julia> init_LG_matrix(2, 4)
4×8 Matrix{Int64}:
 1  2  4  5   7   8  10  11
 2  3  5  6   8   9  11  12
 5  6  8  9  11  12  14  15
 4  5  7  8  10  11  13  14

FiniteElements2dDirichlet.init_EQ_vector_and_m — Method

init_EQ_vector_and_m(Nx, Ny)

Return the vector that defines which terms are important and the dominium dimension.

The size of the vector is (Nx+1)*(Ny+1) and m is (Nx-1)*(Ny-1).

Arguments

Nx::Integer: the number of elements in the x axis.
Ny::Integer: the number of elements in the y axis.

Examples

julia> init_EQ_vector_and_m(2, 5)
([5, 5, 5, 5, 1, 5, 5, 2, 5, 5, 3, 5, 5, 4, 5, 5, 5, 5], 4)

Fe and F Initializers

FiniteElements2dDirichlet.init_Fe_vector! — Method

init_Fe_vector!(f, Xs, Ys, Fe, P, W)

Mutates the given vector Fe (of length 4) to compute the local Fe vector for finite element analysis.

Arguments

f::Function: input function from the equation.
Xs::Vector{Float64}: vector of the element coordinates on the x axis.
Ys::Vector{Float64}: vector of the element coordinates on the y axis.
Fe::Vector{Float64}: vector that is mutated.
P::Vector{Float64}: vector of gauss points for numerical integration.
W::Vector{Float64}: vector of gauss weights for numerical integration.

Examples

using GaussQuadrature
Fe = zeros(4);
P, W = legendre(5);
init_Fe_vector!((x1, x2) -> 64, [0, 0.25, 0.25, 0], [0, 0, 0.25, 0.25], Fe, P, W);
Fe;
# output
4-element Vector{Float64}:
 0.9999999999999993
 0.9999999999999998
 1.0000000000000002
 0.9999999999999998

FiniteElements2dDirichlet.init_F_vector — Method

init_F_vector(f, X_matrix, Y_matrix, m, EQ, LG)

Return a vector $F$ with length m for solving the differential equation.

Arguments

f::Function: input function from the equation.
X_matrix::Matrix{Float64}: x axis coordinates of the mesh generated by the init_mesh function.
Y_matrix::Matrix{Float64}: y axis coordinates of the mesh generated by the init_mesh function.
m::Int: m value generated by the EQ matrix generated by the initEQvectorandm function.
EQ::Matrix{Int}: EQ matrix generated by the initEQvectorandm function.
LG::Matrix{Int}: LG matrix generated by the initLGmatrix function.

Examples

X, Y = init_mesh(4, 3);
LG = init_LG_matrix(4, 3);
EQ, m = init_EQ_vector_and_m(4, 3);
init_F_vector((x1, x2) -> 48, X, Y, m, EQ, LG);
# output
6-element Vector{Float64}:
 3.999999999999999
 3.999999999999999
 3.999999999999999
 3.999999999999999
 3.9999999999999996
 3.999999999999999

Ke and K Initializers

FiniteElements2dDirichlet.init_Ke_matrix! — Method

init_Ke_matrix!(alpha, beta, Xs, Ys, Ke, P, W)

Mutates a given matrix Ke with dimensions 4×4 to become the local K matrix.

Arguments

alpha::Float64: constant α from the equation.
beta::Float64: constant β from the equation.
Xs::Vector{Float64}: vector of the element coordinates on the x axis.
Ys::Vector{Float64}: vector of the element coordinates on the y axis.
Ke::Matrix{Float64}: matrix that suffer the mutation.
P::Vector{Float64}: vector of gauss points for numerical integration.
W::Vector{Float64}: vector of gauss weights for numerical integration.

Examples

using GaussQuadrature
Ke = zeros(4,4);
P, W = legendre(5);
init_Ke_matrix!(6, 0, [0, 0.25, 0.25, 0], [0, 0, 0.25, 0.25], Ke, P, W);
Ke;
# output
4×4 Matrix{Float64}:
  4.0  -1.0  -2.0  -1.0
 -1.0   4.0  -1.0  -2.0
 -2.0  -1.0   4.0  -1.0
 -1.0  -2.0  -1.0   4.0

FiniteElements2dDirichlet.init_K_matrix — Method

init_K_matrix(alpha, beta, X_matrix, Y_matrix, m, EQ, LG)

Return a matrix $K$ with dimensions m×m for solving the differential equation.

Arguments

alpha::Float64: constant α from the equation.
beta::Float64: constant β from the equation.
X_matrix::Matrix{Float64}: x axis coordinates of the mesh generated by the init_mesh function.
Y_matrix::Matrix{Float64}: y axis coordinates of the mesh generated by the init_mesh function.
m::Int: m value generated by the EQ matrix generated by the initEQvectorandm function.
EQ::Vector{Int}: EQ vector generated by the initEQvectorandm function.
LG::Matrix{Int}: LG matrix generated by the initLGmatrix function.

Examples

X, Y = init_mesh(3,3);
EQ, m = init_EQ_vector_and_m(3,3);
LG = init_LG_matrix(3,3);
init_K_matrix(1, 1, X, Y, m, EQ, LG);
# output
4×4 SparseArrays.SparseMatrixCSC{Float64, Int64} with 16 stored entries:
  2.71605   -0.320988  -0.320988  -0.330247
 -0.320988   2.71605   -0.330247  -0.320988
 -0.320988  -0.330247   2.71605   -0.320988
 -0.330247  -0.320988  -0.320988   2.71605

Mesh Initializer

FiniteElements2dDirichlet.init_mesh — Method

init_mesh(Nx, Ny; ns=false, plot=false)

Return the x and y axis mesh for solving the system.

Arguments

Nx::Int: the number of elements in the x axis.
Ny::Int: the number of elements in the y axis.
ns::Bool: indicates if noise should be added to the mesh.
plot::Bool: indicates if the mesh should be plotted.

Examples

jldoctest julia> init_mesh(2, 3) ([0.0 0.0 0.0 0.0; 0.5 0.5 0.5 0.5; 1.0 1.0 1.0 1.0], [0.0 0.3333333333333333 0.6666666666666666 1.0; 0.0 0.3333333333333333 0.6666666666666666 1.0; 0.0 0.3333333333333333 0.6666666666666666 1.0])

System Solver

FiniteElements2dDirichlet.solve_system — Method

solve_system(alpha, beta, f, Nx, Ny; EQLG=false, XY_matrix=false, noise=false)

Return the solution of the differential equation as a vector of the approximate solution function evaluated in the inner knots of the mesh.

Arguments

alpha::Float64: constant α from the equation.
beta::Float64: constant β from the equation.
f::Function: input function from the equation.
Nx::Int: the number of elements in the x axis.
Ny::Int: the number of elements in the y axis.
EQLG::Bool: indicates if the EQ vector and LG matrix should be returned.
XY_matrix::Bool: indicates if the X and Y matrices from the mesh should be returned.
noise::Bool: indicates if noise should be added to the mesh.

Examples

julia> solve_system(1, 1, (x, y) -> (2*pi^2 +1)*sin(pi*x)*sin(pi*y), 2, 3)
2-element Vector{Float64}:
 1.0040408191040564
 1.0040408191040564

Error Convergence

FiniteElements2dDirichlet.error_convergence — Method

error_convergence(lb, ub, alpha, beta, u, f; see_plot=false, ns=false)

Return the errors of a given solution for discretizations with knots ranging from 2^lb to 2^ub. The number of knots increases in powers of 2, starting from 2^lb up to 2^ub.

Arguments

lb::Int: lower-bound for the number of knots as a power of 2.
ub::Int: upper-bound for the number of knots as a power of 2.
alpha::Float64: constant α from the equation.
beta::Float64: constant β from the equation.
u::Function: the analytical solution of the equation.
f::Function: the input function from the equation.
see_plot::Bool: if true, plots the error convergence and saves it as error-convergence.png.
ns::Bool: if true, adds noise to the mesh.

Examples

julia> error_convergence(2, 4, 1, 1, (x,y) -> sin(pi * x) * sin(pi * y), (x,y) -> (2*pi^2 + 1) * sin(pi * x) * sin(pi * y))
3-element Vector{Float64}:
 0.02946436125614044
 0.007348142187147988
 0.001836025960324437