# ISOGEOMETRIC ANALYSIS FOR SINGULARLY PERTURBED PROBLEMS IN 1-D: ERROR ESTIMATES.

1. Introduction. We consider second order singularly perturbed problems (SPPs) in one-dimension, of reaction-convection-diffusion type, whose solution contains boundary layers; see, e.g., [13]. The approximation of the solution to SPPs has received a lot of attention in the last few decades, mainly using finite differences (FDs) and finite elements (FEs) on layer adapted meshes; see, e.g., [9]. Various formulations and results are available in the literature, both theoretical and computational (cf. [9, 13], and the references therein). One method that has not, to our knowledge, been applied to general SPPs is isogeometric analysis (IGA). Since the introduction of IGA by T. R. Hughes et al. [7], the method has been successfully applied to a large number of problem classes. Even though much attention has been given to convection-dominated problems [3], the method has not been applied, as far as we know, to a typical singularly perturbed problem, such as (2.1), (2.2) ahead.

Our goal in this article is to provide the theoretical justification for the application of IGA to SPPs and, in particular, for the approximation of the solution to the boundary value problem (BVP) given by (2.1), (2.2). We use a Galerkin formulation with B-splines as basis functions and select appropriate knot vectors, such that as the polynomial degree increases, the error in the approximation, measured in the energy norm, decays exponentially and independently of the singular perturbation parameters. This is the analog of performing p refinement in the FEM. Even though we focus on the one-dimensional case, our results will serve as the stepping stone to higher dimensions for two reasons: the boundary layer effect is one dimensional (in the direction normal to the boundary) and B-splines in two-dimensions are constructed via tensor products. Hence, this is the first step towards higher dimensions.

The rest of the paper is organized as follows: in Section 2 we describe the model problem, its regularity and the Galerkin formulation we will use. In Section 3 we give a brief introduction to IGA, as described in [3], and present the discrete problem. Section 4 contains our main result of uniform exponential convergence and finally, Section 5 shows the result of a numerical experiment to illustrate the theory.

With I [subset] R an interval with boundary [partial derivative]I and measure |I|, we will denote by [C.sup.k](I) the space of continuous functions on I with continuous derivatives up to order k. We will use the usual Sobolev spaces [W.sup.k,m](I) of functions on I with 0,1,2,..., k generalized derivatives in [L.sup.m] (I), equipped with the norm and seminorm [[parallel] * [parallel].sub.k,m,I] and [| * |.sub.k,m,I], respectively. When m = 2, we will write [H.sup.k] (I) instead of [W.sup.k,2] (I), and for the norm and seminorm, we will write [[parallel]*[parallel].sub.k,I] I and [|*|.sub.k,I], respectively. The usual [L.sup.2](I) inner product will be denoted by [<*,*>.sub.I], with the subscript omitted when there is no confusion. We will also use the space

[H.sub.0.sup.1](I) = {u [member of] [H.sup.1] (I): u[|.sub.[partial derivative]I] = 0} .

The norm of the space [L.sup.[infinity]](I) of essentially bounded functions is denoted by [[parallel] * [parallel].sub.[infinity],I]. Finally, the letter C will denote a generic positive constant, independent of any parameters and possibly having different values in each occurence.

2. The model problem and its regularity. We consider the following model BVP: find u such that

(2.1) -[[epsilon].sub.1]u"(x) + [[epsilon].sub.2]b(x)u'(x) + c(x)u(x) = f(x), x [member of] I = (0, 1),

(2.2) u(0) = u(1) = 0,

where 0 < [[epsilon].sub.1], [[epsilon].sub.2] [less than or equal to] 1 are given parameters that can approach zero and the functions b, c, f are given and sufficiently smooth. We assume that there exist constants [beta],[gamma], [rho], independent of [[epsilon].sub.1], [[epsilon].sub.2], such that [for all]x [member of] I,

(2.3) b(x) [greater than or equal to] [beta] [greater than or equal to] 0, c(x) [greater than or equal to] [gamma] > 0, c(x) - [[[epsilon].sub.2]/[2]]b' (x) [greater than or equal to] [rho] > 0.

The reason for the above assumptions is to ensure existence and uniqueness of a weak solution to (2.1), (2.2).

The structure of the solution to (2.1) depends on the roots of the characteristic equation associated with the differential operator. For this reason, we let [[lambda].sub.0](x), [[lambda].sub.1](x) be the solutions of the characteristic equation and set

[mathematical expression not reproducible]

or equivalently,

(2.4) [mathematical expression not reproducible]

The following holds true [14, 17]:

(2.5) [mathematical expression not reproducible]

The values of [[mu].sub.0], [[mu].sub.1] determine the strength of the boundary layers and since |[[lambda].sub.0] (x)| < |[[lambda].sub.1](x)| the layer at x = 1 is stronger than the layer at x = 0. Essentially, there are three regimes, as shown in Table 1, which is taken from [9]. Figure 2.1 shows the behavior of the solution to (2.1), (2.2), in all three regimes.

The above considerations suggest the following two cases:

1. [[epsilon].sub.1] is large compared to [[epsilon].sub.2]: this is similar to a 'regular perturbation' of reaction-diffusion type. If we consider the limiting case [[epsilon].sub.2] = 0, then we see that there are two boundary layers, one at each endpoint, of width O ([[epsilon].sub.1.sup.1/2] |ln [[epsilon].sub.1.sup.1/2]|). This situation has been studied in the literature (see, e.g., [10]) and will not be considered further in this article.

2. [[epsilon].sub.1] is small compared to [[epsilon].sub.2]: before discussing the different regimes, it is instructive to consider the limiting case [[epsilon].sub.1] = 0. Then there is an exponential layer (of length scale O([[epsilon].sub.2])) at the left endpoint. The homogeneous equation suggests that the different regimes are [[epsilon].sub.1] << [[epsilon].sub.2.sup.2],[[epsilon].sub.1] [approximately equal to] [[epsilon].sub.2.sup.2], [[epsilon].sub.1] >> [[epsilon].sub.2.sup.2].

