Uniformly convergent finite volume difference scheme for 2D convection-dominated problem with discontinuous coefficients

Applied Mathematics and Computation 163 (2005) 645–665 www.elsevier.com/locate/amc

Uniformly convergent finite volume difference scheme for 2D convection-dominated problem with discontinuous coefficients Iliya A. Brayanov Center of Applied Mathematics and Informatics, University of Rousse, Studentska str. 8, 7017 Rousse, Bulgaria

Abstract Two dimensional singularly perturbed convection–diffusion problem with discontinuous coefficients is considered. The problem is discretized using an inverse-monotone finite volume method on Shishkin meshes. We established first-order global pointwise convergence that is uniform with respect to the perturbation parameter. Numerical experiments that support the theoretical results are given. Ó 2004 Elsevier Inc. All rights reserved. Keywords: Convection–diffusion problems; Singular perturbation; Asymptotic analysis; Finite volume methods; Modified upwind approximations; Uniform convergence; Shishkin mesh

1. Introduction Let us consider the following model convection–diffusion problem: Lu ¼ eDu þ a  ru þ bu ¼ f ;

ðx; yÞ 2 X [ Xþ ;

½u C ¼ uð0:5 þ 0; yÞ  uð0:5  0; yÞ ¼ 0;

ujoX ¼ g;


 ou  e ¼ QðyÞ; ox C

E-mail address: [email protected] (I.A. Brayanov). 0096-3003/$ - see front matter Ó 2004 Elsevier Inc. All rights reserved. doi:10.1016/j.amc.2004.04.007

ð1Þ

on C;

ð2Þ

646

I.A. Brayanov / Appl. Math. Comput. 163 (2005) 645–665

where 0 < e 1, X ¼ ð0; 0:5Þ  ð0; 1Þ, Xþ ¼ ð0:5; 1Þ  ð0; 1Þ, X ¼ X [ Xþ , C ¼ fx ¼ 0:5; y 2 ð0; 1Þg, a ¼ ða1 ðx; yÞ; a2 ðx; yÞÞ P ða1 ; a2 Þ > ð0; 0Þ and bðx; yÞ P b > 0. The functions a, b and f could be discontinuous on the segment C. Note that the sign pattern of the convection coefficients is essential for the behavior of the solution, see [7] for more detailed discussion in one dimensional case. The reduced problem is first-order discontinuous hyperbolic problem. As a result boundary and interior layers appear near the outer boundaries of X and Xþ : fx ¼ 0:5; 0 6 y 6 1g, f0 6 x 6 0:5; y ¼ 1g and fx ¼ 1; 0 6 y 6 1g, f0:5 6 x 6 1; y ¼ 1g. This problem may be viewed as modeling a steady-state convection–diffusion process. The importance of accurate numerical methods lies in the fact that it can also be regarded as a simple linear model of the Navier–Stokes flow equations. Boundary and interior layers are normally present in the solutions of problems involving such equations. These layers are thin regions in the domain where the gradient of the solution steepens as the singular perturbation parameter e tends to zero. There is a vast literature dealing with convection–diffusion problems of type (1) with smooth coefficients and smooth right-hand side––see [12,13,16] for survey. The one dimensional problem with discontinuous input data is discussed recently by several authors see [4,7,9,14] and references there. But for two dimensional problems with discontinuous data only handful of papers are available. One of the reason is the complexity of this problem in two dimensions, where boundary, corner and interior layers appear and it’s analytical investigation is still not well done. Singularly perturbed two dimensional reaction–diffusion problem with anisotropic coefficients is considered in [5]. Second order uniformly convergent scheme on Shishkin meshes is derived there. Two dimensional convection–diffusion problem with discontinuous coefficients is considered in [6] where the monotone schemes from [8] are developed further on Shishkin meshes. It was shown there that the scheme are uniformly convergent with respect to the perturbation parameter. Our objective in this paper to derive global pointwise convergence that is uniform with respect to the small parameter e. We construct partially uniform Shishkin meshes in each subdomain X and Xþ . The resulting global mesh is condensed closely to the boundary and interior layers. To derive uniform difference scheme we use an inverse-monotone finite volume discretization on layer adapted meshes. This scheme was introduced by Samarskii [15] and derived later in [8] as a finite volume difference scheme for two dimensional convection–diffusion problem. The key of our analysis is a hybrid stability inequality that shows that the maximal nodal error is bounded by a discrete L1 norm of its truncation error. A result like this is proved by Andreev in [1] and used later in [2,10] to prove uniform convergence in maximum norm.

I.A. Brayanov / Appl. Math. Comput. 163 (2005) 645–665

647

An outline of the paper is as follows. In Section 2 we describe some properties of the differential solution, construct Shishkin decomposition of the solution and derive e-uniform bounds for their derivatives. In Section 3 we construct Shishkin mesh condensed near the boundary and interface layers and obtain finite volume difference scheme. Some a priori bounds of the discrete problem are given in Section 4. The uniform convergence of the constructed scheme is proved is Section 5. In Section 6 we give some numerical experiments that support the theoretical results.

2. Properties of analytical solution In this section we describe some properties of the analytical solution. 2.1. Maximum principle We start with Lemma 1 (Maximum principle). Suppose that a function vðx; yÞ 2 C 2 ðX Þ\ C 2 ðXþ Þ \ CðXÞ, satisfies 8x 2 X [ Xþ ;

vjoX P 0;

Lv P 0;

½v C ¼ 0;

e½v0x C P 0;

y 2 ð0; 1Þ:

Then vðxÞ P 0 for all x 2 X. Proof. Let P0 ðx0 ; y0 Þ be any point at which v attains its minimum value in X. If vðP0 Þ P 0, there is nothing to prove. Suppose therefore that vðP0 Þ < 0. With the above assumption on the boundary values either P0 2 X [ Xþ or P0 2 C. In the first case it is well known fact from the real analysis that v00xx ðP0 Þ P 0, v00yy ðP0 Þ P 0 and v0x ðP0 Þ ¼ v0y ðP0 Þ ¼ 0. Therefore LvðP0 Þ ¼ eðv00xx ðP0 Þ þ v00yy ðP0 ÞÞ þ a1 v00x ðP0 Þ þ a2 v00y ðP0 Þ þ qvðP0 Þ < 0: The last inequality is in contradiction with the assumption Lv P 0. In the second case the argument depends on whether, or not v is differentiable at P0 . If v0x ð0:5 þ 0; y0 Þ 6¼ v0x ð0:5  0; y0 Þ then at least one of these one side derivatives is not zero. Therefore, ½v0x C ðy0 Þ 6¼ 0 and v0x ð0:5 þ 0; y0 Þ P 0, v0x ð0:5  0; y0 Þ 6 0 and it is clear that ½v0x C ðy0 Þ > 0, which is a contradiction. It remains to consider the case v0x ð0:5 þ 0; y0 Þ ¼ v0x ð0:5  0; y0 Þ ¼ 0: Note that there exists 0 < k3 < k such that x3 ¼ 0:5 þ k3 s. Then, by the above considerations

648

I.A. Brayanov / Appl. Math. Comput. 163 (2005) 645–665

v00xx ðx3 ; y0 Þ ¼ v00xx ð0:5 þ k3 s; y0 Þ ¼ u00 ðs; k3 ; 0Þ > 0: Also, v00yy ðx3 ; y0 Þ P 0 and vðx3 ; y0 Þ < 0. Thus Lvðx3 ; y0 Þ < 0;

ð3Þ

which is the required contradiction. We first prove that v00yy P 0 in a neighborhood of P0 on the right of the line x ¼ 0:5. For this we consider the function uðtÞ ¼ vð0:5 þ kt; y0 þ ltÞ, k > 0, jlj 6 1, k, l––real parameters. The function uðtÞ has strong local minimum, i.e. uðtÞ > uð0Þ ¼ vðP0 Þ, uðt; k; lÞ ¼ u0 ð0Þ ¼ 0 for sufficiently small t. By Taylor expansion we find u00 ðhtÞ ¼

2 ðuðtÞ  uð0ÞÞ > 0; t2

ð4Þ

0
for sufficiently small t > 0. Let us fix such a t and the corresponding h. Denoting s ¼ ht, we have u00 ðsÞ ¼ v00xx ð0:5 þ ks; y0 þ lsÞk2 þ 2v00xy ð0:5 þ ks; y0 þ lsÞkl þ v00yy ð0:5 þ ks; y0 þ lsÞl2 : þ

Since v00xx , v00yy , v00xy are bounded on X , we deduce from (4) that there exists k > 0, such that v00yy ð0:5 þ ks; y0 þ lsÞ > 0 for all 0 < k < k ; jlj 6 1: Hence v00yy ð0:5 þ ks; y0 Þ P 0;

0 < k < k :

Now, choose a point P1 ðx1 ; y0 Þ such that 0 < x1 < ks. Since vðP1 Þ > vðP0 Þ we conclude from the mean values theorem that, for some 0 < x2 < x1 v0x ðx2 ; y0 Þ ¼

vðP0 Þ  vðP1 Þ > 0: x1

Also, for some 0 < x3 < x2 v00xx ðx3 ; y0 Þ ¼

v0x ð0:5; y0 Þ  v0x ðx2 ; y0 Þ > 0: x2



Using the maximum principle we obtain that the solution of (1), (2) is uniformly bounded with respect to the small parameter e. Remark 1. Note that the sign pattern of a1 ðx; yÞ is essential in the proof. If a1 ðx; yÞ changes it’s sign and a1 ðx; yÞ > 0 in X , a1 ðx; yÞ < 0 in Xþ , then the proof still holds, but if a1 ðx; yÞ < 0 in X and a1 ðx; yÞ > 0 in Xþ then (3) is not

I.A. Brayanov / Appl. Math. Comput. 163 (2005) 645–665

649

satisfied. In this case for one dimensional problem in [7] was shown that the solution is not e-uniformly bounded. 2.2. Asymptotic analysis Further in the numerical analysis we shall need a decomposition of the solution into regular and singular part. We describe it briefly here in the case of constant coefficients in each subdomain X . In the construction we shall follow the ideas in [11,17]. For more detailed discussion concerning the discontinuous problem we refer [11]. The regular part V is presented as a sum V ¼ U þ Rnþ1 , where  U ðx; yÞ ¼

if ðx; yÞ 2 X ; if ðx; yÞ 2 Xþ ;

U  ðx; yÞ U þ ðx; yÞ

U ¼

n X

ei Ui ;

i¼0

and Ui are solutions of the reduced problems a:rUi ðx; yÞ þ bUi ðx; yÞ ¼ Fi ðx; yÞ;

ðx; yÞ 2 X ;

F0 ðx; yÞ ¼ f ðx; yÞ;

 Fi ðx; yÞ ¼ DUi1 ðx; yÞ;

U0 ð0; yÞ ¼ gð0; yÞ;

U0 ðx; 0Þ ¼ gðx; 0Þ;

Ui ð0; yÞ ¼ Ui ðx; 0Þ ¼ 0;

i P 1;

Uiþ ð0:5; yÞ ¼ Ui ð0:5; yÞ þ Eix ð0; yÞ; U0þ ðx; 0Þ ¼ gðx; 0Þ;

i P 1;

Uiþ ðx; 0Þ ¼ 0;

i P 0; i P 1:

The singular part W is expressed in the form  x  E ðg ; yÞ þ Ey ðx; h Þ þ Exy ðg ; h Þ; W ðx; yÞ ¼ Exþ ðgþ ; yÞ þ Eyþ ðx; hþ Þ þ Exyþ ðgþ ; hþ Þ;

in X ; in Xþ ;

ð5Þ

where Ex ¼

n X

Exy ¼

ei Eix ðg ; yÞ;

i¼0 n X

Ey ¼

n X

ei Eiy ðx; h Þ;

i¼0

ei Eixy ðg ; h Þ

i¼0

and g ¼ ð0:5  xÞ=e, gþ ¼ ð1  xÞ=e, h ¼ ð1  yÞ=e are the stretch variables. The boundary layer terms Eiy are solutions of the problems

650

I.A. Brayanov / Appl. Math. Comput. 163 (2005) 645–665



o2 Eiy ðx; h Þ 2

oðh Þ

 a2

oEiy ðx; h Þ ¼ Fi ðx; h Þ; oh

h 2 ð0; 1Þ;

oE0y ðx; h Þ  bE0y ðx; h Þ; ox y y o2 Ei2 ðx; h Þ oEi1 ðx; h Þ y  bEi1 Fi ðx; h Þ ¼  a ðx; h Þ; 1 2 ox ox Ei ðx; 0Þ ¼ W Ei ðx; 1Þ ¼ 0; i ðxÞ; F1 ðx; h Þ ¼ a1

F0 ¼ 0;

 W 0 ðxÞ ¼ gðx; 1Þ  U0 ðx; 1Þ;

 W i ðxÞ ¼ Ui ðx; 1Þ;

i P 2;

i P 1:

Analogously the boundary layer terms Eix are solutions of the problems 

o2 Eix ðg ; yÞ oðg Þ

2

 a1

oEix ðg ; yÞ ¼ Fi ðg ; yÞ; og

g 2 ð0; 1Þ;

oE0x ðg ; yÞ  bE0x ðg ; yÞ; oy x o2 Ei2 ðg ; yÞ oEx ðg ; yÞ x  bEi1  a2 i1 ðg ; yÞ; i P 2; Fi ¼ 2 oy oy Ei ð1; yÞ ¼ 0; Ei ð0; yÞ ¼ Ui ðyÞ; i P 0; þ þ þ E0 ð0; yÞ ¼ gð1; yÞ  U0 ð1; yÞ; Ei ð0; yÞ ¼ Uiþ ð1; yÞ; i P 1:

F0 ¼ 0;

F1 ðg ; yÞ ¼ a2

The corner layer functions Eixy ðg ; h Þ satisfy  DEixy ðg ; h Þ þ a  rEixy ðg ; h Þ ¼ Fi ðg ; h Þ; F01 ¼ 0;

xy  Fi ðg ; h Þ ¼ bEi1 ðg ; h Þ;

Eixy ð0; h Þ Eixy ðg ; 0Þ

¼

Eiy ð0:5; h Þ;

¼ Eix ðg ; 1Þ;

i P 1;

Eixyþ ð0; hþ Þ Eixy ! 0;

¼ Eiyþ ð1; hþ Þ;

g ; h ! 1:

The functions Ui ðyÞ are chosen so that the function Eix and Eixy to satisfy the conditions oE0x ð0; yÞ oE0xy ð0; h Þ þ ¼ QðyÞ; og og þ  oEix ð0; yÞ oEixy ð0; h Þ oUi1 ð0:5; yÞ oUi1 ð0:5; yÞ  þ ¼ og og ox ox y oEyþ ð0:5; hþ Þ oEi1 ð0:5; h Þ  ; i P 1: þ i1 ox ox The boundary and corner layer function decrease exponentially to zero away from the layers, see [11,17]. Thus we can consider them as multiplied by a cutoff function X , W ¼ XW that makes them zero away from the layers. Then we obtain the following problem for the remainder term Rn :

I.A. Brayanov / Appl. Math. Comput. 163 (2005) 645–665

651

LRn ðx; yÞ ¼ en Fn ðx; yÞ ¼ Oðen Þ; ðx; yÞ 2 X ;   x o2 Enx ðg ; yÞ o2 Eny ðx; h Þ o2 En1 ðg ; yÞ   Fn ðx; yÞ ¼ e DUn þ þ þ oy 2 ox2 oy 2 þ

½Rn C ¼ 0; Wn ðyÞ ¼

y o2 En1 ðx; h Þ oEny ðx; h Þ oEx ðg ; yÞ  a2 n  a 1 2 ox ox oy

 bEnx ðg ; yÞ  bEny ðx; h Þ  bEnxy ðg ; h Þ;
 oRn e ¼ enþ1 Wn ðyÞ ¼ Oðenþ1 Þ; Rn joX ¼ 0: ox C

oUnþ ð0:5; yÞ oUn ð0:5; yÞ oEnyþ ð0:5; hþ Þ oEny ð0:5; h Þ  þ  : ox ox ox ox

Thus we can make the following proposition. Proposition 1. Let the coefficients of problem (1) and (2) are sufficiently smooth and satisfy all necessary compatibility conditions, that guarantee u 2 CðXÞ \ C 4 ðX [ Xþ Þ. Then the solution u can be decomposed into regular part V that satisfies for all k; l-integer, k; l ¼ 0; . . . ; 4; k þ l 6 4 the estimates  kþl   o V ðx; yÞ   þ   ð6Þ  oxk oy l  6 C; ðx; yÞ 2 X [ X and singular part W given by (5) that satisfies  kþl x   o E ðx; yÞ  k    oxk oy l  6 Ce expðm1 ð0:5  xÞ=eÞ;  kþl y   o E ðx; yÞ  l    oxk oy l  6 Ce expðm2 ð1  yÞ=eÞ;   kþl xy  o E ðx; yÞ   6 Cekl expððm1 ð0:5  xÞ þ m2 ð1  yÞÞ=eÞ;    oxk oy l  kþl xþ   o E ðx; yÞ  k    oxk oy l  6 Ce expðm3 ð1  xÞ=eÞ;  kþl yþ   o E ðx; yÞ  l    oxk oy l  6 Ce expðm4 ð1  yÞ=eÞ;  kþl xyþ   o E ðx; yÞ    6 Cekl expððm3 ð1  xÞ þ m4 ð1  yÞÞ=eÞ;   oxk oy l where C, m1 , m2 , m3 , m4 are independent of e positive constants.

ð7Þ

652

I.A. Brayanov / Appl. Math. Comput. 163 (2005) 645–665

3. Numerical solution 3.1. Grid and grid functions It is well known that in singularly perturbed problems cannot be achieved an uniform convergence on uniform meshes. In order to obtain e-uniformly convergent difference scheme we construct a piecewise uniform (Shishkin) mesh  condensed near the boundary and the interior layers, see Fig. 1a. w  ¼ fðxi ; yj Þ; xi ¼ xi1 þ hxi ; yj ¼ yj1 þ hyj ; i ¼ 1; . . . ; 2N ; j ¼ 1; . . . ; M; x x0 ¼ 0; xN ¼ 0:5; x2N ¼ 1; y0 ¼ 0; yM1 ¼ 1g; where 8 2ð0:5dx Þ h1 ¼ N1 1 ; i ¼ 1; . . . ; N2 ; > > > > > < h ¼ 2dx1 ; i ¼ N2 þ 1; . . . ; N ; 2 x N hi ¼ x > h ¼ 2ð0:5d2 Þ ; i ¼ N þ 1; . . . ; 3N ; > 3 > 2 N > > : 2dx2 3N h4 ¼ N ; i ¼ 2 þ 1; . . . ; 2N ; ( 2ð1dy Þ k1 ¼ M ; j ¼ 1; . . . ; M2 ; hyj ¼ y k2 ¼ 2dM ; j ¼ M2 þ 1; . . . ; M; dx1 ¼ minf2e ln N =m1 ; 1=4g;

dx2 ¼ minf2e ln N =m3 ; 1=4g;

dy ¼ minfmaxf2e ln M=m2 ; 2e ln M=m4 g; 1=2g:



Fig. 1. Grid with local refinement on boundary and interior layers. (a) Shishkin mesh in X, –– boundary points, d––interior points in x and xþ , ––interface points and (b) control volume ei;j .



I.A. Brayanov / Appl. Math. Comput. 163 (2005) 645–665

653

For each point ðxi ; yj Þ of w we consider a control volume the rectangle ei;j ¼ eðxi ; yj Þ, see Fig. 1. There are three types of grid points, boundary, interior and interface points.  Let v, w, g are given grid functions of a discrete arguments ðxi ; yj Þ 2 x. Denote gi;j ¼ gðxi ; yj Þ, gN 0;j ¼ gðxN  0; yj Þ. Further we shall use the standard notations hx þ hxiþ1 hx hx  ; i ¼ 1; . . . ; 2N  1;  hx0 ¼ 1 ; hx2N ¼ 2N ; hxi ¼ i 2 2 2 hxN þ1 gN þ0;j þ hxN gN 0;j vi;j  vi1;j vi;j  vi1;j gN ;j ¼ ; vx;i;j ¼ ; vx;i;j ¼ ; hxi hxi 2 hxN vx;i;j  vx;i;j vx;i;j ¼ vx;iþ1;j ; v^x;i;j ¼ vx;iþ1;j ; vx^x;i;j ¼ ; hxi ð8Þ and similarly ones for y. We shall use the following discrete scalar products: ðv; wÞ0;x ¼

2N 1 M1 X X i¼1

ðv; w ;x;x ¼

 hxi  hyj vi;j wi;j ;

j¼1

2N M 1 X X i¼1

hxi  hyj vi;j wi;j ;

ðv; w ;y;x ¼

j¼1

2N 1 X

M X

i¼1

j¼1

hx hyj vi;j wi;j ; i

and the norms kvk0;x ¼ kvk1;x ¼

qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi ðv; vÞ0;x ; N 1 M i 1 X X i¼1

kv j;x;x ¼

 hxi  hyj jvi;j j;

j¼1

qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi ðv; v ;x;x ;

kv j;y;x ¼

qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi ðv; v ;y;x ;

kvk1;x ¼ max jvi;j j: ðxi ;yj Þ2x

3.2. Finite difference approximation The finite difference approximation is derived from the balance equation. Integrating the Eqs. (1) over cell ei;j , that does not intersect the interface C, using the Green’s formula and dividing by  hxi hyj we get Z Z 1 1 W ds ¼ ðf ðx; yÞ  aðx; yÞ  ru  bðx; yÞuðx; yÞÞ dx dy;   hyj oei;j hyj ei;j hxi  hxi  ð9Þ where W ¼ eruðsÞ  n;

654

I.A. Brayanov / Appl. Math. Comput. 163 (2005) 645–665

and n is the unit normal vector to the boundary of eij . After approximation of the integrals in (9), we obtain
 oui1=2;j e ouiþ1=2;j ð1  R1;iþ1;j Þ  ð1 þ R1;i;j Þ  x ox ox hi
 oui;j1=2 e oui;jþ1=2  y ð1  R2;i;jþ1 Þ  ð1 þ R2;i;j Þ þ bi;j ui;j ¼ fi;j ; ð10Þ oy oy hj where q1;i;j ¼ a1 ðxi  hxi =2; yj Þ; R1;i;j ¼ hxi q1;i;j =2e;

q2;i;j ¼ a2 ðxi ; yj  hyj =2Þ;

R2;i;j ¼ hyj q2;i;j =2e:

Now we approximate ouiþ1=2;j  Ux;i;j ; ox

oui1=2;j  Ux;i;j ; ox

oui;jþ1=2  Uy;i;j ; oy

oui;j1=2  Uy ;i;j : oy

To approximate the terms 1  Rm , m ¼ 1; 2 we use that 1  Rm ¼

1 R2m  ; 1 þ Rm 1 þ Rm

1 þ Rm ¼

1 R2m þ 2Rm  ; 1 þ Rm 1 þ Rm

and ignore the last addend R2m =ð1 þ Rm Þ. Thus using (10) at the interior points in w that do not lay on the interface, we obtain the so-called Samarskii’s or modified upwind difference scheme (MUDS), see [8,15]. Lh Ui;j ¼ eðj1 Ux Þ^x;i;j  eðj2 Uy Þ^y ;i;j þ q1;i;j Ux;i;j þ q2;i;j Uy ;i;j þ bi;j Ui;j ¼ fi;j ; ð11Þ where jm;i;j ¼ ð1 þ Rm;i;j Þ1 . þ Let now the rectangle ei;j intersects the interface C. Denote by e i;j , ei;j the parts of ei;j that lay respectively on the left and on the right of the interface. Denote by SC;i;j the intersection of ei;j and C. Using the interface conditions (2), we obtain "Z # Z 1 W ds þ W ds yj xi h h oe =S oeþ =S i;j C;i;j i;j C;i;j "Z Z 1 1 ¼ x y QðyÞ dy þ x y ðf  a  ru  buÞ dx dy hi hj SC;i;j hi hj e i;j # Z þ

ðf  a  ru  buÞ dx dy :

eþ i;j

ð12Þ

I.A. Brayanov / Appl. Math. Comput. 163 (2005) 645–665

655

After approximation of the integrals in (12) we obtain
 oui1=2;j e ouiþ1=2;j ð1  R1;iþ1;j Þ  ð1 þ R1;i;j Þ  x ox ox hi
 oui;j1=2 e oui;jþ1=2 Qj  y ð1  R2;i;jþ1 Þ  ð1 þ R2;i;j Þ þ bi;j ui;j ¼ fi;j þ x : oy oy hi hj ð13Þ 

The definition of f ,  b, R2 that appears here is given in (8). Using (13) similarly as above, we get the following approximation at the interface points: 2;N ;j Uy ;N ;j þ bN ;j UN ;j Lh UN ;j ¼ eðk1 Ux Þ^x;N ;j  eðk2 Uy Þ^y ;N ;j þ q1;N ;j Ux;N ;j þ q Qj ¼ fN ;j þ x : ð14Þ hN Applying the discrete Green’s formula we obtain the following lemma. Lemma 2. For arbitrary mesh functions v, w on x, satisfying v; wjox ¼ 0 we have ðq1 vx ; wÞ0;x ¼ ðv; ðq1 wÞ^x Þ0;x ;

ððj1 vx Þ^x ; wÞ0;x ¼ ðj1 vx ; wx ;x;x ;

ð15Þ

ðq2 vy ; wÞ0;x ¼ ðv; ðq2 wÞ^y Þ0;x ;

ððj2 vy Þ^y ; wÞ0;x ¼ ðj2 vy ; wy ;y;x :

ð16Þ

Setting the boundary conditions U jox ¼ g;

ð17Þ

we obtain the finite difference problem (FDP) (11), (14) and (17). The FDP can be written as a system of linear algebraic equations Au ¼ F ;

ðxi ; yj Þ 2 x;

ð18Þ

where in the right-hand side F we have taken the boundary conditions (17) into account. Using a similar argument as in [8], we can prove the following lemma. Lemma 3. The FDP (18) satisfies the discrete maximum principle and the corresponding matrix A is an M-matrix.

4. A priori estimates In this section we derive some a priori estimates for the solution of (18). The Green’s function G of the operator Lh associated with the mesh node ðnk ; gm Þ 2 x is defined by

656

I.A. Brayanov / Appl. Math. Comput. 163 (2005) 645–665 h;x h;y Lh Gk;m i;j ¼ di;k dj;m

on x;

G¼0

on ox;

h;y where dh;x i;k and dj;m are the discrete analog of the delta-function:   1= hxi if xi ¼ nk ; 1=hyj if yj ¼ gm ; h;x h;y di;k ¼ dj;m ¼ 0 otherwise; 0 otherwise:

From estimates (15) and (16) in Lemma 2 we deduce that the adjoint operator Lh of Lh is Lh; Ui;j ¼ eðj1 Ux Þ^x;i;j  eðj2 Uy Þ^y ;i;j  ðq1 U Þ^x;i;j  ðq2 U Þ^y ;i;j þ bi;j Ui;j ; i 6¼ N ; Lh; UN ;j ¼ eðj1 Ux Þ  eðj2 Uy Þ  ðq1 U Þ  ð q2 U Þ þ bN ;j UN ;j : ^x;N ;j

^y ;N ;j

^x;N ;j

^y ;N ;j

Then G as a function of the second argument satisfies for any ðxi ; yj Þ 2 x h;x h;y Lh; Gk;m i;j ¼ di;k dj;m

on x;

G ¼ 0 on ox:

ð19Þ

Lemma 4. The discrete Green’s function is nonnegative and satisfies max

M 1 X

xi ;yj ;nk 2x

1 y Gk;l i;j hm 6 a1 ;

max

xi ;yj ;gm 2x

m¼1

N 1 X

1 x Gk;m i;j hk 6 a2 :

ð20Þ

k¼1

Proof. The fact that G is nonnegative immediately follows from the inverse monotonicity of the operator Lh . We shall proof the first estimate in (20). The second one is proved analogously. Define k

Gi;j ¼

M 1 X

y Gk;m i;j hm :

m¼1

We multiple (19) by hym and sum over m ¼ 1; . . . ; M  1. Using the positivity of G and a discrete version of mean value theorem, we know there exist mesh ~1 , q ~1 the function functions ~ a1 P a1 and ~ b P b such that with the associated j k Gi;j solves k

k j1 Gi;j Þn^n;k  ð~ q1 Gi;j Þ^n;k þ ~bk Gi;j ¼ dh;x ½e L h; 1 Gi;j ¼ eð~ i;k  Fi;j ; k

0

2N

k ¼ 1; . . . ; 2N  1; Gi;j ¼ 0; Gi;j ¼ 0; where Fi;jk ¼

! e j2;k;1 e j2;k;M1 k;M1 2;k;1 þ y q Gi;j : Gk;1 i;j þ h1 hyM

ð21Þ

I.A. Brayanov / Appl. Math. Comput. 163 (2005) 645–665

657

2;k;j ¼ j2;k;j and q 2;k;j ¼ q2;k;j for k 6¼ N . Note that the coefficients in front Here j of Gki;j in the expression for Fi;jk are nonnegative, therefore Fi;jk is also nonnegative. e 1 denotes the Green’s function of the operator e Let G L h1 adjoint to the h; e operator L 1 in (21) then the solution of (21) can be written in the form 2N 1 X

k e 1;i;k  Gi;j ¼ G

x e 1;s;k F s  G i;j hs :

s¼1

e 1 and F and the estimates for G e 1 (see [3], in continuous The nonnegativity of G case and [4] in discontinuous case) gives k

e 1;i;k 6 a1 : Gi;j 6 G 1



Lemma 4 yields the following a priori bounds. Lemma 5. Let the discrete function vh satisfies Lh vh ¼ f ¼ f1 þ f2 ;

on x; vh jox ¼ 0

ð22Þ

with arbitrary mesh functions f1 , f2 . Then kvh k1;x 6 a1 1

2N 1 X i¼1

 hxi max jf1;i;j j þ a1 2 j¼1;...;M1

M 1 X

hyj

j¼1

max

i¼1;...;2N 1

jf2;i;j j:

ð23Þ

Proof. The solution of (22) admits the representation vhi;j ¼

2N 1 X M1 X k¼1

 hxk  hyl Gk;l i;j ðf1;k;l þ f2;k;l Þ:

l¼1

Therefore jvhi;j j 6

2N 1 X k¼1

 hxk

max jf1;k;l j

l¼1;...;M1

M 1 X

 hyl Gk;l i;j þ

l¼1

M 1 X l¼1

hy l

max

k¼1;...;2N 1

jf2;k;l j

Now from Lemma 4 we get (23). 

5. Uniform convergence Now using Lemma 4 we can prove the following theorem.

2N 1 X k¼1

hx Gk;l k i;j :

658

I.A. Brayanov / Appl. Math. Comput. 163 (2005) 645–665

Theorem 1. Let u be a solution of the differential problem (1) and (2) that satisfies the hypothesis of Proposition 1. Let U be a solution of the discrete problem  þ (18). Then if u 2 C 3 ðX Þ \ C 3 ðX Þ \ CðXÞ then the following estimate holds: kU  uk1;w 6 CðN 1 ln N þ M 1 ln MÞ: 

ð24Þ

þ

If in addition u 2 C 4 ðX Þ \ C 4 ðX Þ \ CðXÞ for sufficiently large M and N , independent of e holds kU  uk1;w 6 CðN 1 þ M 1 Þ;

ð25Þ

where C is some positive constant independent of the mesh and the small parameter e.

Proof. We shall proof the second estimate (25). The first one (24) is proved analogously. Consider the truncation error wi;j ¼ Lh ðU  uÞi;j ¼ ðLu  Lh uÞi;j . We begin with the internal points. The finite difference scheme (11) can be written in the form   1 1 hx q1 Ux L Ui;j ¼ eUx^x;i;j  eUy^y ;i;j þ ðq1;i;j Ux;i;j þ q1;iþ1;j U^x;i;j Þ  2 2 1 þ R1 1 ^x;i;j   1 1 hy q2 Uy þ ðq2;i;j Uy ;i;j þ q2;i;jþ1 U^y ;i;j Þ  þ bi;j Ui;j ¼ fi;j : 2 2 1 þ R1 2 ^y ;i;j h

Then the truncation error can be written in the form

wi;j ¼ gxi;j þ gyi;j ;

gxi;j ¼

3 X k¼1

gxk;i;j ; gyi;j ¼

3 X

gyk;i;j ;

k¼1

where

jgx1;i;j j

   o2 ui;j   ¼ eux^x;i;j  ox2   3   4    o uðx; yj Þ     þ ðhx Þ2 max  o uðx; yj Þ  ; 6 Ce jhxiþ1  hxi j maxx  i    3 x 4 x2Di x2Di ox ox

I.A. Brayanov / Appl. Math. Comput. 163 (2005) 645–665

659

  1 oui;j  x  jg2;i;j j ¼  ðq1;i;j ux;i;j þ q1;iþ1;j ux;i;j Þ  a1;i;j 2 ox   2   3    o uðx; yj Þ     þ ðhx Þ2 max  o uðx; yj Þ  ; 6 C jhxiþ1  hxi j maxx  i    2 x 3 x2Di x2Di ox ox  ( )  2     x x 2   o uðx; yj Þ  h q ð h Þ   1 x x i   : jg3;i;j j ¼  u maxx   6 C min hi ;   1 þ R1 x x2D ox2  e 1

^x;i;j

i

   o2 ui;j  jgy1;i;j j ¼ euy^y ;i;j  oy 2  6 Ce

jgy2;i;j j



 

o hyj j maxy  y2Dj

3

  4 !  o uðxi ; yÞ  uðxi ; yÞ  2 y  ; þ ðhj Þ maxy  oy 3  oy 4  y2Dj

ð26Þ

  1 oui;j   ¼  ðq2;i;j uy ;i;j þ q2;i;jþ1 vy;i;j Þ  a2;i;j 2 oy  6C

jgy3;i;j j

jhyjþ1

jhyjþ1



 

o hyj j maxy  y2Dj

2

  3 !  o uðxi ; yÞ  uðxi ; yÞ  y 2   ;  þ ðhj Þ maxy  oy 2  oy 3  y2Dj

 ( )  2     y 2  o uðxi ; yÞ  hy q2   y ð hj Þ  :  ¼ uy maxy   6 C min hj ;   1 þ R1 oy 2  e y2Dj 2 ^y ;i;j

Here Dxi ¼ ½xi1 ; xiþ1 , Dyj ¼ ½yj1 ; yjþ1 . At the third and sixth estimates in (26) 1 y we have used the fact that the functions ðe= hxi Þ=ð1 þ R1 1;i;j Þ and ðe=hj Þ=ð1 þ R2;i;j Þ are e-uniformly bounded. Similarly at the interface points we have

wN ;j ¼ gxN ;j þ gyN ;j ;

gxN ;j ¼

3 X k¼1

gx k;N ;j þ

3 X

y gxþ k;N ;j ; gN ;j ¼

k¼1

3 X

gyk;N ;j ;

k¼1

2 , R2 instead of a2 , q2 , R2 and a2 , q where gyk;N ;j are the same as above with 

jgx 1;N ;j j

  e  ouN 0;j hxN o2 uN 0;j   ¼ x   ux;N ;j þ ox 2 ox2  hN   2  o3 uðx; yj Þ  ðhx Þ  ; 6 Ce Nx max x h x2D  ox3  N

N

660

I.A. Brayanov / Appl. Math. Comput. 163 (2005) 645–665

  e  ouN þ0;j hxN þ1 o2 uN þ0;j   ¼ x ux;N þ1;j  ox 2 ox2  hN   2  o3 uðx; yj Þ  ðhx Þ   6 Ce Nþ1 max  ox3 ; x hN x2Dxþ N   hxN  ouN 0;j  x jg2;N ;j j ¼ x   q1;N ;j ux;N ;j þ a1;N 0;j ox  2hN    3 !  o2 uðx; yj Þ  ðhxN Þ3  o uðx; yj Þ  ðhxN Þ2     6C  ox2  þ h  ox3  ; x max x max x2Dx x2Dx h N N N N   hxN þ1  ouN þ0;j   q j ¼ u þ a jgxþ 1;N þ0;j 1;N þ1;j x;N ;j 2;N ;j ox  2 hxN   2   3 !  o uðx; yj Þ  ðhxN þ1 Þ3  o uðx; yj Þ  ðhxN þ1 Þ2     ; 6C max  max  þ hx  2 3 xþ   ox ox x2D hxN x2Dxþ N N N

jgxþ 1;N ;j j

ð27Þ

hxN jgx j  q1N ;j ux;N ;j j 3;N ;j j ¼ x 2hN ð1 þ R1 1;N ;j Þ ( )   2  ouðx; yj Þ  hxN ðhxN Þ   6 C min x ; x max  ox ; x2Dx hN e hN N hxN þ1 jgxþ j  q1N þ1;j ux;N ;j j 3;N ;j j ¼ x 2hN ð1 þ R1 Þ 1;N þ1 ( )   2 x  ouðx; yj Þ  hN þ1 ðhxN þ1 Þ   6 C min ; x max  ox :  e hN x2Dxþ hxN N xþ Here Dx N ¼ ½xN 1 ; xN  0 , DN ¼ ½xN þ 0; xN þ1 . At the last two estimates in (27) we have used again the e-uniformly boundness of the function y 1 ðe=hxi Þ=ð1 þ R1 1;i;j Þ and ðe=hj Þ=ð1 þ R2;i;j Þ. Using the fact, that the solution u can be decomposed into regular part v and singular part w ¼ ex þ ey þ exy , u ¼ v þ w, which satisfy the estimates (6) and (7) in Proposition 1, we can write the error in the form sy sxy wi;j ¼ wri;j þ wsx i;j þ wi;j þ wi;j ;

where similarly as above y;r wri;j ¼ ðLv  Lh vÞi;j ¼ gx;r i;j þ gi;j ; y;sx x;sx x h x wsx i;j ¼ ðL e  L e Þi;j ¼ gi;j þ gi;j ;

x;sxy xy h xy wsxy þ gy;sxy i;j ; i;j ¼ ðL e  L e Þi;j ¼ gi;j x;sy y;sy y h y wsy i;j ¼ ðL e  L e Þi;j ¼ gi;j þ gi;j :

Using (26) and (27) and the bounds of the regular part from (6), we obtain the following e-uniform estimates: 1 jgy;r i;j j 6 CM ;

jgx;r N ;j j 6 C;

1 jgx;r i;j j 6 CN

for i 6¼ N ;

I.A. Brayanov / Appl. Math. Comput. 163 (2005) 645–665

661

Therefore 2N 1 X

1  ; hxi jgx;r i;j j 6 CN

i¼1

M 1 X

1  : hyj jgy;r i;j j 6 CM

ð28Þ

j¼1

Now we shall find bounds for wsx i;j . Using again the estimates in (26), (27) and the bounds of the singular part ex in (7) we have 1 jgy;sx i;j j 6 CM :

Therefore M 1 X

1  hyj jgy;sx i;j j 6 CM :

ð29Þ

j¼1  To estimate gx;sx i;j in X we shall consider two cases. x (1) Let d1 ¼ 1=4 then e P m1 ln N =8 and the mesh with respect to x in X becomes uniform. Then the truncation error is 3 x 2 jgx;sx i;j j 6 Ce ðhi Þ expðm1 ð0:5  xiþ1 Þ=eÞ;

i ¼ 1; . . . ; N  1;

2 x 2 2  hxN jgx;sx ln2 N : N ;j j 6 Ce ðhN Þ 6 CN

ð30Þ ð31Þ

Therefore N 1 X

2 x 2 2 N x;sx  ln2 N : hxi jgx;sx i;j j þ hx jgN ;j j 6 Ce ðhi Þ 6 CN

ð32Þ

i¼1

(2) Let dx1 ¼ 2e ln N =m1 . Then for i ¼ 1; . . . ; N  1, i 6¼ N =2  1; N =2 and for i ¼ N the estimates (30) and (31) hold. For i ¼ N =2  1; N =2, we use that L1 ex ¼ e

o2 e x oex þ a 1 ox2 ox

is e-uniformly bounded (see Section 2). Denote by Lh1 Ui;j ¼ eðj1 Ux Þ^x;i;j þ q1;i;j Ux;i;j : Then h x  x hxN =2 jgx;sx N =2;j j ¼ hN =2 jðL1  L1 Þ eN =2;j j !    ouðx; yj Þ  x   þ max juðx; yj Þj 6C  hN =2 þ e max x2DxN =2  ox  x2DxN=2

6 Cð hxN =2 þ expðm1 ð0:5  xN =2þ1 Þ=eÞÞ 6 Cð hx þ N 2 Þ 6 CN 1 : N =2

ð33Þ

662

I.A. Brayanov / Appl. Math. Comput. 163 (2005) 645–665

Analogously 2 x x  Þ hxN =21 jgx;sx N =21;j j 6 CðhN =21 þ expðm1 ð0:5  xN =2 Þ=eÞÞ 6 CðhN =21 þ N

6 CN 1 :

ð34Þ

Now using (30) we have NX =22

3 3 1 x  hxi jgx;sx i;j j 6 Ce ðh1 Þ expðm1 e d1 Þ

1 X

expðm1 e1 ih1 Þ

i¼1

i¼1 3

e3 ðh1 Þ expðm1 e1 h1 Þ expðm1 e1 dx1 Þ 1  expðm1 e1 h1 Þ t3 expðm1 tÞ 6 N 2 max 6 CN 2 : t>0 1  expðm1 tÞ ¼

ð35Þ

Using the inequality 1  expðtÞ P t expðtÞ we get N 1 X

3 3  hxi jgx;sx i;j j 6 Ce ðh2 Þ

expðm1 e1 ih2 Þ

i¼0

i¼N =2þ1 3

¼

1 X

3

e ðh2 Þ 2 6 Ce2 ðh2 Þ expðm1 e1 h2 Þ 1  expðm1 e1 h2 Þ

6 CN 2 ln2 N expðN 2 ln2 N Þ 6 CN 2 ln2 N :

ð36Þ

From (31), (33)–(36) we obtain N 1 X

1 N x;sx  : hxi jgx;sx i;j j þ hx jgN ;j j 6 CN

ð37Þ

i¼1

Analogously 2N 1 X

1 N xþ;sx  : hxi jgx;sx i;j j þ hx jgN ;j j 6 CN

ð38Þ

i¼N þ1

Thus (32), (37) and (38) gives 2N 1 X

1  : hxi jgx;sx i;j j 6 CN

ð39Þ

i¼1

In a similar way we can prove 2N 1 X

1  ; hxi jgx;sy i;j j 6 CN

i¼1 2N 1 X i¼1

M 1 X

1  hyj jgy;sy i;j j 6 CM ;

ð40Þ

i¼1 1  ; hxi jgx;sxy i;j j 6 CN

M 1 X i¼1

1  : hyj jgy;sxy i;j j 6 CM

ð41Þ

I.A. Brayanov / Appl. Math. Comput. 163 (2005) 645–665

663

To prove (41) we use that L exy is e-uniformly bounded. Now if we take in Lemma 5 f1 ¼ gx , f2 ¼ gy , taking into account the estimates (28), (29), (39)–(41) we get (25).  

þ

Remark 2. If u 2 C 3 ðX Þ \ C 3 ðX Þ \ CðXÞ then everywhere in the prove we shall have N 1 ln N and N 1 instead of N 2 ln2 N , N 2 . In this case we obtain (24). Remark 3. In one dimensional case in [4] is constructed modification of Sammarskii’s scheme that gives almost second order OðN 2 ln2 N Þ order of convergence. In this case we need the so-called Shishkin’s decomposition with singular part satisfying Lw ¼ 0. But for 2D problem of type like one considered in this paper this modification is not possible to be constructed.

6. Numerical results Consider the problems (1) and (2) with boundary conditions gðx; yÞ ¼ 1 and coefficients a1 ðx; yÞ ¼ 2 þ 4x; a1 ðx; yÞ ¼ 1 þ x;

a2 ðx; yÞ ¼ 1 þ 2y; a2 ðx; yÞ ¼ 1 þ y;

bðx; yÞ ¼ 2 þ xy; bðx; yÞ ¼ 4  xy;

ðx; yÞ 2 X ; ðx; yÞ 2 Xþ ;

where the right-hand side f ðx; yÞ and the interface function QðyÞ are chosen such that    8 < 1 þ 1expð2ð12xÞ=eÞ  ð1  2xÞ2 1expð3ð1yÞ=eÞ  ð1  yÞ2 in X ;  1expð2=eÞ  1expð3=eÞ  u¼ : 1 þ 1expð2ð1xÞ=eÞ  4ð1  xÞ2 1expð2ð1yÞ=eÞ  ð1  yÞ2 in Xþ 1expð1=eÞ 1expð2=eÞ is the exact solution. This function exhibits typical boundary and interface layer behavior. For our tests we take M ¼ 2N . Table 1 displays the results of our numerical experiments. We observe first-order e-uniform convergence. The convergence rate is taken to be qN ¼ log2 ðkEN k1;w =kE2N k1;w Þ; where kEN k1;w is the maximum error norm error for the corresponding value of N . Fig. 2 shows the approximate solution and the maximal error for N ¼ 16 and e ¼ 0:01. It illustrates very well the boundary and interior layers behavior of the solution. Thus the numerical results support the theoretical ones and show the effectivity of special meshes.

664

I.A. Brayanov / Appl. Math. Comput. 163 (2005) 645–665

Table 1 Error of the solution on Shishkin’s meshes N e

N ¼4

N ¼8

N ¼ 16

N ¼ 32

N ¼ 64

N ¼ 128

e¼1 qN

7.65e)3 1.87

2.09e)3 1.95

5.41e)4 1.97

1.38e)4 1.99

3.48e)5 1.99

8.76e)6 –

e ¼ 102 qN

2.88e)1 0.98

