Parameter-robust numerical method for a system of singularly perturbed initial value problems

Numer Algor (2012) 59:185–195 DOI 10.1007/s11075-011-9483-4 ORIGINAL PAPER

Parameter-robust numerical method for a system of singularly perturbed initial value problems Sunil Kumar · Mukesh Kumar

Received: 3 November 2010 / Accepted: 9 June 2011 / Published online: 8 July 2011 © Springer Science+Business Media, LLC 2011

Abstract In this work we study a system of M(≥ 2) first-order singularly perturbed ordinary differential equations with given initial conditions. The leading term of each equation is multiplied by a distinct small positive parameter, which induces overlapping layers. A maximum principle does not, in general, hold for this system. It is discretized using backward Euler difference scheme for which a general convergence result is derived that allows to establish nodal convergence of O(N −1 ln N) on the Shishkin mesh and O(N −1 ) on the Bakhvalov mesh, where N is the number of mesh intervals and the convergence is robust in all of the parameters. Numerical experiments are performed to support the theoretical results. Keywords Singular perturbation · Initial-value problem · Layer-resolving meshes · Parameter-robust convergence 1 Introduction In this work we study the following system of M(≥ 2) coupled singularly perturbed initial value problems Lu := Eu + Au = f in  = (0, 1], u(0) = η,

S. Kumar Department of Mathematics, Indian Institute of Technology Delhi, Hauz Khas, New Delhi, 110 016, India e-mail: [email protected] M. Kumar (B) Institute of Applied Mathematics and Information Technology-CNR, Via Ferrata 1, Pavia, 27100, Italy e-mail: [email protected]

(1)

186

Numer Algor (2012) 59:185–195

where u = (u1 , . . . , u M )T , f = ( f1 , . . . , f M )T and η = (η1 , . . . , η M )T are column vectors; A = ( aij(x) ) is an M × M matrix, and E = diag(ε1 , . . . , ε M ) is a diagonal matrix with small perturbation parameters 0 < εk  1, k = 1, . . . , M. Such systems of singularly perturbed initial value problems arise in many areas of science and engineering (see [1, 3]). In recent years there has been a growing interest in design and analysis of parameter-robust numerical methods for systems of singularly perturbed differential equations (see [6, 10]). Parameter-robust numerical methods for systems of type (1) are analyzed in [4, 8, 11]. In [8], a system of two equations having small parameter associated with only one equation is considered. The case when all the perturbation parameters are equal is considered in [4]. While, the most difficult case allowing perturbation parameters of different magnitudes associated with each equation (as is the case in the current work) is considered in [11]. All these works impose conditions on the coupling matrix A such as an M-matrix structure so that the operator L and its discretization obey certain maximum principles, which are the essential parts of their analysis. In the current work we relax the assumptions on the coupling matrix A, as we do not impose any M-matrix structure onto A, see (2). As a result, the operator L and its discretization do not, in general, obey maximum principles. We show that problem (1) can be analyzed without using maximum principles for the operator L and its discretization. Thus, we covers a larger class of problems. Furthermore, we derive a general convergence result that allows to establish the robust convergence of the present method on both Shishkin and Bakhvalov meshes, whereas the works in [4, 8, 11] are limited to Shishkin meshes. In the literature two techniques are commonly used to analyze parameterrobust numerical methods for singularly perturbed problems; truncation error and barrier function technique and discrete Green’s function technique (see [6, 10]). The former is simpler than the latter, but so far has been limited to piecewise-uniform Shishkin meshes. While, the latter allows to deal with a larger class of layer-resolving meshes, including piecewise-uniform Shishkin meshes and graded Bakhvalov meshes, in the same framework. A similar weakened assumptions on the coupling matrix A, as given in (2), has recently been considered for systems of singularly perturbed reaction-diffusion twopoint boundary value problems in [7]. The authors used discrete Green’s function technique to establish the robust convergence of the method on both Shishkin and Bakhvalov meshes. In contrast, here we consider a relatively simple truncation error and barrier function technique and show that it can be used to establish the robust convergence of the present method on both Shishkin and Bakhvalov meshes in the same framework. This paper is organized as follows. Assumptions on the continuous problem and the properties of its solution are given in Section 2. In Section 3, the continuous problem is discretized using backward Euler difference scheme for which a general convergence result is derived in Section 4 that allows to establish the robust convergence of O(N −1 ln N) on the Shishkin mesh and

