First-order system least squares finite-elements for singularly perturbed reaction-diffusion equations

We propose a new first-order-system least squares (FOSLS) finite-element discretization for singularly perturbed reaction-diffusion equations. Solutions to such problems feature layer phenomena, and are ubiquitous in many areas of applied mathematics and modelling. There is a long history of the development of specialized numerical schemes for their accurate numerical approximation. We follow a well-established practice of employing a priori layer-adapted meshes, but with a novel finite-element method that yields a symmetric formulation while also inducing a so-called "balanced" norm. We prove continuity and coercivity of the FOSLS weak form, present a suitable piecewise uniform mesh, and report on the results of numerical experiments that demonstrate the accuracy and robustness of the method.



There are no comments yet.



Singularly perturbed reaction-diffusion problems as first order systems

We consider a singularly perturbed reaction diffusion problem as a first...

A convection-diffusion problem with a small variable diffusion coefficient

Consider a singularly perturbed convection-diffusion problem with small,...

A discontinuous least squares finite element method for time-harmonic Maxwell equations

We propose and analyze a discontinuous least squares finite element meth...

Multilevel Schwarz preconditioners for singularly perturbed symmetric reaction-diffusion systems

We present robust and highly parallel multilevel non-overlapping Schwarz...

A Boundary-Layer Preconditioner for Singularly Perturbed Convection Diffusion

Motivated by a wide range of real-world problems whose solutions exhibit...

Optimization of two-level methods for DG discretizations of reaction-diffusion equations

We analyze and optimize two-level methods applied to a symmetric interio...
This week in AI

Get the week's most popular data science and artificial intelligence research sent straight to your inbox every Saturday.

1 Introduction

The numerical solution of singularly perturbed differential equations (SPDEs) is of great interest to numerical analysts, given the importance of these equations in computational modelling, and the challenges they present for classical numerical schemes and the mathematical methods used to analyse them; see [15] for a survey. In this work, we focus on linear second-order reaction-diffusion problems of the form


for , where we assume there exist constants for every . Like all SPDEs, (1) is characterised by a small positive parameter that multiplies the highest derivative. It is “singular” in the sense that the problem is ill-posed if one formally sets . As approaches this limit, the solution typically exhibits layers: regions of rapid change, whose length is determined by . The over-arching goal is to devise methods that resolve these layers, and for which the error (measured in a suitable norm) is independent of . Many classical techniques make the tacit assumption that derivatives of are bounded, which does not hold, uniformly in , for solutions to (1). Numerous specialised methods, usually based around layer-adapted meshes, have been developed with the goal of resolving these layers and the attendant mathematical conundrums. The celebrated piecewise uniform meshes of Shishkin have been particularly successful in this regard; and analyses of finite-difference methods for (1) and its many variants is largely complete [13].

Finite-element methods (FEMs) applied on layer-adapted meshes have also been successfully applied to (1), but their analysis is more problematic. This is highlighted to great effect by Lin and Stynes who demonstrated that the usual energy norm associated with (1) is too weak to adequately express the layers present in the solution [10]. They proposed a first-order FEM (see §2) for which the associated norm is sufficiently strong to capture layers; they coined the term “balanced norm” to describe this.

A flurry of activity on balanced norms was prompted by [10], including the first-order system Petrov-Galerkin (FOSPeG) approach proposed by the authors [1]

, and we refer it its introduction for a survey of the progress up to 2015. Since then, developments have continued apace. Broadly speaking, studies can be classified as one of two types.

  1. Those that give analyses of standard FEMs, but in norms that are not induced by the associated bilinear forms; see, e.g., [16] on sparse grid FEMs, and [12] on -FEMs.

  2. Those that propose new formulations for which the associated norm is naturally “balanced”; see, e.g., the discontinuous Petrov-Galerkin method of Heuer and Karulik [7].

The present study belongs to the second of these classes: we propose a new FEM for which the induced norm is balanced. This method is related to our earlier work [1], but instead uses a weighted least-squares FEM to obtain a symmetric discrete system. In this first-order system least-squares (FOSLS) approach [4, 5], care is taken in choosing the weight, so that the resulting norms are indeed balanced.

