# Robust Numerical Methods for Singularly Perturbed Differential Equations: A Survey Covering 2008–2012

Authors: Hans-Görg Roos
Journal: ISRN Applied Mathematics
 (2012)

---

## Abstract

We present new results in the numerical analysis of singularly perturbed convection-diffusion-reaction problems that have appeared in the last five years. Mainly discussing layer-adapted meshes, we present also a survey on stabilization methods, adaptive methods, and on systems of singularly perturbed equations.


---

## Body

## 1. Introduction

Since the publication of [1] many new ideas in the numerical analysis of singularly perturbed differential equations appeared. Some of them are contained in [2] on layer-adapted meshes, some mentioned in the (very short) overview [3] for the years 2000–2009. In this survey we try to sketch the most important new developments of the last five years, but, of course, the choice is inevitably personal and reflects our own main interests.We mostly discuss linear convection-diffusion-reaction problems in the case of partial differential equations and present results for boundary value problems of ordinary differential equations only in exceptional cases. One exception is Section5 on singularly perturbed systems. In Section 2 we describe mostly stabilization methods including discontinuous Galerkin on shape-regular meshes. Uniform estimates with respect to the perturbation parameter are discussed in Section 3, while adaptive methods based on a posteriori error estimation are the contents of Section 4.


## 2. Semi-Robust Methods

### 2.1. Steady-State Problems

