How to apply a Galerkin finite element method to a linear, one-dimensional boundary value problem

Question

I have the following boundary value problem:

$$-(\alpha u')' + \gamma u = f $$

in $\Omega = (a,b)$ with b.c. $u(a) = u(b) = 0$

and $\alpha > 0, \gamma ≥ 0$ and $f:(a,b) \to \Re $

The weak formulation I got is

$$ \int_a^b \alpha u' v' + \int_a^b \gamma uv = \int_a^b fv~~~~~\forall v\in V=H^1_0(a,b)$$

Then to discretize this, using Galerkin's method

$$ find~~u_h \in V_h:a(u_h,v_h)=F(v_h)~~~\forall v \in V_h \subset V $$

And define $V_h$ as

$$V_h = \Big\{ v_h \in X^1_h:v_h(a)=v_h(b)=0\Big\} $$

So I have a basis $ \phi_i(x) = \cases{ \frac{x-x_{i-1}}{x_i-x_{i-1}},~~~~x_{i-1} ≤ x ≤ x_i \cr \cr \frac{x_{i+1}-x}{x_{i+1}-x_i},~~~~x_i ≤ x ≤ x_{i+1} \cr \cr 0,~~~~~~~~~~~~otherwise}$

Now after all this, how can I build the problem in the form $Au = f$ ?? I'm really confused when trying to do this, can someone explain to me in a clearly and basic way, it's the first time I see FEM. If you have a link to any reference (books, pdfs, videos) that can help would be perfect.

Answer 1

Constructing the matrix $A$ and vector $f$ from the (bi)linear forms $a$ and $F$ is called assembly. The key observation is that since $a$ is bilinear and $V_h$ is finite-dimensional, you can substitute $$ u_h(x) = \sum_{i=1}^N u_i \phi_i(x) $$ (where $u_i$ is the $i$th component of the solution vector $u$ and $N$ is the number of nodes $a=x_1<\dots<x_N=b$) and only need to test with all basis functions $\phi_1,\dots,\phi_N$ instead of all $v_h$. Doing this and pulling the sums out of the bilinear form yields $$ \sum_{i=1}^N a(\phi_i,\phi_j)u_i = F(\phi_j) \qquad \text{for all }j=1,\dots,N. $$ This is a linear system for $u=(u_1,\dots,u_N)^T$ with the matrix $A=(a_{ij})$ and right hand side $F=(F_j)$ where $$ a_{ij} = \int_a^b \alpha\phi_i'\phi_j'\,dx + \int_a^b\gamma \phi_i\phi_j\,dx,\qquad F_j = \int_a^b f \phi_j\,dx. $$ Since the $\phi_j$ are given analytically, you can compute these integrals explicitly for constant $\alpha$ and $\gamma$. If $\alpha$ or $\gamma$ is a function, you would approximate the integrals using Gauss quadrature. Note that the product $\phi_i\phi_j$ is zero unless $i-1 \leq j \leq i +1$, which you can use to reduce the number of entries to compute. For homogeneneous Dirichlet conditions, you can just fix the values of $u_1 = u_N = 0$ (since these are the only functions which are nonzero at the boundary).

The above procedure is called node-based assembly, which should suffice for your one-dimensional model problem. A better alternative in higher dimensions is element-based assembly, where instead of looping over all (non-trivial) pairs of $i$ and $j$, you loop over all elements $[x_i,x_{i+1}]$ and compute all non-zero integrals over this interval. Summing all contributions for each pair $(i,j)$ gives you $a_{ij}$.

This is explained in any standard textbook on finite element methods; one with a particular focus on the algorithmic side (and less on the mathematical theory) is Mark Gockenbach, Understanding and Implementing the Finite Element Method, SIAM, 2006. For an even simpler test problem in 1D, both approaches are explained in Chapter 1.4 of these lecture notes.

Other good (mathematical) references which address implementation are

• Guermond Ern, Theory and Practice of Finite Elements, Springer 2004
• Dietrich Braess, Finite Elements, Cambridge University Press, 2007

Between these three books, you should get a good coverage of material. For self-study, I'd particularly recommend the last book, since it is quite readable and nicely presents the theory without too much technical detail (and there are worked solutions to some of the exercises - available from the link above - so you can check your answers).