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 secondorder reactiondiffusion problems of the form
(1) 
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 illposed if one formally sets . As approaches this limit, the solution typically exhibits layers: regions of rapid change, whose length is determined by . The overarching 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 layeradapted 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 finitedifference methods for (1) and its many variants is largely complete [13].
Finiteelement methods (FEMs) applied on layeradapted 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 firstorder 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 firstorder system PetrovGalerkin (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.

Those that propose new formulations for which the associated norm is naturally “balanced”; see, e.g., the discontinuous PetrovGalerkin 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 leastsquares FEM to obtain a symmetric discrete system. In this firstorder system leastsquares (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 leastsquares 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 firstorder system reformulation of (1), writing the equivalent system as
(2) 
for . Rather than forming a leastsquares finiteelement 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
(3) 
In [10], it is shown that is coercive and continuous with respect to the norm,
(4) 
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 firstorder system approach proposed by Lin and Stynes to include a curl constraint, in the same style as [5], leading to the firstorder system reformulation of (1) as
(5) 
Then, writing
(6) 
leads to the weak form
(7) 
Building on the theory of [10], this form is shown to be coercive and continuous with respect to the balanced norm
(8) 
Furthermore, in [1]
, the authors show that, when discretized using piecewise bilinear finite elements on a tensorproduct Shishkin mesh, this weak form leads to a
parameterrobustdiscretization, with an error estimate independent of the perturbation parameter
.3 Firstorder system least squares finiteelement methods
While theoretical and numerical results in [1] show the effectiveness of the firstorder system PetrovGalerkin approach proposed therein, the nonsymmetric 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 finiteelement 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 leastsquares 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, writingwith the associated norm written as . Slightly reweighting the firstorder system from (5), we have
(9) 
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 .
Proof
For the continuity bound, we note that
Thus, by the CauchySchwarz 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
Thus,
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
(10) 
Then,
for every .
Proof
A direct calculation shows that
Consequently,
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 quasioptimal 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
(11a)  
where each is a layer associated with the edge and each is a layer associated with the corner . There exists a constant such that  
(11b)  
(11c)  
(11d)  
(11e) 
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 wellrepresented 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 wellapproximated 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 wellrepresented 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 wellrepresented 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 twodimensional 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 tensorproduct Shishkin mesh (see, e.g., [1, §3] for more details). To do this, we select a transition point, , and construct a onedimensional mesh with equalsized elements on each of the intervals and . The twodimensional mesh is created as a tensorproduct 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 finiteelement 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 superconvergencetype 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 
(Reduction Rate w.r.t. N)  

/  32  64  128  256  512 
3.086e01  1.921e01 (0.62)  1.137e01 (0.59)  6.531e02 (0.57)  3.680e02 (0.56)  
3.086e01  1.921e01 (0.62)  1.137e01 (0.59)  6.532e02 (0.57)  3.681e02 (0.56)  
3.086e01  1.921e01 (0.62)  1.137e01 (0.59)  6.533e02 (0.57)  3.681e02 (0.56)  
3.086e01  1.921e01 (0.62)  1.137e01 (0.59)  6.533e02 (0.57)  3.681e02 (0.56)  
(Reduction Rate w.r.t. N)  
/  32  64  128  256  512 
6.935e02  1.981e02 (0.29)  6.436e03 (0.32)  2.051e03 (0.32)  6.448e04 (0.31)  
6.945e02  1.983e02 (0.29)  6.444e03 (0.32)  2.053e03 (0.32)  6.455e04 (0.31)  
6.946e02  1.984e02 (0.29)  6.445e03 (0.32)  2.054e03 (0.32)  6.456e04 (0.31)  
6.946e02  1.984e02 (0.29)  6.445e03 (0.32)  2.054e03 (0.32)  6.456e04 (0.31) 
(Reduction Rate w.r.t. N)  

/N  32  64  128  256  512 
9.307e02  3.854e02 (0.41)  1.394e02 (0.36)  4.655e03 (0.33)  1.484e03 (0.32)  
9.306e02  3.854e02 (0.41)  1.394e02 (0.36)  4.656e03 (0.33)  1.485e03 (0.32)  
9.306e02  3.854e02 (0.41)  1.394e02 (0.36)  4.656e03 (0.33)  1.485e03 (0.32)  
9.306e02  3.854e02 (0.41)  1.394e02 (0.36)  4.656e03 (0.33)  1.485e03 (0.32)  
(Reduction Rate w.r.t. N)  
/N  32  64  128  256  512 
1.512e02  1.817e03 (0.12)  2.715e04 (0.15)  3.609e05 (0.13)  4.133e06 (0.11)  
1.518e02  1.823e03 (0.12)  2.730e04 (0.15)  3.622e05 (0.13)  4.145e06 (0.11)  
1.519e02  1.823e03 (0.12)  2.733e04 (0.15)  3.624e05 (0.13)  4.147e06 (0.11)  
1.519e02  1.823e03 (0.12)  2.734e04 (0.15)  3.625e05 (0.13)  4.147e06 (0.11) 
(Reduction Rate w.r.t. N)  

/N  32  64  128  256  512 
2.786e02  7.800e03 (0.28)  1.748e03 (0.22)  3.419e04 (0.20)  6.185e05 (0.18)  
2.785e02  7.800e03 (0.28)  1.749e03 (0.22)  3.420e04 (0.20)  6.187e05 (0.18)  
2.785e02  7.800e03 (0.28)  1.749e03 (0.22)  3.420e04 (0.20)  6.187e05 (0.18)  
2.785e02  7.800e03 (0.28)  1.749e03 (0.22)  3.420e04 (0.20)  6.187e05 (0.18)  
(Reduction Rate w.r.t. N)  
/N  32  64  128  256  512 
4.364e03  6.989e04 (0.16)  9.807e05 (0.14)  1.148e05 (0.12)  1.198e06 (0.10)  
4.370e03  6.993e04 (0.16)  9.812e05 (0.14)  1.149e05 (0.12)  1.199e06 (0.10)  
4.371e03  6.994e04 (0.16)  9.813e05 (0.14)  1.149e05 (0.12)  1.199e06 (0.10)  
4.371e03  6.994e04 (0.16)  9.813e05 (0.14)  1.149e05 (0.12)  1.199e06 (0.10) 
5 Conclusions
In the paper, we propose and analyse a new weightednorm firstorder system least squares methodology tuned for singularly perturbed reactiondiffusion 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 convectiondiffusion equations, and investigating efficient linear solvers for the resulting discretizations.References
 [1] Adler, J.H., MacLachlan, S., Madden, N.: A firstorder system PetrovGalerkin discretisation for a reactiondiffusion 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 firstorder 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 firstorder system least squares. Numerical Linear Algebra with Applications 19, 343–366 (2012)

[4]
Cai, Z., Lazarov, R., Manteuffel, T., McCormick, S.: Firstorder system least squares for secondorder partial differential equations: Part I. SIAM J. Numer. Anal. pp. 1785–1799 (1994)
 [5] Cai, Z., Manteuffel, T., McCormick, S.: Firstorder system least squares for secondorder partial differential equations. II. SIAM J. Numer. Anal. 34(2), 425–454 (1997). https://doi.org/10.1137/S0036142994266066
 [6] De Sterck, H., Manteuffel, T., McCormick, S., Nolting, J., Ruge, J., Tang, L.: Efficiencybased  and refinement strategies for finite element methods. Numer. Linear Algebra Appl. 15(23), 89–114 (2008). https://doi.org/10.1002/nla.567
 [7] Heuer, N., Karkulik, M.: A robust DPG method for singularly perturbed reactiondiffusion problems. SIAM J. Numer. Anal. 55(3), 1218–1242 (2017). https://doi.org/10.1137/15M1041304
 [8] Lee, E., Manteuffel, T.A., Westphal, C.R.: Weightednorm firstorder 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.: Weightednorm firstorder system leastsquares (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 reactiondiffusion problems. SIAM J. Numer. Anal. 50(5), 2729–2743 (2012). https://doi.org/10.1137/110837784
 [11] Liu, F., Madden, N., Stynes, M., Zhou, A.: A twoscale sparse grid method for a singularly perturbed reaction–diffusion problem in two dimensions. IMA J. Numer. Anal. 29(4), 986–1007 (2009). https://doi.org/10.1093/imanum/drn048
 [12] Melenk, J.M., Xenophontos, C.: Robust exponential convergence of FEM in balanced norms for singularly perturbed reactiondiffusion equations. Calcolo 53(1), 105–132 (2016). https://doi.org/10.1007/s100920150139y
 [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). https://doi.org/10.1142/9789814390743
 [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. SpringerVerlag, Berlin, second edn. (2008)
 [16] Russell, S., Stynes, M.: Balancednorm error estimates for sparse grid finite element methods applied to singularly perturbed reactiondiffusion problems. J. Numer. Math. 27(1), 37–55 (2019). https://doi.org/10.1515/jnma20170079
Comments
There are no comments yet.