We consider the two-dimensional convection-diffusion problem(2.1a)Lu:=-εΔu-b·∇u+cu=finΩ,(2.1b)u=0onΓ=∂Ω,whereε is a small positive parameter, b,c are smooth, and f∈L2(Ω). Assuming
(2.2)c+12divb≥α0>0,
the given problem admits a unique solution u∈H01(Ω), for convex Ω we have moreover u∈H2(Ω).Let us introduce theε-weighted H1 norm by
(2.3)∥v∥ε:=ε1/2|v|1+∥v∥0.
Then for the Galerkin finite element method with piecewise linear or bilinear elements one can prove (C denotes a generic constant that is independent of ε and of the mesh)
(2.4)∥u-uh∥ε≤Ch|u|2
on quite general triangulations. However, estimate (2.4) is of no worth: in general, |u|2 tends to infinity for ε→0 due to the presence of layers. The very weak stability properties of standard Galerkin lead to wild nonphysical oscillations in the discrete solution.Therefore, stabilized Galerkin methods should be used. Most stabilized methods modify the standard Galerkin method: find some elementuh from a finite element space Vh such that
(2.5)aG(uh,vh):=ε(∇uh,∇vh)-(b·∇uh-cuh,vh)=(f,vh)∀vh∈Vh.
Instead, a stabilized method satisfies
(2.6)ah(uh,vh):=aG(uh,vh)+ast(uh,vh)=fh(vh),
where ast(·,·) represents a stabilization term. If the discrete bilinear form ah(·,·) is uniformly Vh-elliptic with respect to ε (or satisfies some inf-sup condition) in some norm |||·|||, stronger than the ε-weighted H1 norm, then one can typically prove on a shape-regular mesh for piecewise polynomials of degree k(2.7)|∥u-uh∥|≤C(ε1/2+h1/2)hk|u|k+1.
Although often C is independent of ε and h, the factor |u|k+1, in general, depends strongly on ε. Thus we call these methods “semi-robust.” Often, one can prove local versions of that estimate which show the numerical solution to behave fine away from layer regions.The streamline-diffusion FEM (SDFEM or SUPG-streamline upwind Petrov-Galerkin) is the most popular stabilization method. In the SDFEM we have(2.8)ast(w,v):=∑TδT(-εΔw-b∇w+cw,-b·∇v)T,
moreover
(2.9)fh(vh):=(f,vh)+∑TδT(f,-b·∇vh)T.
The streamline-diffusion parameters δT≥0 are user chosen. The SDFEM is consistent, that is,
(2.10)ah(u,vh)=fh(vh).
Therefore, one can use the orthogonality property ah(u-uh,vh)=0 for all vh∈Vh in the error estimate. Introducing the streamline-diffusion norm
(2.11)|∥v∥|SD2:=∥v∥ε2+∑TδT(b·∇v,b·∇v)T,
we have for bounded δT coercivity of the bilinear form in the SD-norm and the error estimate (2.7).There are many variants of SDFEM, see, for example, the surveys [4, 5]. Here we shall not discuss these variants and other methods closely related to SDFEM like residual free bubbles (RFB) or variational multiscale methods (VMS), see [1, 6–9].The SDFEM solution is not free of oscillations. See Section3.5.2 in [10] for the analysis of the related difference equations and its oscillation properties. Similar behavior can be expected for all other stabilization techniques we are going to discuss, although we are not aware of a rigorous analysis for these methods analogously to that of SDFEM.In 1D one can chooseδ for a problem with constant coefficients in such a way that the discrete solution is exact at the nodes. In two dimensions δT is usually determined by minimizing the factor ε+δT+δT-1hT2 under constraints on δT guaranteeing coercivity. In [10] the authors propose the smart choice
(2.12)δT={hT*2|bT|(1-1PT)ifPT:=|bT|hT*2ε>1,0ifPT≤1,
where |bT| is the Euclidean norm of the wind b at the element center and hT* the element length in the direction of the wind. But for anisotropic meshes the optimal choice of the SD-parameter is not clear. Later we shall see that near characteristic layers (2.12) is not optimal.The choice of theSD-parameter is still the topic of several recent papers. Sometimes the relation between SDFEM and RFB or VMS is used to define new formulae [11], sometimes just the oscillation properties near boundaries are carefully studied [12, 13] or even optimization methods are used a posteriori [14].Because of the observation that solutions computed with SDFEM often possess spurious oscillations in the vicinity of layers, a number of spurious oscillations at layers diminishing (SOLD) methods have been developed. A critical survey of these methods and a number of numerical studies can be found in [15–17].In [18] the authors compare the performance of several numerical schemes: an exponentially fitted finite volume scheme (see Section 3), SDFEM, SDFEM with SOLD, continuous interior penalty (CIP), a discontinuous Galerkin method (dG), and a FEM scheme with flux-correction-transport (FCT). The authors draw the conclusion that a favored method could not be identified and there is still the need to construct better methods than those which are currently available!We now comment on the methods just mentioned, and moreover on LPS (local projection stabilization), which was not discussed in [18].CIP is characterized by the symmetric stabilization term(2.13)ast(uh,vh):=∑EτE(b·[∇uh]E,b·[∇vh]E)E
which adds jumps of the streamline derivative on interior edges E to the Galerkin bilinear form [19]. In [20] it was pointed out that it is extremely useful to incorporate the Dirichlet boundary conditions weakly. Similarly to SDFEM, a nonlinear term of SOLD-type can be added to the CIP method as proposed in [21]. Tests with several SOLD-terms turned out that the results for CIP-SOLD were generally worse than the results for SDFEM-SOLD [18].Much more activity in the last years was directed on LPS, another (inconsistent) symmetric stabilization method. LPS preserves the stability properties of residual-based stabilizations, but LPS is symmetric. Therefore, the analysis of unsteady problems is simpler than for SDFEM (see the next subsection) and applied to optimization problems, discretization, and optimization commute.The most general version of LPS for convection-diffusion problems was described and analyzed in [22]. Let ℳh be a set consisting of a finite number of open subsets M of Ω with certain properties [22], πM a continuous linear projection operator which maps the space L2(M) onto some finite-dimensional subspace DM. Denote by bM a piecewise constant approximation of b and by κM:=id-πM the so-called fluctuation operator. Then, LPS is characterized by
(2.14)ast(uh,vh):=∑M‍τM(κM(bM·∇uh),κM(bM·∇vh))M
(see also [23–26]).Surprisingly, under mild conditions LPS is now stable in a norm containing additionally a streamline-derivative term, namely,(2.15)|∥v∥|LPSD2:=|∥v∥|LPS2+∑M‍τM∥b·∇v∥0,M2
with
(2.16)|∥v∥|LPS2:=∥v∥ε2+ast(v,v).
More precisely, ah(·,·) satisfies for 0≤τM≤ChM2(max(ε,hM))-1 the inf-sup-condition
(2.17)supvh∈Vhah(uh,vh)|∥vh∥|LPSD≥β|∥vh∥|LPSD.
In contrast to the traditional error analysis of LPS [25], there is no need to use special interpolation operators in the error estimation of the convection term. See also [27] for some insight in the stabilization mechanism of Galerkin discretizations of higher order.In [28] the authors compare the one-level and the two-level approach of LPS, see also [29], while in [30] even exponentials are used instead of bubbles to enrich the finite element space in a one-level version of LPS.Let us finally mention the RELP method [31], a residual local projection technique and the attempt to add also a non-linear SOLD term to LPS, see [32].Discontinuous Galerkin becomes more and more popular. Since 2008, for instance, four books on dG were published [33–36].In general, dG methods impose Dirichlet boundary conditions weakly. Surprisingly, only in the last years it was observed that in many cases one should prefer weakly imposed boundary conditions, see [20, 37, 38]. Recently [39] it was proved that the nonsymmetric Nitsche-type method is stable without penalty. Actually, hybrid dG attracts a lot of attention. The key idea of hybrid dG is to introduce additional degrees of freedom at interfaces. For hybrid primal dG see [40], for mixed hybrid [41], for convection-reaction-diffusion problems [42–44]. The so called DPG method combines the method of Bottasso et al. [45, 46] with the concept of “optimal” test functions, see [47–49]. There exist also several proposals to combine continuous and discontinuous Galerkin ideas, for instance, [50, 51].Remark that the numerical experiments in [18] are only performed for the SIP and NIP version of dG. Thus a final assessment of the many existing dG versions in comparison with other stabilization techniques remains open.Algebraic flux correction methods (FEM-FCT) are not mentioned in [1] but achieve in many cases good numerical results. But the methods are nonlinear, and based on the assessment in [18] quite inefficient for linear convection-diffusion problems. The methods start with a standard discretization, like the Galerkin finite element method. Then, the resulting matrices are changed in order to obtain, for instance, a positivity preserving, but still too diffusive scheme. In a second step the diffusion is locally removed using an appropriate antidiffusive term, see [52, 53].


### 2.2. Unsteady Problems

We consider the parabolic initial-boundary value problem(2.18a)ut+Lu=f(x,t)inQ=Ω×(0,T],(2.18b)u(x,0)=u0(x)forx∈Ω,(2.18c)u=0forx∈∂Ω,t∈(0,T].Again, standard Galerkin in space is unsuitable and one combines a stabilization technique in space with some discretization in time with nice stability properties, in the simplest case Euler implicit.For the symmetric stabilization techniques mentioned in Section2.1, there are semi-robust error estimates in norms typically used for parabolic problems in the literature. Reference [54] combines CIP or LPS with the θ-scheme or BDF2 in time, [55] combines LPS with dG in time while dG in space and time is analyzed in [56]. Reference [57] study a two stage IMEX RK scheme, where a symmetric stabilization is used and stabilization and convection are handled explicitly, but diffusion implicitly. For non-linear convection-diffusion problems, dG in space and several discretizations in time, see [58, 59].Streamline diffusion (and other non-symmetric methods) still offer open problems. For consistency reasons, the stabilization term reads(2.19)ast(w,v):=∑T‍δT(wt-εΔw-b∇w+cw,-b·∇v)T,
and the term which corresponds to testing the derivative in time with the streamline derivative causes trouble. If, for instance, implicit Euler is used for the time discretization with step size τ, the standard analysis in the convection-dominated regime leads to
(2.20)δ≈h≈τ.
But it seems not appropriate that the stabilization term disappears if the time step goes to zero. See the detailed discussion in [60]. For semi-discretization in space, optimal error estimates are proved in [60, 61], see also [62].


## 2.1. Steady-State Problems

We consider the two-dimensional convection-diffusion problem(2.1a)Lu:=-εΔu-b·∇u+cu=finΩ,(2.1b)u=0onΓ=∂Ω,whereε is a small positive parameter, b,c are smooth, and f∈L2(Ω). Assuming
(2.2)c+12divb≥α0>0,
the given problem admits a unique solution u∈H01(Ω), for convex Ω we have moreover u∈H2(Ω).Let us introduce theε-weighted H1 norm by
(2.3)∥v∥ε:=ε1/2|v|1+∥v∥0.
Then for the Galerkin finite element method with piecewise linear or bilinear elements one can prove (C denotes a generic constant that is independent of ε and of the mesh)
(2.4)∥u-uh∥ε≤Ch|u|2
on quite general triangulations. However, estimate (2.4) is of no worth: in general, |u|2 tends to infinity for ε→0 due to the presence of layers. The very weak stability properties of standard Galerkin lead to wild nonphysical oscillations in the discrete solution.Therefore, stabilized Galerkin methods should be used. Most stabilized methods modify the standard Galerkin method: find some elementuh from a finite element space Vh such that
(2.5)aG(uh,vh):=ε(∇uh,∇vh)-(b·∇uh-cuh,vh)=(f,vh)∀vh∈Vh.
Instead, a stabilized method satisfies
(2.6)ah(uh,vh):=aG(uh,vh)+ast(uh,vh)=fh(vh),
where ast(·,·) represents a stabilization term. If the discrete bilinear form ah(·,·) is uniformly Vh-elliptic with respect to ε (or satisfies some inf-sup condition) in some norm |||·|||, stronger than the ε-weighted H1 norm, then one can typically prove on a shape-regular mesh for piecewise polynomials of degree k(2.7)|∥u-uh∥|≤C(ε1/2+h1/2)hk|u|k+1.
Although often C is independent of ε and h, the factor |u|k+1, in general, depends strongly on ε. Thus we call these methods “semi-robust.” Often, one can prove local versions of that estimate which show the numerical solution to behave fine away from layer regions.The streamline-diffusion FEM (SDFEM or SUPG-streamline upwind Petrov-Galerkin) is the most popular stabilization method. In the SDFEM we have(2.8)ast(w,v):=∑TδT(-εΔw-b∇w+cw,-b·∇v)T,
moreover
(2.9)fh(vh):=(f,vh)+∑TδT(f,-b·∇vh)T.
The streamline-diffusion parameters δT≥0 are user chosen. The SDFEM is consistent, that is,
(2.10)ah(u,vh)=fh(vh).
Therefore, one can use the orthogonality property ah(u-uh,vh)=0 for all vh∈Vh in the error estimate. Introducing the streamline-diffusion norm
(2.11)|∥v∥|SD2:=∥v∥ε2+∑TδT(b·∇v,b·∇v)T,
we have for bounded δT coercivity of the bilinear form in the SD-norm and the error estimate (2.7).There are many variants of SDFEM, see, for example, the surveys [4, 5]. Here we shall not discuss these variants and other methods closely related to SDFEM like residual free bubbles (RFB) or variational multiscale methods (VMS), see [1, 6–9].The SDFEM solution is not free of oscillations. See Section3.5.2 in [10] for the analysis of the related difference equations and its oscillation properties. Similar behavior can be expected for all other stabilization techniques we are going to discuss, although we are not aware of a rigorous analysis for these methods analogously to that of SDFEM.In 1D one can chooseδ for a problem with constant coefficients in such a way that the discrete solution is exact at the nodes. In two dimensions δT is usually determined by minimizing the factor ε+δT+δT-1hT2 under constraints on δT guaranteeing coercivity. In [10] the authors propose the smart choice
(2.12)δT={hT*2|bT|(1-1PT)ifPT:=|bT|hT*2ε>1,0ifPT≤1,
where |bT| is the Euclidean norm of the wind b at the element center and hT* the element length in the direction of the wind. But for anisotropic meshes the optimal choice of the SD-parameter is not clear. Later we shall see that near characteristic layers (2.12) is not optimal.The choice of theSD-parameter is still the topic of several recent papers. Sometimes the relation between SDFEM and RFB or VMS is used to define new formulae [11], sometimes just the oscillation properties near boundaries are carefully studied [12, 13] or even optimization methods are used a posteriori [14].Because of the observation that solutions computed with SDFEM often possess spurious oscillations in the vicinity of layers, a number of spurious oscillations at layers diminishing (SOLD) methods have been developed. A critical survey of these methods and a number of numerical studies can be found in [15–17].In [18] the authors compare the performance of several numerical schemes: an exponentially fitted finite volume scheme (see Section 3), SDFEM, SDFEM with SOLD, continuous interior penalty (CIP), a discontinuous Galerkin method (dG), and a FEM scheme with flux-correction-transport (FCT). The authors draw the conclusion that a favored method could not be identified and there is still the need to construct better methods than those which are currently available!We now comment on the methods just mentioned, and moreover on LPS (local projection stabilization), which was not discussed in [18].CIP is characterized by the symmetric stabilization term(2.13)ast(uh,vh):=∑EτE(b·[∇uh]E,b·[∇vh]E)E
which adds jumps of the streamline derivative on interior edges E to the Galerkin bilinear form [19]. In [20] it was pointed out that it is extremely useful to incorporate the Dirichlet boundary conditions weakly. Similarly to SDFEM, a nonlinear term of SOLD-type can be added to the CIP method as proposed in [21]. Tests with several SOLD-terms turned out that the results for CIP-SOLD were generally worse than the results for SDFEM-SOLD [18].Much more activity in the last years was directed on LPS, another (inconsistent) symmetric stabilization method. LPS preserves the stability properties of residual-based stabilizations, but LPS is symmetric. Therefore, the analysis of unsteady problems is simpler than for SDFEM (see the next subsection) and applied to optimization problems, discretization, and optimization commute.The most general version of LPS for convection-diffusion problems was described and analyzed in [22]. Let ℳh be a set consisting of a finite number of open subsets M of Ω with certain properties [22], πM a continuous linear projection operator which maps the space L2(M) onto some finite-dimensional subspace DM. Denote by bM a piecewise constant approximation of b and by κM:=id-πM the so-called fluctuation operator. Then, LPS is characterized by
(2.14)ast(uh,vh):=∑M‍τM(κM(bM·∇uh),κM(bM·∇vh))M
(see also [23–26]).Surprisingly, under mild conditions LPS is now stable in a norm containing additionally a streamline-derivative term, namely,(2.15)|∥v∥|LPSD2:=|∥v∥|LPS2+∑M‍τM∥b·∇v∥0,M2
with
(2.16)|∥v∥|LPS2:=∥v∥ε2+ast(v,v).
More precisely, ah(·,·) satisfies for 0≤τM≤ChM2(max(ε,hM))-1 the inf-sup-condition
(2.17)supvh∈Vhah(uh,vh)|∥vh∥|LPSD≥β|∥vh∥|LPSD.
In contrast to the traditional error analysis of LPS [25], there is no need to use special interpolation operators in the error estimation of the convection term. See also [27] for some insight in the stabilization mechanism of Galerkin discretizations of higher order.In [28] the authors compare the one-level and the two-level approach of LPS, see also [29], while in [30] even exponentials are used instead of bubbles to enrich the finite element space in a one-level version of LPS.Let us finally mention the RELP method [31], a residual local projection technique and the attempt to add also a non-linear SOLD term to LPS, see [32].Discontinuous Galerkin becomes more and more popular. Since 2008, for instance, four books on dG were published [33–36].In general, dG methods impose Dirichlet boundary conditions weakly. Surprisingly, only in the last years it was observed that in many cases one should prefer weakly imposed boundary conditions, see [20, 37, 38]. Recently [39] it was proved that the nonsymmetric Nitsche-type method is stable without penalty. Actually, hybrid dG attracts a lot of attention. The key idea of hybrid dG is to introduce additional degrees of freedom at interfaces. For hybrid primal dG see [40], for mixed hybrid [41], for convection-reaction-diffusion problems [42–44]. The so called DPG method combines the method of Bottasso et al. [45, 46] with the concept of “optimal” test functions, see [47–49]. There exist also several proposals to combine continuous and discontinuous Galerkin ideas, for instance, [50, 51].Remark that the numerical experiments in [18] are only performed for the SIP and NIP version of dG. Thus a final assessment of the many existing dG versions in comparison with other stabilization techniques remains open.Algebraic flux correction methods (FEM-FCT) are not mentioned in [1] but achieve in many cases good numerical results. But the methods are nonlinear, and based on the assessment in [18] quite inefficient for linear convection-diffusion problems. The methods start with a standard discretization, like the Galerkin finite element method. Then, the resulting matrices are changed in order to obtain, for instance, a positivity preserving, but still too diffusive scheme. In a second step the diffusion is locally removed using an appropriate antidiffusive term, see [52, 53].


## 2.2. Unsteady Problems

We consider the parabolic initial-boundary value problem(2.18a)ut+Lu=f(x,t)inQ=Ω×(0,T],(2.18b)u(x,0)=u0(x)forx∈Ω,(2.18c)u=0forx∈∂Ω,t∈(0,T].Again, standard Galerkin in space is unsuitable and one combines a stabilization technique in space with some discretization in time with nice stability properties, in the simplest case Euler implicit.For the symmetric stabilization techniques mentioned in Section2.1, there are semi-robust error estimates in norms typically used for parabolic problems in the literature. Reference [54] combines CIP or LPS with the θ-scheme or BDF2 in time, [55] combines LPS with dG in time while dG in space and time is analyzed in [56]. Reference [57] study a two stage IMEX RK scheme, where a symmetric stabilization is used and stabilization and convection are handled explicitly, but diffusion implicitly. For non-linear convection-diffusion problems, dG in space and several discretizations in time, see [58, 59].Streamline diffusion (and other non-symmetric methods) still offer open problems. For consistency reasons, the stabilization term reads(2.19)ast(w,v):=∑T‍δT(wt-εΔw-b∇w+cw,-b·∇v)T,
and the term which corresponds to testing the derivative in time with the streamline derivative causes trouble. If, for instance, implicit Euler is used for the time discretization with step size τ, the standard analysis in the convection-dominated regime leads to
(2.20)δ≈h≈τ.
But it seems not appropriate that the stabilization term disappears if the time step goes to zero. See the detailed discussion in [60]. For semi-discretization in space, optimal error estimates are proved in [60, 61], see also [62].


## 3. Uniformly Convergent Methods

### 3.1. Standard Meshes

Consider the convection-diffusion-reaction problem(3.1)-εΔu-b∇u+cu=finΩ∈R2
with
(3.2)u=0on∂Ω,
where c+(1/2)divb≥α0>0 and Ω is polygonal and convex. Let us discretize (3.1) with standard finite elements on a shape-regular mesh. We ask: is it possible to have uniform convergence with respect to ε in the ε-weighted H1 norm, that is,(3.3)∥u-uh∥ε≤Chrwithsomer>0,
and C independent of ε,h?Surprisingly, the answer is yes. For instance, ifb≡0 and f/c∈H1, one can prove [63] (see also [64] for the LDG least-squares method)
(3.4)∥u-uh∥ε≤Ch1/2.Remark 3.1.
Estimates analogously to (3.4) are also known for the fourth-order problem
(3.5)ε2Δ2u-Δu+cu=f,u=∂u∂n=0on∂Ω.
See [65] for the first result in that direction and [66] for a detailed discussion and further references.
A result like (3.4) is possible because the ε-weighted norm is not balanced with the layer structure, the typical layer term E of the reaction-diffusion problem (b≡0 in (3.1)) gives
(3.6)∥E∥ε=O(ε1/4).
We conjecture that for convection-diffusion problems with b≢0 uniform convergence on standard meshes is also possible if the layers are relatively weak. This can happen for b·n=0 on ∂Ω (∂Ω is characteristic) or for other boundary conditions generating weak layers.
Now let us discuss the case of strong layers and uniform convergence on standard meshes. What is new?
Reference [67] repeat the known fact that in a Petrov-Galerkin scheme it is favorable to use local Green's functions as test functions. But, unfortunately, in 2D those functions cannot be derived in a closed analytical form. Therefore, the authors present a semi-analytical approach.
Optimal test functions play also an important role in the papers [47, 68, 69]. Assuming for the related bilinear form
(3.7)|a(u,v)|≤U∥u∥U∥v∥V
and the inf-sup condition
(3.8)inf∥u∥U=1sup∥v∥V=1a(u,v)≥γ>0,
optimal test functions are constructed as follows: If Un=span(un), define vn=Tun by
(3.9)(Tun,v)V=a(un,v)∀v∈V.
It is remarkable that then the discrete inf-sup condition follows (and, consequently, stability) and, moreover
(3.10)∥u-uh∥E=infwh∈Un∥u-wh∥E,
where ∥u∥E:=sup∥v∥V=1a(u,v) (see also the discussion in [1], pages 88–90). In the next step the test functions are approximated; new is the idea to do that in combination with discontinuous Galerkin methods.
Another old idea is to construct difference schemes in such a way that the coefficients of the scheme satisfy certain exactness conditions, including exponentials. Han and Huang call that “tailored finite point method,” see [70]. The error analysis uses an assumption of the type (for the case of a characteristic layer of width 𝒪(ε1/2))
(3.11)ε1/2h2≤hγ0,
such that the layer term is small in the interior mesh points.
For exponentials on triangles (see Remark 3.104 in [1]), there is still no satisfactorily uniform convergence result in 2D, although it is repeatedly the subject of actual research [71, 72]. Let us finally mention [73] where similar to papers for semiconductor problems [74, 75] an exponential transformation is used to generate a symmetric problem which is approximated by means of a suitable discontinuous Galerkin scheme.


### 3.2. Layer-Adapted Meshes

#### 3.2.1. Solution Decomposition and Meshes Used

In the simplest case of a one-dimensional problem with an exponential boundary layer atx=0 one can estimate
(3.12)|u(i)(x)|≤C(1+ε-iexp(-βxε))fori=0,1,…,q,
where q depends on the regularity of the data of the problem. Surprisingly, (3.12) is equivalent to the existence of a decomposition of the solution u=S+E into a smooth part S and a layer component E with
(3.13)|S(i)(x)|≤C,|E(i)(x)|≤Cε-iexp(-βxε)fori=0,1,…,q,LS=fLE=0.
We call such a decomposition S-decomposition because it was introduced by Shishkin to analyze upwind finite difference schemes.In two dimensions sufficient conditions for the existence of such a decomposition are known in certain special cases and for small values ofq only: for convection-diffusion problems with exponential layers and for some problems with characteristic layers; see [76–79].For the analysis and supercloseness analysis of linear and bilinear finite elements for a convection-diffusion problem with exponential layers, one assumes typically the existence of a solution decomposition(3.14a)u=S+E1+E2+E12
such that
(3.14b)|∂i+jS∂xi∂yj(x,y)|≤C,|∂i+jE1∂xi∂yj(x,y)|≤Cε-ie-β1x/ε,(3.14c)|∂i+jE2∂xi∂yj(x,y)|≤Cε-je-β2y/ε,|∂i+jE12∂xi∂yj(x,y)|≤Cε-(i+j)e-(β1x+β2y)/ε,for all(x,y)∈Ω- and 0≤i+j≤3.The validity of (3.14a), (3.14b), and (3.14c) requires additional compatibility conditions of the data at the corners which are sometimes unrealistic. Therefore, it is reasonable to search for analysis techniques of FEM on layer-adapted meshes that use the pointwise information of (3.14a), (3.14b) and (3.14c) for, say 0≤i+j≤2 only, and weaker L2 information for certain third-order derivatives. We are not going to discuss this issue in detail here and refer the reader to [1]. For the analysis of higher order finite element methods one requires a solution decomposition with analogous information on higher order derivatives.Based on the above solution decomposition we construct layer-adapted meshes. Let us assume that for a one-dimensional problem on the interval[0,1] a layer of type exp(-βx/ε) is located at x=0. As early as 1969, Bakhvalov [80] proposed a special mesh with mesh points xi near x=0 given by
(3.15)q(1-exp(-βxiσε))=ξi:=iN.
The parameter q∈(0,1) determines how many mesh points are used to resolve the layer, while σ>0 controls the spacing within the layer region. Outside the layer an equidistant mesh is used. To be precise, Bakhvalov's mesh is specified by xi=φ(i/N), i=0,1,…,N, where
(3.16)φ(ξ)={χ(ξ):=-σεβlnq-ξqforξ∈[0,τ],χ(τ)+ξ-τ1-τ(1-χ(τ))forξ∈[τ,1].
Here τ is a transition point between the fine and coarse submeshes. Originally Bakhvalov chose τ to ensure that the mesh generating function φ lay in C1(0,1) with φ(1)=1. However the explicit definition
(3.17)τ=γεβ|lnε|(e-βτ/ε=εγ)
is also possible and gives a mesh we shall refer to as a B-type mesh.From the numerical point of view, the choice of the transition point should reflect the smallness of the layer term with respect to the discretization error instead of the smallness with respect toε. Assume the formal order of the method to be σ. Then imposing
(3.18)exp(-βτε)=N-σ
yields the choice τ=(σε/β)lnN for the transition point. We call a mesh an S-type mesh if it is generated by
(3.19)φ(ξ)={σεβφ^(ξ)withφ^(0)=0,φ^(12)=lnNforξ∈[0,12],1-2(1-σεβlnN)(1-ξ)forξ∈[12,1].
In particular, when φ^(ξ)=2(lnN)ξ, the mesh generated is piecewise equidistant. This S-mesh was introduced by Shishkin in 1988. For surveys concerning layer adapted meshes, see [2, 81, 82].The analysis of certain difference methods for one-dimensional problems in [82] shows: if the pointwise error of a particular method on an S-mesh is proportional to (N-1lnN)σ, then on a B-type mesh (and on S-type meshes with certain optimality properties of the function φ^) the error is of order 𝒪(N-σ).For finite element methods the situation is different. So far, except for [83], there are no optimal error estimates for B-type meshes. If a transition point in the sense of Bakhvalov is used and a piecewise constant or locally uniform meshes, then the error weakly depends on ε (see [84]). In some papers, for instance in [85], the different impact of the choice of the transition point in the sense of Bakhvalov or Shishkin is discussed. On a piecewise constant or polynomial graded mesh (i.e., not on the original Bakhvalov mesh) the authors of [85] demonstrate numerical results which indicate the different error behavior on these two types of meshes.The simple structure ofS-type meshes allows error estimates for many methods as we will see in Section 3.2.2. For simplicity, we restrict ourselves often to S-meshes in this paper—in most cases a generalization of the results to S-type meshes is possible.For the two-dimensional problem (2.1a) and (2.1b) in Ω=(0,1)2 with exponential boundary layers let us define
(3.20)λx:=min{q,σεβ1lnN},λy=min{q,σεβ2lnN}.
This definition allows to include the non-singularly perturbed case. Divide the domain Ω as in Figure 1.Figure 1
Subregions ofΩ with exponential layers.The nodes of our rectangular mesh are obtained from the tensor product of a set ofNx points in x-direction and Ny points in y-direction. A one-dimensional Shishkin mesh is characterized by an equidistant mesh size h in [0,λx] and H in [λx,1], in the transition point λx the mesh switches from coarse to fine. For simplicity, let us assume
(3.21)β1=β2,λx=σεβ1lnN,Nx=Ny=N.
Thus
(3.22)h=2λxN,H=2(1-λx)N,xi=yi=ihfori=0,1,…,N2,xi=yi=(1-λx)+H(i-N2)fori=N2+1,…,N.
How to choose σ? For linear or bilinear elements one often chooses σ=2 (or σ=2.5 in supercloseness estimates), for elements based on the polynomial degree p we take σ=p+1.Remark 3.2.
In the caseb1>β1>0, b2≡0 the parabolic boundary layers at y=0 and y=1 are of width 𝒪(ε1/2|lnε|). Therefore, they require a different choice of transition point in y-direction: τy=𝒪(ε1/2lnN), while τx=𝒪(εlnN) remains unchanged; see Figure 2.Figure 2
Subregions ofΩ with exponential and parabolic layers.For more complicated domains the construction of layer-adapted meshes is of course much more involved; see, for instance, [86] for a description of the generation of such meshes. If one wants to use non-matching grids, it can be useful to apply Nitsche-mortaring, see [87].As also mentioned in [1], layer-adapted meshes can also be generated using a monitor function or defined recursively. We do not know any paper which analyses layer-adapted meshes defined via a monitor-function in 2D. Assume recursively (first for a 1D problem)
(3.23)xi:=xi-1+g(ε,H,xi-1).
We call a mesh G-type mesh (after Gartland), if g is of the type
(3.24)g(ε,H,xi-1)=εHexp(xi-1/(με)),
while a D-type mesh (after Duran) for convection-diffusion problems is characterized by
(3.25)x1=μHε;xi:=xi-1+μHxi-1.
Finite element methods for G-type meshes in 2D were analysed already 1997 in [88], for a more recent result for higher order equations but 1D see [89]. Convergence results on D-meshes we shall comment in the next section.


#### 3.2.2. Error Estimates and Supercloseness

For the Galerkin finite element method on Shishkin meshes we know since [90, 91] for convection-diffusion problems with exponential layers and linear or bilinear elements
(3.26)∥u-uN∥ε≤CN-1lnN,
improved for the more general S-type meshes in [92]. Bakhvalov-type meshes are more complicated to analyze, so far we have [93]
(3.27)∥u-uN∥ε≤CN-1Q(ε,N),
with Q≤ln10 for N≥10 and ε≥10-100. On D-meshes and problems with exponential layers it was proved [94]
(3.28)∥u-uN∥ε≤CN-1(ln(1ε))2.D-type meshes are locally uniform (or locally quasi-equidistant), thus it could be that some dependence of the error of ε is necessary. See the detailed discussion in [95].For problems with characteristic or parabolic boundary layers andS-type meshes, see [96]. Remark that optimal error estimates in the L∞ norm for the Galerkin method on layer-adapted meshes in 2D are unknown.In [97] the authors state that a layer-adapted mesh alone cannot stabilize a scheme in all cases, “a common believe of the community of Shishkin type grids.” But nobody in the community is this stupid; moreover, the authors do not understand the nature of the example presented. In a test case of a problem with characteristic layers and c=0 the Galerkin approximation on a Shishkin mesh oscillates strongly depending whether or not the number of mesh points used in the streamline direction is odd or not. The reason: in the continuous problem there is no control on the L2 norm, the continuous problem is already unstable. For an analogous example with c>0 this phenomenon disappears.It was first observed numerically in [98] that both for a Galerkin method and for streamline diffusion the convergence rates in L∞ for linear and bilinear elements on Ω∖Ω11 different significantly: the rates for bilinears are twice the rates for linears! This fact can be explained with superconvergence phenomena (supercloseness) for bilinears and is the reason for us to prefer bilinears in layer regions. Thus we have for bilinear elements on a Shishkin mesh (see [1]) the supercloseness result
(3.29)∥uI-uN∥ε≤C(N-1lnN)2.
Supercloseness allows optimal estimates in L2 and is also an important property concerning postprocessing. Often the so called Lin identities for bilinears, see for example [99], are used to prove supercloseness.In [100] Stynes and Tobiska analyzed the SDFEM for bilinears on an S-mesh. The SD-parameter is chosen by—assuming throughout that ε≤CN-1—
(3.30)δK={N-1ifK∈Ω11,0otherwise.
A detailed analysis shows, for instance, on Ω12 the stabilization parameter δ12 should satisfy δ12≤CεN-2. This value is much smaller than the natural diffusion parameter ε and therefore, switching of the stabilization by setting δ12=0 is reasonable.The final supercloseness result reads(3.31)∥uI-uN∥SD≤C(εN-3/2+N-2ln2N).
However, if bilinear elements are used in the layer region on Ω∖Ω11, but linear elements on the coarse-mesh region Ω11, then only the standard ingredients of the SDFEM analysis are available to estimate the error contribution on Ω11. Fortunately, on Ω11 the layer components are small. One obtains for the method on a hybrid mesh consisting of rectangles and triangles [83]
(3.32)∥uI-uN∥SD≤C(ε1/2N-1+N-3/2).
Remark that recently [101] pointwise error estimates were proved for the SDFEM with bilinears on a Shishkin mesh for problems with exponential layers, for linear elements see [78].For the problem with characteristic boundary layers it is more difficult to tune theSD-parameter. In the region Ω21 (see Figure 2) the recommendation (2.12) gives δ21=𝒪(N-1), but this choice is not appropriate [102]. For bilinears it was shown in [103] that δ21≤Cε-1/4N-2 should be satisfied.Based on (3.31) and the standard postprocessing approach of [104] one can construct a local postprocessing operator P such that PuN approximates u better than uN with respect to ∥·∥ε.Consider a family ofS-meshes 𝒯N where we require N/2 to be even. Then we can build a coarser mesh composed of disjoint macro rectangles M, each comprising four mesh rectangles from 𝒯N, where M belongs to only one of the four domains Ω11,Ω12,Ω21,Ω22. Associate with each macro rectangle M an interpolation operator PM:C(M¯)→Q2(M) defined by standard biquadratic interpolation. As usual, PM can be extended to a continuous global interpolation operator P:C(Ω¯)→WN, where WN is the space of piecewise quadratic elements.Then properties (7.2), (7.3) and (7.4) of [104]—approximation property, stability of P and consistency of P—give
(3.33)∥u-PuN∥ε≤C(εN-3/2+N-2ln2N).
The proof starts from the consistency of P:
(3.34)u-PuN=u-Pu+P(uI-uN).
Then the stability of P is applied in order to give
(3.35)∥u-PuN∥ε≤∥u-Pu∥ε+∥P(uI-uN)∥ε≤∥u-Pu∥ε+C∥uI-uN∥ε.
The terms on the right-hand side will be bounded using the approximation property of P and the supercloseness of the method.However, when linear elements on triangles are used inΩ11, supercloseness results are known only for special triangulations. Moreover, only for isosceles and Friedrichs-Keller triangulations postprocessing procedures are presented in [104].One can also prove supercloseness onD-meshes, see [105] for bilinear elements and convection-diffusion problems with exponential layers. Unfortunately, the final error estimate contains the factor (ln1/ε)5. In [106] the authors study reaction-diffusion problems and prove energy norm estimates for a graded mesh which is independent of the perturbation parameter ε. But for such a mesh one cannot have point-wise convergence uniformly in ε and we conjecture that this result is only possible because the norm is not balanced. For some nonconforming method, see [107].For higher order finite elements one can prove the expected error estimates on layer-adapted meshes if one has enough smoothness of the solution and a corresponding solution decomposition. But what about supercloseness forQp elements?For the Galerkin finite element method this seems to be an open question, see however the numerical experiments in [108]. For SDFEM it was proved in [109]
(3.36)∥πu-uN∥ε≤CN-(p+1/2).
Here instead of the standard interpolant uI the so called vertices-edge-interpolant πu is used. See also [110].Now, let us discuss other stabilization methods on layer-adapted meshes. Often stabilization is only used on the coarse part of the mesh.First, let us discuss the CIP method for a convection-diffusion problem with exponential layers. Using bilinears and Galerkin in the layer region, otherwise CIP stabilization one can prove for a “mixed” interpolantπ*u(3.37)∥uN-π*u∥CIP≤C(N-3/2+N-2ln2N);
see [111] for details of the analysis. The case of characteristic layers is discussed in [112].Next we consider the local-projection stabilization on the coarse mesh region of ourS-mesh for a problem with exponential layers. We consider the version of LPS-schemes characterized by enrichment of the original finite element space and the use of a single mesh rather than the version based on macro meshes. The stabilization term is given by
(3.38)aLPS(uh,vh):=∑T∈Ω11‍τN-1(b·∇uh-π(b·∇uh),b·∇vh-π(b·∇vh)).
For the LPS-scheme let Vh⊂H01(Ω) be the finite element space consisting of bilinears on Ω∖Ω11 and linear elements enriched by a single bubble per element on Ω11. Furthermore, let π project onto the space of piecewise constants on the triangulation of Ω11.For the analysis of the method a special interpolantjhw∈Vh of a given function w is used. On a triangle K it is defined by
(3.39)jhw(Pi)=w(Pi),∀threeverticesPiofK,(w-jhw,q)K=0,forarbitraryconstantsq.
Note, jhw used on Ω-11 and the standard bilinear interpolant wI used in Ω∖Ω11 match continuously. We call its composite π*.The well known results for the bilinear Galerkin method in the layer region and the technique of [25] give
(3.40)∥uN-π*u∥LPS≤C(N-3/2+N-2ln2N).
In [113] the author considers Qp-elements for P≥2 enriched by six additional functions (such that the element contains Pp+1) and local projection stabilization on Ω11 with τ=O(N-2). The resulting error estimates are however not optimal; it is better not to enrich the space on Ω/Ω11 as in [114]. Here it is proved for enriched Qp elements on the coarse mesh
(3.41)∥π*u-uN∥Lp≤CN-(p+1/2).
For problems with characteristic layers, see [115].It is also possible to combine the Galerkin finite element method with bilinears in the layer region with discontinuous Galerkin on the coarse mesh.For NIPG (non-symmetric interior penalty) in the associated dG-norm the supercloseness result(3.42)∥uN-πu∥dG≤C(ε1/2N-1+N-3/2)
was proved in  [116]. Here πu denotes the L2-projection onto P1disc(Ω11) and Q1disc(Ω11), respectively, on Ω11, but the standard bilinear interpolant on Ω∖Ω11.For the so-called local discontinuous Galerkin method LDG see [85, 117]. The authors present numerical studies which indicate superconvergence phenomena as well for the L2 norm of the solution as for the one side flux. A first convergence analysis in 2D for problems with exponential layers and bilinear elements was presented in [118]. Recently, a new type of error analysis for Qp elements on Shishkin meshes resulted in a surprising statement: uniform convergence of order N-(p+1/2) for u in a dG version of the ε-weighted energy norm [119].For edge stabilization and LPS-stabilization with bilinears one can use the same recovery technique as for SDFEM to get improved approximations with respect to∥·∥ε. This is because the finite element spaces used are continuous. However, for dGFEM recovery on Ω11 requires new ingredients. See, for instance, [120].It is possible to apply the so called two-scale finite element method also for singularly perturbed problems on Shishkin meshes. WritingN for the maximum number of mesh intervals in each coordinate direction of a tensor-product mesh for bilinear elements, the combination method simply adds or subtracts solutions that have been computed on N×N, N×N and N×N meshes. It is shown for a convection-diffusion problem with exponential layers [121], that the combination FEM yields (up to a factor lnN) the same order of accuracy as the Galerkin FEM on an N×N mesh. But it requires only 𝒪(N3/2) degrees of freedom compared with the 𝒪(N2) used by the Galerkin FEM. For a reaction-diffusion problem, see [122].Next we comment recent results concerningbalanced norms. Let us first consider reaction-diffusion problems.When linear or bilinear elements are used on a Shishkin mesh, one can prove under certain additional assumptions concerningf for the interpolation error of the Lagrange interpolant uI∈VN on Shishkin meshes
(3.43)∥u-uI∥ε≤C(ε1/4N-1lnN+N-2)
(see [91] or [1]). For a reaction-diffusion problem it follows immediately that the error u-uN also satisfies such an estimate.However, the typical boundary layer functionexp(-x/ε1/2) measured in the norm ∥·∥ε is of order 𝒪(ε1/4). That means: for reaction-diffusion problems, the standard energy norm is not balanced. Consequently, error estimates in this norm are less valuable as for convection diffusion equations where the layers are of the structure exp(-x/ε). Wherefore we first ask the question:Is it possible to prove error estimates in the balanced norm(3.44)∥v∥b:=ε1/4|v|1+∥v∥0?A proof is not trivial because the related bilinear form is not coercive with respect to the balanced norm. New ingredients are needed to prove(3.45)∥u-uN∥b≤C(N-1(lnN)3/2+N-2)
for the Galerkin FEM, see [123]. In [124] the authors introduce a mixed finite element method and prove also uniform convergence in a new balanced norm.For convection-diffusion problems with exponential layers only the energy norm is fine and balanced. If, however, layers of different structure exist an error estimate in a balanced norm is stronger than an estimate in the standard energy norm. Consider, for instance, the boundary value problem(3.46a)-εΔu-bux+cu=finΩ=(0,1)2,(3.46b)u=0on∂Ωwith exponential and characteristic layers. Here a balanced norm is given by(3.47)∥·∥b:=ε∥vx∥02+ε1/2∥vy∥02+γ∥v∥02.
Again the error analysis in that norm is not easy but could be realized for a version of the SDFEM on the corresponding Shishkin mesh, see [125].There exists also an attempt to simulate Shishkin mesh structures and in this way implicitly to stabilize the FEM even in domains with nontrivial geometries [126].In general, in the analysis of discretization methods on layer-adapted meshes relatively strong assumptions concerning the smoothness of the solution are posed. Concerning finite difference methods and reaction-diffusion problems, exceptions are [127, 128]. In [127] the domain is L-shaped, such that one has only u∈C2/3(Ω-). In [128] a Dirichlet and a Neumann condition meet in the corner of a square, thus u∈Cμ(Ω-) with μ<1. For finite element methods and reaction-diffusion problems on polygonal domains, see [86].The first paper on finite elements on layer-adapted meshes for convection-diffusion problems with a corner singularity is [129]. The authors study a problem with characteristic layers(3.48a)-εΔu-bux+cu=finΩ=(0,1)2,(3.48b)u=0on∂Ω,but without compatibility conditions in the corners, such that in each corner one has a singularity of the typer2lnr in the solution. For the Galerkin FEM and SDFEM on a classical Shishkin mesh uniform first order convergence is proved. It remains open whether or not supercloseness holds. In [130] even a problem with a discontinuity in the boundary condition in an inflow corner is studied numerically.For a convection-diffusion problem with exponential boundary layers in a domainΩω, a square sheared into a parallelogram with interior angle ω∈[π/2,π), in [131] the authors study the supercloseness property of the Galerkin-FEM and SDFEM on Shishkin meshes. Because u∉H3, new ingredients are required to prove supercloseness. The final supercloseness estimate takes the form
(3.49)∥uI-uN∥ε≤C(ε1/2(N-1lnN)π/ω-δ+N-2(lnN)2)
with some δ>0. See also [132]. For point-wise error estimates of finite difference methods and non-smooth solutions of convection-diffusion problems, see [133, 134].For unsteady problems of the form(3.50a)ut+Lu=f(x,t)inQ=Ω×(0,T],(3.50b)u(x,0)=u0(x)forx∈Ω,(3.50c)u=0forx∈∂Ω,t∈(0,T],there exist only a few papers from the last years combining some discretization in time with a discretization in space using a layer-adapted mesh.In [135] the authors use implicit Euler in time and a second order difference method on a Shishkin mesh and prove for a problem 1D in space
(3.51)maxi,j|u(xi,tj)-Uij|≤C(τ+(N-1lnN)2).
In [136] the given problem was discretized with the θ-scheme in time and the Galerkin FEM on a Shishkin mesh. In a norm adapted to the nature of the parabolic problem it was proved
(3.52)∥u-U∥≤C(τ+N-1lnN).
Using BDF2 with respect to time the order in τ was improved in [137]. In a recent paper dG in time of arbitrary order is combined with some finite element method in space, see [138].


## 3.1. Standard Meshes

Consider the convection-diffusion-reaction problem(3.1)-εΔu-b∇u+cu=finΩ∈R2
with
(3.2)u=0on∂Ω,
where c+(1/2)divb≥α0>0 and Ω is polygonal and convex. Let us discretize (3.1) with standard finite elements on a shape-regular mesh. We ask: is it possible to have uniform convergence with respect to ε in the ε-weighted H1 norm, that is,(3.3)∥u-uh∥ε≤Chrwithsomer>0,
and C independent of ε,h?Surprisingly, the answer is yes. For instance, ifb≡0 and f/c∈H1, one can prove [63] (see also [64] for the LDG least-squares method)
(3.4)∥u-uh∥ε≤Ch1/2.Remark 3.1.
Estimates analogously to (3.4) are also known for the fourth-order problem
(3.5)ε2Δ2u-Δu+cu=f,u=∂u∂n=0on∂Ω.
See [65] for the first result in that direction and [66] for a detailed discussion and further references.
A result like (3.4) is possible because the ε-weighted norm is not balanced with the layer structure, the typical layer term E of the reaction-diffusion problem (b≡0 in (3.1)) gives
(3.6)∥E∥ε=O(ε1/4).
We conjecture that for convection-diffusion problems with b≢0 uniform convergence on standard meshes is also possible if the layers are relatively weak. This can happen for b·n=0 on ∂Ω (∂Ω is characteristic) or for other boundary conditions generating weak layers.
Now let us discuss the case of strong layers and uniform convergence on standard meshes. What is new?
Reference [67] repeat the known fact that in a Petrov-Galerkin scheme it is favorable to use local Green's functions as test functions. But, unfortunately, in 2D those functions cannot be derived in a closed analytical form. Therefore, the authors present a semi-analytical approach.
Optimal test functions play also an important role in the papers [47, 68, 69]. Assuming for the related bilinear form
(3.7)|a(u,v)|≤U∥u∥U∥v∥V
and the inf-sup condition
(3.8)inf∥u∥U=1sup∥v∥V=1a(u,v)≥γ>0,
optimal test functions are constructed as follows: If Un=span(un), define vn=Tun by
(3.9)(Tun,v)V=a(un,v)∀v∈V.
It is remarkable that then the discrete inf-sup condition follows (and, consequently, stability) and, moreover
(3.10)∥u-uh∥E=infwh∈Un∥u-wh∥E,
where ∥u∥E:=sup∥v∥V=1a(u,v) (see also the discussion in [1], pages 88–90). In the next step the test functions are approximated; new is the idea to do that in combination with discontinuous Galerkin methods.
Another old idea is to construct difference schemes in such a way that the coefficients of the scheme satisfy certain exactness conditions, including exponentials. Han and Huang call that “tailored finite point method,” see [70]. The error analysis uses an assumption of the type (for the case of a characteristic layer of width 𝒪(ε1/2))
(3.11)ε1/2h2≤hγ0,
such that the layer term is small in the interior mesh points.
For exponentials on triangles (see Remark 3.104 in [1]), there is still no satisfactorily uniform convergence result in 2D, although it is repeatedly the subject of actual research [71, 72]. Let us finally mention [73] where similar to papers for semiconductor problems [74, 75] an exponential transformation is used to generate a symmetric problem which is approximated by means of a suitable discontinuous Galerkin scheme.


## 3.2. Layer-Adapted Meshes

### 3.2.1. Solution Decomposition and Meshes Used

In the simplest case of a one-dimensional problem with an exponential boundary layer atx=0 one can estimate
(3.12)|u(i)(x)|≤C(1+ε-iexp(-βxε))fori=0,1,…,q,
where q depends on the regularity of the data of the problem. Surprisingly, (3.12) is equivalent to the existence of a decomposition of the solution u=S+E into a smooth part S and a layer component E with
(3.13)|S(i)(x)|≤C,|E(i)(x)|≤Cε-iexp(-βxε)fori=0,1,…,q,LS=fLE=0.
We call such a decomposition S-decomposition because it was introduced by Shishkin to analyze upwind finite difference schemes.In two dimensions sufficient conditions for the existence of such a decomposition are known in certain special cases and for small values ofq only: for convection-diffusion problems with exponential layers and for some problems with characteristic layers; see [76–79].For the analysis and supercloseness analysis of linear and bilinear finite elements for a convection-diffusion problem with exponential layers, one assumes typically the existence of a solution decomposition(3.14a)u=S+E1+E2+E12
such that
(3.14b)|∂i+jS∂xi∂yj(x,y)|≤C,|∂i+jE1∂xi∂yj(x,y)|≤Cε-ie-β1x/ε,(3.14c)|∂i+jE2∂xi∂yj(x,y)|≤Cε-je-β2y/ε,|∂i+jE12∂xi∂yj(x,y)|≤Cε-(i+j)e-(β1x+β2y)/ε,for all(x,y)∈Ω- and 0≤i+j≤3.The validity of (3.14a), (3.14b), and (3.14c) requires additional compatibility conditions of the data at the corners which are sometimes unrealistic. Therefore, it is reasonable to search for analysis techniques of FEM on layer-adapted meshes that use the pointwise information of (3.14a), (3.14b) and (3.14c) for, say 0≤i+j≤2 only, and weaker L2 information for certain third-order derivatives. We are not going to discuss this issue in detail here and refer the reader to [1]. For the analysis of higher order finite element methods one requires a solution decomposition with analogous information on higher order derivatives.Based on the above solution decomposition we construct layer-adapted meshes. Let us assume that for a one-dimensional problem on the interval[0,1] a layer of type exp(-βx/ε) is located at x=0. As early as 1969, Bakhvalov [80] proposed a special mesh with mesh points xi near x=0 given by
(3.15)q(1-exp(-βxiσε))=ξi:=iN.
The parameter q∈(0,1) determines how many mesh points are used to resolve the layer, while σ>0 controls the spacing within the layer region. Outside the layer an equidistant mesh is used. To be precise, Bakhvalov's mesh is specified by xi=φ(i/N), i=0,1,…,N, where
(3.16)φ(ξ)={χ(ξ):=-σεβlnq-ξqforξ∈[0,τ],χ(τ)+ξ-τ1-τ(1-χ(τ))forξ∈[τ,1].
Here τ is a transition point between the fine and coarse submeshes. Originally Bakhvalov chose τ to ensure that the mesh generating function φ lay in C1(0,1) with φ(1)=1. However the explicit definition
(3.17)τ=γεβ|lnε|(e-βτ/ε=εγ)
is also possible and gives a mesh we shall refer to as a B-type mesh.From the numerical point of view, the choice of the transition point should reflect the smallness of the layer term with respect to the discretization error instead of the smallness with respect toε. Assume the formal order of the method to be σ. Then imposing
(3.18)exp(-βτε)=N-σ
yields the choice τ=(σε/β)lnN for the transition point. We call a mesh an S-type mesh if it is generated by
(3.19)φ(ξ)={σεβφ^(ξ)withφ^(0)=0,φ^(12)=lnNforξ∈[0,12],1-2(1-σεβlnN)(1-ξ)forξ∈[12,1].
In particular, when φ^(ξ)=2(lnN)ξ, the mesh generated is piecewise equidistant. This S-mesh was introduced by Shishkin in 1988. For surveys concerning layer adapted meshes, see [2, 81, 82].The analysis of certain difference methods for one-dimensional problems in [82] shows: if the pointwise error of a particular method on an S-mesh is proportional to (N-1lnN)σ, then on a B-type mesh (and on S-type meshes with certain optimality properties of the function φ^) the error is of order 𝒪(N-σ).For finite element methods the situation is different. So far, except for [83], there are no optimal error estimates for B-type meshes. If a transition point in the sense of Bakhvalov is used and a piecewise constant or locally uniform meshes, then the error weakly depends on ε (see [84]). In some papers, for instance in [85], the different impact of the choice of the transition point in the sense of Bakhvalov or Shishkin is discussed. On a piecewise constant or polynomial graded mesh (i.e., not on the original Bakhvalov mesh) the authors of [85] demonstrate numerical results which indicate the different error behavior on these two types of meshes.The simple structure ofS-type meshes allows error estimates for many methods as we will see in Section 3.2.2. For simplicity, we restrict ourselves often to S-meshes in this paper—in most cases a generalization of the results to S-type meshes is possible.For the two-dimensional problem (2.1a) and (2.1b) in Ω=(0,1)2 with exponential boundary layers let us define
(3.20)λx:=min{q,σεβ1lnN},λy=min{q,σεβ2lnN}.
This definition allows to include the non-singularly perturbed case. Divide the domain Ω as in Figure 1.Figure 1
Subregions ofΩ with exponential layers.The nodes of our rectangular mesh are obtained from the tensor product of a set ofNx points in x-direction and Ny points in y-direction. A one-dimensional Shishkin mesh is characterized by an equidistant mesh size h in [0,λx] and H in [λx,1], in the transition point λx the mesh switches from coarse to fine. For simplicity, let us assume
(3.21)β1=β2,λx=σεβ1lnN,Nx=Ny=N.
Thus
(3.22)h=2λxN,H=2(1-λx)N,xi=yi=ihfori=0,1,…,N2,xi=yi=(1-λx)+H(i-N2)fori=N2+1,…,N.
How to choose σ? For linear or bilinear elements one often chooses σ=2 (or σ=2.5 in supercloseness estimates), for elements based on the polynomial degree p we take σ=p+1.Remark 3.2.
In the caseb1>β1>0, b2≡0 the parabolic boundary layers at y=0 and y=1 are of width 𝒪(ε1/2|lnε|). Therefore, they require a different choice of transition point in y-direction: τy=𝒪(ε1/2lnN), while τx=𝒪(εlnN) remains unchanged; see Figure 2.Figure 2
Subregions ofΩ with exponential and parabolic layers.For more complicated domains the construction of layer-adapted meshes is of course much more involved; see, for instance, [86] for a description of the generation of such meshes. If one wants to use non-matching grids, it can be useful to apply Nitsche-mortaring, see [87].As also mentioned in [1], layer-adapted meshes can also be generated using a monitor function or defined recursively. We do not know any paper which analyses layer-adapted meshes defined via a monitor-function in 2D. Assume recursively (first for a 1D problem)
(3.23)xi:=xi-1+g(ε,H,xi-1).
We call a mesh G-type mesh (after Gartland), if g is of the type
(3.24)g(ε,H,xi-1)=εHexp(xi-1/(με)),
while a D-type mesh (after Duran) for convection-diffusion problems is characterized by
(3.25)x1=μHε;xi:=xi-1+μHxi-1.
Finite element methods for G-type meshes in 2D were analysed already 1997 in [88], for a more recent result for higher order equations but 1D see [89]. Convergence results on D-meshes we shall comment in the next section.


### 3.2.2. Error Estimates and Supercloseness

For the Galerkin finite element method on Shishkin meshes we know since [90, 91] for convection-diffusion problems with exponential layers and linear or bilinear elements
(3.26)∥u-uN∥ε≤CN-1lnN,
improved for the more general S-type meshes in [92]. Bakhvalov-type meshes are more complicated to analyze, so far we have [93]
(3.27)∥u-uN∥ε≤CN-1Q(ε,N),
with Q≤ln10 for N≥10 and ε≥10-100. On D-meshes and problems with exponential layers it was proved [94]
(3.28)∥u-uN∥ε≤CN-1(ln(1ε))2.D-type meshes are locally uniform (or locally quasi-equidistant), thus it could be that some dependence of the error of ε is necessary. See the detailed discussion in [95].For problems with characteristic or parabolic boundary layers andS-type meshes, see [96]. Remark that optimal error estimates in the L∞ norm for the Galerkin method on layer-adapted meshes in 2D are unknown.In [97] the authors state that a layer-adapted mesh alone cannot stabilize a scheme in all cases, “a common believe of the community of Shishkin type grids.” But nobody in the community is this stupid; moreover, the authors do not understand the nature of the example presented. In a test case of a problem with characteristic layers and c=0 the Galerkin approximation on a Shishkin mesh oscillates strongly depending whether or not the number of mesh points used in the streamline direction is odd or not. The reason: in the continuous problem there is no control on the L2 norm, the continuous problem is already unstable. For an analogous example with c>0 this phenomenon disappears.It was first observed numerically in [98] that both for a Galerkin method and for streamline diffusion the convergence rates in L∞ for linear and bilinear elements on Ω∖Ω11 different significantly: the rates for bilinears are twice the rates for linears! This fact can be explained with superconvergence phenomena (supercloseness) for bilinears and is the reason for us to prefer bilinears in layer regions. Thus we have for bilinear elements on a Shishkin mesh (see [1]) the supercloseness result
(3.29)∥uI-uN∥ε≤C(N-1lnN)2.
Supercloseness allows optimal estimates in L2 and is also an important property concerning postprocessing. Often the so called Lin identities for bilinears, see for example [99], are used to prove supercloseness.In [100] Stynes and Tobiska analyzed the SDFEM for bilinears on an S-mesh. The SD-parameter is chosen by—assuming throughout that ε≤CN-1—
(3.30)δK={N-1ifK∈Ω11,0otherwise.
A detailed analysis shows, for instance, on Ω12 the stabilization parameter δ12 should satisfy δ12≤CεN-2. This value is much smaller than the natural diffusion parameter ε and therefore, switching of the stabilization by setting δ12=0 is reasonable.The final supercloseness result reads(3.31)∥uI-uN∥SD≤C(εN-3/2+N-2ln2N).
However, if bilinear elements are used in the layer region on Ω∖Ω11, but linear elements on the coarse-mesh region Ω11, then only the standard ingredients of the SDFEM analysis are available to estimate the error contribution on Ω11. Fortunately, on Ω11 the layer components are small. One obtains for the method on a hybrid mesh consisting of rectangles and triangles [83]
(3.32)∥uI-uN∥SD≤C(ε1/2N-1+N-3/2).
Remark that recently [101] pointwise error estimates were proved for the SDFEM with bilinears on a Shishkin mesh for problems with exponential layers, for linear elements see [78].For the problem with characteristic boundary layers it is more difficult to tune theSD-parameter. In the region Ω21 (see Figure 2) the recommendation (2.12) gives δ21=𝒪(N-1), but this choice is not appropriate [102]. For bilinears it was shown in [103] that δ21≤Cε-1/4N-2 should be satisfied.Based on (3.31) and the standard postprocessing approach of [104] one can construct a local postprocessing operator P such that PuN approximates u better than uN with respect to ∥·∥ε.Consider a family ofS-meshes 𝒯N where we require N/2 to be even. Then we can build a coarser mesh composed of disjoint macro rectangles M, each comprising four mesh rectangles from 𝒯N, where M belongs to only one of the four domains Ω11,Ω12,Ω21,Ω22. Associate with each macro rectangle M an interpolation operator PM:C(M¯)→Q2(M) defined by standard biquadratic interpolation. As usual, PM can be extended to a continuous global interpolation operator P:C(Ω¯)→WN, where WN is the space of piecewise quadratic elements.Then properties (7.2), (7.3) and (7.4) of [104]—approximation property, stability of P and consistency of P—give
(3.33)∥u-PuN∥ε≤C(εN-3/2+N-2ln2N).
The proof starts from the consistency of P:
(3.34)u-PuN=u-Pu+P(uI-uN).
Then the stability of P is applied in order to give
(3.35)∥u-PuN∥ε≤∥u-Pu∥ε+∥P(uI-uN)∥ε≤∥u-Pu∥ε+C∥uI-uN∥ε.
The terms on the right-hand side will be bounded using the approximation property of P and the supercloseness of the method.However, when linear elements on triangles are used inΩ11, supercloseness results are known only for special triangulations. Moreover, only for isosceles and Friedrichs-Keller triangulations postprocessing procedures are presented in [104].One can also prove supercloseness onD-meshes, see [105] for bilinear elements and convection-diffusion problems with exponential layers. Unfortunately, the final error estimate contains the factor (ln1/ε)5. In [106] the authors study reaction-diffusion problems and prove energy norm estimates for a graded mesh which is independent of the perturbation parameter ε. But for such a mesh one cannot have point-wise convergence uniformly in ε and we conjecture that this result is only possible because the norm is not balanced. For some nonconforming method, see [107].For higher order finite elements one can prove the expected error estimates on layer-adapted meshes if one has enough smoothness of the solution and a corresponding solution decomposition. But what about supercloseness forQp elements?For the Galerkin finite element method this seems to be an open question, see however the numerical experiments in [108]. For SDFEM it was proved in [109]
(3.36)∥πu-uN∥ε≤CN-(p+1/2).
Here instead of the standard interpolant uI the so called vertices-edge-interpolant πu is used. See also [110].Now, let us discuss other stabilization methods on layer-adapted meshes. Often stabilization is only used on the coarse part of the mesh.First, let us discuss the CIP method for a convection-diffusion problem with exponential layers. Using bilinears and Galerkin in the layer region, otherwise CIP stabilization one can prove for a “mixed” interpolantπ*u(3.37)∥uN-π*u∥CIP≤C(N-3/2+N-2ln2N);
see [111] for details of the analysis. The case of characteristic layers is discussed in [112].Next we consider the local-projection stabilization on the coarse mesh region of ourS-mesh for a problem with exponential layers. We consider the version of LPS-schemes characterized by enrichment of the original finite element space and the use of a single mesh rather than the version based on macro meshes. The stabilization term is given by
(3.38)aLPS(uh,vh):=∑T∈Ω11‍τN-1(b·∇uh-π(b·∇uh),b·∇vh-π(b·∇vh)).
For the LPS-scheme let Vh⊂H01(Ω) be the finite element space consisting of bilinears on Ω∖Ω11 and linear elements enriched by a single bubble per element on Ω11. Furthermore, let π project onto the space of piecewise constants on the triangulation of Ω11.For the analysis of the method a special interpolantjhw∈Vh of a given function w is used. On a triangle K it is defined by
(3.39)jhw(Pi)=w(Pi),∀threeverticesPiofK,(w-jhw,q)K=0,forarbitraryconstantsq.
Note, jhw used on Ω-11 and the standard bilinear interpolant wI used in Ω∖Ω11 match continuously. We call its composite π*.The well known results for the bilinear Galerkin method in the layer region and the technique of [25] give
(3.40)∥uN-π*u∥LPS≤C(N-3/2+N-2ln2N).
In [113] the author considers Qp-elements for P≥2 enriched by six additional functions (such that the element contains Pp+1) and local projection stabilization on Ω11 with τ=O(N-2). The resulting error estimates are however not optimal; it is better not to enrich the space on Ω/Ω11 as in [114]. Here it is proved for enriched Qp elements on the coarse mesh
(3.41)∥π*u-uN∥Lp≤CN-(p+1/2).
For problems with characteristic layers, see [115].It is also possible to combine the Galerkin finite element method with bilinears in the layer region with discontinuous Galerkin on the coarse mesh.For NIPG (non-symmetric interior penalty) in the associated dG-norm the supercloseness result(3.42)∥uN-πu∥dG≤C(ε1/2N-1+N-3/2)
was proved in  [116]. Here πu denotes the L2-projection onto P1disc(Ω11) and Q1disc(Ω11), respectively, on Ω11, but the standard bilinear interpolant on Ω∖Ω11.For the so-called local discontinuous Galerkin method LDG see [85, 117]. The authors present numerical studies which indicate superconvergence phenomena as well for the L2 norm of the solution as for the one side flux. A first convergence analysis in 2D for problems with exponential layers and bilinear elements was presented in [118]. Recently, a new type of error analysis for Qp elements on Shishkin meshes resulted in a surprising statement: uniform convergence of order N-(p+1/2) for u in a dG version of the ε-weighted energy norm [119].For edge stabilization and LPS-stabilization with bilinears one can use the same recovery technique as for SDFEM to get improved approximations with respect to∥·∥ε. This is because the finite element spaces used are continuous. However, for dGFEM recovery on Ω11 requires new ingredients. See, for instance, [120].It is possible to apply the so called two-scale finite element method also for singularly perturbed problems on Shishkin meshes. WritingN for the maximum number of mesh intervals in each coordinate direction of a tensor-product mesh for bilinear elements, the combination method simply adds or subtracts solutions that have been computed on N×N, N×N and N×N meshes. It is shown for a convection-diffusion problem with exponential layers [121], that the combination FEM yields (up to a factor lnN) the same order of accuracy as the Galerkin FEM on an N×N mesh. But it requires only 𝒪(N3/2) degrees of freedom compared with the 𝒪(N2) used by the Galerkin FEM. For a reaction-diffusion problem, see [122].Next we comment recent results concerningbalanced norms. Let us first consider reaction-diffusion problems.When linear or bilinear elements are used on a Shishkin mesh, one can prove under certain additional assumptions concerningf for the interpolation error of the Lagrange interpolant uI∈VN on Shishkin meshes
(3.43)∥u-uI∥ε≤C(ε1/4N-1lnN+N-2)
(see [91] or [1]). For a reaction-diffusion problem it follows immediately that the error u-uN also satisfies such an estimate.However, the typical boundary layer functionexp(-x/ε1/2) measured in the norm ∥·∥ε is of order 𝒪(ε1/4). That means: for reaction-diffusion problems, the standard energy norm is not balanced. Consequently, error estimates in this norm are less valuable as for convection diffusion equations where the layers are of the structure exp(-x/ε). Wherefore we first ask the question:Is it possible to prove error estimates in the balanced norm(3.44)∥v∥b:=ε1/4|v|1+∥v∥0?A proof is not trivial because the related bilinear form is not coercive with respect to the balanced norm. New ingredients are needed to prove(3.45)∥u-uN∥b≤C(N-1(lnN)3/2+N-2)
for the Galerkin FEM, see [123]. In [124] the authors introduce a mixed finite element method and prove also uniform convergence in a new balanced norm.For convection-diffusion problems with exponential layers only the energy norm is fine and balanced. If, however, layers of different structure exist an error estimate in a balanced norm is stronger than an estimate in the standard energy norm. Consider, for instance, the boundary value problem(3.46a)-εΔu-bux+cu=finΩ=(0,1)2,(3.46b)u=0on∂Ωwith exponential and characteristic layers. Here a balanced norm is given by(3.47)∥·∥b:=ε∥vx∥02+ε1/2∥vy∥02+γ∥v∥02.
Again the error analysis in that norm is not easy but could be realized for a version of the SDFEM on the corresponding Shishkin mesh, see [125].There exists also an attempt to simulate Shishkin mesh structures and in this way implicitly to stabilize the FEM even in domains with nontrivial geometries [126].In general, in the analysis of discretization methods on layer-adapted meshes relatively strong assumptions concerning the smoothness of the solution are posed. Concerning finite difference methods and reaction-diffusion problems, exceptions are [127, 128]. In [127] the domain is L-shaped, such that one has only u∈C2/3(Ω-). In [128] a Dirichlet and a Neumann condition meet in the corner of a square, thus u∈Cμ(Ω-) with μ<1. For finite element methods and reaction-diffusion problems on polygonal domains, see [86].The first paper on finite elements on layer-adapted meshes for convection-diffusion problems with a corner singularity is [129]. The authors study a problem with characteristic layers(3.48a)-εΔu-bux+cu=finΩ=(0,1)2,(3.48b)u=0on∂Ω,but without compatibility conditions in the corners, such that in each corner one has a singularity of the typer2lnr in the solution. For the Galerkin FEM and SDFEM on a classical Shishkin mesh uniform first order convergence is proved. It remains open whether or not supercloseness holds. In [130] even a problem with a discontinuity in the boundary condition in an inflow corner is studied numerically.For a convection-diffusion problem with exponential boundary layers in a domainΩω, a square sheared into a parallelogram with interior angle ω∈[π/2,π), in [131] the authors study the supercloseness property of the Galerkin-FEM and SDFEM on Shishkin meshes. Because u∉H3, new ingredients are required to prove supercloseness. The final supercloseness estimate takes the form
(3.49)∥uI-uN∥ε≤C(ε1/2(N-1lnN)π/ω-δ+N-2(lnN)2)
with some δ>0. See also [132]. For point-wise error estimates of finite difference methods and non-smooth solutions of convection-diffusion problems, see [133, 134].For unsteady problems of the form(3.50a)ut+Lu=f(x,t)inQ=Ω×(0,T],(3.50b)u(x,0)=u0(x)forx∈Ω,(3.50c)u=0forx∈∂Ω,t∈(0,T],there exist only a few papers from the last years combining some discretization in time with a discretization in space using a layer-adapted mesh.In [135] the authors use implicit Euler in time and a second order difference method on a Shishkin mesh and prove for a problem 1D in space
(3.51)maxi,j|u(xi,tj)-Uij|≤C(τ+(N-1lnN)2).
In [136] the given problem was discretized with the θ-scheme in time and the Galerkin FEM on a Shishkin mesh. In a norm adapted to the nature of the parabolic problem it was proved
(3.52)∥u-U∥≤C(τ+N-1lnN).
Using BDF2 with respect to time the order in τ was improved in [137]. In a recent paper dG in time of arbitrary order is combined with some finite element method in space, see [138].


## 3.2.1. Solution Decomposition and Meshes Used

In the simplest case of a one-dimensional problem with an exponential boundary layer atx=0 one can estimate
(3.12)|u(i)(x)|≤C(1+ε-iexp(-βxε))fori=0,1,…,q,
where q depends on the regularity of the data of the problem. Surprisingly, (3.12) is equivalent to the existence of a decomposition of the solution u=S+E into a smooth part S and a layer component E with
(3.13)|S(i)(x)|≤C,|E(i)(x)|≤Cε-iexp(-βxε)fori=0,1,…,q,LS=fLE=0.
We call such a decomposition S-decomposition because it was introduced by Shishkin to analyze upwind finite difference schemes.In two dimensions sufficient conditions for the existence of such a decomposition are known in certain special cases and for small values ofq only: for convection-diffusion problems with exponential layers and for some problems with characteristic layers; see [76–79].For the analysis and supercloseness analysis of linear and bilinear finite elements for a convection-diffusion problem with exponential layers, one assumes typically the existence of a solution decomposition(3.14a)u=S+E1+E2+E12
such that
(3.14b)|∂i+jS∂xi∂yj(x,y)|≤C,|∂i+jE1∂xi∂yj(x,y)|≤Cε-ie-β1x/ε,(3.14c)|∂i+jE2∂xi∂yj(x,y)|≤Cε-je-β2y/ε,|∂i+jE12∂xi∂yj(x,y)|≤Cε-(i+j)e-(β1x+β2y)/ε,for all(x,y)∈Ω- and 0≤i+j≤3.The validity of (3.14a), (3.14b), and (3.14c) requires additional compatibility conditions of the data at the corners which are sometimes unrealistic. Therefore, it is reasonable to search for analysis techniques of FEM on layer-adapted meshes that use the pointwise information of (3.14a), (3.14b) and (3.14c) for, say 0≤i+j≤2 only, and weaker L2 information for certain third-order derivatives. We are not going to discuss this issue in detail here and refer the reader to [1]. For the analysis of higher order finite element methods one requires a solution decomposition with analogous information on higher order derivatives.Based on the above solution decomposition we construct layer-adapted meshes. Let us assume that for a one-dimensional problem on the interval[0,1] a layer of type exp(-βx/ε) is located at x=0. As early as 1969, Bakhvalov [80] proposed a special mesh with mesh points xi near x=0 given by
(3.15)q(1-exp(-βxiσε))=ξi:=iN.
The parameter q∈(0,1) determines how many mesh points are used to resolve the layer, while σ>0 controls the spacing within the layer region. Outside the layer an equidistant mesh is used. To be precise, Bakhvalov's mesh is specified by xi=φ(i/N), i=0,1,…,N, where
(3.16)φ(ξ)={χ(ξ):=-σεβlnq-ξqforξ∈[0,τ],χ(τ)+ξ-τ1-τ(1-χ(τ))forξ∈[τ,1].
Here τ is a transition point between the fine and coarse submeshes. Originally Bakhvalov chose τ to ensure that the mesh generating function φ lay in C1(0,1) with φ(1)=1. However the explicit definition
(3.17)τ=γεβ|lnε|(e-βτ/ε=εγ)
is also possible and gives a mesh we shall refer to as a B-type mesh.From the numerical point of view, the choice of the transition point should reflect the smallness of the layer term with respect to the discretization error instead of the smallness with respect toε. Assume the formal order of the method to be σ. Then imposing
(3.18)exp(-βτε)=N-σ
yields the choice τ=(σε/β)lnN for the transition point. We call a mesh an S-type mesh if it is generated by
(3.19)φ(ξ)={σεβφ^(ξ)withφ^(0)=0,φ^(12)=lnNforξ∈[0,12],1-2(1-σεβlnN)(1-ξ)forξ∈[12,1].
In particular, when φ^(ξ)=2(lnN)ξ, the mesh generated is piecewise equidistant. This S-mesh was introduced by Shishkin in 1988. For surveys concerning layer adapted meshes, see [2, 81, 82].The analysis of certain difference methods for one-dimensional problems in [82] shows: if the pointwise error of a particular method on an S-mesh is proportional to (N-1lnN)σ, then on a B-type mesh (and on S-type meshes with certain optimality properties of the function φ^) the error is of order 𝒪(N-σ).For finite element methods the situation is different. So far, except for [83], there are no optimal error estimates for B-type meshes. If a transition point in the sense of Bakhvalov is used and a piecewise constant or locally uniform meshes, then the error weakly depends on ε (see [84]). In some papers, for instance in [85], the different impact of the choice of the transition point in the sense of Bakhvalov or Shishkin is discussed. On a piecewise constant or polynomial graded mesh (i.e., not on the original Bakhvalov mesh) the authors of [85] demonstrate numerical results which indicate the different error behavior on these two types of meshes.The simple structure ofS-type meshes allows error estimates for many methods as we will see in Section 3.2.2. For simplicity, we restrict ourselves often to S-meshes in this paper—in most cases a generalization of the results to S-type meshes is possible.For the two-dimensional problem (2.1a) and (2.1b) in Ω=(0,1)2 with exponential boundary layers let us define
(3.20)λx:=min{q,σεβ1lnN},λy=min{q,σεβ2lnN}.
This definition allows to include the non-singularly perturbed case. Divide the domain Ω as in Figure 1.Figure 1
Subregions ofΩ with exponential layers.The nodes of our rectangular mesh are obtained from the tensor product of a set ofNx points in x-direction and Ny points in y-direction. A one-dimensional Shishkin mesh is characterized by an equidistant mesh size h in [0,λx] and H in [λx,1], in the transition point λx the mesh switches from coarse to fine. For simplicity, let us assume
(3.21)β1=β2,λx=σεβ1lnN,Nx=Ny=N.
Thus
(3.22)h=2λxN,H=2(1-λx)N,xi=yi=ihfori=0,1,…,N2,xi=yi=(1-λx)+H(i-N2)fori=N2+1,…,N.
How to choose σ? For linear or bilinear elements one often chooses σ=2 (or σ=2.5 in supercloseness estimates), for elements based on the polynomial degree p we take σ=p+1.Remark 3.2.
In the caseb1>β1>0, b2≡0 the parabolic boundary layers at y=0 and y=1 are of width 𝒪(ε1/2|lnε|). Therefore, they require a different choice of transition point in y-direction: τy=𝒪(ε1/2lnN), while τx=𝒪(εlnN) remains unchanged; see Figure 2.Figure 2
Subregions ofΩ with exponential and parabolic layers.For more complicated domains the construction of layer-adapted meshes is of course much more involved; see, for instance, [86] for a description of the generation of such meshes. If one wants to use non-matching grids, it can be useful to apply Nitsche-mortaring, see [87].As also mentioned in [1], layer-adapted meshes can also be generated using a monitor function or defined recursively. We do not know any paper which analyses layer-adapted meshes defined via a monitor-function in 2D. Assume recursively (first for a 1D problem)
(3.23)xi:=xi-1+g(ε,H,xi-1).
We call a mesh G-type mesh (after Gartland), if g is of the type
(3.24)g(ε,H,xi-1)=εHexp(xi-1/(με)),
while a D-type mesh (after Duran) for convection-diffusion problems is characterized by
(3.25)x1=μHε;xi:=xi-1+μHxi-1.
Finite element methods for G-type meshes in 2D were analysed already 1997 in [88], for a more recent result for higher order equations but 1D see [89]. Convergence results on D-meshes we shall comment in the next section.


## 3.2.2. Error Estimates and Supercloseness

For the Galerkin finite element method on Shishkin meshes we know since [90, 91] for convection-diffusion problems with exponential layers and linear or bilinear elements
(3.26)∥u-uN∥ε≤CN-1lnN,
improved for the more general S-type meshes in [92]. Bakhvalov-type meshes are more complicated to analyze, so far we have [93]
(3.27)∥u-uN∥ε≤CN-1Q(ε,N),
with Q≤ln10 for N≥10 and ε≥10-100. On D-meshes and problems with exponential layers it was proved [94]
(3.28)∥u-uN∥ε≤CN-1(ln(1ε))2.D-type meshes are locally uniform (or locally quasi-equidistant), thus it could be that some dependence of the error of ε is necessary. See the detailed discussion in [95].For problems with characteristic or parabolic boundary layers andS-type meshes, see [96]. Remark that optimal error estimates in the L∞ norm for the Galerkin method on layer-adapted meshes in 2D are unknown.In [97] the authors state that a layer-adapted mesh alone cannot stabilize a scheme in all cases, “a common believe of the community of Shishkin type grids.” But nobody in the community is this stupid; moreover, the authors do not understand the nature of the example presented. In a test case of a problem with characteristic layers and c=0 the Galerkin approximation on a Shishkin mesh oscillates strongly depending whether or not the number of mesh points used in the streamline direction is odd or not. The reason: in the continuous problem there is no control on the L2 norm, the continuous problem is already unstable. For an analogous example with c>0 this phenomenon disappears.It was first observed numerically in [98] that both for a Galerkin method and for streamline diffusion the convergence rates in L∞ for linear and bilinear elements on Ω∖Ω11 different significantly: the rates for bilinears are twice the rates for linears! This fact can be explained with superconvergence phenomena (supercloseness) for bilinears and is the reason for us to prefer bilinears in layer regions. Thus we have for bilinear elements on a Shishkin mesh (see [1]) the supercloseness result
(3.29)∥uI-uN∥ε≤C(N-1lnN)2.
Supercloseness allows optimal estimates in L2 and is also an important property concerning postprocessing. Often the so called Lin identities for bilinears, see for example [99], are used to prove supercloseness.In [100] Stynes and Tobiska analyzed the SDFEM for bilinears on an S-mesh. The SD-parameter is chosen by—assuming throughout that ε≤CN-1—
(3.30)δK={N-1ifK∈Ω11,0otherwise.
A detailed analysis shows, for instance, on Ω12 the stabilization parameter δ12 should satisfy δ12≤CεN-2. This value is much smaller than the natural diffusion parameter ε and therefore, switching of the stabilization by setting δ12=0 is reasonable.The final supercloseness result reads(3.31)∥uI-uN∥SD≤C(εN-3/2+N-2ln2N).
However, if bilinear elements are used in the layer region on Ω∖Ω11, but linear elements on the coarse-mesh region Ω11, then only the standard ingredients of the SDFEM analysis are available to estimate the error contribution on Ω11. Fortunately, on Ω11 the layer components are small. One obtains for the method on a hybrid mesh consisting of rectangles and triangles [83]
(3.32)∥uI-uN∥SD≤C(ε1/2N-1+N-3/2).
Remark that recently [101] pointwise error estimates were proved for the SDFEM with bilinears on a Shishkin mesh for problems with exponential layers, for linear elements see [78].For the problem with characteristic boundary layers it is more difficult to tune theSD-parameter. In the region Ω21 (see Figure 2) the recommendation (2.12) gives δ21=𝒪(N-1), but this choice is not appropriate [102]. For bilinears it was shown in [103] that δ21≤Cε-1/4N-2 should be satisfied.Based on (3.31) and the standard postprocessing approach of [104] one can construct a local postprocessing operator P such that PuN approximates u better than uN with respect to ∥·∥ε.Consider a family ofS-meshes 𝒯N where we require N/2 to be even. Then we can build a coarser mesh composed of disjoint macro rectangles M, each comprising four mesh rectangles from 𝒯N, where M belongs to only one of the four domains Ω11,Ω12,Ω21,Ω22. Associate with each macro rectangle M an interpolation operator PM:C(M¯)→Q2(M) defined by standard biquadratic interpolation. As usual, PM can be extended to a continuous global interpolation operator P:C(Ω¯)→WN, where WN is the space of piecewise quadratic elements.Then properties (7.2), (7.3) and (7.4) of [104]—approximation property, stability of P and consistency of P—give
(3.33)∥u-PuN∥ε≤C(εN-3/2+N-2ln2N).
The proof starts from the consistency of P:
(3.34)u-PuN=u-Pu+P(uI-uN).
Then the stability of P is applied in order to give
(3.35)∥u-PuN∥ε≤∥u-Pu∥ε+∥P(uI-uN)∥ε≤∥u-Pu∥ε+C∥uI-uN∥ε.
The terms on the right-hand side will be bounded using the approximation property of P and the supercloseness of the method.However, when linear elements on triangles are used inΩ11, supercloseness results are known only for special triangulations. Moreover, only for isosceles and Friedrichs-Keller triangulations postprocessing procedures are presented in [104].One can also prove supercloseness onD-meshes, see [105] for bilinear elements and convection-diffusion problems with exponential layers. Unfortunately, the final error estimate contains the factor (ln1/ε)5. In [106] the authors study reaction-diffusion problems and prove energy norm estimates for a graded mesh which is independent of the perturbation parameter ε. But for such a mesh one cannot have point-wise convergence uniformly in ε and we conjecture that this result is only possible because the norm is not balanced. For some nonconforming method, see [107].For higher order finite elements one can prove the expected error estimates on layer-adapted meshes if one has enough smoothness of the solution and a corresponding solution decomposition. But what about supercloseness forQp elements?For the Galerkin finite element method this seems to be an open question, see however the numerical experiments in [108]. For SDFEM it was proved in [109]
(3.36)∥πu-uN∥ε≤CN-(p+1/2).
Here instead of the standard interpolant uI the so called vertices-edge-interpolant πu is used. See also [110].Now, let us discuss other stabilization methods on layer-adapted meshes. Often stabilization is only used on the coarse part of the mesh.First, let us discuss the CIP method for a convection-diffusion problem with exponential layers. Using bilinears and Galerkin in the layer region, otherwise CIP stabilization one can prove for a “mixed” interpolantπ*u(3.37)∥uN-π*u∥CIP≤C(N-3/2+N-2ln2N);
see [111] for details of the analysis. The case of characteristic layers is discussed in [112].Next we consider the local-projection stabilization on the coarse mesh region of ourS-mesh for a problem with exponential layers. We consider the version of LPS-schemes characterized by enrichment of the original finite element space and the use of a single mesh rather than the version based on macro meshes. The stabilization term is given by
(3.38)aLPS(uh,vh):=∑T∈Ω11‍τN-1(b·∇uh-π(b·∇uh),b·∇vh-π(b·∇vh)).
For the LPS-scheme let Vh⊂H01(Ω) be the finite element space consisting of bilinears on Ω∖Ω11 and linear elements enriched by a single bubble per element on Ω11. Furthermore, let π project onto the space of piecewise constants on the triangulation of Ω11.For the analysis of the method a special interpolantjhw∈Vh of a given function w is used. On a triangle K it is defined by
(3.39)jhw(Pi)=w(Pi),∀threeverticesPiofK,(w-jhw,q)K=0,forarbitraryconstantsq.
Note, jhw used on Ω-11 and the standard bilinear interpolant wI used in Ω∖Ω11 match continuously. We call its composite π*.The well known results for the bilinear Galerkin method in the layer region and the technique of [25] give
(3.40)∥uN-π*u∥LPS≤C(N-3/2+N-2ln2N).
In [113] the author considers Qp-elements for P≥2 enriched by six additional functions (such that the element contains Pp+1) and local projection stabilization on Ω11 with τ=O(N-2). The resulting error estimates are however not optimal; it is better not to enrich the space on Ω/Ω11 as in [114]. Here it is proved for enriched Qp elements on the coarse mesh
(3.41)∥π*u-uN∥Lp≤CN-(p+1/2).
For problems with characteristic layers, see [115].It is also possible to combine the Galerkin finite element method with bilinears in the layer region with discontinuous Galerkin on the coarse mesh.For NIPG (non-symmetric interior penalty) in the associated dG-norm the supercloseness result(3.42)∥uN-πu∥dG≤C(ε1/2N-1+N-3/2)
was proved in  [116]. Here πu denotes the L2-projection onto P1disc(Ω11) and Q1disc(Ω11), respectively, on Ω11, but the standard bilinear interpolant on Ω∖Ω11.For the so-called local discontinuous Galerkin method LDG see [85, 117]. The authors present numerical studies which indicate superconvergence phenomena as well for the L2 norm of the solution as for the one side flux. A first convergence analysis in 2D for problems with exponential layers and bilinear elements was presented in [118]. Recently, a new type of error analysis for Qp elements on Shishkin meshes resulted in a surprising statement: uniform convergence of order N-(p+1/2) for u in a dG version of the ε-weighted energy norm [119].For edge stabilization and LPS-stabilization with bilinears one can use the same recovery technique as for SDFEM to get improved approximations with respect to∥·∥ε. This is because the finite element spaces used are continuous. However, for dGFEM recovery on Ω11 requires new ingredients. See, for instance, [120].It is possible to apply the so called two-scale finite element method also for singularly perturbed problems on Shishkin meshes. WritingN for the maximum number of mesh intervals in each coordinate direction of a tensor-product mesh for bilinear elements, the combination method simply adds or subtracts solutions that have been computed on N×N, N×N and N×N meshes. It is shown for a convection-diffusion problem with exponential layers [121], that the combination FEM yields (up to a factor lnN) the same order of accuracy as the Galerkin FEM on an N×N mesh. But it requires only 𝒪(N3/2) degrees of freedom compared with the 𝒪(N2) used by the Galerkin FEM. For a reaction-diffusion problem, see [122].Next we comment recent results concerningbalanced norms. Let us first consider reaction-diffusion problems.When linear or bilinear elements are used on a Shishkin mesh, one can prove under certain additional assumptions concerningf for the interpolation error of the Lagrange interpolant uI∈VN on Shishkin meshes
(3.43)∥u-uI∥ε≤C(ε1/4N-1lnN+N-2)
(see [91] or [1]). For a reaction-diffusion problem it follows immediately that the error u-uN also satisfies such an estimate.However, the typical boundary layer functionexp(-x/ε1/2) measured in the norm ∥·∥ε is of order 𝒪(ε1/4). That means: for reaction-diffusion problems, the standard energy norm is not balanced. Consequently, error estimates in this norm are less valuable as for convection diffusion equations where the layers are of the structure exp(-x/ε). Wherefore we first ask the question:Is it possible to prove error estimates in the balanced norm(3.44)∥v∥b:=ε1/4|v|1+∥v∥0?A proof is not trivial because the related bilinear form is not coercive with respect to the balanced norm. New ingredients are needed to prove(3.45)∥u-uN∥b≤C(N-1(lnN)3/2+N-2)
for the Galerkin FEM, see [123]. In [124] the authors introduce a mixed finite element method and prove also uniform convergence in a new balanced norm.For convection-diffusion problems with exponential layers only the energy norm is fine and balanced. If, however, layers of different structure exist an error estimate in a balanced norm is stronger than an estimate in the standard energy norm. Consider, for instance, the boundary value problem(3.46a)-εΔu-bux+cu=finΩ=(0,1)2,(3.46b)u=0on∂Ωwith exponential and characteristic layers. Here a balanced norm is given by(3.47)∥·∥b:=ε∥vx∥02+ε1/2∥vy∥02+γ∥v∥02.
Again the error analysis in that norm is not easy but could be realized for a version of the SDFEM on the corresponding Shishkin mesh, see [125].There exists also an attempt to simulate Shishkin mesh structures and in this way implicitly to stabilize the FEM even in domains with nontrivial geometries [126].In general, in the analysis of discretization methods on layer-adapted meshes relatively strong assumptions concerning the smoothness of the solution are posed. Concerning finite difference methods and reaction-diffusion problems, exceptions are [127, 128]. In [127] the domain is L-shaped, such that one has only u∈C2/3(Ω-). In [128] a Dirichlet and a Neumann condition meet in the corner of a square, thus u∈Cμ(Ω-) with μ<1. For finite element methods and reaction-diffusion problems on polygonal domains, see [86].The first paper on finite elements on layer-adapted meshes for convection-diffusion problems with a corner singularity is [129]. The authors study a problem with characteristic layers(3.48a)-εΔu-bux+cu=finΩ=(0,1)2,(3.48b)u=0on∂Ω,but without compatibility conditions in the corners, such that in each corner one has a singularity of the typer2lnr in the solution. For the Galerkin FEM and SDFEM on a classical Shishkin mesh uniform first order convergence is proved. It remains open whether or not supercloseness holds. In [130] even a problem with a discontinuity in the boundary condition in an inflow corner is studied numerically.For a convection-diffusion problem with exponential boundary layers in a domainΩω, a square sheared into a parallelogram with interior angle ω∈[π/2,π), in [131] the authors study the supercloseness property of the Galerkin-FEM and SDFEM on Shishkin meshes. Because u∉H3, new ingredients are required to prove supercloseness. The final supercloseness estimate takes the form
(3.49)∥uI-uN∥ε≤C(ε1/2(N-1lnN)π/ω-δ+N-2(lnN)2)
with some δ>0. See also [132]. For point-wise error estimates of finite difference methods and non-smooth solutions of convection-diffusion problems, see [133, 134].For unsteady problems of the form(3.50a)ut+Lu=f(x,t)inQ=Ω×(0,T],(3.50b)u(x,0)=u0(x)forx∈Ω,(3.50c)u=0forx∈∂Ω,t∈(0,T],there exist only a few papers from the last years combining some discretization in time with a discretization in space using a layer-adapted mesh.In [135] the authors use implicit Euler in time and a second order difference method on a Shishkin mesh and prove for a problem 1D in space
(3.51)maxi,j|u(xi,tj)-Uij|≤C(τ+(N-1lnN)2).
In [136] the given problem was discretized with the θ-scheme in time and the Galerkin FEM on a Shishkin mesh. In a norm adapted to the nature of the parabolic problem it was proved
(3.52)∥u-U∥≤C(τ+N-1lnN).
Using BDF2 with respect to time the order in τ was improved in [137]. In a recent paper dG in time of arbitrary order is combined with some finite element method in space, see [138].


## 4. Adaptive Methods

Adaptive methods refine the mesh (h-method) or change locally the polynomial degree (p-method) based on an a posteriori error estimatorη. η should be locally computable based on the actual numerical solution uh and input data. Moreover, preferably η should be equivalent to the error in some norm:
(4.1)dlη≤∥u-uh∥≤duη
(remark that in the DWR method one tries to control some functional instead of the norm). If for a singularly perturbed problem the constants dl,du are independent of ε, we call the estimatorrobust. If for convection-diffusion problems and some norms du is independent of ε but dl slightly depends on ε, we call the estimatorsemi-robust with respect to that norm.First we discuss error estimators with respect to an energy or some related norm. Based on [139] in the last years for convection- diffusion problems instead of using the norm ∥·∥ε robustness of several estimators with respect to the norm
(4.2)|∥v∥|2:=∥v∥ε2+supφ(b·∇v,φ)∥φ∥ε
was proved, see [140–145]. Unfortunately, this norm is not computable.Relatively new, not discussed in [1] and actually used in many papers, is the idea of flux reconstruction in H(div,Ω). Consider first, for simplicity a reaction-diffusion problem following [146]. The derivation of the error estimator starts from
(4.3)a(u-uh,v)=(f-cuh-∇·σh,v)-(ε∇uh+σh,∇v).
Now σh∈H(div,Ω) is chosen to approximate the numerical flux -ε∇uh. Additionally, σh shall satisfy
(4.4)(∇·σh+cuh,1)K=(f,1)K∀K∈𝒯h.
For details of the computation of the recovered flux σh see [43].The combination of the first term of (4.3) with (4.4) yields the residual part of the new estimator (with the weights mK from [139])
(4.5)ηK,res:=mK∥f-cuh-∇·σh∥K.
The second term of (4.3) yields after some manipulation (the direct application of Cauchy-Schwarz yields a non-robust estimator) a more complicated diffusive flux estimator ηDF. The resulting full estimator of [146] is robust and equivalent to Verfürth's residual estimator for reaction-diffusion problems.For convection-diffusion problems the derivation of estimators using flux reconstruction works similarly. Consider a convection-diffusion problem of the form(4.6a)-εΔu+∇·(bu)+cu=finΩ=(0,1)2,(4.6b)u=0on∂Ω.Nowσh has to approximate -ε∇uh+buh and should again satisfy (4.4). See [147] for details. For a similar approach concerning fully computable a posteriori error estimators for stabilized FEM (even in 3D) see [148].One can also use a discontinuous numerical approximationuh as typical for dG methods. Than an additional reconstruction of the potential uh by a function in H1(Ω) is introduced, see [149, 150]. For the separation of the error into an element residual error and flux errors with respect to goal-oriented a posteriori error estimates (DWR method), see [151]. Nonlinear unsteady problems are discussed in [152].There exist some attempts to derive pointwise a posteriori error estimates in the singularly-perturbed case. Using the Green's function of the continuous problem the starting point is, in general, the representation(4.7)(u-uh)(x)=a(u-uh,G)=(f,G)-a(uh,G)
or (in the strong form using distributions)
(4.8)(u-uh)(x)=∫‍L(u-uh)G.
In [153] the author studies in 1D higher order FEM based on (4.7) and gets estimates using L1 information on G. The estimator contains discrete derivatives of the numerical solution uh. Kopteva [154] studied the finite difference method in 2D based on (4.8).It is open to use the information on the Green's function provided in [155, 156] to derive pointwise a posteriori error estimates for finite elements in 2D. For time-dependent problems, but one-dimensional in space, see [157].An error estimator provides the informationwhere the mesh should be refined. But, in general, we do not have additional information on the directional behavior of the error in order to know how the mesh should be refined. An exception is [158], where the author tries to detect anisotropies, see also the recent [159]. There seems to be no new theoretical perception concerning metric-based algorithms for mesh generation in the singularly perturbed case.It is still mostly standard to derive error estimators based on the assumption that the mesh is shape-regular and locally uniform. For anisotropic meshes, there seems to be so far nothing better than the theory developed by Kunert [160] ten years ago. For some survey on anisotropic refinement methods in FEM, see [161].


## 5. Singularly Perturbed Systems

In 2009, Linss and Stynes [162] presented a survey on the numerical solution of singularly perturbed systems. In this Section we only comment on some recent results not contained in [162] which sparkle that survey.First we study systems of reaction-diffusion equations of the form(5.1)-E2u′′+Au=fin(0,1),u(0)=u(1)=0,
where E=diag(ε1,ε2,…,εl) with ε1≤ε2≤⋯≤εl. If the matrix A satisfies certain conditions, the asymptotic behavior of such a system is well understood. Assume that A has positive diagonal entries, moreover the matrix Γ defined by
(5.2)γii=1,γij=-∥aijaii∥∞fori≠j
satisfies Γ-1≥0. Then, in [163] the existence of a solution decomposition is proved. Other authors assume that A is an M-matrix or that A is point-wise positive definite. See [162, Theorem 2.2] for establishing a connection between positive definiteness and the property Γ-1≥0. In [164] a full asymptotic expansion is derived for positive definite A in the case of two equations, including information on analytic regularity.Systems of convection-diffusion problems are more delicate to handle. Consider first weakly coupled systems of the form(5.3)Lu:=-Eu′′-diag(b)u'+Au=f,u(0)=u(1)=0,
assuming |bi|≥βi>0andΓ~-1≥0with(5.4)γ~ii=1,γ~ij=-min(∥aijaii∥∞,∥aijbi∥∞)fori≠j.
Then it was shown in [165] for ν=0,1(5.5)|uk(ν)(x)|≤C{1+εk-νe-βk(1-x)/εifbk<0,1+εk-νe-βkx/εifbk>0.
When only first order derivatives are considered, there is no strong interaction between the layers of different components unlike the reaction-diffusion case.But, consider for example a set of two equations withb1>0 and b2<0 for ε1=ε2. Then the layer at x=1 in the first component generates a weak layer at x=1 in the second component; the situation at x=0 is analogous. Under certain conditions [166], one can prove the existence of the following solution decomposition for ν≤2:
(5.6)u1=S1+E10+E11,u2=S2+E20+E21
with
(5.7)∥S1(ν)∥0,∥S2(ν)∥0≤C,|E10(ν)(x)|≤Cε1-νe-α(x/ε),|E11(ν)(x)|≤Cε-νe-α((1-x)/ε),|E20(ν)(x)|≤Cε-νe-α(x/ε),|E21(ν)(x)|≤Cε1-νe-α((1-x)/ε).
Here α is some positive parameter. This observation is important for control problems governed by convection-diffusion equations:
(5.8)miny,qJ(y,q):=miny,q(12∥y-y0∥02+λ2∥q∥02)
subject to
(5.9)Ly:=-εy′′+by'+cy=f+qin(0,1),y(0)=y(1)=0.
For strongly coupled systems of convection-diffusion equations full layer-interaction takes place. Consider the system of two equations
(5.10)Lu:=-εu′′-Bu'+Au=f,u(0)=u(1)=0
assuming (V1)
B is symmetric.(V2)
A+1/2B'  is positive semidefinite.(V3)
The eigenvalues ofB satisfy |λ1,2|>α>0  for all x.If both eigenvalues ofB are positive, both solution components do have overlapping layers at x=0, the reduced solution solves an initial value problem [167]. But if the eigenvalues do have a different sign, both solution components do have layers at x=0andx=1; we have full layer interaction. It is remarkable that the reduced solution, in general, does not satisfy any of the given boundary conditions [168].Even more complicated are strongly coupled systems with several small parameters of the form(5.11)Lu:=-Eu′′-Bu'+Au=f,u(0)=u(1)=0.
Some a priori estimates are to find in [169], information on the layer structure in [164, 170].



---

*Source: 379547-2012-12-06.xml*