The remainder of the paper is outlined as follows. Section 2 gives a brief discussion on balanced norms, where the Lin and Stynes and FOSPeG methods are summarized. In Section 3, we discuss the weighted least-squares approach and provide the necessary analysis, which applies in one, two and three dimensions. In Section 4, we focus on the particular case of ; we present a suitable Shishkin mesh of the problem, and present numerical results that support our findings. Some concluding remarks are given in Section 5.

2 Balanced norms

In [10], Lin and Stynes propose a first-order system reformulation of (1), writing the equivalent system as


for . Rather than forming a least-squares finite-element discretization as in [4, 5], they choose to close the system in a nonsymmetric manner, defining and

then writing the solution of (1) as that of the weak form


In [10], it is shown that is coercive and continuous with respect to the norm,


which is shown to be a balanced norm for the problem, in the sense that all the components in (4) have the same order of magnitude.

In [1], the authors augmented the first-order system approach proposed by Lin and Stynes to include a curl constraint, in the same style as [5], leading to the first-order system reformulation of (1) as


Then, writing


leads to the weak form


Building on the theory of [10], this form is shown to be coercive and continuous with respect to the balanced norm


Furthermore, in [1]

, the authors show that, when discretized using piecewise bilinear finite elements on a tensor-product Shishkin mesh, this weak form leads to a


discretization, with an error estimate independent of the perturbation parameter


3 First-order system least squares finite-element methods

While theoretical and numerical results in [1] show the effectiveness of the first-order system Petrov-Galerkin approach proposed therein, the non-symmetric nature of the weak form also has disadvantages. Primary among these is that the weak form no longer can be used as an accurate and reliable error indicator, contrary to the common practice for FOSLS finite-element approaches [2, 3, 4, 5, 6]. Standard techniques to symmetrize the weak form in (7

) fail, however, either sacrificing the balanced nature of the norm (and, thus, any guarantee of parameter robustness of the resulting discretization) or coercivity or continuity of the weak form (destroying standard error estimates). Here, we propose a FOSLS approach for the problem in (

1), made possible by considering a weighted norm with spatially varying weight function. Weighted least-squares formulations have been used for a wide variety of problems including those with singularities due to the domain [8, 9].

To this end, we define the weighted inner product on both scalar and vector

spaces, writing

with the associated norm written as . Slightly reweighting the first-order system from (5), we have


and pose the weighted FOSLS weak form as

This form leads to a natural weighted product norm given by

As shown below, under a reasonable assumption on the weight function, , the FOSLS weak form is coercive and continuous with respect to this norm.

Theorem 3.1

Let be given such that there exists for which

for every , and let be given. Then,

for all .


For the continuity bound, we note that

Thus, by the Cauchy-Schwarz and triangle inequalities, we have

For the coercivity bound, we note

Now consider

where we use the fact that on the boundary in the integration by parts step. Note that

and, consequently, that


Finally, consider

By our assumption on ,

and, so,

This gives

A natural question, in light of this result, is whether a suitable choice of exists. We now give a concrete construction of one such family of functions, , for which the assumption above is satisfied. This family is constructed for the case of with boundary layers along each boundary adjacent to the origin (i.e., where for some ). The extension to boundary layers along all boundary faces is straightforward from the construction.

Theorem 3.2

Let be given, and define . Take



for every .


A direct calculation shows that


Note, however, that

This gives

Substituting in the chosen value for gives the stated result.

The final question to be resolved is whether as given in (10) is a “good” choice, in the sense of whether quasi-optimal approximation in the resulting norm is expected to give a good approximation to the layer structure in a typical solution. We consider the case of , the unit square. Following Lemmas 1.1 and 1.2 of [11], we require that the problem data satisfy the assumptions of [11, §2.1], specifically that and that vanishes at the corners of the domain. Denoting the four edges of the domain by , , numbered clockwise with the edge as , and the four corners of the domain by , , numbered clockwise with the origin as , we have the following result.

Lemma 1 ([11, Lemmas 1.1 and 1.2])