(a) In the regime [[epsilon].sub.1] << [[epsilon].sub.2.sup.2], we have [[mu].sub.0] = O([[epsilon].sub.2.sup.-1]) and [[mu].sub.1] = O([[epsilon].sub.2][[epsilon].sub.1.sup.-1]). Hence, [[mu].sub.1] is much larger than [[mu].sub.0] and the boundary layer in the vicinity of x = 1 is stronger. Consequently, there is a layer of width O([[epsilon].sub.2]) at the left endpoint (the one that arose from the analysis of the case [[epsilon].sub.1] = 0) and additionally, there is another layer at the right endpoint, of width O([[epsilon].sub.1]/[[epsilon].sub.2]).

(b) In the regime [[epsilon].sub.1] [approximately equal to] [[epsilon].sub.2.sup.2], there are layers at both endpoints of width O ([[epsilon].sub.1.sup.1/2]).

(c) In the regime [[epsilon].sub.2.sup.2] << [[epsilon].sub.1] << 1, there are layers at both endpoints of width O ([[epsilon].sub.1.sup.1/2]).

It was shown in [9, p. 46] (see also [14, Lemma 2.2], which is the result we quote below) that under the assumptions b, c, f [member of] [C.sup.q](I) for some q [greater than or equal to] 1 and q [[parallel]b'[parallel].sub.[infinity,I] [[epsilon].sub.2] [less than or equal to] C(1 - l) for some C, l [member of] (0, 1), the solution u to (2.1), (2.2) can be decomposed into a smooth part [E.sub.S], a boundary layer part at the left endpoint [E.sub.L], and a boundary layer part at the right endpoint [E.sub.R], viz.

u = [E.sub.S] + [E.sub.L] + [E.sub.R],

with

[mathematical expression not reproducible]

for all x [member of] [bar.I] and for n = 0, 1, 2,...,q. If one wants to approximate the solution using a fixed order method (e.g., the h version of the Finite Element Method (FEM)), then this regularity result is sufficient for proving convergence; for a higher order method a more refined regularity result is needed for the smooth part. To this end, we assume that b, c, f satisfy, [for all] n = 0, 1, 2, ...,

(2.6) [mathematical expression not reproducible]

for some positive constants C, [[gamma].sub.f], [[gamma].sub.c], [[gamma].sub.b], independent of [[epsilon].sub.1], [[epsilon].sub.2]. In terms of classical differentiability, we have that the solution to (2.1), (2.2) satisfies (see, e.g., [9, p. 39])

(2.7) [[parallel]u[parallel].sub.[infinity],I] [less than or equal to] C.

The following lemma provides an estimate for u'.

LEMMA 2.1. Let u be the solution of (2.1), (2.2) and assume that (2.3), (2.6) hold. Then there exists a positive constant C, such that

[mathematical expression not reproducible]

Proof. The proof follows [11]. Let

[mathematical expression not reproducible]

and note that A(1) = 0 and A'(x) = -[[[epsilon].sub.2]/[[epsilon].sub.1]]b(x). Then multiplying (2.1) by [e.sup.A(x)] and integrating from x to 1 gives

[mathematical expression not reproducible]

Multiplying by [e.sup.-A(x)] yields

(2.8) [mathematical expression not reproducible]

Integrating from 0 to 1, we further get

(2.9) [mathematical expression not reproducible]

Since we wish to first estimate u'(1), we need upper and lower bounds for [??]dx. From (2.3) we have

(2.10) [mathematical expression not reproducible]

Similarly,

(2.11) [mathematical expression not reproducible]

Also, to estimate the remaining term in (2.9), we consider

[mathematical expression not reproducible]

for some [zeta] between t and x. Hence,

[mathematical expression not reproducible]

Using (2.9)-(2.11), we get

[mathematical expression not reproducible]

Inserting this bound in (2.8) gives

[mathematical expression not reproducible]

as desired.

Using an inductive argument we are able to prove the following.

THEOREM 2.2. Let u be the solution of (2.1), (2.2). Then there exist positive constants K, C, independent of [[epsilon].sub.1], [[epsilon].sub.2], and u, such that for n = 0,1, 2,...

(2.12) [mathematical expression not reproducible]

Proof. The proof is by induction on n and follows from [10]. Equation (2.7) and Lemma 2.1 give the result for n = 0,1, so we assume it holds for 0 [less than or equal to] ? [less than or equal to] n + 1 and show that it holds for n + 2. Differentiating n times (2.1) gives

[mathematical expression not reproducible]

By the induction hypothesis we have

[mathematical expression not reproducible]

Using the estimates below (which follow by standard considerations)

[mathematical expression not reproducible]

[mathematical expression not reproducible]

we obtain

[mathematical expression not reproducible]

Choose the constant K > max{1, [[gamma].sub.f], [[gamma].sub.b], [[gamma].sub.c]} such that the expression in brackets above is bounded by 1, and we have

[mathematical expression not reproducible]

Dividing by [[epsilon].sub.1], gives the desired result.

For simplicity, we will focus on the case

[[epsilon].sub.1.sup.2] << [[epsilon].sub.2].

For the remainder of the article we will make the following assumption:

Assumption 1: Assume there exist positive constants C, [bar.K], [K.sub.1], l, [K.sub.2], [delta] > 0, independent of [[epsilon].sub.1], [[epsilon].sub.2], such that the solution u of (2.1), (2.2) can be decomposed into a smooth part [u.sub.S], two boundary layers at each endpoint [u.sub.BL.sup.[+ or -]], and a remainder [u.sub.R], viz.

(2.13) u = [u.sub.S] + [u.sub.BL.sup.-] + [u.sub.BL.sup.+] + [u.sub.R],

with the following estimates: for every n [member of] [N.sub.0] there holds

(2.14) [mathematical expression not reproducible]

(2.15) [mathematical expression not reproducible]

(2.16) [mathematical expression not reproducible]

for all x [member of] [bar.I].

We provide in the Appendix, for the convenience of the reader, a proof of Assumption 1 in the case of constant coefficients, in order to illustrate the procedure for obtaining such results. The case of variable coefficients is considerably more tedious (cf. [16]).

3. Discretization using isogeometric analysis. Isogeometric analysis may be combined with anumber of formulations; we use Galerkin's approach, i.e., we multiply (2.1) by a suitable test function, integrate by parts and use the boundary conditions (2.2). The resulting variational formulation reads: find u [member of] [H.sub.0.sup.1](I) such that

(3.1) B (u, v) = [<f, v>.sub.I] [for all]v [member of] [H.sub.0.sup.1](I),

where

(3.2) B (u, v) = [[epsilon].sub.1] [<u', v'>.sub.I] + [[epsilon].sub.2] [<bu', v>.sub.I] + [<cu, v>.sub.I].

The bilinear form B (*, *) given by (3.2) is coercive (due to (2.3)) with respect to the energy norm

(3.3) [[parallel]v[parallel].sub.E,I.sup.2] := [[epsilon].sub.1] [|v|.sub.1,I.sup.2] + [rho] [[parallel]v[parallel].sub.0,I.sup.2],

i.e.,

(3.4) B (v, v) [greater than or equal to] [[parallel]v[parallel].sub.E,I.sup.2] [for all]v [member of] [H.sub.0.sup.1] (I).

Next, we restrict our attention to a finite dimensional subspace [V.sub.N] [subset] [H.sub.0.sup.1] (I), that will be selected shortly, and obtain the discrete version of (3.1) as: find [u.sub.N] [member of] [V.sub.N] such that

(3.5) B ([u.sub.N],v) = [<f,v>.sub.I] [for all]v [member of] [V.sub.N].

There holds [9]

B ([u.sub.N] - u, v) = 0 [for all]v [member of] [V.sub.N].

In order to define the space [V.sub.N], we first review the concept of IGA. In this article we use B-splines as basis functions and follow [3] closely. To this end let [XI] = {[[xi].sub.1], [[xi].sub.2], ..., [[xi].sub.n+p+1]} be a knot vector, where [[xi].sub.i] [member of] R is the ith knot, i = 1,2,..., N + p + 1, p is the polynomial order, and N is the number of basis functions. The numbers in [XI] are non-decreasing and may be repeated, in which case we are talking about a non-uniform knot vector. If the first and last knot values appear p + 1 times, the knot vector is called open (see [3] for more details). With a knot vector [XI] in hand, the B-spline basis functions are defined recursively, starting with piecewise constants (p = 0):

[mathematical expression not reproducible]

For p =1,2,..., they are defined by the Cox-de Boor recursion formula [5, 6]:

[mathematical expression not reproducible]

We also mention the recursive formula for obtaining the derivative of a B-spline [3]:

[mathematical expression not reproducible]

We will be considering open knot vectors, having [[xi].sub.1],..., [[xi].sub.m] distinct knots, each with multiplicity ri. Then

[mathematical expression not reproducible]

and there holds [??] [r.sub.i] = N+p+1. Since we are using open knots, we have [r.sub.1] = [r.sub.m] = p+1. The regularity of the B-spline at each knot [[xi].sub.i] is determined by [r.sub.i], in that the B-spline has p - [r.sub.i] continuous derivatives at [[xi].sub.i]. For this reason, we define [k.sub.i] = p - [r.sub.i] + 1 as a measure of the regularity at the knot [[xi].sub.i] and set k = [[k.sub.1],..., [k.sub.m]]. Note that [k.sub.1] = [k.sub.m] = 0 in the case of an open knot vector.

B-splines form a partition of unity and they span the space of continuous piecewise polynomials of degreep on the subdivision {[[xi].sub.1],..., [[xi].sub.m]}.

Each basis function is positive and has support in [[[xi].sub.i], [[xi].sub.i+p+1]]. In the sections that follow, we will approximate the solution to the BVP under consideration, using the space

(3.6) [S.sub.k.sup.p] = span[{[B.sub.k,p]}.sub.k=1.sup.N]

of dimension

N = dim ([S,sub.k.sup.p]) = mp - [m.summation over (i=1)][k.sub.i]

We point out that we are using a uniform polynomial degree p, while we allow for the regularity at each knot to (possibly) vary. A more general approach would be to allow p to vary as well. We will refer to N as the number of degrees of freedom, DOF.

Returning to our problem, the space [V.sub.N] in (3.5) is chosen as [V.sub.N] = [S.sub.k.sup.p], given by (3.6). Thus, we may write the approximate solution as

[u.sub.N] = [N.summation over (k=0)][[alpha].sub.k][B.sub.k,p],

with [alpha] = [[[[alpha].sub.1],..., [[alpha].sub.N]].sup.T] unknown coefficients, and subsitute in (3.5) to obtain the linear system of equations

(3.7) [mathematical expression not reproducible]

where

[mathematical expression not reproducible]

for i, j = 1,..., N. The linear system (3.7) has a unique solution since the matrix M is invertible (due to (3.4) and the linear independence of the basis functions [B.sub.j,p]).

If [[epsilon].sub.1] and [[epsilon].sub.2] are large, then no boundary layers are present and approximating u may be done using a fixed mesh (of say one element) and increasing p. On the other hand, if [[epsilon].sub.1] and [[epsilon].sub.2] are small, then classical techniques fail and the mesh must be chosen carefully. The challenge lies in approximating the typical boundary layer function [e.sup.-x/[epsilon]]. In the context of FDs and FEs, the mesh points must depend on [epsilon], as is well documented in the literature under the name layer-adapted meshes [9]. A similar observation holds for IGA, in the sense that the knot vector must depend on [epsilon]. This was illustrated in [8] through numerical experiments and in the next section we establish it mathematically.

4. Error estimates. We begin by citing the main tool for the analysis (see Corollary 2 in [4]).

PROPOSITION 4.1. Given the subdivision {0 = [[xi].sub.1],..., [[xi].sub.m] = 1} of the reference domain I = (0, 1), let [I.sub.i] = ([[xi].sub.i],[[xi].sub.i+1]), [h.sub.i] = [[xi].sub.i+1] - [[xi].sub.i],i = 1,...,m - 1,k,pnon-negative integers with p [greater than or equal to] 2k - 1 and [u.sup.(k)] [member of] [H.sup.s] (I) for some 0 [less than or equal to] s [less than or equal to] [kappa]=p - k + 1. Then there exists a quasi-interpolation operator [[pi].sub.p,k] such that, for i = 1,...,m - 1, j = 1,...,k - 1,

[mathematical expression not reproducible]

In [2], the above was specialized for maximal smoothness by selecting p = 2q - 1, for some q [greater than or equal to] 1, and then k = [kappa] = q, so the estimate of Proposition 4.1 becomes

(4.1) [mathematical expression not reproducible]

for 0 < s [less than or equal to] q, j = 1, ..., q - 1.

In view of the above, from now on we will be using the symbol q to denote the polynomial degree.

The knot vector described below, is 'inspired' by the so-called spectral boundary layer mesh used in the hp-FEM for such problems (see [12] and the references therein).

DEFINITION 4.2. Let [[mu].sub.0], [[mu].sub.1] be given by (2.4) and let [lambda] [greater than or equal to] 1 be a user specified parameter. Then, if [lambda]q[[mu].sub.1.sup.-1] [greater than or equal to] 1/2,

(4.2) [mathematical expression not reproducible]

and, if [lambda]q[[mu].sub.0.sup.-1] < 1/2,

(4.3) [mathematical expression not reproducible]

where q is the polynomial degree.

The following auxiliary result will be used in the sequel.

LEMMA 4.3. There exists [tau] [member of] [1/2, 2/3] such that for every q [greater than or equal to] 2, 3,4,..., there holds s := [tau]q [member of] N and

[mathematical expression not reproducible]

where

(4.4) [mathematical expression not reproducible]

Proof. We first show [tau]q [member of] N, for some [tau] [member of] [1/2,2/2]. If q is even, then we simply take [tau] = 1/2. If q is odd, i.e., q = 2n + 1, n [member of] N, then we choose [tau] = (n + 1)/(2n + 1).

Now, with [tau]q [member of] N, we have (q [+ or -] [tau]q)! = G(q [+ or -] [tau]q + 1) and, as q [right arrow] [infinity] [1],

[mathematical expression not reproducible]

THEOREM 4.4. Let u be the solution to (2.1), (2.2) and let [[pi].sub.2q-1,q]u [member of] [S.sub.q.sup.2q-1] be its quasi-interpolant, constructed using the knot vector of Definition 4.2 and satisfying (4.1). Then there exist positive constants [sigma], C, independent of [[epsilon].sub.1], [[epsilon].sub.2], such that as q [right arrow] [infinity] there holds

[[parallel]u - [[pi].sub.2q-1,q]u[parallel].sub.E,I] [less than or equal to] [Ce.sup.-[sigma]q].

Proof The proof is separated into two cases. In Case 1, we make use of the classical differentiability result (2.12) and construct a quasi-interpolant of u on the entire interval I with the desired properties. In Case 2, we use the decomposition (2.13) and construct quasi-interpolants for each piece separately, as follows: for the smooth part, the quasi-interpolant is constructed on I, while for the boundary layers, it has support only in the layer regions. No quasi-interpolant is constructed for the remainder, since it is already exponentially small. We next give the details.

Case 1: [lambda]q[[mu].sub.1.sup.-1] [greater than or equal to] 1/2 or equivalently 2[lambda]q [greater than or equal to] [[epsilon].sub.1.sup.-1]. In this case, the knot vector is given by (4.2). We make use of (2.12) and (4.1) to obtain a quasi-interpolant [[pi].sub.2q-1,q]u on I, such that

[mathematical expression not reproducible]

We choose s = [tau]q, with [tau] [member of] [1/2,2/3] as asserted by Lemma 4.3, and we note that, from (2.5) and the fact that in this case there holds 2[lambda]q [greater than or equal to] [[epsilon].sub.1.sup.-1], we have

[mathematical expression not reproducible]

since 2[lambda] [greater than or equal to] [tau] + 1. Then, with the aid of Lemma 4.3, we get

[mathematical expression not reproducible]

with [[sigma].sub.1] given by (4.4). By Stirling s formula [mathematical expression not reproducible], hence

[mathematical expression not reproducible]

Selecting [lambda] such that [(2eK [lambda]).sup.2[tau]] [(eK/2).sup.2] < 1, gives

[mathematical expression not reproducible]

where C > 0 is independent of q [member of] N and [[epsilon].sub.1]. Reducing the value of [[sigma].sub.1] to [??] [member of] (0, [[sigma].sub.1]), there exists a positive constant C, independent of q and [[epsilon].sub.1], such that

(4.5) [mathematical expression not reproducible]

Case 2: [lambda]q[[mu].sub.0.sup.-1] < 1/2 or equivalently, 2q[[epsilon].sub.1] < 1/[lambda]. In this case, the solution is decomposed as in (2.13) with the bounds (2.14)-(2.16) being valid. We will approximate each term in (2.13) separately, using Proposition 4.1 to construct appropriate quasi-interpolants over the three intervals

(4.6) [I.sub.1] = (0, q[[mu].sub.0.sup.-1]), [I.sub.2] = (q[[mu].sub.0.sup.-1], 1 - q[[mu].sub.1.sup.-1]), [I.sub.3] = (q[[mu].sub.1.sup.-1], 1).

Throughout the proof, we choose

[mathematical expression not reproducible] (4.7)

with K, [K.sub.1], [K.sub.2] the constants in (2.14)-(2.16). So we have, with the obvious notation,

u - [[pi].sub.2q-1,q]u = [u.sub.S] + [u.sub.BL.sup.-] + [u.sub.B.sup.+] + [u.sub.R] - [[pi].sub.2q-1,q][u.sub.S] - [[pi].sub.2q-1,q][u.sub.BL.sup.-] - [[pi].sub.2q-1,q][u.sub.BL.sup.+],

hence

|u - [[pi].sub.2q-1,q]u| [less than or equal to] |[u.sub.S] - [[pi].sub.2q-1,q][u.sub.S]| + |[u.sub.BL.sup.-] - [[pi].sub.2q-1,q][u.sub.BL.sup.-]| + |[u.sub.B.sup.+] - [[pi].sub.2q-1,q][u.sub.BL.sup.+]| + |[u.sub.R]|.

For the smooth part, Proposition 4.1 ensures that [[pi].sub.2q-1,q][u.sub.S], is such that

[mathematical expression not reproducible]

Choosing s = [??]q with [??] [member of] [1/2,2/3] as in Lemma 4.3, and using Lemma 1 of [2], we get

[mathematical expression not reproducible]

by (4.7). Therefore, there exist constants C, [[sigma].sub.2] > 0 such that for every [[epsilon].sub.1], [[epsilon].sub.2] [member of] (0, 1], and for every q [member of] N, there holds

(4.8) [mathematical expression not reproducible]

For the left boundary layer [u.sub.BL.sup.-] we obtain from Proposition 4.1, a quasi-interpolant [[pi].sub.2q-1,q][u.sub.BL.sup.-] on the interval [I.sub.1] = (0, [lambda]q[[mu].sub.0.sup.-1]), such that

[mathematical expression not reproducible]

Using (2.5), we obtain

[mathematical expression not reproducible]

where (4.7) was used. Choosing s = [bar.[tau]]q, with [bar.[tau]] [member of] [1/2,2/3] as in Lemma 4.3, we get (for [[sigma].sub.3] given by (4.4) with [tau] replaced by [bar.[tau]]),

[mathematical expression not reproducible]

Since, by (4.7) and Stirling s formula, [mathematical expression not reproducible], we have

[mathematical expression not reproducible]

with C > 0 independent of [[epsilon].sub.1] and q. On I\[I.sub.1], [u.sub.BL.sup.-] is already exponentially small, i.e., by (2.15), we have

[mathematical expression not reproducible]

hence it will not be approximated. Thus, there exist constants C,[[bar.[sigma]].sub.3] > 0 such that for every [[epsilon].sub.1], [[epsilon].sub.2] [member of] (0, 1], and for every q [member of] N, there holds

(4.9) [mathematical expression not reproducible]

where again (2.5) was used. For the [L.sup.2] error, we have

[mathematical expression not reproducible]

Choosing s as a non-integer multiple of q, and repeating the same steps as for the [H.sup.1]-seminorm, we see that there exist constants C, [[sigma].sub.4] > 0 such that for every [[epsilon].sub.1], [[epsilon].sub.2] [member of] (0, 1], and for every q [member of] N, there holds

[mathematical expression not reproducible]

On I\[I.sub.1], [u.sub.BL.sup.-] is already exponentially small, hence there exist constants C, [[sigma].sub.5] > 0 such that for every [[epsilon].sub.1], [[epsilon].sub.2] [member of] (0, 1], and for every q [member of] N, there holds

(4.10) [mathematical expression not reproducible]

We follow the same steps for the right boundary layer [u.sub.BL.sup.+]: Proposition 4.1 gives [[pi].sub.2q-1,q][u.sub.BL.sup.+] on [I.sub.3], such that

[mathematical expression not reproducible]

Using (2.5), we obtain

[mathematical expression not reproducible]

Choosing s = [??]q, with [??] [member of] [1/2,2/3] as in Lemma 4.3, and following the same steps as in the approximation of the left boundary layer (using (4.7)), we arrive at

[mathematical expression not reproducible]

with [[sigma].sub.6] given by (4.4) with [tau] replaced by [??]. On I\[I.sub.3], [u.sub.BL.sup.+] is already exponentially small, i.e., by (2.15) we have

[mathematical expression not reproducible]

hence it will not be approximated. Thus

(4.11) [mathematical expression not reproducible]

for some positive constant [??], where (2.5) was used once more. For the [L.sup.2] error, we have

[mathematical expression not reproducible]

where we used (2.5) and (4.7). The remaining steps are the same as above, so they are omitted. We arrive at

(4.12) [mathematical expression not reproducible]

for some positive constant [[sigma].sub.7]. Finally,

(4.13) [mathematical expression not reproducible]

for some positive constant [??].

It remains to consider the remainder, which by (2.16) is exponentially small:

(4.14) [mathematical expression not reproducible]

since [lambda]p[[mu].sub.0.sup.-1] < 1/2 (hence [lambda]p[[epsilon].sub.2] < 1/2), where [zeta] > 0 is a constant independent of [[epsilon].sub.1], [[epsilon].sub.2].

The proof is completed by using the definition of the energy norm (3.3), along with (2.16), (4.5), (4.8)-(4.14), to get

[mathematical expression not reproducible]

for some positive consant [sigma].

We next estimate the difference between the IGA solution [u.sub.N] and the quasi-interpolant [[pi].sub.2q-1,q]u.

LEMMA 4.5. Let [u.sub.N] [member of] [S.sub.q.sup.2q-1] be the approximation to u, the solution of (2.1), (2.2), based on the knot vector of Defintion 4.2, and let [[pi].sub.2q-1,q] be the approximation operator of Proposition 4.1. Then there exist constants C, [sigma] > 0, independent of [[epsilon].sub.1], [[epsilon].sub.2], such that [for all]q [member of] N,

[mathematical expression not reproducible]

Proof Set [xi] := [[pi].sub.2q-1,q]u - [u.sub.N]. Then, by coercivity of the bilinear form [B.sub.[epsilon] (eq. (3.4)), there holds

[mathematical expression not reproducible]

where we also used Galerkin orthogonality. Hence,

[mathematical expression not reproducible]

The first and last term may be estimated using Cauchy-Schwarz:

[mathematical expression not reproducible]

For the second term, we will consider the two ranges of q separately: in the asymptotic range of q, i.e., [lambda]q[[mu].sub.1.sup.-1] [greater than or equal to] 1/2 or equivalently [lambda]q[[epsilon].sub.1] [greater than or equal to] 1/2, we have

[mathematical expression not reproducible]

In the pre-asymptotic range of q, i.e., [lambda]q[[mu].sub.0.sup.-1] < 1/2, we first use integration by parts to obtain

[mathematical expression not reproducible]

Next, we consider the three subintervals given by (4.6): on the first subinterval we have

[mathematical expression not reproducible]

where we used an inverse inequality; see, e.g., [15, Thm. 3.91]. Thus

[mathematical expression not reproducible]

Similarly, on the second subinterval we have

[mathematical expression not reproducible]

Finally, on the third subinterval,

[mathematical expression not reproducible]

Therefore, by (4.12), we get

[mathematical expression not reproducible]

and

[mathematical expression not reproducible]

which completes the proof.

We conclude with our main result.

THEOREM 4.6. Let u be the solution to (2.1), (2.2) and let [u.sub.N] [member of] [S.sub.q.sup.2q-1] be its approximation, based on the knot vector of Defintion 4.2. Then there exist positive constants C, [sigma], independent of [[epsilon].sub.1], [[epsilon].sub.2], such that, [for all]q [member of] N,

[[parallel]u -[u.sub.N][parallel].sub.E,I] [less than or equal to] [Ce.sup.-[sigma]q].

Proof. We begin with the triangle inequality:

[mathematical expression not reproducible]

where [[pi].sub.2q-1,q] is the approximation operator of Theorem 4.4. The first term is handled by Lemma 4.4 and the second by Lemma 4.5.

5. Numerical example. We will be considering the following BVP and we refer to [8] for more numerical experiments. Find u such that

[mathematical expression not reproducible]

An exact solution is not available, so we use as a reference solution the one computed with the highest polynomial degree, denoted by [u.sub.REF]. Instead of the error in the energy norm, we will be measuring

[mathematical expression not reproducible]

where [{[x.sub.k]}.sub.k=1.sup.n] [member of] I are points in (0,1), chosen uniformly in the layer region and outside--we use n = 400 in each region for our computations.

Figure 5.1 shows the Error vs. the number of degrees of freedom, DOF, in a semi-log scale. The curves indicate exponential convergence and the fact that they coincide indicates robustness. The two curves that are not on top of the rest, correspond to even smaller errors and could be due to the fact that we are using a reference solution.

6. Conclusions. In this article we performed the numerical analysis of IGA for one-dimensional reaction-convection-diffusion problems with two small parameters. We established that if the knot vector is chosen appropriately and depending on the singular perturbation parameters, then p-refinement yields robust, exponential rates of convergence, in the energy norm. We also presented one numerical example agreeing with, and even extending the theory.

This was a step towards the study of two-dimensional SPPs, especially fourth order SPPs, for which there are very few available methods for general two-dimensional domains.

Appendix A. Here we establish eqs. (2.14)-(2.16), in the case of constant coefficients, i.e., b(x) = b > 0, c(x) = c > 0. In [16], the non-constant coefficient case is considered; we present here the much simpler case of constant coefficients in order to illustrate the procedure for obtaining such regularity estimates, for the benefit of the reader.

Recall that we are focusing on the case [[epsilon].sub.1.sup.2] << [[epsilon].sub.2], hence we anticipate a layer of width O([[epsilon].sub.2]) at the left endpoint and a layer of width O([[epsilon].sub.1]/[[epsilon].sub.2]) at the right endpoint. To deal with this we define the stretched variables x = x/[[epsilon].sub.2] and x = (1 - x)[[epsilon].sub.2]/[[epsilon].sub.1] and make the formal ansatz

(A.1) [mathematical expression not reproducible]

with [u.sub.i,j], [??], [??] to be determined. Substituting (A.1) into (2.1), separating the slow (i.e., x) and fast (i.e., x, x) variables, and equating like powers of [[epsilon].sub.1] and [[epsilon].sub.2], we get

(A.2) [mathematical expression not reproducible]

(A.3) [mathematical expression not reproducible]

(A.4) [mathematical expression not reproducible]

The last two equations are supplemented with the following boundary conditions (in order for (2.2) to be satisfied) for all i, j [greater than or equal to] 0:

(A.5) [mathematical expression not reproducible]

Note that, by (A.3) and the fact that b, c > 0, we automatically have [mathematical expression not reproducible].

Next, we would like to describe the regularity of the functions [u.sub.i,j], [??], [??], defined by (A.2)-(A.5) above. We begin with [u.sub.i,j], and we have the following.

LEMMA A.1. Let uij be defined by (A.2) and assume (2.6) holds. Then there exist positive constants C, K and a complex neighborhood G of [bar.I] such that the complex extension ofu (denoted again by u) satisfies

[mathematical expression not reproducible]

Proof. The proof is by induction on i. The case i = 0 holds trivially, so assume the result holds for i and establish it for i + 1. Let [tau] [member of] (0, 1) and let K > 0 be a constant so that ([[2]/[K.sup.2]] + [[1]/[K]]) [less than or equal to] C. We have by (A.2), the induction hypothesis with [G.sub.(1-[tau])[delta]] [contains] [G.sub.[delta]], and Cauchys Integral Theorem,

[mathematical expression not reproducible]

Choose[tau] = 1/(i + 1). Then

[mathematical expression not reproducible]

so by the choice of K the expression in brackets is bounded and this completes the proof.

LEMMA A.2. Let [u.sub.i,j] be defined by (A.2) and assume (2.6) holds. Then there exist positive constants C, [K.sub.1], [K.sub.2] such that

[mathematical expression not reproducible]

Proof. This follows immediately from Lemma A.1 and Cauchy's Integral Theorem for derivatives:

[mathematical expression not reproducible]

with [K.sub.1] = e, [K.sub.2] = K/[delta].

The following auxiliary lemma, which is an analog of Lemma 7.3.6 in [10], will be used in the proof of Lemma A.4.

LEMMA A.3. Let [lambda], [gamma] [member of] C with Re([lambda]) > 0,Re([gamma]) > 0. Let F be an entire function satisfying, for some i,j [member of] [N.sub.0] and C > 0,

[mathematical expression not reproducible]

Let [alpha] [member of] C and let v : (0, [infinity]) [right arrow] C, be the solution of the problem

v' + [lambda]v = F on (0, [infinity]), v(0) = [alpha].

Then v can be extended to an entire function (denoted again by v), which satisfies

[mathematical expression not reproducible]

Proof. Using an integrating factor, we find

[mathematical expression not reproducible]

from which we have

[mathematical expression not reproducible]

where we used the assumption on F. The result follows.

LEMMA A.4. The functions [??], which satisfy (A3), are entire and there exist positive constants C, K, [??], depending only on the data, such that

(A.6) [mathematical expression not reproducible]

Proof. The proof is by induction. We first note that from (A.3), we may calculate

[mathematical expression not reproducible]

Thus, using Lemma A.2 to bound the term |[u.sub.i,0](0)|, we get

[mathematical expression not reproducible]

where [??] = [K.sub.2], [beta] = c/b, hence the result holds for j = 0 and for all i [greater than or equal to] 0. We then proceed by induction on j. We assume (A.6) holds for j [greater than or equal to] 1 (and for all i [greater than or equal to] 0) and show it for j + 1. We note that, by (A3), [??] satisfies, [for all] i, j [greater than or equal to] 0,

[mathematical expression not reproducible]

By the induction hypothesis and Cauchy's Integral Theorem for Derivatives (we take as the contour the unit circle centered at z), we get

[mathematical expression not reproducible]

Lemma A3 is then applicable (with [lambda] = [[c]/[b]], [gamma] = [??], F = [[1]/[b]]([??])" and [alpha] = -[u.sub.i,j+1](0)) and with the aid of Lemma A.2 (to bound |[alpha]|) we obtain

[mathematical expression not reproducible]

where we used the fact that [??] = [K.sub.2]. Hence, the induction is complete and this concludes the proof.

COROLLARY A.5. The functions [??], which satisfy (A3), are entire and there exist positive constants C, K, [K.sub.1], [??], depending only on the data, such that, [for all] n [member of] [N.sub.0] and x [greater than or equal to] 0, there holds

[mathematical expression not reproducible]

Proof. Cauchy's Integral Theorem for Derivatives gives

[mathematical expression not reproducible]

The proof is completed by observing that

(A.7) [mathematical expression not reproducible]

REMARK A.6. An analogous result may be proven for the functions [??], which satisfy (A.4), (A.5):

(A.8) [mathematical expression not reproducible]

Indeed, from (A.4) we find

[mathematical expression not reproducible]

and an inductive argument, like in the proof of Lemma A.4, can be used to establish (A.8). Next, we define, for some M [member of] Z,

(A.9) [mathematical expression not reproducible]

(A.10) [mathematical expression not reproducible]

(A.11) [mathematical expression not reproducible]

and we have the following decomposition

(A.12) [mathematical expression not reproducible]

THEOREM A.7. Assume that (2.6), (2.3) hold. Then there exist positive constants C,[K.sub.1],[K.sub.2], K,K,[??],[??], [delta], independent of [[epsilon].sub.1], [[epsilon].sub.2], such that the solution u of (2.1), (2.2) can be decomposed as in (A.12), with

(A.13) [mathematical expression not reproducible]

(A.14) [mathematical expression not reproducible]

(A.15) [mathematical expression not reproducible]

(A.16) [mathematical expression not reproducible]

provided [epsilon]2eM max {[K.sub.2],[??], [??]} < 1 and [[[epsilon].sub.1]/[[epsilon].sub.2.sup.2]]eM max{[??],[??]} < 1.

Proof. We first show (A.13): from (A.9) and Lemma A.2 we have

[mathematical expression not reproducible]

since both sums are convergent geometric series due to the assumptions [[epsilon].sub.2]M [K.sub.2] < 1 and [[epsilon].sub.1]/[[epsilon].sub.2.sup.2] < 1.

Next, we show (A.14): by (A.10) and Lemma A.4, we have

[mathematical expression not reproducible]

Now, [(i +j).sup.i+j] [less than or equal to] [e.sup.i][i.sup.i][e.sup.j][j.sup.j] (cf. (A.7)) hence we get

[mathematical expression not reproducible]

since both sums are convergent geometric series, due to the assumptions [??]eM[[epsilon].sub.2] < 1, [[[epsilon].sub.1]/[[epsilon].sub.2.sup.2]]eM < 1. The result follows.

Similarly, we show (A.15): by (A.11) and (A.8),

[mathematical expression not reproducible]

It remains to show (A.16). To this end, note that

[mathematical expression not reproducible]

By (A.8),

[mathematical expression not reproducible]

for some positive [delta], independent of [[epsilon].sub.1], [[epsilon].sub.2] and bounded away from 0. Similarly,

[mathematical expression not reproducible]

Combining the two results, we have

[mathematical expression not reproducible]

Now, let L:= -[[epsilon].sub.1[d.sup.2]/[dx.sup.2] + [[epsilon].sub.2]b[[d]/[dx]] + c Id, with Id the identity operator, and consider

[mathematical expression not reproducible]

with [u.sub.i,j] satisfying (A.2). After some calculations, we find

[mathematical expression not reproducible]

hence

[mathematical expression not reproducible]

Using Lemma A.2, we further obtain

[mathematical expression not reproducible]

since the finite sum can be bounded by a converging geometric series.

We also consider the operator L in the stretched variable x:

[mathematical expression not reproducible]

and we find, after some calculations,

[mathematical expression not reproducible]

where (A.3) was used. Hence, using (A.6), we have

[mathematical expression not reproducible]

Similarly, in the stretched variable x we have

[mathematical expression not reproducible]

and, with the help of (A.4),

[mathematical expression not reproducible]

thus, by following the exact same steps as above,

[mathematical expression not reproducible]

Therefore,

[mathematical expression not reproducible]

Under the assumptions of the theorem, we have shown that the remainder [r.sub.M] has exponentially small values at the endpoints of I, and [Lr.sub.M] is uniformly bounded by an arbitrarily small quantity on I. By, e.g., stability (see [9]) we have the desired result.

REFERENCES

[1] G. E. ANDREWS, R. ASKEY, AND R. ROY, Special Functions, Cambridge University Press, Cambridge, 1999.

[2] A. BUFFA, G. SANGALLI, AND C. SCHWAB, Exponential convergence of the hp version of isogeometric analysis in 1D, in Spectral and High Order Methods for Partial DifferentialEequations--ICOSAHOM 2012, M. Azaiez, H. El Fekih, J. S. Hesthaven, eds, Lect. Notes Comput. Sci. Eng., 95, Springer, Cham (2014), pp. 191-203.

[3] J. A. COTTRELL, T. R. HUGHES, AND Y. BASILEVS, Isogeometric Analysis: Toward Integration of CAD and FEA, Wiley, Chichester, 2009.

[4] L. BEIRAO DA VEIGA, A. BUFFA, J. RIVAS, AND G. SANGALLI, Some estimates for h-p-k refinement in isogeometric analysis, Numer. Math., 118 (2011), pp. 271-305.

[5] M. G. COX, The numerical evaluation of B-splines, National Physics Laboratory Report DNAC 4, Teddington, 1971.

[6] C. DE BOOR, On calculation with B-splines, J. Approx. Theory, 6 (1972), pp. 50-62.

[7] T. R. HUGHES, J. A. COTTRELL, AND Y. BASILEVS, Isogeometric analysis: CAD, finite elements, NURBS, exact geometry and mesh refinement, Comput. Meth. Appl. Mech. Eng., 194 (2005), pp. 4125-4195.

[8] K. LIOTATI AND C. XENOPHONTOS, Isogeometric analysis for singularly perturbed problems in 1-D: a numerical study, Preprint on arXiv, 2018, http://arxiv.org/abs/1901.01949. To appear in the Proceedings of BAIL 2018.

[9] T. LINSS, Layer-Adapted Meshes for Reaction-Convection-Diffusion Problems, Springer, Berlin, 2010.

[10] J. M. MELENK, On the robust exponential convergence of hp finite element methods for problems with boundary layers, IMA J. Numer. Anal., 17 (1997), pp. 577-601.

[11] J. M. MELENK AND C. SCHWAB, An hp finite element method for convection-diffusion problems in one-dimension, IMA J. Numer. Anal., 19 (1999), pp. 425-453.

[12] J. M. MELENK AND C. XENOPHONTOS, Robust exponential convergence of hp-FEM in balanced norms for singularly perturbed reaction-diffusion equations, Calcolo, 53 (2016), pp. 105-132.

[13] H.-G. ROOS, M. STYNES, AND L. TOBISKA, Robust Numerical Methods for Singularly Perturbed Differential Equations, Springer, Berlin, 2008.

[14] H.-G. ROOS AND Z. UZELAC, The SDFEM for a convection-diffusion problem with two small parameters, Comput. Methods Appl. Math., 3 (2003), pp. 443-458.

[15] C. SCHWAB, p- and hp-Finite Element Methods, Oxford University Press, New York, 1998.

[16] I. SYKOPETRITOU AND C. XENOPHONTOS, Analytic regularity for a singularly perturbed reaction-convectiondiffusion boundary value problem with two small parameters, Preprint on arXiv, 2019. http://arxiv.org/abs/1901.09397 Article to appear in the Mediterranean Journal of Mathematics.

[17] L. TEOFANOV AND H.-G. ROOS, An elliptic singularly perturbed problem with two parameters I: Solution decompostition, J. Comput. Appl. Math., 206 (2007), pp. 1802-1097.

CHRISTOS XENOPHONTOS ([dagger]) AND IRENE SYKOPETRITOU ([dagger])

(*) Received October 2, 2018. Accepted November 8, 2019. Published online on January 16, 2020. Recommended by Torsten LinB.

([dagger]) Department of Mathematics and Statistics, University of Cyprus, P.O. Box 20537, 1678 Nicosia, Cyprus (xenophontos@ucy.ac.cy).

DOI: 10.1553/etna_vol52s1
TABLE 2.1
Different regimes based on the relationship between [[epsilon].sub.1]
and [[epsilon].sub.2].

convection-diffusion       [[epsilon].sub.1] << [[epsilon].sub.2] =1
convection-reaction-diffusion  [[epsilon].sub.1] << [[epsilon]
.sub.2.sup.2] << 1
reaction-diffusion       [[epsilon].sub.2.sup.2] << [[epsilon]
.sub.1] << 1

[[mu].sub.0]

convection-diffusion       1
convection-reaction-diffusion  [[epsilon].sub.2.sup.-1]
reaction-diffusion       [[epsilon].sub.1.sup.-1/2]

[[mu].sub.1]

convection-diffusion       [[epsilon].sub.1.sup.-1]
convection-reaction-diffusion  [[epsilon].sub.2]/[[epsilon].sub.1]
reaction-diffusion       [[epsilon].sub.1.sup.-1/2]
COPYRIGHT 2020 Institute of Computational Mathematics
No portion of this article can be reproduced without the express written permission from the copyright holder.