COMPUTATIONAL METHODS IN APPLIED MATHEMATICS, Vol.3(2003), No.3, pp.443-458 © 2003 Editorial Board of the Journal "Computational Methods in Applied Mathematics" Joint Co. Ltd.

THE SDFEM FOR A CONVECTION-DIFFUSION PROBLEM WITH TWO SMALL PARAMETERS

HANS-GORG ROOS

Technische Universität Dresden, Germany

ZORICA UZELAC 1 Faculty of Technical Sciences, Novi Sad, Yugoslavia Dedicated to John J.H.Miller on the occasion of his 65th birthday.

Abstract — A singularly perturbed convection-diffusion problem with two small parameters is considered. The problem is solved using the streamline-diffusion finite element method on a Shishkin mesh. We prove that the method is convergent independently of the perturbation parameters. Numerical experiments support these theoretical results.

2000 Mathematics Subject Classification: 34E15, 65L10, 65L12.

Keywords: convection-diffusion problems, singular perturbation, streamline-diffusion method, two small parameters, Shishkin-type mesh.

1. Introduction

We consider the singularly perturbed convection-diffusion problem

Lu := -eiu" + £2b(x)u' + c(x)u = f (x) x e I = (0,1), u(0) = 0, u(1) = 0, (1.1)

where b, c, f are smooth functions, b(x) ^ b0 > 0, c(x) ^ c0 > 0 and 0 < £1>2 << 1 are both small perturbation parameters. Although physical problems frequently involve the solution of boundary-value problems with many small parameters, not much is known about robust numerical methods for solving two parameter problems ( see O'Malley [8] for some physical problems involving two small parameters and for details concerning asymptotic expansions).

Shishkin and Titov [12] proved, for an exponentially fitted difference scheme on an equidistant mesh, for some v < 2/5, that

|u(xi) - Ui\ ^ Chv,

where C is independent of £ and h. Instead of using exponential fitting, the use of layer-adapted meshes becomes more and more popular. Vulanovic [14] considers Bakhvalov and

1 This paper was written during a visit to the Technical University of Dresden in September-October 2002 supported by DAAD grant and partly supported by Ministry of Science of the Republic of Serbia under grant no 1840.

Shishkin meshes but assumes e2 = ef with p > 0. Then the problem (1.1) behaves similarly as the standard reaction-diffusion problem with e2 = 0. Two-parameter problems are considered also in Li [2] and in Madden and Stynes [7], but are of different type than the one considered here.

Recently, simple upwinding on a properly chosen Shishkin mesh was analyzed independently in [5,9]. Both papers are based on similar estimates for a priori bounds for derivatives. These estimates are also the basis of our analysis, and because [5] and the present paper were written almost in parallel, we present the estimates in Section 2 ( it turns out that Lemma 2.1 and Lemma 2.2 are almost identical with Lemma 3.1 in [9]).

Simple upwinding results in a first order scheme. It is well known that for a convection-diffusion problem with e2 = 1 the streamline diffusion finite element method (SDFEM) can generate a second order scheme [13]. Therefore in the present paper we shall analyze the SDFEM on a properly chosen Shishkin mesh for the two-parameter problem and prove almost second-order pointwise convergence uniformly with respect to the parameters e1; e2.

Our analysis is absolutely different from the analysis in [13]. While Stynes and To-biska use the consistency error of the difference scheme generated by the SDFEM, our approach uses mainly finite element tools. This technique was already presented in [10] for a convection-diffusion problem with one small parameter and a point source.

In the literature so far convection-duffusion problems (e2 = 1) and reaction-diffusion problems (e2 = 0) are handled separately. Our approach offers a unified treatment of problem (1.1) for all possible classes of subproblems.

Throughout the paper C will represent a constant independent of e1,e2, and of the mesh. For a continuous function g, the value at a point xi will be denoted by gi. Furthermore , || • 11^ denotes the maximum norm, while || • \ \Ll denotes the L1 norm.

2. Some properties of the continuous problem

The solution of the reduced problem cu0 = f, in general, does not satisfy the boundary conditions, so there exist boundary layers both at x = 0 and at x =1. To describe the layers let us introduce the characteristic equation

—e1A2(x) + e2b(x)\(x) + c(x) = 0.

It has two real solutions A0(x) < 0 and A1(x) > 0. Because |Ao| < A1, the layer at x = 1 is stronger than the layer at x = 0. Let

fi0 = — max A0(x) and fi1 = min A1(x).

£€[0,1] ^€[0,1]

Lemma 2.1. For any 0 < p < 1 we have up to the a certain order q that depends on the smoothness of the data

|u(k)(x)| ^ C{1 + n0e-p№* + e-™1 (1-£)} for 0 ^ k ^ q. (2.1)

Proof. We shall prove (2.1) by induction with respect to k. For k = 0 the estimate follows immediately from the maximum principle.

(Lu)' = Lu' + e2b'u' + c'u,

then by differentiating the original equation k times, we get

Lku(k) := -£lu(k+2) + e2bu(k+l) + cku(k) = gk, for k =1,... (2.2)

c0 = c, g0 = f

Ck = c + k£2b, gk = gk-1 - c'k-1u(k-l), for k =1,.... The (k + 1)th derivative satisfies the equation

Lk+iu(k+1) = gk+i

|gk+i| ^ C{1 + ¡4e-pil0x + ¡ike-p^(1-x)}.

In the first step of the proof we assume that we have already proved the inequalities

|u(k+1)(0)| ^ C1ik0+1 and |u(k+1)(1)| ^ Cik+1. (2.3)

We shall bound u(k+1) using (2.2), (2.3) and the following barrier function

$ = C (1 + i'k+1wo,p(x) + ik1+1W1,p (x)),

wo,p(x) = e-pil°x and w^x) = e-p^l(1-x).

It is sufficient to show that Lk+1w0p ^ a0w0>p and Lk+1w1p ^ a1w1>p with some a0,a1 > 0. Since

Lk+1w0 ,p = {-£1p2iI - £2Pl0b + c + O£)}w0,p ^ {p(-£1l0 - £2bl0) + c + O(£2)}w0,p