The solution of (1) can be decomposed as

where each is a layer associated with the edge and each is a layer associated with the corner . There exists a constant such that

with analogous bounds for , , , and .

Thus, as a “stereotypical” solution of (1) in the case where boundary layers only form along the edges and of , we can consider

Next, we check if is “balanced”, not only in the sense of all terms having the same order, but in addition that each component in the stereotypical solution above is well-represented in the norm. This means the norm can be bounded from above and below by -independent values, so that it is not seen as being well-approximated by zero in the norm (unless truly vanishingly small), nor that the norm blows up as . For this case, (10) simplifies as

and the checks rely on two direct calculations:

With this, assuming that is over a nontrivial fraction of the domain, we conclude that

because of the separable nature of the calculation. Thus, the regular part of the solution is well-represented in the norm.

For the layer term, we write and calculate from the above that

Noting that all derivatives of this term with respect to are zero and that , we compute

Again, this shows that the layer term is well-represented in the norm. Similar calculations show the same to be true for the layer and corner terms in the stereotypical solution.

4 Numerical Results

To test the above approach, we consider a two-dimensional problem with constant posed on the unit square. We construct a problem whose solution mimics the stereotypical solution discussed above, with two edge layers and one corner layer. Specifically, we choose so that the solution is

We note that this has somewhat more complex layer behaviour than the stereotypical solution, but still obeys the bounds of Lemma 1. Also, the solution is constructed so as to obey the homogeneous Dirichlet boundary conditions. For numerical stability, we rescale the equations by defining and making corresponding changes in weights to preserve the balanced nature of the norm. With this, we pick to match the powers of in the weighting terms of both and in , equivalent to taking above.

We discretize the test problem on a tensor-product Shishkin mesh (see, e.g., [1, §3] for more details). To do this, we select a transition point, , and construct a one-dimensional mesh with equal-sized elements on each of the intervals and . The two-dimensional mesh is created as a tensor-product of this mesh with itself, with rectangular (quadrilateral) elements. For the choice of , we slightly modify the standard choice from the literature (see, for example, [1, 10, 11]) to account for both the layer functions present in the solution decomposition and in the definition of in (10). As such, we take

where is the degree of the polynomial space ( for bilinear elements, for biquadratic, and for bicubic), so that this factor matches the expected rate of convergence of the approximation, while the terms decrease appropriately as does, but increase (corresponding to increasing layer width) with decreases in or . In the results that follow, we take , implying . All numerical results were computed using Firedrake [14] for the discretization and a direct solver for the resulting linear systems.

Table 1 shows the expected reduction rates in errors with respect to the mesh parameter, , if we were to have standard estimates of approximation error in the -norm on the Shishkin meshes considered here. Tables 2, 3 and 4 show the measured errors (relative to the manufactured solution) for the bilinear, biquadratic, and bicubic discretizations, respectively. Expected behaviour for the bilinear case is a reduction like for (where represents the manufactured solution, and its gradient) and like for the discrete maximum norm of the error,

, which is measured at the nodes of the mesh corresponding to the finite-element degrees of freedom. These are both expected to be raised by one power in the biquadratic case, and a further one power for bicubics. In Tables

2, 3, and 4, we see convergence behaviour comparable to these rates, with the exception of the results for the discrete maximum norm in Table 3. These seem to show a superconvergence-type phenomenon, although we have no explanation for this observation at present.

