1-D Poisson Equation with the Finite Element Method

Why start FEM with Poisson?

In previous Numerics Lab examples, we explored 1-D linear convection using Finite Differences (FDM), Finite Volumes (FVM), and Discontinuous Galerkin (DG) methods. While those are excellent for learning about waves and advection, the Finite Element Method (FEM) is best introduced with an elliptic boundary-value problem.

The Poisson equation is mathematically symmetric, positive-definite, and perfectly suited for the standard continuous Galerkin formulation. This makes it an ideal starting point to learn the core mechanics of FEM before tackling more complex, convection-dominated equations that require stabilization.

Strong form of the problem

The classical "strong form" of the 1-D Poisson equation relates the second spatial derivative of a field variable $u$ to a given source term $f(x)$:

We will solve this on the one-dimensional domain $x \in [0, 1]$.

Boundary conditions

To obtain a unique solution, we must specify boundary conditions. We will apply homogeneous Dirichlet boundary conditions, which fix the solution to zero at both ends of our domain:

Manufactured analytic solution

To verify our numerical implementation, we use the method of manufactured solutions. We propose an exact, continuous solution that perfectly satisfies our boundary conditions:

Plugging this into the strong form gives us the corresponding source term that we must feed into our solver:

From strong form to weak form

A crucial distinction of FEM is that it starts from a weak form rather than directly approximating derivatives (like FDM) or flux balances (like FVM).

We multiply the strong form by an arbitrary test function $v(x)$ (which is zero at the boundaries) and integrate over the domain. Using integration by parts, we shift one spatial derivative from our unknown $u$ onto the test function $v$. The boundary terms vanish due to our boundary conditions, resulting in the weak form:

The problem is now to find a function $u_h$ in a finite element space such that this integral equality holds for all test functions $v_h$ in that same space.

Linear Lagrange basis functions

We partition the domain into a uniform mesh of elements. Over each element, we approximate the solution as a piecewise-linear function using linear Lagrange basis functions. In a standard local coordinate system $\xi \in [-1, 1]$, these functions are simply straight lines (often called "hat" functions when connected across elements).

Element stiffness matrix

By substituting our basis functions into the weak form, the continuous integrals turn into discrete matrices. For a 1-D element of length $h$, the local element stiffness matrix $K_e$ is given by:

Assembly into the global system

A core mechanic of FEM is assembly. We loop over all individual elements and add their local stiffness matrices $K_e$ and load vectors $F_e$ into a large global stiffness matrix $K$ and global load vector $F$. The overlap of continuous basis functions at shared nodes naturally connects the adjacent elements together into one coupled system.

Applying Dirichlet boundary conditions

Once assembled, the global matrix is singular. We must enforce our homogeneous Dirichlet boundary conditions by fixing the values at the endpoints to zero. A simple and robust way to achieve this is to remove the rows and columns corresponding to the boundary nodes, solving a reduced linear system strictly for the interior unknown values.

Solving the linear system

After applying the boundary conditions, we are left with a symmetric, positive-definite linear system. Solving this system gives us the discrete nodal values of our numerical solution. Note: In this foundational example, we use standard dense linear algebra, but production FEM solvers leverage the matrix's sparsity for efficiency.

Error and convergence

To quantify accuracy, we measure the error between our discrete piecewise-linear solution and the exact continuous sine wave. By refining the mesh (decreasing $h$) and repeating the process, we observe that the $L_2$ error drops quadratically. This confirms that linear elements achieve optimal $\mathcal{O}(h^2)$ convergence for the Poisson equation.

Relationship to FDM, FVM, and DG examples

To understand where FEM fits in our previous examples:

• FDM: Discretizes the strong form directly, evaluating derivatives at grid points.
• FVM: Discretizes the integral conservation law, focusing on fluxes crossing cell faces.
• DG: Uses a weak formulation on independent elements, coupled by numerical interface fluxes.
• FEM: Uses a global weak formulation where elements are strongly coupled by shared nodes, requiring a global matrix solve.

What comes next: diffusion, advection-diffusion, and stabilized convection

With the weak form, element integration, matrix assembly, and boundary conditions established via this Poisson foundation, we are ready to tackle more complex physics. The natural next steps include adding transient diffusion, introducing a velocity field for advection-diffusion, and eventually utilizing specialized stabilization techniques to handle purely convective problems within the continuous FEM framework.