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

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.

-
Are $\alpha$ and $\gamma$ constants or functions? Are they given analytically, or as vectors of point/cell values? –  Christian Clason Nov 27 '12 at 12:09
They are constants. –  BRabbit27 Nov 27 '12 at 12:20

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

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).

-
So in the sum $\sum_{i=1}^N a(\phi_i,\phi_j)u_i = F(\phi_j)$ the term $u_i$ is it a function? And what does the $\phi_j$ are? I mean what's the difference between $\phi_i$ and $\phi_j$. I'll check the book you say maybe with one worked solution I can understand better. –  BRabbit27 Nov 27 '12 at 13:07
$u_i$ is the $i$th component of the solution vector $u$ (in your notation). Similarly, $\phi_i$ is exactly as given in your question; $\phi_j$ is defined in the same way, but with the node $x_j$ in place of $x_i$. Both $i$ and $j$ run over all possible values from $1$ to $N$. –  Christian Clason Nov 27 '12 at 13:17
Thanks !! Finally I think I got it !! –  BRabbit27 Nov 27 '12 at 14:20
In general, $u_i$ is the coefficient of the ith basis function. Since lagrange basis functions (like the one you are using) have a value of 1 at the ith node, then each $u_i$ corresponds to the value of the function at the ith node. It is possible to use basis functions whose values are not unity at the ith node, in which case the coefficient $u_i$ is no longer the function value. Instead, we simply interpret the solution as a linear combination of the chosen basis functions. –  Paul Nov 27 '12 at 14:41