Understanding Finite element method

Question

Suppose we have Poisson in 1D: $u'' = f(t)$ where $0<t<1$ and $u(0)=0$ and $u(1)=1$

We approximate the solution by $U(t) \approx \sum_{i=1}^n x_i \phi_i(t) $ where $\phi_i(t)$ are some basis functions. Once we put this into our equation we get some residual $r(x_i,t) = U''(t) - f(t) $. The idea is to multiply by some ${\bf weight }$ function and solve

$$ \int\limits_0^1 r(x_i,t) w(t) d t $$

In the Galerkin method, we take $w(t) = \phi_i(t)$ and use integration by parts to derive a system of equations where we find $x_1,...,x_n$. The question is

Suppose the hat function is used as the shape function (piecewise linear shape function) and the domain is decomposed into subintervals of lenght $h=0.2$ Compute the element stiffness matrix ${\bf K}$

If I understood correctly, we use

$$ \phi_i(x) = \begin{cases} \frac{ x - x_{i-1} }{h} & x_{i-1} \leq x \leq x_i \\ \frac{x_{i+1}-x}{h} & x_i \leq x \leq x_{i+1} \end{cases} $$

as our basis functions. We are given that

$$ K_{ij} = \sum_{e} K_{ij}^e $$

where $K_{ij}^e$ is the element stiffness matrix for element $\Omega_e$. In our case, we have $5$ elements as $h=0.2$. For example, for the first element: $[x_0,x_1]$ we have

$$ K_{ij}^1 = \int\limits_{x_0}^{x_1} ( \phi_i' \phi_j' + \phi_i \phi_j) dx $$

Is this how finite element method works?

Answer 1

Let's take $n$ test functions $\phi_i\in H_0^1$ that are null at the domain boundaries, here $\phi_i(0)=\phi_i(1)=0$ for all $i=1,\dots,n$. You can always write:

$$ \int_0^1 u''(x) \phi_i(x) dx = \int_0^1 f(x) \phi_i(x) dx $$ Using integration by parts with $(u'(x) \phi_i(x))'=u''(x) \phi_i(x) + u'(x) \phi'_i(x)$ (or more generally some Green formula) you can write $$ \int_0^1 u''(x)\phi_i(x) dx = u'(1) \phi_i(1) - u'(0) \phi_i(0) -\int_0^1 u'(x)\phi'_i(x) dx = -\int_0^1 u'(x)\phi'_i(x) dx $$

Now to manage the Dirichlet boundary conditions: $u(0)=0$ and $u(1)=1$ you introduce any function $\tilde{u}(x)$ which satisfies $\tilde{u}(0)=0$ and $\tilde{u}(1)=1$ (a simple choice is $\tilde{u}(x)=x$ for instance). By linearity you can search a solution of the form $u(x)=u_0(x)+ \tilde{u}(x)$ where $u_0(0)=0$ and $u_0(1)=0$ (homogeneous Dirichlet conditions).

By construction $u_0$ is null at boundaries hence you can also expand it using the same basis functions $\phi_i\in H_0^1$. $$ u_0(x)=\sum\limits_{j=1}^n \alpha_j \phi_j(x) $$ Back to our previous equation $$ -\int_0^1 u'(x)\phi'_i(x) dx = \int_0^1 f(x)\phi_i(x)dx $$ we introduce $u(x)=u_0(x)+ \tilde{u}(x)$ which gives: $$ -\int_0^1 u_0'(x)\phi'_i(x) dx = \int_0^1 f(x)\phi_i(x)dx+\int_0^1 \tilde{u}'(x)\phi'_i(x)dx $$ $$ \sum_{j=1}^n\left(-\int_0^1 \phi'_j(x)\phi'_i(x) dx\right)\alpha_j = \int_0^1 f(x)\phi_i(x)dx+\int_0^1 \tilde{u}'(x)\phi'_i(x)dx $$ At the end you get a symmetric definite positive linear system of $n$ equations: $$ A\alpha = F $$ where $i=1,\dots,n$ and $j=1,\dots,n$, $$A_{i,j}=-\int_0^1 \phi'_j(x) \phi'_i(x) dx$$ and $i=1,\dots,n$ $$F_i=\int_0^1 f(x)\phi_i(x)dx+\int_0^1 \tilde{u}'(x)\phi'_i(x)dx$$

You have a lot of freedom to chose the basis function $\phi_i$. However it is generally convenient to split the spacial domain into simplices (the mesh) and use continuous piece-wise linear functions on these elements (what we call Lagrange elements). This is the basis functions you used in your question. And yes this approach is what we call the Finite Element Method.

Some clarifications concerning $\phi_i\in H_0^1$

• for Dirichlet conditions we chose $\phi_i\in H_0^1$, which are null at the boundary, to cancel the term: $u'(1) \phi_i(1) - u'(0) \phi_i(0)$. Reason? we have absolutly no idea of the value of $u'$ at the boundary, so it is helpful to be able to cancel this term.

• if instead of Dirichlet conditions we have a problem with Neumann conditions where you impose $u'(0)$ and $u'(1)$, we would have taken $\phi_i \in H^1$ (that can take any value at the boundaries). The situation is even simpler here because we do not have to introduce $\tilde{u}$ anymore and the Neuman conditions are introduced in the term $u'(1) \phi_i(1) - u'(0) \phi_i(0)$ which is then moved in the right hand side of the equation.

See for instance: Advanced Finite Element Methods for further details