Oscillations in singularly perturbed reaction-diffusion problems with finite elements

Question

When FEM-discretizing and solving a reaction-diffusion problem, e.g., $$ - \varepsilon \Delta u + u = 1 \text{ on } \Omega\\ u = 0 \text{ on } \partial\Omega $$ with $0 < \varepsilon \ll 1$ (singular perturbation), the solution of the discrete problem will typically exhibit oscillatory layers close to the boundary. With $\Omega = (0,1)$, $\varepsilon=10^{-5}$ and linear finite elements, the solution $u_h$ looks like

I see there's a lot of literature out there for such unwanted effects when they are caused by convection (e.g., upwind discretizations), but when it comes to reaction, people seem to focus on refined meshes (Shishkin, Bakhvalov).

Are there discretizatons that avoid such oscillations, i.e., that preserve monotonicity? What else may be useful in this context?

Answer 1

In the case you show, the solution has a boundary layer. If you can't resolve it because your mesh is too coarse, then for all practical matters the solution is discontinuous to the numerical scheme.

Now, if you just apply a standard discretization to this problem, the discrete solution is the result of applying a linear projection operator to the exact solution, projecting onto a finite dimensional space. This is really no different than, for example, taking the first $N$ terms of a Fourier series. But there you know what's going to happen when you apply it to a discontinuous function: you get Gibbs phenomenon with over- and undershots. The situation here is really no different and it will happen with any linear scheme.

There are essentially two approaches you can follow to working around this problem: (i) You add artificial diffusion, i.e., you add a diffusion term that is proportional to the mesh size to some power; this increases the size of the boundary layer to a width that can be represented by the mesh, but it goes back to $\varepsilon$ as $h\rightarrow 0$. (ii) You use a nonlinear projection scheme; the hyperbolic literature has any number of such schemes, e.g., shock capturing, slope limiting, etc.

Answer 2

TL;DR: Your options are limited 1) go brute force adaptive for accurate and expensive solution 2) use numerical diffusion for a less accurate but stable solution or (my favorite) 3) leverage the fact that this is a singular perturbation problem and solve two inexpensive inner/outer problems and let matched asymptotics do its magic!

---

If you really must obtain a uniform numerical solution to the problem, there is really not much you could do beyond adaptive mesh refinement. You are facing a singular perturbation problem that develops a boundary layer of thickness $\delta = \mathcal{O}(\sqrt{\epsilon})$ near the boundary. This boundary layer separates the inner and outer solutions.

To be more precise, consider the left boundary, where $x = \mathcal{O}(\delta)$. In this region, it is possible to rescale the equation, by using stretched coordinate $\eta = x/\delta$, for the inner solution:

$$ -\Delta u_i + u_i = 1 $$

subjected to the boundary conditions $u(0) = 0$ and $u_i(\eta\rightarrow\infty) = u_o(x\rightarrow0)$. For the outer solution, $u_o$, one assumes that $x = \mathcal{O}(1)$ and thus the dominant balance is between $u$ and $1$ which simply gives $u_0 = 1$. With the outer solution at hand, you can now solve the regularized inner solution with ease -- in this case even analytically.

This is in fact the technique that was (and still is) very popular for solving laminar boundary layer problems in fluid mechanics back in the day. In fact if you look at the Navier-Stokes equations, at high Reynolds numbers, you are effectively facing a singular perturbation problem, which like the one you mentioned here, develops a boundary layer (fun fact: the terms "boundary layer" in perturbation analysis actually comes from the fluid boundary layer problem I just described).

So, back in the day where computational resources were very limited, the typical thing to do was to separate the outer and inner problems. For the outer, one would first solve the potential flow (analogous to $u_0 = 1$ here) and for the inner, well the boundary layer equations! All this without the needing fancy adaptive methods if you were to tackle the whole problem at once and in a single computational domain.