^ {p(-£1\l(x) + £2b(x)\0(x)) + c(x) + O(£2)}w0,p = {c(x)(1 - p) + O(£2)}w0,p ,

for any 0 < p < 1, we have the desired property if £2 is small enough. Similarly, we estimate Lk+1 w1 p.

In the next step we shall prove (2.3) starting from

Lku(k) = gk and |u(k)(0)| ^ Cik, u(k)(1)| ^ Cik.

Setting

u(k) = y + y0 with y0 = u(k) (0)(1 - x), we obtain for y the partially homogeneous problem

Lky = gk - £2by0 - cky0, y(0) = 0, y(1) = uk (1).

Now we bound y by a barrier function of the type

$ = C1(1 + e-p^°x - 2e-d^0x)ik + C2(e-p^l(1-x) - e-p^)ik,

where 9 > 0 is large enough. It is easy to check that $(0) = 0 and $(1) ^ y(1) if C2 is sufficiently large. Since Lke-d^°x ^ 0 and Lke-p^°x ^ 0,Lke-pw(1-x) ^ 0, if 9 and p are chosen adequately, $ majorizes y if C1 is sufficiently large. Then the definition of the derivative gives

y (0)| ^ C^k+1 + C^1+1e-Wl ^ Ctf?1-

Similarly we prove the second estimate in (2.3). □

From (2.1) we obtain immediately the existence of a Shishkin decomposition [4, p. 892]: Lemma 2.2. The solution u has the representation u = S + E0 + E1, where

S(k)(x)| ^ C for 0 ^ k ^ g

|E0k)(x)| ^ Cvke-p^°x for 0 ^ k ^ g,

|E(k)(x)| ^ C^1 e-Wl(1-x) for 0 ^ k ^ g.

Let us for a moment assume b and c to be constant. Then the two extremal layer situations are characterized by

(i) e2 << e1

which implies ¡i0 ~ ¡i1 ~ \Jc/e1, and we have layers similar to the reaction-diffusion case

e2 = 0;

( ii ) e 1 << e 22

which implies ¡i0 ~ c/(be2), but fi1 ~ (be2)/e1 is much larger than ¡i0, so the layer at x = 1 is much stronger than the layer at x = 0.

In the intermediate case, when neither (i) nor (ii) hold, ¡i0 and fi1 move in the interval given by the values from (i) and (ii).

Fig. 1 below presents the exact solutions, for various e1 and e2, for the problem considered in [3]:

—e1u" + e2u' + u = cos nx x E (0,1), u(0) = u(1) = 0.

Figure 1. The exact solution for £\ = 0.001, e2 = 0.01 and £\ = 0.001, e2 = 0.1

The exact solution is

u(x) = a cos nx + b sin nx + AeXox + Be-Xl(1-x),

£in2 + 1 £2^ . 1 + e-Xl

b -;-~--T^, A = — a-

4n2 + (£rn2 + 1)2' £\n2 + {£n2 + 1)2' 1 - eXo-Xl'