0.60 0.58 0.57 0.56
0.36 0.34 0.33 0.32
0.22 0.20 0.19 0.18
0.13 0.12 0.11 0.10
Table 1: Expected error reduction rates on a Shishkin mesh.
(Reduction Rate w.r.t. N)
/ 32 64 128 256 512
3.086e-01 1.921e-01 (0.62) 1.137e-01 (0.59) 6.531e-02 (0.57) 3.680e-02 (0.56)
3.086e-01 1.921e-01 (0.62) 1.137e-01 (0.59) 6.532e-02 (0.57) 3.681e-02 (0.56)
3.086e-01 1.921e-01 (0.62) 1.137e-01 (0.59) 6.533e-02 (0.57) 3.681e-02 (0.56)
3.086e-01 1.921e-01 (0.62) 1.137e-01 (0.59) 6.533e-02 (0.57) 3.681e-02 (0.56)
(Reduction Rate w.r.t. N)
/ 32 64 128 256 512
6.935e-02 1.981e-02 (0.29) 6.436e-03 (0.32) 2.051e-03 (0.32) 6.448e-04 (0.31)
6.945e-02 1.983e-02 (0.29) 6.444e-03 (0.32) 2.053e-03 (0.32) 6.455e-04 (0.31)
6.946e-02 1.984e-02 (0.29) 6.445e-03 (0.32) 2.054e-03 (0.32) 6.456e-04 (0.31)
6.946e-02 1.984e-02 (0.29) 6.445e-03 (0.32) 2.054e-03 (0.32) 6.456e-04 (0.31)
Table 2: -weighted norm and discrete max norm errors for model problem with bilinear discretization.
(Reduction Rate w.r.t. N)
/N 32 64 128 256 512
9.307e-02 3.854e-02 (0.41) 1.394e-02 (0.36) 4.655e-03 (0.33) 1.484e-03 (0.32)
9.306e-02 3.854e-02 (0.41) 1.394e-02 (0.36) 4.656e-03 (0.33) 1.485e-03 (0.32)
9.306e-02 3.854e-02 (0.41) 1.394e-02 (0.36) 4.656e-03 (0.33) 1.485e-03 (0.32)
9.306e-02 3.854e-02 (0.41) 1.394e-02 (0.36) 4.656e-03 (0.33) 1.485e-03 (0.32)
(Reduction Rate w.r.t. N)
/N 32 64 128 256 512
1.512e-02 1.817e-03 (0.12) 2.715e-04 (0.15) 3.609e-05 (0.13) 4.133e-06 (0.11)
1.518e-02 1.823e-03 (0.12) 2.730e-04 (0.15) 3.622e-05 (0.13) 4.145e-06 (0.11)
1.519e-02 1.823e-03 (0.12) 2.733e-04 (0.15) 3.624e-05 (0.13) 4.147e-06 (0.11)
1.519e-02 1.823e-03 (0.12) 2.734e-04 (0.15) 3.625e-05 (0.13) 4.147e-06 (0.11)
Table 3: -weighted norm and discrete max norm errors for model problem with biquadratic discretization.
(Reduction Rate w.r.t. N)
/N 32 64 128 256 512
2.786e-02 7.800e-03 (0.28) 1.748e-03 (0.22) 3.419e-04 (0.20) 6.185e-05 (0.18)
2.785e-02 7.800e-03 (0.28) 1.749e-03 (0.22) 3.420e-04 (0.20) 6.187e-05 (0.18)
2.785e-02 7.800e-03 (0.28) 1.749e-03 (0.22) 3.420e-04 (0.20) 6.187e-05 (0.18)
2.785e-02 7.800e-03 (0.28) 1.749e-03 (0.22) 3.420e-04 (0.20) 6.187e-05 (0.18)
(Reduction Rate w.r.t. N)
/N 32 64 128 256 512
4.364e-03 6.989e-04 (0.16) 9.807e-05 (0.14) 1.148e-05 (0.12) 1.198e-06 (0.10)
4.370e-03 6.993e-04 (0.16) 9.812e-05 (0.14) 1.149e-05 (0.12) 1.199e-06 (0.10)
4.371e-03 6.994e-04 (0.16) 9.813e-05 (0.14) 1.149e-05 (0.12) 1.199e-06 (0.10)
4.371e-03 6.994e-04 (0.16) 9.813e-05 (0.14) 1.149e-05 (0.12) 1.199e-06 (0.10)
Table 4: -weighted norm and discrete max norm errors for model problem with bicubic discretization.

5 Conclusions

In the paper, we propose and analyse a new weighted-norm first-order system least squares methodology tuned for singularly perturbed reaction-diffusion equations that lead to boundary layers. The analysis includes a standard ellipticity result for the FOSLS formulation in a weighted norm, and shows that this norm is suitably weighted to be considered a “balanced norm” for the problem. Numerical results confirm the effectiveness of the method. Future work includes completing the error analysis by proving the necessary interpolation error estimates, with respect to

