On the uniform convergence of a finite difference scheme for time ...
Applied Mathematics and Computation 216 (2010) 1478–1488
Contents lists available at ScienceDirect
Applied Mathematics and Computation journal homepage: www.elsevier.com/locate/amc
On the uniform convergence of a finite difference scheme for time dependent singularly perturbed reaction-diffusion problems C. Clavero *, J.L. Gracia Department of Applied Mathematics, University of Zaragoza, 50018 Zaragoza, Spain
a r t i c l e
i n f o
a b s t r a c t
Keywords: Parabolic reaction–diffusion problems Semidiscrete problems Asymptotic behavior Uniform convergence Special nonuniform mesh
In this work we are interested in the numerical approximation of 1D parabolic singularly perturbed problems of reaction–diffusion type. To approximate the multiscale solution of this problem we use a numerical scheme combining the classical backward Euler method and central differencing. The scheme is defined on some special meshes which are the tensor product of a uniform mesh in time and a special mesh in space, condensing the mesh points in the boundary layer regions. In this paper three different meshes of Shishkin, Bahkvalov and Vulanovic type are used, proving the uniform convergence with respect to the diffusion parameter. The analysis of the uniform convergence is based on a new study of the asymptotic behavior of the solution of the semidiscrete problems, which are obtained after the time discretization by the Euler method. Some numerical results are showed corroborating in practice the theoretical results on the uniform convergence and the order of the method. Ó 2010 Elsevier Inc. All rights reserved.
1. Introduction In this paper we consider 1D initial-boundary value parabolic reaction–diffusion singularly perturbed problems given by
8 ðx; tÞ 2 Q ¼ X ð0; T ð0; 1Þ ð0; T; > < ut þ Lx;e u ¼ f ðx; tÞ; uðx; 0Þ ¼ 0; x 2 X; > : uð0; tÞ ¼ uð1; tÞ ¼ 0; t 2 ð0; T;
ð1Þ
where the spatial differential operator is defined as
Lx;e u euxx þ bðx; tÞu:
ð2Þ
We assume that 0 < e 6 1 and it can be arbitrarily small, the reaction term satisfies b P b > 0 for all x 2 Q (the forthcoming analysis can be easily extended to the case b P b P 0 using a standard change of variable), and the data are sufficiently smooth functions b; f 2 C ð4;2Þ ðQ Þ. Moreover we assume compatibility conditions between data of problem (1) in order that the exact solution u 2 C ð4;2Þ ðQ Þ. Sufficient conditions for this regularity are given by
f ð0; 0Þ ¼ f ð1; 0Þ ¼ 0; ðLx;e f Þð0; 0Þ ¼ ft ð0; 0Þ;
ðLx;e f Þð1; 0Þ ¼ ft ð1; 0Þ:
* Corresponding author. E-mail addresses: [email protected] (C. Clavero), [email protected] (J.L. Gracia). 0096-3003/$ - see front matter Ó 2010 Elsevier Inc. All rights reserved. doi:10.1016/j.amc.2010.02.050
C. Clavero, J.L. Gracia / Applied Mathematics and Computation 216 (2010) 1478–1488
1479
pffiffiffi It is well-known (see [6],[12]) that the solution of (1) has a boundary layer at x ¼ 0; 1 of a width O ej ln ej . In the next result we remember the asymptotic behavior with respect to e of that solution showing the presence of two boundary layers. Theorem 1 (See [6]). The solution of (1) satisfies
juðk;mÞ ðx; tÞj 6 C 1 þ ek=2 Be ðxÞ ;
0 6 k þ 2m 6 4;
ð3Þ
where the function Be ðxÞ is given by
Be ðxÞ ¼ e
pffiffiffiffiffi
b=ex
þ e
pffiffiffiffiffi
b=eð1xÞ
:
This behavior motivates the use of uniformly convergent methods (see [11,12]) to obtain efficient numerical approximations for any value of e. There are many papers in the literature dealing with the numerical approximation of the solution for time dependent singularly perturbed reaction–diffusion problems. In [2,8,10] 1D problems like (1) are considered, in [1,4] the case of 2D parabolic problems is analyzed and coupled linear systems are solved in [5,9]. Different techniques have been used to prove the uniform convergence of the numerical method defined for each problem. In some papers the totally discrete method is directly considered and the analysis of the convergence is based on appropriate bounds for the truncation error and the discrete maximum principle for the totally discrete operator (see [9]). Nevertheless, this technique cannot be extended to higher order time approximation than the Euler method due to the lack of a discrete maximum principle for the totally discrete method. To elude this difficulty, in other papers (see [1–4]) the analysis of the uniform convergence is made in two steps. In the first one, the problem is discretized in time and appropriate estimates for the corresponding local error are proved. In the second step, the set of one dimensional problems resulting from the time discretization are discretized in space and it is necessary to prove the uniform convergence of the corresponding spatial discretization. So far, and motivated by the classical analysis of the convergence for initial value problems, in the previous papers this analysis is based on the uniform convergence of the numerical scheme which discretizes some auxiliary problems arising in the definition of the method. To deduce that uniform convergence, the asymptotic behavior of the solution of these auxiliary problems was needed. The proof of the uniform convergence of the totally discrete schemes combining the two steps, uses a recursive argument which needs that the powers of a discrete transition operator arising in the analysis be uniformly bounded with respect to e. To our knowledge this question is only completely solved in the case that the backward Euler method is considered and it is an open question for higher order time approximations. The idea of this paper is to show a new strategy which can be used to overcome the main difficulties associated to the second technique. With this purpose, the uniform convergence of the numerical method considered to discretize the singularly perturbed problems resulting from the time discretization, will be proved without any use of auxiliary problems. This new proof requires to know the asymptotic behavior of the solution of these problems, which is based on an inductive argument and the well-known behavior of steady singularly perturbed reaction–diffusion problems (see [11,12]). In this new analysis, it is not necessary any recursive argument and therefore the stability of any transition operator, making simpler the analysis of the uniform convergence. The paper is structured as follows. In Section 2 we remember that the backward Euler method defined on an equidistant mesh is a first order uniformly convergent method. Moreover, the asymptotic behavior of the solution of the semidiscrete problems resulting after the discretization in time is proved, which will be crucial for the posterior analysis of the uniform convergence of the numerical scheme. In Section 3 central differencing is used to discretize the previous semidiscrete problems. The scheme is constructed on different nonuniform meshes; concretely we use a piecewise uniform Shishkin mesh (see [12]), a Bahkvalov mesh (see [7]) and a mesh based on a generating function of polynomial type (see [13]). It is proved that the scheme has almost second order uniform convergence on the Shishkin and the Vulanovic meshes and it has second order on the Bahkvalov mesh. In Section 4 some numerical results for two different examples are given illustrating in practice the efficiency of the method and corroborating the order of uniform convergence theoretically proved. Henceforth, C denotes a generic positive constant independent of the diffusion parameter e and also of the discretization parameters N and s. We always use the (pointwise) maximum norm k kD , where D is a closed and bounded set. 2. The time semidiscretization: uniform convergence and asymptotic behavior The first step to obtain the totally discrete method is based on the time discretization of the continuous problem (1). To simplify we discretize in ½0; T using a uniform mesh, which is defined by
M ¼ ft k ¼ ks; 0 6 k 6 M; s ¼ T=Mg: x M we consider the classical backward Euler method given by Then, on x
8 zðx; 0Þ ¼ 0; > : zð0; tn Þ ¼ zð1; t n Þ ¼ 0:
1 6 n 6 M;
ð4Þ
1480
C. Clavero, J.L. Gracia / Applied Mathematics and Computation 216 (2010) 1478–1488
Estimates for the global error are deduced from appropriate bounds of the local error, where the auxiliary problems
ðI þ sLx;e Þ^zðx; tn Þ ¼ sf ðx; tn Þ þ uðx; t n1 Þ; ^zð0; t n Þ ¼ ^zð1; t n Þ ¼ 0
1 6 n 6 M;
ð5Þ
are introduced to define the local error. Lemma 1 (See [4]). The local error associated to the method (4), defined by En ¼ uðx; t n Þ ^zðx; tn Þ, satisfies
kEn kX 6 C s2 ;
1 6 n 6 M:
Using that kðI þ sLx;e Þ1 kX 6 1=ð1 þ sbÞ, the uniform convergence of the backward Euler method trivially follows. Theorem 2 (See [4]). The global error associated to the method (4), defined by en ¼ uðx; tn Þ zðx; tn Þ, satisfies
ken kX 6 C s;
1 6 n 6 M;
ð6Þ
and therefore the backward Euler method is a first order uniformly convergent scheme. Corollary 1. For the exact solution of problems (4) it holds
kzðx; tn Þ zðx; t n1 ÞkX 6 C s;
1 6 n 6 M:
ð7Þ
Proof. We have
zðx; t n Þ zðx; tn1 Þ ¼ ðzðx; t n Þ uðx; tn ÞÞ þ ðuðx; tn Þ uðx; t n1 ÞÞ þ ðuðx; tn1 Þ zðx; t n1 ÞÞ: Then, using (3), (6), and uðx; tn Þ uðx; t n1 Þ ¼ sut ðx; nÞ, with n 2 ðtn1 ; t n Þ, the result trivially follows. h Previous to the construction of the numerical method and to study its uniform convergence, we give the asymptotic behavior, with respect to e, of the exact solution of the semidiscrete problems (4). In the next section, the bounds proved for the solution and their derivatives will be used in the analysis of the error of the spatial discretization. With this aim we will use the following technical results. Lemma 2. Let uðx; tÞ be the solution of (1) and ^zðx; tn Þ the solution of (5). Then, it holds
kðuðx; tn Þ ^zðx; t n ÞÞ ðuðx; t n1 Þ ^zðx; tn1 ÞÞkX 6 C s3 ;
1 6 n 6 M;
and therefore, from the uniform stability of the Euler method, it follows
kðuðx; tn Þ zðx; t n ÞÞ ðuðx; t n1 Þ zðx; tn1 ÞÞkX 6 C s2 ;
1 6 n 6 M:
Proof. It is straightforward to obtain
ðI þ sLx;e Þ½ðuðx; t n Þ ^zðx; t n ÞÞ ðuðx; t n1 Þ ^zðx; tn1 ÞÞ ¼ uðx; tn Þ 2uðx; t n1 Þ þ uðx; t n2 Þ þ sðut ðx; tn1 Þ ut ðx; t n ÞÞ ¼ Oðs3 Þ; using in the last equality Taylor expansions and the bounds for the time partial derivatives given in (3). Also, in the boundary we have
ðuð0; t n Þ ^zð0; tn ÞÞ ðuð0; tn1 Þ ^zð0; tn1 ÞÞ ¼ 0; ðuð1; tn Þ ^zð1; tn ÞÞ ðuð1; t n1 Þ ^zð1; tn1 ÞÞ ¼ 0: Then, the maximum principle for the operator ðI þ sLx;e Þ proves the result. h Lemma 3. Let zðx; t n2 Þ; zðx; t n1 Þ and zðx; tn Þ be the solution of problems (4) at the time levels n 2; n 1; and n, respectively. Then, it holds
zðx; tn Þ 2zðx; tn1 Þ þ zðx; t n2 Þ 6 C: s2 Proof. We have
zðx; t n Þ 2zðx; t n1 Þ þ zðx; t n2 Þ ¼ ðzðx; tn Þ uðx; t n ÞÞ þ ðuðx; t n Þ 2uðx; t n1 Þ þ uðx; t n2 ÞÞ þ ð2uðx; t n1 Þ 2zðx; tn1 ÞÞ þ ðzðx; tn2 Þ uðx; t n2 ÞÞ: Then using Taylor expansions, Theorem 1 and Lemma 2 the result trivially follows. h
ð8Þ
1481
C. Clavero, J.L. Gracia / Applied Mathematics and Computation 216 (2010) 1478–1488
Lemma 4. Let zðx; tn1 Þ and zðx; tn Þ be the solution of problems (4) at the time levels n 1, and n respectively. Then, it holds
00 z ðx; tn Þ z00 ðx; t n1 Þ 6 C e1 ; s
1 6 n 6 M:
ð9Þ
Proof. The semidiscrete problem at level n is given by
zðx; tn Þ zðx; tn1 Þ
s
ez00 ðx; t n Þ þ bðx; t n Þzðx; tn Þ ¼ f ðx; t n Þ;
ð10Þ
and at level n 1 is given by
zðx; tn1 Þ zðx; t n2 Þ
s
ezxx ðx; t n1 Þ þ bðx; tn1 Þzðx; t n1 Þ ¼ f ðx; t n1 Þ:
ð11Þ
Then, subtracting Eqs. (10) and (11), we obtain
zðx; tn Þ 2zðx; t n1 Þ þ zðx; t n2 Þ
s2 ¼e
z00 ðx; t n Þ z00 ðx; tn1 Þ
s
þ bðx; t n Þ
zðx; tn Þ zðx; tn1 Þ
s
f ðx; t n Þ f ðx; tn1 Þ
s
þ zðx; t n1 Þ
bðx; t n Þ bðx; t n1 Þ
s
:
Using an inductive argument and the maximum principle for problem (4) can be easily deduced that kzðx; tn ÞkX 6 C; n ¼ 0; . . . ; M. From the boundness of zðx; tn Þ, (7), (8), jf ðx; tn Þ f ðx; tn1 Þj 6 C s and jbðx; t n Þ bðx; t n1 Þj 6 C s, the result trivially follows. h Now we are ready to obtain the main result of this section, showing the behavior of the solution of the semidiscrete problems (4). Theorem 3. Let zðx; tn Þ be the solution of problem (4) at the time level n. Then, it holds
ðkÞ z ðx; t n Þ 6 C 1 þ ek=2 Be ðxÞ ;
0 6 k 6 4; 1 6 n 6 M:
ð12Þ
Proof. The proof proceeds by induction. Let n ¼ 1 be. The semidiscrete problem is given by
(
I
zðx; t 1 Þ ¼ f ðx; t 1 Þ; zð0; t1 Þ ¼ zð1; t 1 Þ ¼ 0: s þ Lx;e
ð13Þ
Then, using that jf ðx; t1 Þj 6 C and the maximum principle for the operator sI þ Lx;e , it follows that jzðx; t 1 Þj 6 C s. Moreover, problem (13) can be written as
(
ez00 ðx; t 1 Þ þ bðx; t1 Þzðx; t 1 Þ ¼ f ðx; t 1 Þ zðx;ts 1 Þ h1 ðx; sÞ; zð0; t1 Þ ¼ zð1; t 1 Þ ¼ 0;
ð14Þ
with jh1 ðx; sÞj 6 C. Then, following [11] it is possible to prove jz0 ð0; t 1 Þj 6 C e1=2 and jz0 ð1; t1 Þj 6 C e1=2 : Differentiating once w.r.t. x the differential equation in (13), we obtain
I
s
þ Lx;e zx ðx; t1 Þ ¼ fx ðx; t 1 Þ bx ðx; t1 Þzðx; t 1 Þ h2 ðx; sÞ;
with jh2 ðx; sÞj 6 C. Taking the barrier function
w1 ðxÞ ¼ C 1 þ e1=2 Be ðxÞ ; I
ð15Þ
I
we have s þ Lx;e w1 ðxÞ P C; and from the maximum principle for s þ Lx;e , it follows
jz0 ðx; t 1 Þj 6 C 1 þ e1=2 Be ðxÞ :
Moreover, from (14) we directly obtain jz00 ð0; t 1 Þj 6 C e1 and jz00 ð1; t 1 Þj 6 C e1 , and differentiating twice w.r.t. x the differential equation in (13), we obtain
I
s
þ Lx;e z00 ðx; t1 Þ ¼ fxx ðx; t 1 Þ bxx ðx; t 1 Þzðx; t 1 Þ 2bx ðx; t 1 Þz0 ðx; t 1 Þ h3 ðx; sÞ;
with jh3 ðx; sÞj 6 ð1 þ e1=2 Be ðxÞÞ. Taking now the barrier function
w2 ðxÞ ¼ C 1 þ e1 Be ðxÞ ;
ð16Þ
in the same way as before we can prove jz ðx; t1 Þj 6 Cð1 þ e Be ðxÞÞ: For the third and the fourth order derivatives, using (9) we know that 00
1
1482
C. Clavero, J.L. Gracia / Applied Mathematics and Computation 216 (2010) 1478–1488
00 z ðx; t 1 Þ z00 ðx; 0Þ 6 C e1 ; s and then a similar argument than the previous one permit to obtain the bounds
ðkÞ z ðx; t 1 Þ 6 C 1 þ ek=2 Be ðxÞ ;
k ¼ 3; 4;
which is the required result. We assume now the induction hypothesis
ðkÞ z ðx; t m Þ 6 C 1 þ ek=2 Be ðxÞ ;
k ¼ 0; 1; 2; 3; 4; 1 6 m 6 n;
and we will prove that similar bounds hold also for the time level n þ 1. We know that
(
I
s þ Lx;e
zðx; t nþ1 Þ ¼ f ðx; t nþ1 Þ þ zðx;ts n Þ g 1 ðx; sÞ;
ð17Þ
zð0; tnþ1 Þ ¼ zð1; t nþ1 Þ ¼ 0;
with jg 1 ðx; sÞj 6 C=s. From the maximum principle we deduce jzðx; t nþ1 Þj 6 C. Moreover, problem (17) can be written as
(
nþ1 Þ ez00 ðx; t nþ1 Þ þ bðx; tnþ1 Þzðx; t nþ1 Þ ¼ f ðx; t nþ1 Þ þ zðx;tn Þzðx;t g 2 ðx; sÞ; s zð0; tnþ1 Þ ¼ zð1; t nþ1 Þ ¼ 0;
ð18Þ
with jg 2 ðx; sÞj 6 C from (7). Using again the technique given in [11] we obtain jz0 ð0; t nþ1 Þj 6 C e1=2 ; jz0 ð1; t nþ1 Þj 6 C e1=2 ; and jz00 ð0; t nþ1 Þj 6 C e1 ; jz00 ð1; tnþ1 Þj 6 C e1 directly follows from (18). Differentiating w.r.t. x the differential equation in (17), we obtain
I
s
z0 ðx; tn Þ þ Lx;e z0 ðx; tnþ1 Þ ¼ fx ðx; tnþ1 Þ bx ðx; tnþ1 Þzðx; t nþ1 Þ þ g 3 ðx; sÞ:
s
From the hypothesis of induction, it holds jg 3 ðx; sÞj 6 s ð1 þ e Be ðxÞÞ: Taking again the barrier functions w1 ðxÞ and w2 ðxÞ given in (15) and (16) respectively, and using the maximum principle for sI þ Lx;e we deduce C
jz0 ðx; tnþ1 Þj 6 C 1 þ e1=2 Be ðxÞ ;
1=2
jz00 ðx; tnþ1 Þj 6 C 1 þ e1 Be ðxÞ :
For the third and fourth order derivatives we proceed as follows. Differentiating twice w.r.t. x the differential equation in (18), we obtain
8 z00 ðx;t n Þz00 ðx;t nþ1 Þ ð4Þ 00 > < ez ðx; t nþ1 Þ þ bðx; t nþ1 Þz ðx; tnþ1 Þ ¼ fxx ðx; tnþ1 Þ þ s 0 bxx ðx; tnþ1 Þzðx; tnþ1 Þ 2bx ðx; t nþ1 Þz ðx; tnþ1 Þ g 4 ðx; sÞ; > : 00 z ð0; t nþ1 Þ; z00 ð1; t nþ1 Þ given:
ð19Þ
From the bound for jz0 ðx; t nþ1 Þj and using (7) we have that jg 4 ðx; sÞj 6 C e1 , and following [11] we deduce jz000 ð0; t nþ1 Þj 6 C e3=2 ; jz000 ð1; tnþ1 Þj 6 C e3=2 . Moreover, from (19) it directly follows jzð4Þ ð0; tnþ1 Þj 6 C e2 ; jzð4Þ ð1; t nþ1 Þj 6 C e2 : Now, differentiating three times w.r.t. x the differential equation in (17), we obtain
I
s
z000 ðx; t n Þ þ Lx;e z000 ðx; t nþ1 Þ ¼ fxxx ðx; t nþ1 Þ bxxx ðx; t nþ1 Þzðx; t nþ1 Þ 3bxx ðx; tnþ1 Þz0 ðx; t nþ1 Þ 3bx ðx; t nþ1 Þz00 ðx; t nþ1 Þ þ
s
g 5 ðx; sÞ: From the hypothesis of induction and the bounds for jzðkÞ ðx; t nþ1 Þj; k ¼ 0; 1; 2; it follows jg 5ðx; sÞj 6 Cs ð1 þ e3=2 Be ðxÞÞ: Then, using the barrier function w3 ðxÞ ¼ Cð1 þ e3=2 Be ðxÞÞ, and the maximum principle for sI þ Lx;e we obtain
jz000 ðx; tnþ1 Þj 6 C 1 þ e3=2 Be ðxÞ : In a completely similar way, it can be proved that jzð4Þ ðx; tnþ1 Þj 6 C 1 þ e2 Be ðxÞ , which is the required result.
h
For the analysis of the uniform convergence of the totally discrete method, we write problem (4) in vector form LM x;e z ¼ f, where
z ¼ ðzðx; t1 Þ; . . . ; zðx; tM ÞÞT ;
LM x;e ¼
h
i1 h iM T LM ; . . . ; LM ; x;e x;e
with
h in zðx; t Þ zðx; t Þ n n1 LM þ Lx;e zðx; tn Þ ¼ f ðx; t n Þ; x;e z
s
1 6 n 6 M:
The initial condition is zðx; 0Þ ¼ 0; x 2 ð0; 1Þ; and the boundary conditions are zð0; tn Þ ¼ zð1; t n Þ ¼ 0; n ¼ 0; . . . ; M:
ð20Þ
C. Clavero, J.L. Gracia / Applied Mathematics and Computation 216 (2010) 1478–1488
1483
3. The totally discrete method: uniform convergence M XN To discretize in space we consider a nonuniform mesh XN ¼ f0 ¼ x0 < < xN ¼ 1g, and we denote by Q M;N ¼ x M;N M;N the corresponding grid for the x; t-variables, and Q ¼Q \ Q. On XN we use central differences to approximate problem (20), which is given by
h
in U n U n1 j j LM;N U þ LNx;e U nj ¼ f ðxj ; t n Þ; e
s
j
1 6 n 6 M; 0 < j < N;
ð21Þ
where
LNx;e U nj ed2x U nj þ bðxj ; t n ÞU nj ; and d2x is the classical central difference approximation on a nonuniform mesh for the second order derivative. The initial condition is given by U 0j ¼ 0; j ¼ 0; . . . ; N, and the boundary conditions are U n0 ¼ U nN ¼ 0; 1 6 n 6 M. The following discrete maximum principle holds for this scheme. Lemma 5. Let LM;N be the discrete operator given in (21) and U ¼ U0 ; . . . ; UM g is a vectorial function with Un ¼ fU n0 ; . . . ; U nN g. If e n M;N 0 n n , then U P 0: U P 0; U 0 P 0; U N P 0; 1 6 n 6 M, and ½LM;N e Uj P 0 in Q From Lemma 5 it follows that the numerical scheme is uniformly stable in the maximum norm, and also the following comparison principle holds. n n M;N Corollary 2. If V, W are two vectorial functions such that V0 P W0 ; V n0 P W n0 ; V nN P W nN ; 1 6 n 6 M, and ½LM;N e Vj P ½Le Wj M;N in Q , then V P W.
To prove the uniform convergence of the numerical solution of the scheme (21), three different special nonuniform meshes, adapted to singularly perturbed reaction–diffusion problems, are used. The Shishkin mesh [11,12]. It is a piecewise uniform mesh, depending on two transition points which are defined by means of the transition parameter
r ¼ min
pffiffiffi 1 ; r0 e ln N ; 4
ð22Þ
where r0 is a constant to be fixed later. A uniform mesh is placed in ½0; r; ½r; 1 r, and ½1 r; 1, such that x0 ¼ 0; xN=4 ¼ r; x3N=4 ¼ 1 r, and xN ¼ 1 and therefore the mesh points are given by
8 j ¼ 0; 1; . . . ; N=4; > < 4jr=N; xj ¼ r þ 2ðj N=4Þð1 2rÞ=N; j ¼ N=4 þ 1; . . . ; 3N=4; > : j ¼ 3N=4 þ 1; . . . ; N: 1 r þ 4ðj 3N=4Þr=N;
ð23Þ
The Vulanovic mesh [13]. This mesh is a generalized Shishkin mesh constructed by using a suitable generating function @, which also depends on two transition points. To simplify, we use the same parameter r given in (22). The grid points are defined by xj ¼ @ðj=NÞ; j ¼ 0; 1; . . . ; N=2, with @ 2 C 2 ½0; 1=2 and
(
4r=N;
@ðzÞ ¼
z 2 ½0; 1=4; 3
pðz 1=4Þ þ 4rðz 1=4Þ þ r; z 2 ½1=4; 1=2:
ð24Þ
The coefficient p is such that @ð1=2Þ ¼ 1=2 and the mesh is symmetric with respect to the point x ¼ 1=2. Note that in ½0; r and ½1 r; 1 the mesh points are the same than in the Shishkin mesh. Otherwise, in ½r; 1 r it is nonuniform but the step sizes satisfy
jhjþ1 hj j 6 CN2 ;
j ¼ N=4; . . . ; 3N=4:
ð25Þ
The Bahkvalov mesh [7]. This mesh is based on user–chosen mesh parameters r0 > 0 and q 2 ð0; 1=2Þ. To simplify, we take pffiffiffi pffiffiffi q ¼ 1=4. The mesh points are xj ¼ j=N; j ¼ 0; 1; . . . ; N if r0 e < q, while when r0 e P q they are given by
xj ¼
uðj=NÞ; for j 6 N=2; 1 uððN jÞ=NÞ; for j > N=2;
ð26Þ
where
(
uðzÞ ¼
pffiffiffi
vðzÞ r0 e lnð1 4zÞ; z 2 ½0; m; pðzÞ vðmÞ þ v0 ðmÞðz mÞ; z 2 ½m; 1=2;
and the transition point 0
m is the solution of the nonlinear equation
ð1 2mÞv ðmÞ ¼ 1 2vðmÞ:
ð27Þ
1484
C. Clavero, J.L. Gracia / Applied Mathematics and Computation 216 (2010) 1478–1488
This defines the mesh in [0,1/2] and the mesh is again symmetric with respect to the central point x ¼ 1=2. Note that in ½0; m and ½1 m; 1 the mesh is graded and in ½m; 1 m the mesh is uniform. Theorem 4. Let Un be the numerical solution of (21) and zðx; t n Þ be the solution of (4) both at time level t n . We assume that the pffiffiffi constant r0 used to define the Shishkin, Vulanovic and Bahkvalov meshes satisfies r0 P 2= b. Then, if the method is constructed on the Shishkin or the Vulanovic meshes, the error satisfies
n
U j zðxj ; t n Þ
XN
6 CðN1 ln NÞ2 ;
ð28Þ
and therefore the space discretization is an almost second order uniformly convergent scheme. Otherwise, if the method is constructed on the Bahkvalov mesh the error satisfies
n
U j zðxj ; t n Þ
XN
6 CN2 ;
ð29Þ
and therefore the space discretization is a second order uniformly convergent scheme. Proof. We begin with the Shishkin mesh. Using similar ideas as in [4], we can obtain that the local error satisfies
( 1 h 2 in M;N L ðU zÞ ¼ Lx;e LN zðxj ; tn Þ 6 CN e þ CN ; if xj ¼ r; 1 r; x;e e 2 1 j otherwise: CðN ln NÞ ; Note that from the uniform stability of the scheme we cannot deduce the almost second order of uniform convergence. Nevertheless, using a standard barrier function in the context of singularly perturbed problems given by ½Wnj ¼ CðN 2 hðxj Þ ln N þ ðN 1 ln NÞ2 Þ; where
8x > : 1x r
ifx 2 ½0; r; if x 2 ½r; 1 r; ; ifx 2 ½1 r; 1;
h in M;N n L ðU zÞ ; and then from Lemma 2 the result follows. it holds LM;N ½W P e j e j In second place, we consider the Vulanovic mesh. Using that in ð0; rÞ [ ð1 r; 1Þ the mesh is the same than for the Shishkin mesh, and in ½r; 1 r the step sizes satisfy (25), it is not difficult to deduce that
( 2 h in M;N ifxj 2 ½r; 1 r; L ðU zÞ 6 CN ; e j CðN1 ln NÞ2 ; otherwise; and therefore the result follows from the uniform stability of the scheme. Finally, for the Bahkvalov mesh the same analysis than this one given in [7] (see Appendix, p. 2094) permits to deduce that the local error satisfies
h in M;N L ðU zÞ 6 CN 2 ; e j and again the result follows from the uniform stability of the scheme. h Theorem 5. Let U be the numerical solution of (21), u the solution of (1), and constructed on the Shishkin or the Vulanovic meshes, the error satisfies
n
U j uðxj ; tn Þ
Q M;N
6C
pffiffiffi
r0 such that r0 P 2= b. Then, if the method is
s þ ðN1 ln NÞ2 :
ð30Þ
Otherwise, if the method is constructed on the Bahkvalov mesh, the error satisfies
n
U j uðxj ; tn Þ
Q M;N
6 Cðs þ N2 Þ:
ð31Þ
Proof. The results follows from Theorems 2 and 4. h Remark 1. Note that in Theorem 5 there is not any relation between the discretization parameter N and s in contrast with previous results in [1,5], where the theoretical restriction N q 6 C s with 0 < q < 1 appears. This fact accords with our numerical experiments, where that restriction between the two discretization parameters has been never needed. Moreover, the proof of the uniform convergence here is considerably simpler than those ones in previous papers because now we have not needed any recursive argument invoking the uniform stability of any transition operator.
C. Clavero, J.L. Gracia / Applied Mathematics and Computation 216 (2010) 1478–1488
1485
4. Numerical results The first example that we consider is given by
8 1þx2 t > < ut euxx þ 2 u ¼ e 1 þ sinðpxÞ; ðx; tÞ 2 ð0; 1Þ ð0; 1; uðx; 0Þ ¼ 0; x 2 X; > : uð0; tÞ ¼ uð1; tÞ ¼ 0; t 2 ð0; 1;
ð32Þ
which solution is unknown. In Fig. 1 we show the numerical solution of this problem for a particular value of e ¼ 104 at two different times. From it we clearly see the boundary layers at the end points of the spatial domain. To estimate the numerical errors we use a variant of the double mesh principle and therefore at times t n ¼ ns for each point xj ; j ¼ 0; 1; . . . ; N the error is approximated by using
e;2N;s Dej;n;N;s ¼ U ej;n;N;s U j;n 4 ;
e;2N;s
where U ej;n;N;s is the numerical solution obtained using a constant time step s and ðN þ 1Þ points in the spatial mesh and U j;n 4 is computed using 4s as time step and ð2N þ 1Þ points in the spatial mesh but with the same transition parameters than in the original mesh. For each fixed value of e, the maximum global errors are estimated by
De;N;s ¼ max Dej;n;N;s ; j;n
and therefore, in standard way, the numerical orders of convergence are given by
q¼
logðDe;N;s =De;2N;s=4 Þ : log 2
From these values we obtain the e-uniform errors and the e-uniform orders of convergence by
DN;s ¼ max De;N;s ; e
quni ¼
logðDN;s =D2N;s=4 Þ : log 2
ð33Þ
Tables 1 and 2 display the results on the Shishkin, Vulanovic and Bahkvalov meshes. In all cases we have taken the constant r0 ¼ 2. Note that the results on Shishkin and Vulanovic meshes are identical and for this reason they appear together in Table 1. The maximum errors appear in the boundary layer regions where both meshes are exactly the same, showing the almost second order uniform convergence. Table 2 shows that the method on the Bahkvalov mesh has second order convergence according with the theoretical results. The second test problem that we consider has a solution with a really small variation in time (see [1]); the reason to take this example is that we wish that the error associated to the spatial discretization stage dominates in the global error. In the previous example to see the influence of the spatial discretization in the global error we have needed to divide the time step
1.4
t=0.5 t=1
1.2 1 0.8 0.6 0.4 0.2 0
0
0.1
0.2
0.3
0.4
0.5
0.6
0.7
Fig. 1. Numerical solution at different times of (32) for
0.8
0.9
1
e ¼ 104 with N ¼ 32.
1486
C. Clavero, J.L. Gracia / Applied Mathematics and Computation 216 (2010) 1478–1488
Table 1 Maximum errors and uniform orders of convergence on the Shishkin and Vulanovic meshes for problem (32).
e
N=32 s ¼ 0:1
N=64
N=128
N=256
N=512
N=1024
s ¼ 0:1=4
s ¼ 0:1=42
s ¼ 0:1=43
s ¼ 0:1=44
s ¼ 0:1=45
26
0.179E1 1.983
0.453E2 1.996
0.113E2 1.999
0.284E3 2.000
0.710E4 2.000
0.177E4
28
0.221E1 1.978
0.562E2 1.994
0.141E2 1.999
0.353E3 2.000
0.883E4 2.000
0.221E4
210
0.246E1 1.980
0.624E2 1.994
0.157E2 1.998
0.392E3 2.000
0.980E4 2.000
0.245E4
212
0.263E1 1.986
0.664E2 1.996
0.166E2 1.757
0.493E3 1.592
0.163E3 1.642
0.524E4
214
0.272E1 1.985
0.687E2 1.997
0.172E2 1.810
0.490E3 1.591
0.163E3 1.641
0.522E4
216
0.277E1 1.986
0.698E2 1.997
0.175E2 1.838
0.489E3 1.591
0.162E3 1.640
0.521E4
218
0.279E1 1.986
0.704E2 1.997
0.176E2 1.852
0.489E3 1.591
0.162E3 1.640
0.521E4
220
0.280E1 1.987
0.707E2 1.997
0.177E2 1.859
0.489E3 1.591
0.162E3 1.640
0.521E4
222
0.281E1 1.987 ...
0.709E2 1.997 ...
0.178E2 1.863 ...
0.489E3 1.591 ...
0.162E3 1.640 ...
0.520E4
230
0.282E1 1.987
0.711E2 1.997
0.178E2 1.866
0.488E3 1.591
0.162E3 1.640
0.520E4
DN;s quni
0.282E1
0.711E2
0.178E2
0.493E3
0.163E3
0.524E4
1.987
1.997
1.854
1.592
1.642
...
...
Table 2 Maximum errors and uniform orders of convergence on the Bahkvalov mesh for problem (32).
e
N=32 s ¼ 0:1
N=64
N=128
N=256
N=512
N=1024
s ¼ 0:1=4
s ¼ 0:1=42
s ¼ 0:1=43
s ¼ 0:1=44
s ¼ 0:1=45
26
0.179E1 1.983
0.453E2 1.996
0.113E2 1.999
0.284E3 2.000
0.710E4 2.000
0.177E4
28
0.223E01 1.981
0.565E02 1.995
0.142E02 1.999
0.355E03 2.000
0.887E04 2.000
0.222E04
210
0.248E01 1.975
0.632E02 1.995
0.158E02 1.998
0.397E03 2.000
0.992E04 2.000
0.248E04
212
0.265E01 1.986
0.668E02 1.994
0.168E02 1.999
0.420E03 2.000
0.105E03 2.000
0.262E04
214
0.273E01 1.988
0.689E02 1.997
0.173E02 1.999
0.432E03 2.000
0.108E03 2.000
0.270E04
216
0.277E01 1.989
0.698E02 1.996
0.175E02 1.999
0.438E03 2.000
0.110E03 2.000
0.274E04
218
0.279E01 1.987
0.704E02 1.996
0.177E02 1.999
0.442E03 2.000
0.110E03 2.000
0.276E04
220
0.280E01 1.986
0.708E02 1.997
0.177E02 1.999
0.443E03 2.000
0.111E03 2.000
0.277E04
222
0.281E01 1.985 ...
0.709E02 1.997 ...
0.178E02 1.999 ...
0.444E03 2.000 ...
0.111E03 2.000 ...
0.278E04
30
2
0.282E01 1.987
0.711E02 1.997
0.178E02 1.999
0.445E03 2.000
0.111E03 2.000
0.278E04
DN;s quni
0.282E01
0.711E02
0.178E02
0.445E03
0.111E03
0.278E04
1.987
1.997
1.999
2.000
2.000
...
...
size by four and therefore the computational cost of the method is higher. In this new example the time step size is only divided by two. The problem is defined as
1487
C. Clavero, J.L. Gracia / Applied Mathematics and Computation 216 (2010) 1478–1488
8 > < ut euxx þ ð10 þ 0:1 sinðpxÞÞu ¼ f ðx; tÞ; ðx; tÞ 2 ð0; 1Þ ð0; 1; uðx; 0Þ ¼ 0; x 2 X; > : t 2 ð0; 1; uð0; tÞ ¼ g 1 ðtÞ; uð1; tÞ ¼ g 2 ðtÞ;
ð34Þ
where f ; g 1 and g 2 are such that the exact solution is given by
pffiffiffiffiffiffiffiffiffiffiffi uðx; tÞ ¼ expð 10=exÞð1 expðt=100ÞÞ: Note that the boundary conditions are not homogeneous, but a simple change of variable transforms the example (34) in a problem of type (1). Now the errors at t n ¼ ns; xj ; j ¼ 0; 1; . . . ; N and the orders of convergence are calculated by
Eej;n;N;s ¼ U ej;n;N;s uðxj ; t n Þ; the maximum global errors and the numerical orders of convergence are given by
Ee;N;s ¼ max Eej;n;N;s ; j;n
p¼
logðEe;N;s =Ee;2N;s=2 Þ : log 2
From these values the e-uniform errors and the e-uniform orders of convergence are obtained by
EN;s ¼ max Ee;N;s ; e
puni ¼
logðEN;s =E2N;s=2 Þ : log 2
ð35Þ
Again we take r0 ¼ 2, but we have obtained similar results for smaller values of this parameter. Also, the maximum errors on the Shishkin and Vulanovic meshes are identical for problem (34). Tables 3 and 4 display the results on the Shishkin/Vula-
Table 3 Maximum errors and uniform orders of convergence on the Shishkin and Vulanovic meshes for problem (34).
e
N=32 s ¼ 0:1
N=64
N=128
N=256
N=512
N=1024
s ¼ 0:1=2
s ¼ 0:1=22
s ¼ 0:1=23
s ¼ 0:1=24
s ¼ 0:1=25
26
0.836E04 1.937
0.218E04 1.958
0.561E05 2.002
0.140E05 2.011
0.348E06 2.023
0.856E07
28
0.267E03 1.673
0.836E04 1.936
0.219E04 1.957
0.563E05 1.999
0.141E05 2.005
0.351E06
210
0.382E03 0.517
0.267E03 1.673
0.837E04 1.936
0.219E04 1.956
0.563E05 1.997
0.141E05
212
0.382E03 0.454
0.279E03 1.183
0.123E03 1.547
0.421E04 1.626
0.136E04 1.687
0.423E05
...
...
...
...
...
...
...
230
0.382E03 0.454
0.279E03 1.183
0.123E03 1.547
0.421E04 1.626
0.136E04 1.687
0.423E05
EN;s puni
0.382E03
0.279E03
0.123E03
0.421E04
0.136E04
0.423E05
0.454
1.183
1.547
1.626
1.687
Table 4 Maximum errors and uniform orders of convergence on the Bahkvalov mesh for problem (34).
e
N=32 s ¼ 0:1
N=64
N=128
N=256
N=512
N=1024
s ¼ 0:1=2
s ¼ 0:1=22
s ¼ 0:1=23
s ¼ 0:1=24
s ¼ 0:1=25
26
0.836E04 1.937
0.218E04 1.958
0.561E05 2.002
0.140E05 2.011
0.348E06 2.023
0.856E07
28
0.429E04 1.901
0.115E04 1.990
0.289E05 2.003
0.721E06 2.022
0.178E06 2.046
0.430E07
210
0.428E04 1.900
0.115E04 1.990
0.289E05 2.004
0.721E06 2.022
0.177E06 2.046
0.430E07
212
0.428E04 1.900 ...
0.115E04 1.990 ...
0.289E05 2.004 ...
0.721E06 2.022 ...
0.177E06 2.046 ...
0.430E07
2
0.428E04 1.900
0.115E04 1.990
0.289E05 2.004
0.721E06 2.022
0.177E06 2.046
0.430E07
EN;s puni
0.836E04
0.218E04
0.561E05
0.140E05
0.348E06
0.856E07
1.937
1.958
2.002
2.011
2.023
... 30
...
1488
C. Clavero, J.L. Gracia / Applied Mathematics and Computation 216 (2010) 1478–1488
novic and Bahkvalov meshes respectively. From these tables we observe again that the method has almost second order on the Shishkin and Vulanovic meshes and second order on the Bahkvalov mesh, in agreement with the theoretical results. Acknowledgments This research was partially supported by the project MEC/FEDER MTM2007-63204 and the Diputación General de Aragón. References [1] B. Bujanda, C. Clavero, J.L. Gracia, J.C. Jorge, A high order uniformly convergent alternating direction scheme for time dependent reaction–diffusion singularly perturbed problems, Num. Math. 107 (2007) 1–25. [2] C. Clavero, J.L. Gracia, High order methods for elliptic and time dependent reaction–diffusion singularly perturbed problems, Appl. Math. Comp. 168 (2005) 1109–1127. [3] C. Clavero, J.L. Gracia, J.C. Jorge, A uniformly convergent alternating direction HODIE finite difference scheme for 2D time dependent convection– diffusion problems, IMA J. Numer. Anal. 26 (2006) 155–172. [4] C. Clavero, J.C. Jorge, F. Lisbona, G.I. Shishkin, An alternating direction scheme on a nonuniform mesh for reaction–diffusion parabolic problem, IMA J. Numer. Anal. 20 (2000) 263–280. [5] J.L. Gracia, F. Lisbona, A uniformly convergent scheme for a system of reaction–diffusion equations, J. Comput. Appl. Math. 206 (2007) 1–16. [6] P.W. Hemker, G.I. Shishkin, L.P. Shiskina, e-uniform schemes with high-order time accuracy for parabolic singular perturbation problems, IMA J. Numer. Anal. 20 (2000) 99–121. [7] B. Kellogg, T. Linß, M. Stynes, A finite difference method on layer-adapted meshes for an elliptic reaction–diffusion system in two dimensions, Math. Comp. 77 (2008) 2085–2096. [8] T. Linß, Layer-adapted meshes and FEM for time-dependent singularly perturbed reaction-diffusion problems, Int. J. Computing Science and Math. 1 (2007) 259–270. [9] T. Linß, A robust method for a time-dependent system of coupled singularly perturbed reaction–diffusion problems, preprint MATH-NM-03-2008, Institut für Numerische Mathematik, TU Dresden. [10] T. Linß, N. Madden, Parameter uniform approximations for time-dependent reaction-diffusion problems, Numer. Methods Partial Differ. Equ. 23 (2007) 1290–1300. [11] J.J.H. Miller, E. O’Riordan, G.I. Shishkin, Fitted Numerical Methods for Singular Perturbation Problems, World-Scientific, Singapore, 1996. [12] H.-G. Roos, M. Stynes, L. Tobiska, Robust numerical methods for singularly perturbed differential equations, Springer Series in Computational Mathematics 24, Springer-Verlag, Berlin, 2008. [13] R. Vulanovic, A high order scheme for quasilinear boundary value problems with two small parameters, Computing 67 (1986) 287–303.
Contents lists available at ScienceDirect
Applied Mathematics and Computation journal homepage: www.elsevier.com/locate/amc
On the uniform convergence of a finite difference scheme for time dependent singularly perturbed reaction-diffusion problems C. Clavero *, J.L. Gracia Department of Applied Mathematics, University of Zaragoza, 50018 Zaragoza, Spain
a r t i c l e
i n f o
a b s t r a c t
Keywords: Parabolic reaction–diffusion problems Semidiscrete problems Asymptotic behavior Uniform convergence Special nonuniform mesh
In this work we are interested in the numerical approximation of 1D parabolic singularly perturbed problems of reaction–diffusion type. To approximate the multiscale solution of this problem we use a numerical scheme combining the classical backward Euler method and central differencing. The scheme is defined on some special meshes which are the tensor product of a uniform mesh in time and a special mesh in space, condensing the mesh points in the boundary layer regions. In this paper three different meshes of Shishkin, Bahkvalov and Vulanovic type are used, proving the uniform convergence with respect to the diffusion parameter. The analysis of the uniform convergence is based on a new study of the asymptotic behavior of the solution of the semidiscrete problems, which are obtained after the time discretization by the Euler method. Some numerical results are showed corroborating in practice the theoretical results on the uniform convergence and the order of the method. Ó 2010 Elsevier Inc. All rights reserved.
1. Introduction In this paper we consider 1D initial-boundary value parabolic reaction–diffusion singularly perturbed problems given by
8 ðx; tÞ 2 Q ¼ X ð0; T ð0; 1Þ ð0; T; > < ut þ Lx;e u ¼ f ðx; tÞ; uðx; 0Þ ¼ 0; x 2 X; > : uð0; tÞ ¼ uð1; tÞ ¼ 0; t 2 ð0; T;
ð1Þ
where the spatial differential operator is defined as
Lx;e u euxx þ bðx; tÞu:
ð2Þ
We assume that 0 < e 6 1 and it can be arbitrarily small, the reaction term satisfies b P b > 0 for all x 2 Q (the forthcoming analysis can be easily extended to the case b P b P 0 using a standard change of variable), and the data are sufficiently smooth functions b; f 2 C ð4;2Þ ðQ Þ. Moreover we assume compatibility conditions between data of problem (1) in order that the exact solution u 2 C ð4;2Þ ðQ Þ. Sufficient conditions for this regularity are given by
f ð0; 0Þ ¼ f ð1; 0Þ ¼ 0; ðLx;e f Þð0; 0Þ ¼ ft ð0; 0Þ;
ðLx;e f Þð1; 0Þ ¼ ft ð1; 0Þ:
* Corresponding author. E-mail addresses: [email protected] (C. Clavero), [email protected] (J.L. Gracia). 0096-3003/$ - see front matter Ó 2010 Elsevier Inc. All rights reserved. doi:10.1016/j.amc.2010.02.050
C. Clavero, J.L. Gracia / Applied Mathematics and Computation 216 (2010) 1478–1488
1479
pffiffiffi It is well-known (see [6],[12]) that the solution of (1) has a boundary layer at x ¼ 0; 1 of a width O ej ln ej . In the next result we remember the asymptotic behavior with respect to e of that solution showing the presence of two boundary layers. Theorem 1 (See [6]). The solution of (1) satisfies
juðk;mÞ ðx; tÞj 6 C 1 þ ek=2 Be ðxÞ ;
0 6 k þ 2m 6 4;
ð3Þ
where the function Be ðxÞ is given by
Be ðxÞ ¼ e
pffiffiffiffiffi
b=ex
þ e
pffiffiffiffiffi
b=eð1xÞ
:
This behavior motivates the use of uniformly convergent methods (see [11,12]) to obtain efficient numerical approximations for any value of e. There are many papers in the literature dealing with the numerical approximation of the solution for time dependent singularly perturbed reaction–diffusion problems. In [2,8,10] 1D problems like (1) are considered, in [1,4] the case of 2D parabolic problems is analyzed and coupled linear systems are solved in [5,9]. Different techniques have been used to prove the uniform convergence of the numerical method defined for each problem. In some papers the totally discrete method is directly considered and the analysis of the convergence is based on appropriate bounds for the truncation error and the discrete maximum principle for the totally discrete operator (see [9]). Nevertheless, this technique cannot be extended to higher order time approximation than the Euler method due to the lack of a discrete maximum principle for the totally discrete method. To elude this difficulty, in other papers (see [1–4]) the analysis of the uniform convergence is made in two steps. In the first one, the problem is discretized in time and appropriate estimates for the corresponding local error are proved. In the second step, the set of one dimensional problems resulting from the time discretization are discretized in space and it is necessary to prove the uniform convergence of the corresponding spatial discretization. So far, and motivated by the classical analysis of the convergence for initial value problems, in the previous papers this analysis is based on the uniform convergence of the numerical scheme which discretizes some auxiliary problems arising in the definition of the method. To deduce that uniform convergence, the asymptotic behavior of the solution of these auxiliary problems was needed. The proof of the uniform convergence of the totally discrete schemes combining the two steps, uses a recursive argument which needs that the powers of a discrete transition operator arising in the analysis be uniformly bounded with respect to e. To our knowledge this question is only completely solved in the case that the backward Euler method is considered and it is an open question for higher order time approximations. The idea of this paper is to show a new strategy which can be used to overcome the main difficulties associated to the second technique. With this purpose, the uniform convergence of the numerical method considered to discretize the singularly perturbed problems resulting from the time discretization, will be proved without any use of auxiliary problems. This new proof requires to know the asymptotic behavior of the solution of these problems, which is based on an inductive argument and the well-known behavior of steady singularly perturbed reaction–diffusion problems (see [11,12]). In this new analysis, it is not necessary any recursive argument and therefore the stability of any transition operator, making simpler the analysis of the uniform convergence. The paper is structured as follows. In Section 2 we remember that the backward Euler method defined on an equidistant mesh is a first order uniformly convergent method. Moreover, the asymptotic behavior of the solution of the semidiscrete problems resulting after the discretization in time is proved, which will be crucial for the posterior analysis of the uniform convergence of the numerical scheme. In Section 3 central differencing is used to discretize the previous semidiscrete problems. The scheme is constructed on different nonuniform meshes; concretely we use a piecewise uniform Shishkin mesh (see [12]), a Bahkvalov mesh (see [7]) and a mesh based on a generating function of polynomial type (see [13]). It is proved that the scheme has almost second order uniform convergence on the Shishkin and the Vulanovic meshes and it has second order on the Bahkvalov mesh. In Section 4 some numerical results for two different examples are given illustrating in practice the efficiency of the method and corroborating the order of uniform convergence theoretically proved. Henceforth, C denotes a generic positive constant independent of the diffusion parameter e and also of the discretization parameters N and s. We always use the (pointwise) maximum norm k kD , where D is a closed and bounded set. 2. The time semidiscretization: uniform convergence and asymptotic behavior The first step to obtain the totally discrete method is based on the time discretization of the continuous problem (1). To simplify we discretize in ½0; T using a uniform mesh, which is defined by
M ¼ ft k ¼ ks; 0 6 k 6 M; s ¼ T=Mg: x M we consider the classical backward Euler method given by Then, on x
8 zðx; 0Þ ¼ 0; > : zð0; tn Þ ¼ zð1; t n Þ ¼ 0:
1 6 n 6 M;
ð4Þ
1480
C. Clavero, J.L. Gracia / Applied Mathematics and Computation 216 (2010) 1478–1488
Estimates for the global error are deduced from appropriate bounds of the local error, where the auxiliary problems
ðI þ sLx;e Þ^zðx; tn Þ ¼ sf ðx; tn Þ þ uðx; t n1 Þ; ^zð0; t n Þ ¼ ^zð1; t n Þ ¼ 0
1 6 n 6 M;
ð5Þ
are introduced to define the local error. Lemma 1 (See [4]). The local error associated to the method (4), defined by En ¼ uðx; t n Þ ^zðx; tn Þ, satisfies
kEn kX 6 C s2 ;
1 6 n 6 M:
Using that kðI þ sLx;e Þ1 kX 6 1=ð1 þ sbÞ, the uniform convergence of the backward Euler method trivially follows. Theorem 2 (See [4]). The global error associated to the method (4), defined by en ¼ uðx; tn Þ zðx; tn Þ, satisfies
ken kX 6 C s;
1 6 n 6 M;
ð6Þ
and therefore the backward Euler method is a first order uniformly convergent scheme. Corollary 1. For the exact solution of problems (4) it holds
kzðx; tn Þ zðx; t n1 ÞkX 6 C s;
1 6 n 6 M:
ð7Þ
Proof. We have
zðx; t n Þ zðx; tn1 Þ ¼ ðzðx; t n Þ uðx; tn ÞÞ þ ðuðx; tn Þ uðx; t n1 ÞÞ þ ðuðx; tn1 Þ zðx; t n1 ÞÞ: Then, using (3), (6), and uðx; tn Þ uðx; t n1 Þ ¼ sut ðx; nÞ, with n 2 ðtn1 ; t n Þ, the result trivially follows. h Previous to the construction of the numerical method and to study its uniform convergence, we give the asymptotic behavior, with respect to e, of the exact solution of the semidiscrete problems (4). In the next section, the bounds proved for the solution and their derivatives will be used in the analysis of the error of the spatial discretization. With this aim we will use the following technical results. Lemma 2. Let uðx; tÞ be the solution of (1) and ^zðx; tn Þ the solution of (5). Then, it holds
kðuðx; tn Þ ^zðx; t n ÞÞ ðuðx; t n1 Þ ^zðx; tn1 ÞÞkX 6 C s3 ;
1 6 n 6 M;
and therefore, from the uniform stability of the Euler method, it follows
kðuðx; tn Þ zðx; t n ÞÞ ðuðx; t n1 Þ zðx; tn1 ÞÞkX 6 C s2 ;
1 6 n 6 M:
Proof. It is straightforward to obtain
ðI þ sLx;e Þ½ðuðx; t n Þ ^zðx; t n ÞÞ ðuðx; t n1 Þ ^zðx; tn1 ÞÞ ¼ uðx; tn Þ 2uðx; t n1 Þ þ uðx; t n2 Þ þ sðut ðx; tn1 Þ ut ðx; t n ÞÞ ¼ Oðs3 Þ; using in the last equality Taylor expansions and the bounds for the time partial derivatives given in (3). Also, in the boundary we have
ðuð0; t n Þ ^zð0; tn ÞÞ ðuð0; tn1 Þ ^zð0; tn1 ÞÞ ¼ 0; ðuð1; tn Þ ^zð1; tn ÞÞ ðuð1; t n1 Þ ^zð1; tn1 ÞÞ ¼ 0: Then, the maximum principle for the operator ðI þ sLx;e Þ proves the result. h Lemma 3. Let zðx; t n2 Þ; zðx; t n1 Þ and zðx; tn Þ be the solution of problems (4) at the time levels n 2; n 1; and n, respectively. Then, it holds
zðx; tn Þ 2zðx; tn1 Þ þ zðx; t n2 Þ 6 C: s2 Proof. We have
zðx; t n Þ 2zðx; t n1 Þ þ zðx; t n2 Þ ¼ ðzðx; tn Þ uðx; t n ÞÞ þ ðuðx; t n Þ 2uðx; t n1 Þ þ uðx; t n2 ÞÞ þ ð2uðx; t n1 Þ 2zðx; tn1 ÞÞ þ ðzðx; tn2 Þ uðx; t n2 ÞÞ: Then using Taylor expansions, Theorem 1 and Lemma 2 the result trivially follows. h
ð8Þ
1481
C. Clavero, J.L. Gracia / Applied Mathematics and Computation 216 (2010) 1478–1488
Lemma 4. Let zðx; tn1 Þ and zðx; tn Þ be the solution of problems (4) at the time levels n 1, and n respectively. Then, it holds
00 z ðx; tn Þ z00 ðx; t n1 Þ 6 C e1 ; s
1 6 n 6 M:
ð9Þ
Proof. The semidiscrete problem at level n is given by
zðx; tn Þ zðx; tn1 Þ
s
ez00 ðx; t n Þ þ bðx; t n Þzðx; tn Þ ¼ f ðx; t n Þ;
ð10Þ
and at level n 1 is given by
zðx; tn1 Þ zðx; t n2 Þ
s
ezxx ðx; t n1 Þ þ bðx; tn1 Þzðx; t n1 Þ ¼ f ðx; t n1 Þ:
ð11Þ
Then, subtracting Eqs. (10) and (11), we obtain
zðx; tn Þ 2zðx; t n1 Þ þ zðx; t n2 Þ
s2 ¼e
z00 ðx; t n Þ z00 ðx; tn1 Þ
s
þ bðx; t n Þ
zðx; tn Þ zðx; tn1 Þ
s
f ðx; t n Þ f ðx; tn1 Þ
s
þ zðx; t n1 Þ
bðx; t n Þ bðx; t n1 Þ
s
:
Using an inductive argument and the maximum principle for problem (4) can be easily deduced that kzðx; tn ÞkX 6 C; n ¼ 0; . . . ; M. From the boundness of zðx; tn Þ, (7), (8), jf ðx; tn Þ f ðx; tn1 Þj 6 C s and jbðx; t n Þ bðx; t n1 Þj 6 C s, the result trivially follows. h Now we are ready to obtain the main result of this section, showing the behavior of the solution of the semidiscrete problems (4). Theorem 3. Let zðx; tn Þ be the solution of problem (4) at the time level n. Then, it holds
ðkÞ z ðx; t n Þ 6 C 1 þ ek=2 Be ðxÞ ;
0 6 k 6 4; 1 6 n 6 M:
ð12Þ
Proof. The proof proceeds by induction. Let n ¼ 1 be. The semidiscrete problem is given by
(
I
zðx; t 1 Þ ¼ f ðx; t 1 Þ; zð0; t1 Þ ¼ zð1; t 1 Þ ¼ 0: s þ Lx;e
ð13Þ
Then, using that jf ðx; t1 Þj 6 C and the maximum principle for the operator sI þ Lx;e , it follows that jzðx; t 1 Þj 6 C s. Moreover, problem (13) can be written as
(
ez00 ðx; t 1 Þ þ bðx; t1 Þzðx; t 1 Þ ¼ f ðx; t 1 Þ zðx;ts 1 Þ h1 ðx; sÞ; zð0; t1 Þ ¼ zð1; t 1 Þ ¼ 0;
ð14Þ
with jh1 ðx; sÞj 6 C. Then, following [11] it is possible to prove jz0 ð0; t 1 Þj 6 C e1=2 and jz0 ð1; t1 Þj 6 C e1=2 : Differentiating once w.r.t. x the differential equation in (13), we obtain
I
s
þ Lx;e zx ðx; t1 Þ ¼ fx ðx; t 1 Þ bx ðx; t1 Þzðx; t 1 Þ h2 ðx; sÞ;
with jh2 ðx; sÞj 6 C. Taking the barrier function
w1 ðxÞ ¼ C 1 þ e1=2 Be ðxÞ ; I
ð15Þ
I
we have s þ Lx;e w1 ðxÞ P C; and from the maximum principle for s þ Lx;e , it follows
jz0 ðx; t 1 Þj 6 C 1 þ e1=2 Be ðxÞ :
Moreover, from (14) we directly obtain jz00 ð0; t 1 Þj 6 C e1 and jz00 ð1; t 1 Þj 6 C e1 , and differentiating twice w.r.t. x the differential equation in (13), we obtain
I
s
þ Lx;e z00 ðx; t1 Þ ¼ fxx ðx; t 1 Þ bxx ðx; t 1 Þzðx; t 1 Þ 2bx ðx; t 1 Þz0 ðx; t 1 Þ h3 ðx; sÞ;
with jh3 ðx; sÞj 6 ð1 þ e1=2 Be ðxÞÞ. Taking now the barrier function
w2 ðxÞ ¼ C 1 þ e1 Be ðxÞ ;
ð16Þ
in the same way as before we can prove jz ðx; t1 Þj 6 Cð1 þ e Be ðxÞÞ: For the third and the fourth order derivatives, using (9) we know that 00
1
1482
C. Clavero, J.L. Gracia / Applied Mathematics and Computation 216 (2010) 1478–1488
00 z ðx; t 1 Þ z00 ðx; 0Þ 6 C e1 ; s and then a similar argument than the previous one permit to obtain the bounds
ðkÞ z ðx; t 1 Þ 6 C 1 þ ek=2 Be ðxÞ ;
k ¼ 3; 4;
which is the required result. We assume now the induction hypothesis
ðkÞ z ðx; t m Þ 6 C 1 þ ek=2 Be ðxÞ ;
k ¼ 0; 1; 2; 3; 4; 1 6 m 6 n;
and we will prove that similar bounds hold also for the time level n þ 1. We know that
(
I
s þ Lx;e
zðx; t nþ1 Þ ¼ f ðx; t nþ1 Þ þ zðx;ts n Þ g 1 ðx; sÞ;
ð17Þ
zð0; tnþ1 Þ ¼ zð1; t nþ1 Þ ¼ 0;
with jg 1 ðx; sÞj 6 C=s. From the maximum principle we deduce jzðx; t nþ1 Þj 6 C. Moreover, problem (17) can be written as
(
nþ1 Þ ez00 ðx; t nþ1 Þ þ bðx; tnþ1 Þzðx; t nþ1 Þ ¼ f ðx; t nþ1 Þ þ zðx;tn Þzðx;t g 2 ðx; sÞ; s zð0; tnþ1 Þ ¼ zð1; t nþ1 Þ ¼ 0;
ð18Þ
with jg 2 ðx; sÞj 6 C from (7). Using again the technique given in [11] we obtain jz0 ð0; t nþ1 Þj 6 C e1=2 ; jz0 ð1; t nþ1 Þj 6 C e1=2 ; and jz00 ð0; t nþ1 Þj 6 C e1 ; jz00 ð1; tnþ1 Þj 6 C e1 directly follows from (18). Differentiating w.r.t. x the differential equation in (17), we obtain
I
s
z0 ðx; tn Þ þ Lx;e z0 ðx; tnþ1 Þ ¼ fx ðx; tnþ1 Þ bx ðx; tnþ1 Þzðx; t nþ1 Þ þ g 3 ðx; sÞ:
s
From the hypothesis of induction, it holds jg 3 ðx; sÞj 6 s ð1 þ e Be ðxÞÞ: Taking again the barrier functions w1 ðxÞ and w2 ðxÞ given in (15) and (16) respectively, and using the maximum principle for sI þ Lx;e we deduce C
jz0 ðx; tnþ1 Þj 6 C 1 þ e1=2 Be ðxÞ ;
1=2
jz00 ðx; tnþ1 Þj 6 C 1 þ e1 Be ðxÞ :
For the third and fourth order derivatives we proceed as follows. Differentiating twice w.r.t. x the differential equation in (18), we obtain
8 z00 ðx;t n Þz00 ðx;t nþ1 Þ ð4Þ 00 > < ez ðx; t nþ1 Þ þ bðx; t nþ1 Þz ðx; tnþ1 Þ ¼ fxx ðx; tnþ1 Þ þ s 0 bxx ðx; tnþ1 Þzðx; tnþ1 Þ 2bx ðx; t nþ1 Þz ðx; tnþ1 Þ g 4 ðx; sÞ; > : 00 z ð0; t nþ1 Þ; z00 ð1; t nþ1 Þ given:
ð19Þ
From the bound for jz0 ðx; t nþ1 Þj and using (7) we have that jg 4 ðx; sÞj 6 C e1 , and following [11] we deduce jz000 ð0; t nþ1 Þj 6 C e3=2 ; jz000 ð1; tnþ1 Þj 6 C e3=2 . Moreover, from (19) it directly follows jzð4Þ ð0; tnþ1 Þj 6 C e2 ; jzð4Þ ð1; t nþ1 Þj 6 C e2 : Now, differentiating three times w.r.t. x the differential equation in (17), we obtain
I
s
z000 ðx; t n Þ þ Lx;e z000 ðx; t nþ1 Þ ¼ fxxx ðx; t nþ1 Þ bxxx ðx; t nþ1 Þzðx; t nþ1 Þ 3bxx ðx; tnþ1 Þz0 ðx; t nþ1 Þ 3bx ðx; t nþ1 Þz00 ðx; t nþ1 Þ þ
s
g 5 ðx; sÞ: From the hypothesis of induction and the bounds for jzðkÞ ðx; t nþ1 Þj; k ¼ 0; 1; 2; it follows jg 5ðx; sÞj 6 Cs ð1 þ e3=2 Be ðxÞÞ: Then, using the barrier function w3 ðxÞ ¼ Cð1 þ e3=2 Be ðxÞÞ, and the maximum principle for sI þ Lx;e we obtain
jz000 ðx; tnþ1 Þj 6 C 1 þ e3=2 Be ðxÞ : In a completely similar way, it can be proved that jzð4Þ ðx; tnþ1 Þj 6 C 1 þ e2 Be ðxÞ , which is the required result.
h
For the analysis of the uniform convergence of the totally discrete method, we write problem (4) in vector form LM x;e z ¼ f, where
z ¼ ðzðx; t1 Þ; . . . ; zðx; tM ÞÞT ;
LM x;e ¼
h
i1 h iM T LM ; . . . ; LM ; x;e x;e
with
h in zðx; t Þ zðx; t Þ n n1 LM þ Lx;e zðx; tn Þ ¼ f ðx; t n Þ; x;e z
s
1 6 n 6 M:
The initial condition is zðx; 0Þ ¼ 0; x 2 ð0; 1Þ; and the boundary conditions are zð0; tn Þ ¼ zð1; t n Þ ¼ 0; n ¼ 0; . . . ; M:
ð20Þ
C. Clavero, J.L. Gracia / Applied Mathematics and Computation 216 (2010) 1478–1488
1483
3. The totally discrete method: uniform convergence M XN To discretize in space we consider a nonuniform mesh XN ¼ f0 ¼ x0 < < xN ¼ 1g, and we denote by Q M;N ¼ x M;N M;N the corresponding grid for the x; t-variables, and Q ¼Q \ Q. On XN we use central differences to approximate problem (20), which is given by
h
in U n U n1 j j LM;N U þ LNx;e U nj ¼ f ðxj ; t n Þ; e
s
j
1 6 n 6 M; 0 < j < N;
ð21Þ
where
LNx;e U nj ed2x U nj þ bðxj ; t n ÞU nj ; and d2x is the classical central difference approximation on a nonuniform mesh for the second order derivative. The initial condition is given by U 0j ¼ 0; j ¼ 0; . . . ; N, and the boundary conditions are U n0 ¼ U nN ¼ 0; 1 6 n 6 M. The following discrete maximum principle holds for this scheme. Lemma 5. Let LM;N be the discrete operator given in (21) and U ¼ U0 ; . . . ; UM g is a vectorial function with Un ¼ fU n0 ; . . . ; U nN g. If e n M;N 0 n n , then U P 0: U P 0; U 0 P 0; U N P 0; 1 6 n 6 M, and ½LM;N e Uj P 0 in Q From Lemma 5 it follows that the numerical scheme is uniformly stable in the maximum norm, and also the following comparison principle holds. n n M;N Corollary 2. If V, W are two vectorial functions such that V0 P W0 ; V n0 P W n0 ; V nN P W nN ; 1 6 n 6 M, and ½LM;N e Vj P ½Le Wj M;N in Q , then V P W.
To prove the uniform convergence of the numerical solution of the scheme (21), three different special nonuniform meshes, adapted to singularly perturbed reaction–diffusion problems, are used. The Shishkin mesh [11,12]. It is a piecewise uniform mesh, depending on two transition points which are defined by means of the transition parameter
r ¼ min
pffiffiffi 1 ; r0 e ln N ; 4
ð22Þ
where r0 is a constant to be fixed later. A uniform mesh is placed in ½0; r; ½r; 1 r, and ½1 r; 1, such that x0 ¼ 0; xN=4 ¼ r; x3N=4 ¼ 1 r, and xN ¼ 1 and therefore the mesh points are given by
8 j ¼ 0; 1; . . . ; N=4; > < 4jr=N; xj ¼ r þ 2ðj N=4Þð1 2rÞ=N; j ¼ N=4 þ 1; . . . ; 3N=4; > : j ¼ 3N=4 þ 1; . . . ; N: 1 r þ 4ðj 3N=4Þr=N;
ð23Þ
The Vulanovic mesh [13]. This mesh is a generalized Shishkin mesh constructed by using a suitable generating function @, which also depends on two transition points. To simplify, we use the same parameter r given in (22). The grid points are defined by xj ¼ @ðj=NÞ; j ¼ 0; 1; . . . ; N=2, with @ 2 C 2 ½0; 1=2 and
(
4r=N;
@ðzÞ ¼
z 2 ½0; 1=4; 3
pðz 1=4Þ þ 4rðz 1=4Þ þ r; z 2 ½1=4; 1=2:
ð24Þ
The coefficient p is such that @ð1=2Þ ¼ 1=2 and the mesh is symmetric with respect to the point x ¼ 1=2. Note that in ½0; r and ½1 r; 1 the mesh points are the same than in the Shishkin mesh. Otherwise, in ½r; 1 r it is nonuniform but the step sizes satisfy
jhjþ1 hj j 6 CN2 ;
j ¼ N=4; . . . ; 3N=4:
ð25Þ
The Bahkvalov mesh [7]. This mesh is based on user–chosen mesh parameters r0 > 0 and q 2 ð0; 1=2Þ. To simplify, we take pffiffiffi pffiffiffi q ¼ 1=4. The mesh points are xj ¼ j=N; j ¼ 0; 1; . . . ; N if r0 e < q, while when r0 e P q they are given by
xj ¼
uðj=NÞ; for j 6 N=2; 1 uððN jÞ=NÞ; for j > N=2;
ð26Þ
where
(
uðzÞ ¼
pffiffiffi
vðzÞ r0 e lnð1 4zÞ; z 2 ½0; m; pðzÞ vðmÞ þ v0 ðmÞðz mÞ; z 2 ½m; 1=2;
and the transition point 0
m is the solution of the nonlinear equation
ð1 2mÞv ðmÞ ¼ 1 2vðmÞ:
ð27Þ
1484
C. Clavero, J.L. Gracia / Applied Mathematics and Computation 216 (2010) 1478–1488
This defines the mesh in [0,1/2] and the mesh is again symmetric with respect to the central point x ¼ 1=2. Note that in ½0; m and ½1 m; 1 the mesh is graded and in ½m; 1 m the mesh is uniform. Theorem 4. Let Un be the numerical solution of (21) and zðx; t n Þ be the solution of (4) both at time level t n . We assume that the pffiffiffi constant r0 used to define the Shishkin, Vulanovic and Bahkvalov meshes satisfies r0 P 2= b. Then, if the method is constructed on the Shishkin or the Vulanovic meshes, the error satisfies
n
U j zðxj ; t n Þ
XN
6 CðN1 ln NÞ2 ;
ð28Þ
and therefore the space discretization is an almost second order uniformly convergent scheme. Otherwise, if the method is constructed on the Bahkvalov mesh the error satisfies
n
U j zðxj ; t n Þ
XN
6 CN2 ;
ð29Þ
and therefore the space discretization is a second order uniformly convergent scheme. Proof. We begin with the Shishkin mesh. Using similar ideas as in [4], we can obtain that the local error satisfies
( 1 h 2 in M;N L ðU zÞ ¼ Lx;e LN zðxj ; tn Þ 6 CN e þ CN ; if xj ¼ r; 1 r; x;e e 2 1 j otherwise: CðN ln NÞ ; Note that from the uniform stability of the scheme we cannot deduce the almost second order of uniform convergence. Nevertheless, using a standard barrier function in the context of singularly perturbed problems given by ½Wnj ¼ CðN 2 hðxj Þ ln N þ ðN 1 ln NÞ2 Þ; where
8x > : 1x r
ifx 2 ½0; r; if x 2 ½r; 1 r; ; ifx 2 ½1 r; 1;
h in M;N n L ðU zÞ ; and then from Lemma 2 the result follows. it holds LM;N ½W P e j e j In second place, we consider the Vulanovic mesh. Using that in ð0; rÞ [ ð1 r; 1Þ the mesh is the same than for the Shishkin mesh, and in ½r; 1 r the step sizes satisfy (25), it is not difficult to deduce that
( 2 h in M;N ifxj 2 ½r; 1 r; L ðU zÞ 6 CN ; e j CðN1 ln NÞ2 ; otherwise; and therefore the result follows from the uniform stability of the scheme. Finally, for the Bahkvalov mesh the same analysis than this one given in [7] (see Appendix, p. 2094) permits to deduce that the local error satisfies
h in M;N L ðU zÞ 6 CN 2 ; e j and again the result follows from the uniform stability of the scheme. h Theorem 5. Let U be the numerical solution of (21), u the solution of (1), and constructed on the Shishkin or the Vulanovic meshes, the error satisfies
n
U j uðxj ; tn Þ
Q M;N
6C
pffiffiffi
r0 such that r0 P 2= b. Then, if the method is
s þ ðN1 ln NÞ2 :
ð30Þ
Otherwise, if the method is constructed on the Bahkvalov mesh, the error satisfies
n
U j uðxj ; tn Þ
Q M;N
6 Cðs þ N2 Þ:
ð31Þ
Proof. The results follows from Theorems 2 and 4. h Remark 1. Note that in Theorem 5 there is not any relation between the discretization parameter N and s in contrast with previous results in [1,5], where the theoretical restriction N q 6 C s with 0 < q < 1 appears. This fact accords with our numerical experiments, where that restriction between the two discretization parameters has been never needed. Moreover, the proof of the uniform convergence here is considerably simpler than those ones in previous papers because now we have not needed any recursive argument invoking the uniform stability of any transition operator.
C. Clavero, J.L. Gracia / Applied Mathematics and Computation 216 (2010) 1478–1488
1485
4. Numerical results The first example that we consider is given by
8 1þx2 t > < ut euxx þ 2 u ¼ e 1 þ sinðpxÞ; ðx; tÞ 2 ð0; 1Þ ð0; 1; uðx; 0Þ ¼ 0; x 2 X; > : uð0; tÞ ¼ uð1; tÞ ¼ 0; t 2 ð0; 1;
ð32Þ
which solution is unknown. In Fig. 1 we show the numerical solution of this problem for a particular value of e ¼ 104 at two different times. From it we clearly see the boundary layers at the end points of the spatial domain. To estimate the numerical errors we use a variant of the double mesh principle and therefore at times t n ¼ ns for each point xj ; j ¼ 0; 1; . . . ; N the error is approximated by using
e;2N;s Dej;n;N;s ¼ U ej;n;N;s U j;n 4 ;
e;2N;s
where U ej;n;N;s is the numerical solution obtained using a constant time step s and ðN þ 1Þ points in the spatial mesh and U j;n 4 is computed using 4s as time step and ð2N þ 1Þ points in the spatial mesh but with the same transition parameters than in the original mesh. For each fixed value of e, the maximum global errors are estimated by
De;N;s ¼ max Dej;n;N;s ; j;n
and therefore, in standard way, the numerical orders of convergence are given by
q¼
logðDe;N;s =De;2N;s=4 Þ : log 2
From these values we obtain the e-uniform errors and the e-uniform orders of convergence by
DN;s ¼ max De;N;s ; e
quni ¼
logðDN;s =D2N;s=4 Þ : log 2
ð33Þ
Tables 1 and 2 display the results on the Shishkin, Vulanovic and Bahkvalov meshes. In all cases we have taken the constant r0 ¼ 2. Note that the results on Shishkin and Vulanovic meshes are identical and for this reason they appear together in Table 1. The maximum errors appear in the boundary layer regions where both meshes are exactly the same, showing the almost second order uniform convergence. Table 2 shows that the method on the Bahkvalov mesh has second order convergence according with the theoretical results. The second test problem that we consider has a solution with a really small variation in time (see [1]); the reason to take this example is that we wish that the error associated to the spatial discretization stage dominates in the global error. In the previous example to see the influence of the spatial discretization in the global error we have needed to divide the time step
1.4
t=0.5 t=1
1.2 1 0.8 0.6 0.4 0.2 0
0
0.1
0.2
0.3
0.4
0.5
0.6
0.7
Fig. 1. Numerical solution at different times of (32) for
0.8
0.9
1
e ¼ 104 with N ¼ 32.
1486
C. Clavero, J.L. Gracia / Applied Mathematics and Computation 216 (2010) 1478–1488
Table 1 Maximum errors and uniform orders of convergence on the Shishkin and Vulanovic meshes for problem (32).
e
N=32 s ¼ 0:1
N=64
N=128
N=256
N=512
N=1024
s ¼ 0:1=4
s ¼ 0:1=42
s ¼ 0:1=43
s ¼ 0:1=44
s ¼ 0:1=45
26
0.179E1 1.983
0.453E2 1.996
0.113E2 1.999
0.284E3 2.000
0.710E4 2.000
0.177E4
28
0.221E1 1.978
0.562E2 1.994
0.141E2 1.999
0.353E3 2.000
0.883E4 2.000
0.221E4
210
0.246E1 1.980
0.624E2 1.994
0.157E2 1.998
0.392E3 2.000
0.980E4 2.000
0.245E4
212
0.263E1 1.986
0.664E2 1.996
0.166E2 1.757
0.493E3 1.592
0.163E3 1.642
0.524E4
214
0.272E1 1.985
0.687E2 1.997
0.172E2 1.810
0.490E3 1.591
0.163E3 1.641
0.522E4
216
0.277E1 1.986
0.698E2 1.997
0.175E2 1.838
0.489E3 1.591
0.162E3 1.640
0.521E4
218
0.279E1 1.986
0.704E2 1.997
0.176E2 1.852
0.489E3 1.591
0.162E3 1.640
0.521E4
220
0.280E1 1.987
0.707E2 1.997
0.177E2 1.859
0.489E3 1.591
0.162E3 1.640
0.521E4
222
0.281E1 1.987 ...
0.709E2 1.997 ...
0.178E2 1.863 ...
0.489E3 1.591 ...
0.162E3 1.640 ...
0.520E4
230
0.282E1 1.987
0.711E2 1.997
0.178E2 1.866
0.488E3 1.591
0.162E3 1.640
0.520E4
DN;s quni
0.282E1
0.711E2
0.178E2
0.493E3
0.163E3
0.524E4
1.987
1.997
1.854
1.592
1.642
...
...
Table 2 Maximum errors and uniform orders of convergence on the Bahkvalov mesh for problem (32).
e
N=32 s ¼ 0:1
N=64
N=128
N=256
N=512
N=1024
s ¼ 0:1=4
s ¼ 0:1=42
s ¼ 0:1=43
s ¼ 0:1=44
s ¼ 0:1=45
26
0.179E1 1.983
0.453E2 1.996
0.113E2 1.999
0.284E3 2.000
0.710E4 2.000
0.177E4
28
0.223E01 1.981
0.565E02 1.995
0.142E02 1.999
0.355E03 2.000
0.887E04 2.000
0.222E04
210
0.248E01 1.975
0.632E02 1.995
0.158E02 1.998
0.397E03 2.000
0.992E04 2.000
0.248E04
212
0.265E01 1.986
0.668E02 1.994
0.168E02 1.999
0.420E03 2.000
0.105E03 2.000
0.262E04
214
0.273E01 1.988
0.689E02 1.997
0.173E02 1.999
0.432E03 2.000
0.108E03 2.000
0.270E04
216
0.277E01 1.989
0.698E02 1.996
0.175E02 1.999
0.438E03 2.000
0.110E03 2.000
0.274E04
218
0.279E01 1.987
0.704E02 1.996
0.177E02 1.999
0.442E03 2.000
0.110E03 2.000
0.276E04
220
0.280E01 1.986
0.708E02 1.997
0.177E02 1.999
0.443E03 2.000
0.111E03 2.000
0.277E04
222
0.281E01 1.985 ...
0.709E02 1.997 ...
0.178E02 1.999 ...
0.444E03 2.000 ...
0.111E03 2.000 ...
0.278E04
30
2
0.282E01 1.987
0.711E02 1.997
0.178E02 1.999
0.445E03 2.000
0.111E03 2.000
0.278E04
DN;s quni
0.282E01
0.711E02
0.178E02
0.445E03
0.111E03
0.278E04
1.987
1.997
1.999
2.000
2.000
...
...
size by four and therefore the computational cost of the method is higher. In this new example the time step size is only divided by two. The problem is defined as
1487
C. Clavero, J.L. Gracia / Applied Mathematics and Computation 216 (2010) 1478–1488
8 > < ut euxx þ ð10 þ 0:1 sinðpxÞÞu ¼ f ðx; tÞ; ðx; tÞ 2 ð0; 1Þ ð0; 1; uðx; 0Þ ¼ 0; x 2 X; > : t 2 ð0; 1; uð0; tÞ ¼ g 1 ðtÞ; uð1; tÞ ¼ g 2 ðtÞ;
ð34Þ
where f ; g 1 and g 2 are such that the exact solution is given by
pffiffiffiffiffiffiffiffiffiffiffi uðx; tÞ ¼ expð 10=exÞð1 expðt=100ÞÞ: Note that the boundary conditions are not homogeneous, but a simple change of variable transforms the example (34) in a problem of type (1). Now the errors at t n ¼ ns; xj ; j ¼ 0; 1; . . . ; N and the orders of convergence are calculated by
Eej;n;N;s ¼ U ej;n;N;s uðxj ; t n Þ; the maximum global errors and the numerical orders of convergence are given by
Ee;N;s ¼ max Eej;n;N;s ; j;n
p¼
logðEe;N;s =Ee;2N;s=2 Þ : log 2
From these values the e-uniform errors and the e-uniform orders of convergence are obtained by
EN;s ¼ max Ee;N;s ; e
puni ¼
logðEN;s =E2N;s=2 Þ : log 2
ð35Þ
Again we take r0 ¼ 2, but we have obtained similar results for smaller values of this parameter. Also, the maximum errors on the Shishkin and Vulanovic meshes are identical for problem (34). Tables 3 and 4 display the results on the Shishkin/Vula-
Table 3 Maximum errors and uniform orders of convergence on the Shishkin and Vulanovic meshes for problem (34).
e
N=32 s ¼ 0:1
N=64
N=128
N=256
N=512
N=1024
s ¼ 0:1=2
s ¼ 0:1=22
s ¼ 0:1=23
s ¼ 0:1=24
s ¼ 0:1=25
26
0.836E04 1.937
0.218E04 1.958
0.561E05 2.002
0.140E05 2.011
0.348E06 2.023
0.856E07
28
0.267E03 1.673
0.836E04 1.936
0.219E04 1.957
0.563E05 1.999
0.141E05 2.005
0.351E06
210
0.382E03 0.517
0.267E03 1.673
0.837E04 1.936
0.219E04 1.956
0.563E05 1.997
0.141E05
212
0.382E03 0.454
0.279E03 1.183
0.123E03 1.547
0.421E04 1.626
0.136E04 1.687
0.423E05
...
...
...
...
...
...
...
230
0.382E03 0.454
0.279E03 1.183
0.123E03 1.547
0.421E04 1.626
0.136E04 1.687
0.423E05
EN;s puni
0.382E03
0.279E03
0.123E03
0.421E04
0.136E04
0.423E05
0.454
1.183
1.547
1.626
1.687
Table 4 Maximum errors and uniform orders of convergence on the Bahkvalov mesh for problem (34).
e
N=32 s ¼ 0:1
N=64
N=128
N=256
N=512
N=1024
s ¼ 0:1=2
s ¼ 0:1=22
s ¼ 0:1=23
s ¼ 0:1=24
s ¼ 0:1=25
26
0.836E04 1.937
0.218E04 1.958
0.561E05 2.002
0.140E05 2.011
0.348E06 2.023
0.856E07
28
0.429E04 1.901
0.115E04 1.990
0.289E05 2.003
0.721E06 2.022
0.178E06 2.046
0.430E07
210
0.428E04 1.900
0.115E04 1.990
0.289E05 2.004
0.721E06 2.022
0.177E06 2.046
0.430E07
212
0.428E04 1.900 ...
0.115E04 1.990 ...
0.289E05 2.004 ...
0.721E06 2.022 ...
0.177E06 2.046 ...
0.430E07
2
0.428E04 1.900
0.115E04 1.990
0.289E05 2.004
0.721E06 2.022
0.177E06 2.046
0.430E07
EN;s puni
0.836E04
0.218E04
0.561E05
0.140E05
0.348E06
0.856E07
1.937
1.958
2.002
2.011
2.023
... 30
...
1488
C. Clavero, J.L. Gracia / Applied Mathematics and Computation 216 (2010) 1478–1488
novic and Bahkvalov meshes respectively. From these tables we observe again that the method has almost second order on the Shishkin and Vulanovic meshes and second order on the Bahkvalov mesh, in agreement with the theoretical results. Acknowledgments This research was partially supported by the project MEC/FEDER MTM2007-63204 and the Diputación General de Aragón. References [1] B. Bujanda, C. Clavero, J.L. Gracia, J.C. Jorge, A high order uniformly convergent alternating direction scheme for time dependent reaction–diffusion singularly perturbed problems, Num. Math. 107 (2007) 1–25. [2] C. Clavero, J.L. Gracia, High order methods for elliptic and time dependent reaction–diffusion singularly perturbed problems, Appl. Math. Comp. 168 (2005) 1109–1127. [3] C. Clavero, J.L. Gracia, J.C. Jorge, A uniformly convergent alternating direction HODIE finite difference scheme for 2D time dependent convection– diffusion problems, IMA J. Numer. Anal. 26 (2006) 155–172. [4] C. Clavero, J.C. Jorge, F. Lisbona, G.I. Shishkin, An alternating direction scheme on a nonuniform mesh for reaction–diffusion parabolic problem, IMA J. Numer. Anal. 20 (2000) 263–280. [5] J.L. Gracia, F. Lisbona, A uniformly convergent scheme for a system of reaction–diffusion equations, J. Comput. Appl. Math. 206 (2007) 1–16. [6] P.W. Hemker, G.I. Shishkin, L.P. Shiskina, e-uniform schemes with high-order time accuracy for parabolic singular perturbation problems, IMA J. Numer. Anal. 20 (2000) 99–121. [7] B. Kellogg, T. Linß, M. Stynes, A finite difference method on layer-adapted meshes for an elliptic reaction–diffusion system in two dimensions, Math. Comp. 77 (2008) 2085–2096. [8] T. Linß, Layer-adapted meshes and FEM for time-dependent singularly perturbed reaction-diffusion problems, Int. J. Computing Science and Math. 1 (2007) 259–270. [9] T. Linß, A robust method for a time-dependent system of coupled singularly perturbed reaction–diffusion problems, preprint MATH-NM-03-2008, Institut für Numerische Mathematik, TU Dresden. [10] T. Linß, N. Madden, Parameter uniform approximations for time-dependent reaction-diffusion problems, Numer. Methods Partial Differ. Equ. 23 (2007) 1290–1300. [11] J.J.H. Miller, E. O’Riordan, G.I. Shishkin, Fitted Numerical Methods for Singular Perturbation Problems, World-Scientific, Singapore, 1996. [12] H.-G. Roos, M. Stynes, L. Tobiska, Robust numerical methods for singularly perturbed differential equations, Springer Series in Computational Mathematics 24, Springer-Verlag, Berlin, 2008. [13] R. Vulanovic, A high order scheme for quasilinear boundary value problems with two small parameters, Computing 67 (1986) 287–303.