D 1 + ex° £2 ^ v/£fT4£7

B — a-:-r~ 1 An 1 — -.

1 - eA°"Ai' 0 ;1 2£1

3. The SDFEM with linear elements 3.1. The scheme

For problem (1.1), the standard weak formulation is: find u E V such that

aG(u,v) :— (£1u',v') + (£2bu! + cu,v) — (f,v) Vv E V,

where V — H^(Q) denotes the usual Sobolev space and (•, •) is the scalar product on L2(Q). On a given mesh

0 — x0 < x1 < • • • < xN-1 < xN — 1,

with the mesh step hi — xi — xi-1, i — 1,..., N, we set hi — (hi + hi+1)/2, i — 1,..., N — 1, and for a given mesh function {u^} we define the forward and backward difference operators D+ and D- by

ui+1 — ui n J^.- ui — ui-1

D^ui :— —-- and D ui :— ---.

hi+1 hi

Let Vh C V be the space of continuous linear finite elements with the basis functions

x — xi 1

for x E [xi-1,xi\,

fa -x for x E [xi,xi+1],

k 0 for x E [xi-1,xi+1].

The Galerkin approach is characterized by: find uh E Vh such that

ao(uh, vh) — (f, vh) Vvh E Vh.

Based on the decomposition u — S + E0 + E1 given in Lemma 2.2, we take p — 1/2 and introduce a corresponding S-mesh. Namely, the mesh is equidistant on Q0 :— [0,a0], :— [a0, 1 — a1] and Q1 :— [1 — a1,1]. We choose

a0 — — ln N, a1 — — ln N Ho H1

(such that e~^°x/2\x=a0 — N-T, for example). Since ¡i0 and /i1 are large and /i1 ^ ¡i0, the mesh is fine on Q0 and Q1 and coarse on Qg, with mesh step sizes

h„ — o (N-lhN), hl — o (N-^nN) and H — h, — o (N-1),

respectively.

If u1 denotes the linear interpolant of u in Vh, simple modifications of the well-known techniques of [6] result in:

Lemma 3.1. The interpolation error satisfies

\\u — u1 ll^n, ^ CN-2, \\u — u1 |UMn9 ^ C(N-1 lnN)2. The Galerkin method generates a difference scheme which can be written in the form tn D+ui — D-uif ui+1 — ui ui — u—A

Lgui := —e1-f--+ e2 Pi+1-f--Pi-1-f- + Ti-1^-1+Yiui+Yi+1ui+1

hi V hi hi J

=h1 f*i'

where 0i denotes again the basis function of Vh with 0i(xj) = 6ij, and

Xi Xi+1

Pi-1 = J b0i-10i < 0, Pi+1 = J b0i+10i > 0,

xi - 1 xi

Xi Xi+1 Xi+1

Yi-1 = t~ c0i-10i, Yi = t c0i0i, Yi+1 = ^r c0i+10i > 0. hi J ^ hi J 1 hi J

When c = const we have Yi-1 = hi/(6hi),Yi = 2/3, Ym = hi+1/(6hi). If we wish to have an M-matrix, we have the conditions

e2$+1 + hi Yi+1 ^t— and Y—hi ^ e2(—A-1) + 7-, hi+1 hi

which for a lumped scheme (Yi—1 = Yi+1 = 0 and YiT™ = Yi-1 + Yi + Ym) reduce to

S2ßi+1 ^

Because these conditions are too restrictive, one has to stabilize the FEM to generate M-matrices.

For the SDFEM we modify the Galerkin bilinear operator aG by introducing lumping for the c term:

, v) := e1 (u',v')+ e2(bu',v) '

The streamline-diffusion method is characterized by the stabilization

aGL(u,v) := ei(u',v') + e2(bu',v) + ^ hkCkUkVk.

a st (u,v) := ^^ / ök(-e1u" + e2bu' + cu)v'

k=1Xfc-i

with some elementwise-defined parameter 6k which will be determined later.

Now we again choose Vh C V to be the space of piecewise linear elements on the given grid. Introducing a(u,v) = aoL(u,v) + aST(u,v), we define the SDFEM discretization: find uh E Vh such that

a(uh, vh) = fh(vh) yvh E Vh

fh(vh):=(f,vh) + Y, / h H.

k-1xk-i

The SDFEM generates a finite difference scheme, which can be written in the form

LNui :— —£1(D+ui — D-ui) + aiD+u + PiD'u + Yiu — fh(Vi), (3.1)

ai — hi+1 J (£2bpi+1pi + bi+1 £2by'i+1w'i + 5i+1cj^i+1^'i), xi

Pi — hi J (£2b^'iVi + $i£2byivi + 8icpi-1pi_1),

Xi Xi+1

Yi — hiCi + Si J ctpi + 5i+1 j ctpi.

xi-l xi

The choice of the parameter Si is determined by the structure of the coefficient matrix of scheme (3.1), as we wish to have an M-matrix. Namely, if the local mesh size is small enough, standard Galerkin can be applied,i.e., it is possible to choose Si — 0, i — 1,...,N. This leads to the condition

ai-1^-0 ^ £1 i.e., £2hibi ^ £1 for i — 1,...,N,

where bi — (1/hi) JX b^i-1. Otherwise, we choose Si from the condition ai-1 — 0, i — 1, ... , N. With the notations

xi xi Ci—h j cw, bi—h j b,

xi—1 xi 1

this condition is equivalent to

Finally we obtain

Si ( + «I = £2bi

0 if -2hibi ^ £1;

-2bihi . (3-2)

otherwise.

£2bi + Ci hi

We can see that

Si = O(hi) if hi is smaller than £2,

Si = O(£2 ) if £2 is smaller than hi.

In the case of £2 = 1 and £1 = £, we choose Si = 0 on the fine mesh. If we check the condition £2hibi ^ £1; we see that it is satisfied on Q1 because fi1 ^ (£2b0)/£1, i.e., £1 ^ £2bo/^1- On , however, the situation is different.

Summarizing, we have

on Hg:

on Hn:

Si = 0;

£2 Hbi

£2bi + CiH

£ 2 hn b i

if £2Hbi ^ £i, otherwise;

if £2hnbi ^ £i,

otherwise.

£2bi + Cihn

Since hn = O(N-1 ln N/^n) and 1/^n = O(£2 + £^2), the stabilization on Hn is necessary if £1 is extremely small, roughly speaking, in comparison with £2N-1 ln N. The corresponding stabilization parameter is of order O(hn).

3.2. Pointwise error analysis

For estimating the pointwise error we use an approach similar to [10,11] based on a discrete Green's function. Let us define a function Gj <E Vh by

a(^i,Gj) = Sij for all i = 1,...,N - 1.

In our estimates we need bounds for the L^ and W1,1 norm of Gj (the L1 norm of its derivative). Since Gj is piecewise linear, the continuous and the discrete versions of these norms are equivalent. For the discrete norms Andreev [1] proved estimates for an upwind finite difference operator in the case of £2 = 1. In the Appendix we shall prove similar estimates for the two-parameter problem.

It is a well known fact that the scheme generated by SDFEM can be written as a generalized upwind scheme (see [13] and especially [11], Section 5). We do not want to repeat the details from [11], but it is not difficult to see that scheme (3.1) with (3.2) is equivalent to

Lh Uh = fh

(LhU)i := -(pi+1^+Ui - PiD-Ui) + £2b* i T i-1 + c*Ui. hi hi

Here all coefficients are the result of rewriting (3.1) with (3.2) in this form, especially

* = + Si- _ Si+1 C

ci ci + y ci J Ci+1.

Do we have c* ^ c0 > 0 if ci ^ cn > 0? It is not difficult to see that we are in trouble only if Si = 0 and Si+1 = 0, but this happens only in the first subinterval of Hg with follows the transition point an. On that interval we have

Si+1 = - if £1 ^ £2HCi.

£2bi + CiH

Do we need c* ^ cQ > 0? The proof in the Appendix shows that the strong positivity is only necessary if e"^ ^ ks1. Let us assume for a moment that b = const. and c = const. Then

c. > c(i -1 )=<> H/2

H (e2b + cH/2y e2b + cH/2' Because e^/n, ^ e1 ^ e2Hb implies e2 ^ nbH, it is true that

H/2 1/2

' ^ ^—r- > 0.

£2 b + cH/2 Kb2 + c/2 The nonconstant-coefficient case can be handled similarly.

The discrete Green's function Gj now satisfies (see the Appendix for some details and estimates for the Green's function)

LhG = Sh.

Using the linear operator P : V ^ Vh with

we have the following representation for the error at the mesh point Xj :

(u - uh)(xj) = P(u1 - u)\x=xj + (Pu - uh)\x=xj. (3.3)

Furthermore,

P(u1 - u)\x=xj = a(u - u1, Gj)

and the consistency error Pu - uh vanishes for a consistent finite element method.

Based on the error representation (3.3) we want to estimate the pointwise error of our SDFEM on a S-mesh. Since the stabilization itself is consistent, the consistency error (Pu - uh)\x=xj reduces to the error induced by lumping of the c-term. That error part can be estimated as in [11] and it is of order (N-1 ln N)2. To estimate the other error part we start from

.I — „IV ^jv\ , I „ („. „.I\

a(u - u1, Gj) = ei((u - u1 )', (Gj)') + £ hkck(u - u1 )(xk)Gj(xk)

+ e2(b(u - u1 )', Gj) + &k(-eiu" + e2b(u - u1 )' + c(u - u1 ), (Gj)')

The first two terms are equal to zero. To estimate the remaining terms we want to apply the results concerning the discrete Green's function from the Appendix.

Now we estimate several error contributions. First, using integration by parts

\e2(b(u - u1 )', Gj)\ ^ C\\u - u1 |U|Gj\\Ll + C\e2(u - u1, (Gj)')\.

Furthermore,

\e2(u - u1, (Gj)')\ ^ Ce2\\u - u1 \\^\(Gj)'\\li ^ C . f2 (N-1 lnN)

V e2 + e1

Next we estimate the contribution of the stabilization on Qg. Let us recall that

£2 + H

Sg = O ( 2 Tr ) if £1 ^ C£2H.

(i) (u", (Gj)'). For the smooth part we get

C£10g j.— ^ ch2, \£1 ^ C£2H, Sg = o(H), J:2— ^ 1 £22 + £1 £22 + £1

For the layer part we get ( only for E1 because that layer is stronger)

C£1Ôg\\E'1 IU1 ||(Gj)'|U ^ C£1Ôg^1HTH-1\\Gj|U This error contribution is small enough because £1v1/^£2 + £1 is bounded.

(ii) Z£20k (b(u — u1 )', (Gj)').

Since (Gj)" = 0, (u — u1 )xj = 0, integration by parts yields

\T,£2Sk(b(u — u1 )', (Gj)')| ^ C£2ÔgWu — u1 |U|(Gj)'WL1 ^ CSgN-2.

(iii) SSk(c(u — u1 ), (Gj)'). Using Sg = O(£2), we get

|£Sk(c(u — u1 ), (Gj)')\ ^ C£2 , 21 Hu — u1|U ^ CN-2.

£22 + £1

Finally, we consider the contributions of the stabilization on Q0. Let us recall that

So = O(ho) if £1 ^ C£2ho,

'N-1 ln N\ 1 Vo J V0

ho = O - , - = O £2 + W£22 + £1

(i') Z£Mu'', (Gj)').

For the smooth part there is no problem. For the layer part E2 (E1 is extremely small on Q0) using £1 ^ C£2h0,h0 — 0(N-1 ln N)/y0, we get

\E£1So(u", (Gj)')\ ^ £1Sov2||(Gj)'|li ^ C£2(hoVo)2^== ^ C (N-1 ln N)2.

£22 + £1

The two other error terms on Oo, which correspond to (ii) and (iii), can be estimated as on

Theorem 3.1. The pointwise error of the SDFEM on a S-mesh satisfies the following estimate:

Hu — uhH™ < C (N-1 ln N)2 .

4. Appendix: The discrete Green's function

Denoting the upwind difference operator on an arbitrary mesh by Lh :

(LhUi) := - j- (D+Ui - D-u) + biS2 hi

Ui - Ui-1 h i

(assuming bi > 0,ci > 0 as above), we define the discrete Green's function G(xi,^j) for Lh for fixed , by

LhG(Xi,Cj) = Sh(xi,0) G(0,Cj) = G(1,0) = 0,

Sh(xi,ij) = < j

hT if i =j,

0 if i = j.

Setting W = \J¿2 + £1, we shall prove

\\G(Xi,tj)\u ^ C.

To do this, let us introduce the barrier function w defined by

11(1 - fhk)~ for j>i, k=i

i -1 n(1+fhk)~ for j<i.

(Lhw)i ^ j-, hi

if the constant C in (4.2) is chosen adequately. We use the abbreviations

^ n - ^ n

jo = y < 0, j1 = y > 0

and introduce

Wi+1 wi wi 1

(1 - hohi+1) 1 1

(1 + jnhi)-1

and we obtain (Lhwi)i

1 - johi+1

C ( £1

W\ hi \1 - johi+1 1 + j1 hi

1 + j1 hi

+ £2bihiT

From [3] we know already that Lhw ^ 0 at the mesh points {x1,..., xi-1} and {xi+1, ...,xN-1}. It remains to verify

hi(1 + j1hi)

Case A: ei ^ K£1

In this case we neglect the b-term and use the existence of a constant C1 such that

£2 + W

fi1 ^ 2C1-— 21 as well as — ¡i0 ^ 2|1.

£1 1 1

and we obtain

^ -=-, - ^

1 — Hohi+1 1 + v1h1+ i-h 1+ v1hi

(T \ ( £1(V1 — Vo) , h i ^

(LhWi)i > hA W (1 + Vh) +Ci W)'

h VW(1 + Ii1hi) W If hi > pW for some p > 0, we have the desired estimate. Otherwise, we have

£1(11 — 10) > £1(11 — I0) 1

W (1+ i\hi) W 1 + i{phi

Since £1(-1 — -0)/W is bounded away from zero and i1W is bounded if ¿2 ^ K£1, it follows that (4.3) holds true if C is chosen adequately.

Case B: £22 > k£1

In this case we neglect the c-term and get

(Lhw)i > C - 1 _ \£1jl1 — £1 — -,1 1hi + £2bi ^1hi W hi(1+-1hi) [ 1 — -0hi+1

C 1 ) /_ - N -1hi + -0hi+1 . 7-7

' £1(1-1 — -0) — £1 -0—-—T--+ £2bi

W hi(1 + j-1hi) \ 1 — i-ohi+1

-1hi + -ohi+1 = -1hi(1 — -o hi+1) + -ohi+1 (1 + -1hi)

we obtain

(Lhw)i ^ C - 1 _ <£1 (-1 — -o) + (£2bi — £1 -o)-1hi — £1 -o-ohi+1(1 + -1hi)

W hi(1 + -1hi){ 1 — -ohi+1

■ £10-1 — -o) + £2bi — £1 -o - h

C ) W + W -1hi £1-o-ohi+1

hi I 1+-1hi W(1 — -ohi+1)

This expression has the structure

(Lhw)i > W i D+Er1 — Fi W { 1 + J

with \F\ < E, \F\ < D and \F\ < —£1 -o/W.

Since D — F is bounded away from zero, it remains to show that E — F is bounded away from zero also. We have the necessary bound since 4 ^ K£1 and E — F > £2bi/W.

Remark 4.1. It is a simple exercise to replace the difference operator by the more general operator

(Lhu)i := -h {pi+\D+Ui — PiD-Hi) + £2% Ui j^1-1 + c*Ui, hi hi

assuming pi ^ p > 0, and again b* > 0,c* > 0. (SDFEM is related to a difference operator of that form)

Remark 4.2. The estimate remains true if we consider the standard upwind difference operator with the difference quotient (ui — ui-1)/hi.

Let us summarize the properties of the discrete Green's function for the two parameter problem:

Lemma 4.1. The discrete Green's function G satisfies

\\G\\Ll ^ C; \\G\\„ ^ J—; \\G'\\Ll ^ Jl--

V£2 + £1 V £2 + £1

The L 1 estimate can be obtained by summing up the equations for G at all meshpoints. The W1,1 bound follows from the L^ bound from the same arguments used by Andreev [1, p. 925].

5. Numerical experiments

In this section we present some numerical results for scheme (3.1), applied to the following example.

Example 5.1.

— £1u" + £2u' + u = cos nx x E (0,1), u(0) = u(1) = 0. The exact solution is

u(x) = a cos nx + b sin nx + AeXox + Be-Xl(1-x),

£1n2 + 1 £2n

£2n2 + (£1n2 + 1)2 ' £2n2 + (£1n2 + 1)2 '

A a 1 + e-Xl B a 1 + eX° A £2 T V4 + 4£1

A — — a-:-r—, B — a-:-r—, An 1 — -.

1 - eA°-Ai 1 - eA°-Ai 0,1 2£i

Table 1. Errors EN at the mesh points and convergence rates OrdN for N = 25

£2 £1

1 10-2 10-4 10-6 10-8 10-10 10-12

1 2.272e-05 2.124e-04 1.513e-03 6.686e-02 9.962e-02 1.031e-01 1.024e-01 E N

2.004 2.002 2.028 3.782 -0.033 -0.013 -0.096 OrdN

10-1 1.189e-05 1.316e-04 8.723e-04 1.173e-02 1.585e-01 2.067e-01 2.215e-01 E N

2.004 2.003 2.006 2.130 0.411 0.0897 0.094 OrdN

10~2 6.416e-06 7.698e-05 4.598e-04 4.171e-03 2.822e-01 3.223e-01 3.384e-01 E N

2.002 2.002 2.002 2.028 3.843 0.097 0.088 OrdN

10-3 3.850e-06 4.882e-05 3.673e-04 1.124e-03 2.330e-02 2.968e-01 2.981e-01 E N

2.002 2.001 1.998 1.996 2.136 0.028 0.100 OrdN

10-4 2.682e-06 3.538e-05 3.068e-04 1.268e-03 3.597e-03 1.772e-01 1.882e-01 E N

2.002 2.000 2.000 1.996 1.973 3.195 -0.147 OrdN

10~5 2.147e-06 2.900e-05 2.691e-04 1.260e-03 4.243e-03 1.392e-02 1.352e-01 E N

2.001 2.001 1.999 1.995 1.980 2.029 0.497 OrdN

10~6 1.887e-06 2.588e-05 2.486e-04 1.221e-03 4.255e-03 1.467e-02 4.869e-02 E N

1.996 1.999 1.999 1.995 1.955 1.878 1.948 OrdN

10-v 1.770e-06 2.441e-05 2.373e-04 1.189e-03 4.216e-03 1.524e-02 4.609e-02 E N

2.001 2.000 1.995 1.991 1.974 1.912 1.670 OrdN

10-8 1.711e-06 2.367e-05 2.316e-04 1.169e-03 4.154e-03 1.505e-02 4.266e-02 E N

2.001 2.000 1.994 1.988 1.978 1.912 1.506 OrdN

10-9 1.682e-06 2.329e-05 2.290e-04 1.158e-03 4.112e-03 1.480e-02 3.760e-02 E N

2.000 2.000 1.995 1.987 1.978 1.912 1.344 OrdN

10-10 1.667e-06 2.311e-05 2.277e-04 1.153e-03 4.088e-03 1.463e-02 3.494e-02 E N

2.000 2.000 1.996 1.986 1.978 1.913 1.264 OrdN

10-11 1.656e-06 2.296e-05 2.267e-04 1.150e-03 4.068e-03 1.449e-02 3.290e-02 E N

2.000 2.000 1.996 1.987 1.978 1.913 1.204 OrdN

10-12 1.653e-06 2.293e-05 2.265e-04 1.149e-03 4.063e-03 1.445e-02 3.238e-02 E N

2.000 2.000 1.996 1.987 1.978 1.913 1.189 OrdN

10-13 1.652e-06 2.292e-05 2.264e-04 1.148e-03 4.062e-03 1.444e-02 3.221e-02 E N

2.000 2.000 1.996 1.987 1.978 1.913 1.184 OrdN

10-14 1.652e-06 2.292e-05 2.264e-04 1.148e-03 4.061e-03 1.444e-02 3.221e-02 E N

2.000 2.000 1.996 1.987 1.978 1.913 1.184 OrdN

In Table 1 and Table 2 we present maximum pointwise errors EN for N = 25 and N = 210, respectively, which are given by

EN = max |u(xi) — ui|.

Assuming the convergence order (N-1)r, for some r, we use the double-mesh method to compute the experimental rates of convergence OrdN. For every fixed £1 and £2, we computed the rate of convergence from

_ N ln EN — ln E2N

Ord =---.

The results clearly show robust uniform convergence of order two.

Table 2. Errors EN at the mesh points and convergence rates OrdN for N = 210

£2 £1

1 10-2 10-4 10-6 10-8 10-10 10-12

1 6.768e-08 7.393e-07 5.712e-06 6.753e-05 1.502e-04 1.407e-01 1.063e-01 EN

4.996 2.577 2.037 2.003 1.700 11.58 -0.002 OrdN

10-1 3.532e-08 4.236e-07 3.283e-06 4.074e-05 3.732e-04 3.813e-04 1.969e-01 EN

4.674 3.404 2.145 2.005 1.701 1.686 -0.443 OrdN

10~2 1.982e-08 2.504e-07 1.595e-06 1.595e-05 2.843e-04 7.012e-04 3.467e-01 EN

3.566 3.223 2.674 2.019 2.000 1.693 10.58 OrdN

10~3 1.268e-08 1.680e-07 1.338e-06 4.249e-06 8.206e-05 8.452e-04 9.594e-04 EN

2.905 2.683 2.351 2.167 2.002 1.698 1.691 OrdN

10-4 9.414e-09 1.281e-07 1.156e-06 4.907e-06 1.433e-05 3.720e-04 9.376e-04 EN

2.499 2.372 2.171 2.059 2.015 2.000 1.699 OrdN

10~5 7.889e-09 1.086e-07 1.033e-06 4.916e-06 1.693e-05 5.448e-05 6.448e-04 EN

2.264 2.196 2.086 2.026 2.005 2.001 1.706 OrdN

10~6 7.158e-09 9.900e-08 9.630e-07 4.782e-06 1.724e-05 6.425e-05 2.149e-04 EN

2.136 2.101 2.044 2.013 2.002 2.000 1.999 OrdN

10-v 6.801e-09 9.426e-08 9.258e-07 4.672e-06 1.693e-05 6.529e-05 2.530e-04 EN

2.069 2.052 2.022 2.006 2.001 1.999 1.998 OrdN

10-8 6.625e-09 9.190e-08 9.067e-07 4.606e-06 1.663e-05 6.404e-05 2.570e-04 EN

2.035 2.026 2.011 2.003 2.000 1.999 1.999 Ordun

10-9 6.538e-09 9.072e-08 8.971e-07 4.570e-06 1.644e-05 6.285e-05 2.516e-04 EN

2.018 2.013 2.006 2.002 2.000 2.000 1.997 OrdN

10-10 6.494e-09 9.013e-08 8.922e-07 4.551e-06 1.634e-05 6.207e-05 2.470e-04 EN

2.009 2.007 2.003 2.001 2.000 1.999 1.998 OrdN

10-11 6.462e-09 8.969e-08 8.885e-07 4.537e-06 1.625e-05 6.144e-05 2.423e-04 EN

2.002 2.002 2.001 2.000 2.000 1.999 1.998 OrdN

10-12 6.454e-09 8.958e-08 8.876e-07 4.533e-06 1.623e-05 6.127e-05 2.410e-04 EN

2.001 2.000 2.000 2.000 2.000 1.999 1.998 Ordun

10-13 6.451e-09 8.955e-08 8.873e-07 4.532e-06 1.623e-05 6.122e-05 2.406e-04 EN

2.000 2.000 2.000 2.000 2.000 1.999 1.998 OrdN

10-14 6.451e-09 8.955e-08 8.873e-07 4.532e-06 1.623e-05 6.121e-05 2.406e-04 En

2.000 2.000 2.000 2.000 2.000 1.999 1.998 OrdN

References

[1] V. B. Andreev, The Green function and a priori estimates of solutions of monotone three-point singularly perturbed finite-difference schemes, Differ. Equations, 37 (2001), pp. 923-933.

[2] J. Li, A robust finite element method for a singularly perturbed elliptic problem with two small parameters, Comput. Math. Appl., 36 (1998), pp. 91-110.

[3] T. Linfi, Untersuchungen zur numerischen Behandlung singular gestorter Randwert- und Anfangsrandwertaufgaben auf Shishkin-Gittern, Diploma thesis, Department of Math., Technical University of Dresden, 1996.

[4] T. Linfi, The necessity of Shishkin decompositions, Appl. Math. Lett., 14 (2001), pp. 891-896.

[5] T. Linfi and H.-G. Roos, Analysis of a finite difference scheme for a singularly perturbed problem with two small parameters, Report MATH-NU-18, TU Dresden, 2002.

[6] T. Linfi, Layer adapted meshes for convection-diffusion problems, Comput. Meth. Appl. Mech. Eng., 192 (2003), pp. 1061-1105.

[7] N. Madden and M. Stynes, A uniformly convergent method for a coupled system of two singularly perturbed linear reaction-diffusion problems, Preprint No. 1, School of Math., University of Cork, 2002.

[8] R. E. O'Malley, Two-parameter singular perturbation problems for second-order equations, J. Math. Mech., 16 (1967), pp. 1143-1164.

[9] E. O' Riordan, M.L. Pickett, and G.I. Shishkin, Singularly perturbed problems modelling reaction-convection-diffusion processes, Report MS-02-13, Dublin City University, 2002.

[10] H.-G. Roos and H. Zarin, The streamline-diffusion method for a convection-diffusion problem with a point source, J. Comp. Appl. Math., 150 (2003), pp. 109-128.

[11] H.-G. Roos and H. Zarin, A second order scheme for singularly perturbed differential equations with discontinuous source term, East-West J. Numer. Math., 10 (2002), pp. 275-289.

[12] G. I. Shishkin and V. A. Titov, A difference scheme for a differential equation with two small parameters at the derivatives, Numer. Meth. Contin. Medium Mech., 7 (1976), No. 2, pp. 145-155.

[13] M. Stynes and L. Tobiska, A finite difference analysis of a streamline diffusion method on a Shishkin mesh, Numer. Algorithms, 18 (1998), pp. 337-360.

[14] R. Vulanovic, A higher-order scheme for quasilinear boundary-value problems with two small parameters, Computing, 67 ( 2001), pp. 287-303.

Received 26 Nov. 2002 Revised 02 Jun. 2003