, investigating the observed superconvergence properties, generalizing the theory to convection-diffusion equations, and investigating efficient linear solvers for the resulting discretizations.


  • [1] Adler, J.H., MacLachlan, S., Madden, N.: A first-order system Petrov-Galerkin discretisation for a reaction-diffusion problem on a fitted mesh. IMA J. Numer. Anal. 36(3), 1281–1309 (2016)
  • [2] Berndt, M., Manteuffel, T.A., McCormick, S.F.: Local error estimates and adaptive refinement for first-order system least squares (FOSLS). Electron. Trans. Numer. Anal. 6, 35–43 (1997)
  • [3] Brezina, M., Garcia, J., Manteuffel, T., McCormick, S., Ruge, J., Tang, L.: Parallel adaptive mesh refinement for first-order system least squares. Numerical Linear Algebra with Applications 19, 343–366 (2012)
  • [4]

    Cai, Z., Lazarov, R., Manteuffel, T., McCormick, S.: First-order system least squares for second-order partial differential equations: Part I. SIAM J. Numer. Anal. pp. 1785–1799 (1994)

  • [5] Cai, Z., Manteuffel, T., McCormick, S.: First-order system least squares for second-order partial differential equations. II. SIAM J. Numer. Anal. 34(2), 425–454 (1997).
  • [6] De Sterck, H., Manteuffel, T., McCormick, S., Nolting, J., Ruge, J., Tang, L.: Efficiency-based - and -refinement strategies for finite element methods. Numer. Linear Algebra Appl. 15(2-3), 89–114 (2008).
  • [7] Heuer, N., Karkulik, M.: A robust DPG method for singularly perturbed reaction-diffusion problems. SIAM J. Numer. Anal. 55(3), 1218–1242 (2017).
  • [8] Lee, E., Manteuffel, T.A., Westphal, C.R.: Weighted-norm first-order system least squares (FOSLS) for problems with corner singularities. SIAM J. Numer. Anal. 44(5), 1974–1996 (2006)
  • [9] Lee, E., Manteuffel, T.A., Westphal, C.R.: Weighted-norm first-order system least-squares (FOSLS) for div/curl systems with three dimensional edge singularities. SIAM J. Numer. Anal. 46(3), 1619–1639 (2008)
  • [10] Lin, R., Stynes, M.: A balanced finite element method for singularly perturbed reaction-diffusion problems. SIAM J. Numer. Anal. 50(5), 2729–2743 (2012).
  • [11] Liu, F., Madden, N., Stynes, M., Zhou, A.: A two-scale sparse grid method for a singularly perturbed reaction–diffusion problem in two dimensions. IMA J. Numer. Anal. 29(4), 986–1007 (2009).
  • [12] Melenk, J.M., Xenophontos, C.: Robust exponential convergence of -FEM in balanced norms for singularly perturbed reaction-diffusion equations. Calcolo 53(1), 105–132 (2016).
  • [13] Miller, J.J.H., O’Riordan, E., Shishkin, G.I.: Fitted numerical methods for singular perturbation problems. World Scientific Publishing Co. Pte. Ltd., Hackensack, NJ, revised edn. (2012).
  • [14] Rathgeber, F., Ham, D.A., Mitchell, L., Lange, M., Luporini, F., McRae, A.T., Bercea, G.T., Markall, G.R., Kelly, P.H.: Firedrake: automating the finite element method by composing abstractions. ACM Transactions on Mathematical Software (TOMS) 43(3),  24 (2017)
  • [15] Roos, H.G., Stynes, M., Tobiska, L.: Robust numerical methods for singularly perturbed differential equations, Springer Series in Computational Mathematics, vol. 24. Springer-Verlag, Berlin, second edn. (2008)
  • [16] Russell, S., Stynes, M.: Balanced-norm error estimates for sparse grid finite element methods applied to singularly perturbed reaction-diffusion problems. J. Numer. Math. 27(1), 37–55 (2019).