1.46e)1 1.08

6.90e)2 1.19

3.03e)2 1.31

1.22e)2 1.43

4.52e)2 –

e ¼ 104 qN

2.97e)1 0.93

1.55e)1 0.93

8.14e)2 0.97

4.16e)2 0.98

2.10e)2 1.00

1.05e)2 –

e ¼ 106 qN

2.97e)1 0.93

1.55e)1 0.93

8.16e)2 0.97

4.18e)2 0.98

2.12e)2 0.99

1.07e)2 –

e ¼ 108 qN

2.97e)1 0.93

1.5e)1 0.93

8.16e)2 0.97

4.18e)2 0.98

2.12e)2 0.99

1.07e)2 –

approximate solution

maximal error

2

0.07

1.8

0.06 0.05

1.6

0.04

1.4

0.03 1.2 0.02 1 0.01 0.8 0

0 0

0.2

0.2

0.4

0.4

0.6

0

0.8

0.5 1

1

0.6

0

0.8

0.5 1

1

Fig. 2. Approximate solution and error on Shishkin mesh.

References [1] V.B. Andreev, Anisotropic estimates for the discrete Green’s function of a monotone difference operator in the singularly perturbed case and its application, Preprint TU Dresden, Math-NM06-2000, 2000.

I.A. Brayanov / Appl. Math. Comput. 163 (2005) 645–665