Numer Algor (2012) 59:185–195

187

O(N −1 ) on the Bakhvalov mesh. Finally, results of numerical experiments are given in Section 5. Notations Throughout the paper, C denotes a generic positive constant that is independent of all perturbation parameters and the discretization parameter. Similarly, C = (C, C, . . . , C)T . Define v ≤ w if vi ≤ wi , 1 ≤ i ≤ M, and |v| = (|v1 |, . . . , |v M |)T . For any function g ∈ C(), define g j = g(x j); if g ∈ C() M then g j = g(x j) = (g1, j, . . . , g M, j)T . For any function g ∈ C(), define g = maxx∈ |g(x)|; if g ∈ C() M , define g = maxi=1,...,M gi . The analogous disN crete maximum norm on the mesh  is denoted by ||.|| N . 2 The continuous problem We consider problem (1) assuming that aii (x) > 0, x ∈ , i = 1, . . . , M,

(2a)

  M   aik    < ξ < 1, i = 1, . . . , M, a 

(2b)

k=1,k=i

ii

for some constant ξ. The coupling matrix A has entries in C1 (), and the column vector f ∈ C1 () M . Without loss of generality we consider the following ordering ε1 ≤ ε2 ≤ · · · ≤ ε M .

(3)

Under assumption (2), the operator L does not, in general, obey a maximum principle. Using an idea similar to [7], we shall prove that it is maximum norm stable. For this purpose, we need the following lemma. Lemma 2.1 Let y ∈ C1 () ∩ C() satisfy μy + b y = f in , with μ > 0 and b > 0 on . Then     f   ||y|| ≤ max   , |y(0)| . b Proof Introducing T v := μv  + b v, we use contradiction argument to prove that the operator T satisfies the maximum principle, viz., if v(0) ≥ 0 and T v ≥ [9]). Then, using this maximum principle with the 0 in  then v ≥ 0 in  (see  f ± barrier function ϕ = max  b  , |y(0)| ± y, we get the desired estimate.

   a  Let ϒ = (γij) be a M × M matrix with entries γii = 1 and γij = −  aijii , for

i = j. From (2b), it is easy to see that ϒ is inverse monotone; that is, ϒ is invertible and ϒ −1 > 0.

188

Numer Algor (2012) 59:185–195

Lemma 2.2 (Stability of continuous problem) Let u be the solution of (1) and let the coupling matrix A satisfy (2). Then    M   fk   , |ηk | , i = 1, . . . , M. ui  ≤ ( ϒ −1 )i k max  (4) a  kk k=1

Proof We split the solution of (1) into two parts u = p + q, where the ith components of p and q are the solutions of εi pi + aii pi = fi in , pi (0) = ηi (0), and εi qi + aii qi = −

M 

aik uk in , qi (0) = 0.

k=1,k=i

For each i, application of Lemma 2.1 gives      M   fi   aik      uk . || pi || ≤ max   , |ηi (0)| and ||qi || ≤ a  aii ii k=1,k=i

Hence   M   aik    uk , i = 1, . . . , M. ui  ≤  pi  + qi  ≤  pi  + a  ii k=1,k=i

Consequently ui  −

     M   aik      uk  ≤ max  fi  , |ηi (0)| , i = 1, . . . , M. a  a  ii ii

k=1,k=i

Invoking the inverse monotonicity of ϒ, the desired result follows.



Remark 2.3 From Lemma 2.2 and standard arguments (see [5]), problem (1) has a unique solution u ∈ C2 () M . Introducing the layer functions Bi (x) := e−xα/εi , i = 1, . . . , M, where α := (1 − ξ ) min min aii (x), i=1,...,M x∈

(5)

we analyze the behavior of the solution of (1) and its derivatives by decomposing it into two parts u = v + w, where v and w are respectively the regular and layer parts of u. Theorem 2.4 Suppose that the coupling matrix A satisf ies (2). Then the solution u of (1) can be decomposed as u = v + w, where v and w are the solutions of Lv = f in , v(0) = A(0)−1 f(0)

(6a)

Numer Algor (2012) 59:185–195

189

and Lw = 0 in , w(0) = u(0) − v(0)

(6b)

vi(k)  ≤ C, k = 0, 1, εi vi(2)  ≤ C, i = 1, . . . , M.

(7)

satisfying

|wi(k) (x)| ≤ C

M 

ε −k B (x), k = 0, 1; i = 1, . . . , M.

(8)

=i

εi |wi(2) (x)| ≤ C

M 

ε −1 B (x), i = 1, . . . , M.

(9)

=1

Proof First, we obtain bounds on the regular part v and its derivatives. From Lemma 2.2, it follows that v ≤ C. From (6a), note that εi vi (0) = 0, i = 1, . . . , M. Differentiating Lv = f , we have Lv  = f  − A v in , v  (0) = 0.

(10)

Now, use Lemma 2.2 and the bound on v to get v   ≤ C. Finally, use (10) along with bounds on v and v  to get εi vi(2)  ≤ C, i = 1, . . . , M. Next, we obtain bounds on the layer part w and its derivatives. Using w = B M (x)w, we transform (6b) to the problem εi wi + (aii − αεi /ε M )wi = −

M 

aik w k in , wi (0) = (ηi − vi (0)).

(11)

k=1,k=i

Here wi is the ith component of w. Using (3) and (5), we see that aii − αεi /ε M ≥ aii − α > 0 in . Applying Lemma 2.1 to (11), we get   M   aik    (12) wi  −  a − α  .wi  ≤ C, i = 1, . . . , M. ii k=1,k=i

   aij  Let  = ( ij) be an M × M matrix with entries ii = 1 and ij = −  aii −α ,

for i = j. Following arguments in [7], it holds   M   aik     a − α  < 1, k = 1, . . . , M. ii

(13)

k=1,k=i

Thus  is an inverse monotone matrix. Hence, from (12), we get w ≤ C. Consequently |wi (x)| ≤ CB M (x), x ∈ , i = 1, . . . , M.

(14)

190

Numer Algor (2012) 59:185–195

Using (6b) and (14), we get |wi (x)| ≤ Cεi−1 B M (x), x ∈ , i = 1, . . . , M.

(15)

We need to sharpen the bound (15) for the first M − 1 components. With this in mind, we consider first M − 1 equations of the system (6b) w = z in , w w :=  L E w  + A i (0) = ηi − vi (0), i = 1, . . . , M − 1,  are matrices obtained by deleting the last row and the last where  E and A  are zi = −aiM w M , column from E and A respectively. The components of z, w w i = wi , i = 1, . . . , M − 1. Similar to the decomposition (6), we construct a  = r + d, where r and d are the solution to decomposition w  −1 z(0)  = z in , r(0) = A(0) Lr

(16)

 = 0 in , di (0) = ηi (0) − vi (0) − ri (0), i = 1 . . . , M − 1. Ld

(17)

and

Observe that d is the solution of a problem similar to (6b), but with ε M eliminated from the system. Thus, previous arguments can be used to obtain |di(k) (x)| ≤ Cεi−k B M−1 (x), x ∈ , i = 1, . . . , M − 1, k = 0, 1.

(18)

 − d along with bounds (14) and (18) to get |ri (x)| ≤ Now, use r = w  =z CB M (x), x ∈ , i = 1, . . . , M − 1. In order to bound r  , differentiate Lr  r in , r  (0) = 0.   = z − A Lr

(19)

 r| ≤ Cε−1 B M (x). Now, Using estimates for w M , w M and r, we see that |z − A M −1  using the transformation r = ε M B M (x)r and imitating the arguments that we used to prove (14), we get r ≤ C. Hence, |ri (x)| ≤ Cε−1 M B M (x), x ∈ , i = 1, . . . , M − 1. Collecting the estimates for r and d, we have |wi(k) (x)| ≤ C(εi−k B M−1 (x) + ε−k M B M (x)), x ∈ , i = 1, . . . , M − 1. Now, the desired estimate (8) follows from an induction argument. Finally, by differentiating (6b) and using (8), we obtain (9). 

3 The discrete problem N

Let  := {x j} Nj=0 be a general discretization mesh with points 0 = x0 < x1 < · · · < x N = 1, such that x j = x j−1 + h j, j = 1, . . . , N. The set of interior mesh points is  N := {x j} Nj=1 . For a mesh function Y, define the backward difference operator: (D− Y) j =

Y j − Y j−1 , j = 1, . . . , N. hj

Numer Algor (2012) 59:185–195

191

Letting U ∈ (R N+1 ) M , we define the backward Euler finite difference scheme as follows: (L N U) j = f j,

j = 1, . . . , N, U 0 = η,

(20)



T

T N where f j = f1, j, . . . , f M, j , and (L N U) j = (L1N U) j, . . . , (L M U) j with (LiN U) j := εi (D− U i ) j +

M 

aik, jU k, j.

k=1

Our next aim is to prove the maximum-norm stability of the discrete operator L N . For this purpose, we need the following discrete version of Lemma 2.1. Lemma 3.1 Let Y ∈ R N+1 satisfy μ(D− Y) j + b jY j = f j, j = 1, . . . , N. Then Y N

   f  ≤ max  b 

 N

, |Y0 | .

Proof Using arguments similar to Lemma 2.1, the proof follows.



Lemma 3.2 (Stability of discrete problem) Suppose that the coupling matrix A satisf ies (2). Then the solution U of (20) satisf ies U i  N ≤

M  k=1

  fk −1 (ϒ)ik max  a

kk

   

N

 , |ηk | , i = 1, . . . , M.

Proof Using an idea similar to Lemma 2.2 and appealing to Lemma 3.1, the proof follows. 

An immediate consequence of the above lemma is the uniqueness of the solution to the discrete problem (20).

4 Convergence analysis In this section we first derive a general convergence result for the finite difference scheme (20). This general result we then use to establish the robust convergence of the finite difference scheme (20) on two special layer-resolving meshes.

192

Numer Algor (2012) 59:185–195

Theorem 4.1 Let u be the exact solution of (1) and that U the solution of the f inite difference scheme (20). Then ||u − U|| N ≤ Cϑ( N ),

 xj M  −1 N where ϑ( ) := max ε B (s) ds. 1+ j=1,...,N

x j−1

(21)

=1

Proof Let Z = u − U denote the error of the finite difference scheme (20). Introducing the operator (LiN Y) j := εi (D− Y) j + aii, jY j, we split the error into two parts Z = p + q with (LiN pi ) j = (LiN Z) j, j = 1, . . . , N, pi,0 = 0, and (LiN qi ) j = −

M 

aik, j Z k, j, j = 1, . . . , N, qi,0 = 0.

k=1,k=i

For each i, application of the triangle inequality and Lemma 3.2 gives Z i  N ≤  pi  N + qi  N ≤  pi  N +

  M   aik    Z k  N a  

k=1,k=i

ii

Consequently Z i 

N

  M   aik    Z k  N ≤  pi  N , i = 1, . . . , M. − a    ii

(22)

k=1,k=i

Invoking the inverse monotonicity of ϒ, we get ui − U i  N ≤ C pi  N , i = 1, . . . , M.

(23)

By Lemma 3.1, we have || pi || N ≤ C max εi |(D− ui − ui ) j|, i = 1, . . . , M, j=1,...,N

since (LiN Z) j = (LiN u) j − (LiN U) j = (LiN u) j − fi, j = εi (D− ui − ui ) j. The Taylor expansion with the integral form of the remainder yield  xj −  εi |(D ui − ui ) j| ≤ εi |ui (s)|ds, i = 1, . . . , M. x j−1

Then, by Theorem 2.4, we have −

εi |(D ui −

ui ) j|

 ≤C

xj x j−1

1+

M  =1

ε −1 B (s)

ds, i = 1, . . . , M.

Numer Algor (2012) 59:185–195

193

Hence || pi || N ≤ Cϑ( N ), i = 1, . . . , M. Combining (23) and (24), we get the desired result.

(24) 

The general convergence result of Theorem 4.1 gives the key to construct a number of special layer-resolving meshes for the discretization of (1), but we restrict ourselve to standard Bakhvalov and Shishkin meshes. Bakhvalov meshes for the discretization of (1) can be constructed by equidistributing the function   κ1 κ M −αs/τ ε M e M Ba (s) := max 1, e−αs/τ ε1 , . . . , ε1 εM with positive user chosen constants κi and τ ; that is, the mesh points x j are chosen such that  xj  j 1 M Ba (s)ds = M Ba (s)ds N 0 0 For τ ≥ 1, we have ϑ( N ) ≤ CN −1 , cf. Linss [6]. Using Theorem 4.1, we deduce that ||u − U|| N ≤ CN −1 . Shishkin meshes for the discretization of (1) are constructed as follows. Let τ > 0 be the mesh parameter. Setting σ0 = 0 and σ M+1 = 1, the transition points σ are defined by   σ +1 τ ε σ = min , ln N , = M, . . . , 1. +1 α Assuming that N is divisible by M + 1, the mesh is constructed by dividing each of the intervals [σ , σ +1 ], = 0, . . . , M, into N/(M + 1) equal subintervals. For the mesh constructed in this way, we have ϑ( N ) ≤ C(N −τ + N −1 ln N), cf. Linss [6]. Then, Theorem 4.1 yields ||u − U|| N ≤ CN −1 ln N, if τ ≥ 1. 5 Numerical results In this section numerical results are presented which confirm the theoretical order of convergence established in the previous section. Two test examples are considered in which the operator L and its discretization may or may not satisfy maximum principles.

194

Numer Algor (2012) 59:185–195

Example 5.1 ε1 u1 + (2 + x)u1 − (1 + x/2)u2 = 5x + 1/2, u1 (0) = 1, ε2 u2 − (1 + x)u1 + (2 + x)u2 = xex , u2 (0) = 1. Example 5.2 ε1 u1 + 4u1 + u2 + u3 = x,

u1 (0) = 0,

ε2 u2 − u1 + (4 + x)u2 + u3 = 1,

u2 (0) = 0,

ε3 u3 + 2u1 − u2 + (5 + x)u3 = 1 + x2 ,

u3 (0) = 0.

For the numerical tests we take α = .99, τ = 1, and κi = τ/α, in the definition of the meshes. As the solutions of Examples 5.1 and 5.2 are not known, we estimate errors in the computed solution and rates of convergence using a variant of the double mesh principle (see [2]). To this end, we compute N not only approximate solution U N on mesh  , but also another approximate  2N on the mesh that contains the mesh points of  N and their solution U midpoints. We then compute     2N  EεN := U N − U  N and E N := max EεN , 

Sε

where the values of small parameters ε1 and ε2 in Example 5.1 belong to the set Sε = {(ε1 , ε2 ) | ε2 = 20 , 2−2 , . . . , 2−30 , ε1 = ε2 , 2−2 ε2 , . . . , 2−40 }, and of small parameters ε1 , ε2 and ε3 in Example 5.2 belong to the set  Sε = (ε1 , ε2 , ε3 ) | ε3 = 20 , 2−2 , . . . , 2−30 , ε2 = ε3 , 2−2 ε3 , . . . , 2−40 ,  ε1 = ε2 , 2−2 ε2 , . . . , 2−60 .

Table 1 Numerical results for Example 5.1

N 3 × 26

3 × 27 3 × 28 3 × 29 3 × 210 3 × 211

Shishkin mesh

Bakhvalov mesh

EN

N

EN

N

8.85E-03 5.17E-03 2.96E-03 1.65E-03 9.06E-04 4.92E-04

0.78 0.80 0.84 0.86 0.88 –

2.74E-03 1.39E-03 6.98E-04 3.50E-04 1.75E-04 8.72E-05

0.98 0.99 1.00 1.00 1.00 –

Numer Algor (2012) 59:185–195 Table 2 Numerical results for Example 5.2

195 N

Shishkin mesh

4 × 26

4 × 27 4 × 28 4 × 29 4 × 210 4 × 211

Bakhvalov mesh

EN

N

EN

N

5.68E-03 3.40E-03 1.96E-03 1.09E-03 6.02E-04 3.26E-04

0.74 0.79 0.85 0.86 0.89 –

1.58E-03 8.06E-04 4.09E-04 2.06E-04 1.03E-04 5.14E-05

0.97 0.98 0.99 1.00 1.00 –

From these values, we calculate the parameter-robust rates of convergence using the formula ln(E N /E2N ) . ln 2 Tables 1 and 2 represent the results of the numerical experiments for Examples 5.1 and 5.2, respectively. From these tables we see that the numerical results are well in accordance with the theoretical ones. N =

Acknowledgements The authors gratefully acknowledge the valuable comments and suggestions from the anonymous referees.

References 1. Athanasios, A.C.: Approximation of Large-Scale Dynamical Systems. SIAM, Philadelphia (2005) 2. Farrell, P.A., Hegarty, A.F., Miller, J.J.H., O’Riordan, E., Shishkin, G.I.: Robust Computational Techniques for Boundary Layers. Chapman & Hall/CRC, Boca Raton (2000) 3. Gaji´c, Z., Lim, M.-T.: Optimal Control of Singularly Perturbed Linear Systems and Applications. Marcel Dekker, New York (2001) 4. Hemavathi, S., Bhuvaneswari, T., Valarmathi, S., Miller, J.: A parameter uniform numerical method for a system of singularly perturbed ordinary differential equations. Appl. Math. Comput. 191, 1–11 (2007) 5. Ladyzhenskaya, O.A., Ural’tseva, N.N.: Linear and Quasilinear Elliptic Equations. Academic Press, New York (1968) 6. Linss, T.: Layer-adapted meshes for reaction-convection-diffusion problems. In: Lecture Notes in Mathematics, vol. 1985. Springer, Berlin (2010) 7. Linss, T., Madden, N.: Layer-adapted meshes for a linear system of coupled singularly perturbed reaction-diffusion problems. IMA J. Numer. Anal. 29, 109–125 (2009) 8. Meenakshi, P.M., Valarmathi, S., Miller, J.: Solving a partially singularly perturbed initial value problem on shishkin meshes. Appl. Math. Comput. 215, 3170–3180 (2010) 9. Protter, M.H., Weinberger, H.F.: Maximum Principles in Differential Equations. PrenticeHall, Englewood Cliffs (1967) 10. Roos, H.-G., Stynes, M., Tobiska, L.: Robust numerical methods for singularly perturbed differential equations. In: Springer Series in Computational Mathematics, 2nd edn. Springer, Berlin (2008) 11. Valarmathi, S., Miller, J.: A parameter-uniform finite difference method for singularly perturbed linear dynamical systems. Int. J. Numer. Anal. Model. 7, 535–548 (2010)