665

[2] V.B. Andreev, Green’s function and uniform convergence of monotone difference schemes for a singularly perturbed convection–diffusion equation on Shishkin’s mesh in one and two dimensions, DCU Maths. Dept. Preprint Series, MS-01-15, 2001. [3] V.B. Andreev, I.A. Savin, On the uniform convergence of the monotone Samarski’s scheme and its modification, Comput. Math. Math. Phys. 35 (5) (1995) 739–752. [4] I.A. Braianov, L.G. Vulkov, Uniform in a small parameter convergence of Samarskii’s monotone scheme and its modification for the convection–diffusion equation with a concentrated source, Comput. Math. Math. Phys. 40 (4) (2000) 534–550. [5] I.A. Braianov, L.G. Vulkov, Numerical solution of a reaction–diffusion elliptic interface problem with strong anisotropy, Computing 71 (2) (2003) 153–173. [6] I.A. Brayanov, L.G. Vulkov, Finite volume difference methods for convection-dominated problems with interface. Springer-Verlag, Lecture Notes in Comp. Sci. 2907 (2004) 429–437. [7] P.A. Farrel, A.F. Hegarty, J.J.H. Miller, E. O’Riordan, G.I. Shishkin, Global maximum norm parameter-uniform numerical method for a singularly perturbed convection–diffusion problem with discontinuous convection coefficient, Preprint MS-01-06, School of Mathematical Sciences, Dublin City University, 2001. [8] R.D. Lazarov, I.D. Michev, P.S. Vassilevski, Finite volume methods for convection–diffusion problems, SIAM J. Numer. Anal. 33 (1) (1996) 31–55. [9] T. Linss, Finite difference schemes for convection–diffusion problems with a concentrated source and discontinuous convection field, Comput. Methods Appl. Math. 2 (1) (2002) 41–49. [10] T. Linss, Uniform pointwise convergence of an upwind finite volume method on layer-adapted meshes, ZAMM 82 (4) (2002) 247–254. [11] T. Linss, M. Stynes, Asymptotic analysis and Shishkin-type decomposition for an elliptic convection–diffusion problem, J. Math. Anal. Appl. 261 (2001) 604–632. [12] J.J.H. Miller, E. O’Riordan, G.I. Shishkin, Fitted Numerical Methods for Singular Perturbation Problems, World Scientific, Singapore, 1996. [13] H.-G. Roos, M. Stynes, L. Tobiska, Numerical Methods for Singularly Perturbed Differential Equations, Springer-Verlag, Berlin, 1996. [14] H.-G. Roos, H. Zarin, The streamline-diffusion method for a convection–diffusion problem with a point source, J. Comput. Appl. Math. 150 (2003) 109–128. [15] A.A. Samarskii, Theory of Difference Schemes, Nauka, Moscow, 1983 (in Russian). [16] G.I. Shishkin, Grid approximations of singularly perturbed elliptic and parabolic equations, Ural Branch, Russian Academy of Sciences, Ekaterinburg, 1992 (in Russian). [17] A.B. Vasil’eva, V.F. Butuzov, Asymptotic Methods in Singular Perturbation Theory, Vysshaya Shkola, Moscow, 1996 (in Russian).