- Email: [email protected]
High order parameter uniform numerical method for singular perturbation problems
Applied Mathematics and Computation 188 (2007) 720–733 www.elsevier.com/locate/amc
High order parameter uniform numerical method for singular perturbation problems Kailash C. Patidar Department of Mathematics and Applied Mathematics, University of the Western Cape, Private Bag X17 Bellville 7535, South Africa
Abstract In this paper, we systematically describe how to derive a method of higher order for the numerical solution of singularly perturbed ordinary differential equations. First we apply this idea to derive a fourth-order method for a self-adjoint singularly perturbed two point boundary value problem. This method is uniformly convergent on a piecewise uniform mesh of Shishkin type. After we have developed and analyzed a fourth-order method, we explain with appropriate details, how can one obtain the methods of order higher than four which looks straightforward but has not been seen in the literature so far. Besides these, the fourth-order e-uniformity in the theoretical estimate has been justified by some numerical experiments. 2006 Elsevier Inc. All rights reserved. Keywords: Singular perturbation problems; Ordinary differential equations; Boundary value problems; Finite difference methods; Shishkin Mesh
1. Introduction One question always strikes in every numerical analyst’s mind: ‘‘Why we need higher-order methods?’’ Not only because they converge faster than the lower order methods but also that they provide highly accurate results in fewer number of steps. Whenever the high order accuracy is/was expected, many researchers try/ tried in various different ways to design the higher-order methods. Out of the main two approaches to obtain higher-order methods; the first one is to use the methods based on the idea of acceleration of convergence whereas the other, considered in this paper, is to use the higher-order methods directly. We explore the idea of deriving some higher order schemes (specially those which uses variable step sizes because the one with the uniform step size can be derived from these trivially) by considering the following self-adjoint singularly perturbed two point boundary value problem 0
Ly eðaðxÞy 0 Þ þ bðxÞy ¼ f ðxÞ;
x 2 ½0; 1;
yð0Þ ¼ g0 ;
E-mail address: [email protected] 0096-3003/$ - see front matter 2006 Elsevier Inc. All rights reserved. doi:10.1016/j.amc.2006.10.040
yð1Þ ¼ g1 ;
ð1:1Þ
K.C. Patidar / Applied Mathematics and Computation 188 (2007) 720–733
721
where g0, g1 are given constants and 0 < e 6 1. Further, f(x), a(x) and b(x) are sufficiently smooth functions satisfying the conditions aðxÞ P a > 0;
bðxÞ P b > 0:
Under these conditions, the operator L admits the following maximum principle [5]: Lemma 1.1. Let L be the operator as in (1.1) such that B0y(0) y(0) = g0, B1y(1) y(1) = g1. Suppose /(x) be any smooth function satisfying B0/(0) P 0, B1/(1) P 0 and let L/(x) P 0, " 0 < x < 1 then /(x) P 0, " 0 6 x 6 1. Proof. The proof is by contradiction. Let x* be such that /(x*) = minx2[0,1]/(x) and assume that /(x*) < 0. Clearly, x* 62 {0, 1} and therefore / 0 (x*) = 0 and /00 (x*) P 0. Further, 0
L/ðx Þ ¼ eðaðx Þ/0 ðx ÞÞ þ bðx Þ/ðx Þ < 0; which is a contradiction. It follows that /(x*) P 0 and thus /(x) P 0 "x 2 [0, 1]. This maximum principle implies that the solution is unique and it also exists (as for linear problems, the existence of the solution is implied by its uniqueness). Further, we show below using this principle that the solution of (1.1) is bounded. h Lemma 1.2. Let y(x) be the solution of the problem (1.1), then we have kyk 6 b1 kf k þ maxðjg0 j; jg1 jÞ: Proof. We construct two barrier functions P± defined by P ðxÞ ¼ b1 kf k þ maxðjg0 j; jg1 jÞ yðxÞ: Then we have P ð0Þ ¼ b1 kf k þ maxðjg0 j; jg1 jÞ yð0Þ ¼ b1 kf k þ maxðjg0 j; jg1 jÞ g0 P 0; P ð1Þ ¼ b1 kf k þ maxðjg0 j; jg1 jÞ yð1Þ
¼ b1 kf k þ maxðjg0 j; jg1 jÞ g1 P 0; and we have 0 0
LP ðxÞ ¼ eðaðxÞðP ðxÞÞ Þ þ bðxÞP ðxÞ ¼ bðxÞðb1 kf k þ maxðjg0 j; jg1 jÞÞ LyðxÞ ¼ bðxÞ½b1 kf k þ maxðjg0 j; jg1 jÞ f ðxÞ P ðkf k f ðxÞÞ þ bðxÞ maxðjg0 j; jg1 jÞ P 0;
since kf k P f ðxÞ:
Therefore, by the maximum principle (Lemma 1.1), we obtain P±(x) P 0, for all x 2 [0, 1], which gives the required estimate. h Parameter dependent differential equations occurs naturally in various fields of science and engineering. Singular perturbation problems (Spps) belong to the class of such problems in which a very small positive parameter is multiplied to the highest order derivative term in the differential equation. Most of the conventional methods fail when this small parameter approaches to zero. Numerous books have appeared and are still appearing continuously in this direction (see, e.g., [4,6,10,11,17–20,22,23,27]) which hold the fascination for researchers in this direction. The list of research works on spps is quite long and the interested readers may refer to author’s recently appeared two large survey articles [8,9].
722
K.C. Patidar / Applied Mathematics and Computation 188 (2007) 720–733
Since the solution of the singular perturbation problems possesses boundary and/or interior layer(s), problems have been experienced with the classical methods when applied to the problems of the type (1.1). Except some special cases, most of the classical methods require severe restrictions on the step size to preserve the stability when e is very small. To avoid such restrictions, there are two ways to design the schemes which have small truncation errors inside these layers. The first one is to choose a difference formula (on a uniform mesh) which reflects the behavior of the solution in these layers. This strategy falls under the class of fitted operator methods. The second one is to choose a fine mesh there which forms the class of fitted mesh methods. Both the strategies have their own merits and demerits. On one hand, fitted operator methods work nicely on uniform mesh but are difficult to extend for higher dimensional problems. On the other hand, the fitted mesh methods require the knowledge of the ‘location and exact width’ of the layer(s) but are easy to extend for higher dimensional and nonlinear problems. Sometimes (but not always) the combination of the two methods is useful. A very good discussion using one or both of the above strategies can be found in Miller et al. [17]. Other notable articles for the problems of the type (1.1) are those of Boglaev [3], Miller [16], Niijima [21], O’Riordan and Stynes [24], Roos [26], Schatz and Wahlbin [30], Stojanovic´ [32,33], etc. The fitted operator methods were first suggested by Allen and Southwell [1] for solving the problem of viscous fluid flow past a cylinder. One type of exponentially fitted methods considered in Liniger and Willoughby [13] which is a special class of h-method given in Lambert [12] has been considered by Doolan et al. [5]. The details about the construction of the fitting factor in the above methods has been discussed in [5]. However, it is to be noted that the highest possible order of convergence in most of the above methods was only one (and two in some special cases). Recently, in [25] we have derived some e-uniform second-order and fourth-order (not e-uniform) methods just to fill the gap mentioned above. However, since the method presented in [25] is still not fourth-order e-uniform, we went through the first strategy and obtained, as we expected, the e-uniform results using the variable mesh version of the method (SNFDM) considered in [25]. Towards the use of the fitted mesh methods, tremendous works have been done so far. Whether one need to use a special scheme on a given mesh or a special mesh is to be used for a given scheme, can be answered by taking into account an excellent quote (by Stynes [35]): Miller has moved from [15] the question ‘‘what scheme should one use on a given mesh?’’ to [6] ‘‘what mesh should one use with a given scheme?’’ Hitherto, various mesh selection strategies (either graded or piecewise uniform meshes) have been provided in the literature, however, the most successful and widely used ones are those of Bakhvalov-type and Shishkin-type meshes. Bakhvalov [2] was the first to use an a priori mesh to solve a singular perturbation problem. His mesh was generalized by Vulanovic´ [37]. A mesh of this kind is generated by a suitable smooth mesh generating function which in the layer corresponds to the inverse of the function describing the layer behavior of the continuous solution. Bakhvalov meshes have been successfully applied to many different types of singular perturbation problem, see, e.g., [38–42]. A new impetus to the a priori mesh approach was given by Shishkin [31], who discovered that simple piecewise equidistant meshes can also produce e-uniform results. A mesh of this kind uses a fine equidistant mesh inside each layer, and a coarse, again equidistant, mesh on the outside. The Shishkin meshes are discussed in detail in [27,34,36] as well as in [28], where Roos surveyed results on layer-adapted meshes. It is shown in [42] that though Shishkin meshes are simpler than Bakhvalov meshes, the latter produce much better numerical results. The reason is that Bakhvalov meshes are better adopted to the layer structure. Comparisons of these two meshes made in [14,29,42,43] have shown that the standard three-point difference schemes produce much better results on Bakhvalov meshes than on Shishkin meshes. However, it is also well known that a higher-order scheme is much easier to construct on an equidistant mesh than on a fully nonequidistant one. The paper [44] shows that the real advantage of Shishkin meshes over Bakhvalov meshes is born out when applying more complicated higher order schemes. Since we are also working on higher order schemes, all of these comparisons suggest us therefore to use a Shishkin type mesh for the problem under consideration. To the best of our knowledge, the results presented in this paper are the first, fourth-order results for a self-adjoint spp (which is in conservation form). The systematic way of constructing the methods of any order higher than four (for such problems) is also a plus to this paper. We would like to mention at this stage that higher order uniform methods do exist for some other type of problems, e.g., semi-linear spps, see, [7,36,42].
K.C. Patidar / Applied Mathematics and Computation 188 (2007) 720–733
723
The rest of the paper is organized as follows. In Section 2, we reduce Eq. (1.1) to the normal form. Shishkin mesh (suitable for problem in this paper) has been described in Section 3. Section 4 deals with the design of a fourth-order fitted mesh finite difference scheme. Convergence analysis of this scheme has been given in Section 5. In Section 6, we explain, how can one get the schemes of order higher than four. Two numerical examples demonstrating the e-uniform convergence of the proposed method along with some comparisons with the results obtained via standard method have been given in Section 7. Finally, Section 8 deals with the conclusion and future directions. Notations and terminology Most of the notations and symbols we employ are fairly standard. However, to avoid notational complexities, we omit subscripted e from ye, Le and so on. 2. Reduction to the normal form Re-writing Eq. (1.1) as y 00 þ P ðxÞy 0 þ QðxÞy ¼ RðxÞ;
ð2:1Þ
where P ðxÞ ¼ Let
a0 ðxÞ ; aðxÞ
QðxÞ ¼
bðxÞ eaðxÞ
and
RðxÞ ¼
f ðxÞ : eaðxÞ
Z 1 x U ðxÞ :¼ exp P ðfÞ df : 2 0
Using yðxÞ ¼ U ðxÞV ðxÞ;
ð2:2Þ
we transform (2.1) into the normal form, that is, e ðxÞV ¼ RðxÞ; e V 00 þ W
ð2:3Þ
where e ðxÞ ¼ QðxÞ 1 P 0 ðxÞ 1 ðP ðxÞÞ2 ; W 2 4 Z x 1 e RðxÞ ¼ RðxÞ exp P ðfÞ df ; 2 0 with V ð0Þ ¼
yð0Þ ¼: a0 ; U ð0Þ
V ð1Þ ¼
yð1Þ ¼: a1 ; a0 ; a1 2 R: U ð1Þ
Multiplying (2.3) throughout by e, we get e LV eV 00 þ W ðxÞV ¼ ZðxÞ; V ð0Þ ¼ a0 ;
V ð1Þ ¼ a1 ;
ð2:4Þ
where e ðxÞ and W ðxÞ ¼ e W
e ZðxÞ ¼ e RðxÞ:
It is to be noted that the operator e L also satisfies the continuous maximum principle: Lemma 2.1. Let w(x) be any sufficiently smooth function such that w(0) P 0 and w(1) P 0. Then e LwðxÞ P 0, " x 2 (0, 1) implies that w(x) P 0, " x 2 [0, 1].
724
K.C. Patidar / Applied Mathematics and Computation 188 (2007) 720–733
Proof. Analogous to the proof of Lemma 1.1. To solve problem (2.4), we design the fitted mesh fourth-order method in Section 4. h 3. Description of the mesh Let n be a positive integer and denote by d the width of the boundary layer. Let the interval [0, 1] be divided into three sub-intervals: ½0; 1 :¼ ½0; d [ ½d; 1 d [ ½1 d; 1: The piecewise uniform mesh (of Shishkin type) in these three sub-intervals is designed as follows. Assuming that n = 2m with m P 5, the intervals (0, d) and (1 d, 1) are each divided into n/4 equal mesh elements, while the interval (d, 1 d) is divided into n/2 equal mesh elements. This insures that there are 1 + n/4 equidistant grid points into the intervals [0, d] and [1 d, 1] whereas one has n/2 1 equidistant grid points in (d, 1 d). We define the parameter d by n pffiffiffiffiffiffiffi o d ¼ min 1=4; 4 e=b ln ðn=16Þ : ð3:1Þ where 0 < b 6 W(x), " x 2 [0, 1]. Let j0 ¼ n=4;
xj0 ¼ d;
xnj0 ¼ 1 d
and ½0; 1 :¼ 0 ¼ x0 < x1 < < xj0 < < xnj0 < < xn ¼ 1; with xj = hj hj1, where the mesh spacing is given by 4dn1 ; j ¼ 1; . . . ; j0 ; n j0 þ 1; . . . ; n; ð3:2Þ hj ¼ 2ð1 2dÞn1 ; j ¼ j0 þ 1; . . . ; n j0 : pffiffiffiffiffiffiffi If d = 1/4, i.e., 1=4 < 4 e=b lnðn=16Þ then n1 is very small relative to e. This is very unlikely in practice and in such a case the method can be analyzed using the standard techniques. We therefore assume that pffiffiffiffiffiffiffi d¼4 e=b ln ðn=16Þ: ð3:3Þ 4. Fitted mesh finite difference scheme We use the notations Vj = V(xj), Wj = W(xj) and Zj = Z(xj) and denote the approximations of the Vj at the grid points xj by the unknowns mj. Taking the Taylor series expansion for m around xj, and neglecting the terms containing fifth and higher order terms, we get the following approximations for mj+1 and mj1: mjþ1 mj þ hjþ1 m0j þ
h2jþ1 00 h3jþ1 000 h4jþ1 ðivÞ m þ m þ m 2! j 3! j 4! j
ð4:1Þ
and mj1 mj hj m0j þ
h2j 00 h3j 000 h4j ðivÞ m m þ mj : 2! j 3! j 4!
ð4:2Þ
Adding Eqs. (4.1) and (4.2), we obtain mj1 2mj þ mjþ1 ðhjþ1 hj Þm0j þ
h2jþ1 þ h2j 00 h3jþ1 h3j 000 h4jþ1 þ h4j ðivÞ mj þ mj þ mj : 2! 3! 4!
ð4:3Þ
Neglecting third and fourth-order term in Eqs. (4.1) and (4.2), and solving the resulting system for m0j and m00j , we get
K.C. Patidar / Applied Mathematics and Computation 188 (2007) 720–733
725
m0j
h i 1 h2j mjþ1 ðh2j h2jþ1 Þmj h2jþ1 mj1 hj hjþ1 ðhj þ hjþ1 Þ
ð4:4Þ
m00j
2 ½hj mjþ1 ðhj þ hjþ1 Þmj þ hjþ1 mj1 hj hjþ1 ðhj þ hjþ1 Þ
ð4:5Þ
and
The relations 0
00 m000 j ¼ ðmj Þ
and
ðivÞ
mj
¼ ðm00j Þ
00
then gives the approximations for h i 1 m000 h2j m00jþ1 ðh2j h2jþ1 Þm00j h2jþ1 m00j1 j hj hjþ1 ðhj þ hjþ1 Þ
ð4:6Þ
ð4:7Þ
and ðivÞ
mj
2 ½hj m00jþ1 ðhj þ hjþ1 Þm00j þ hjþ1 m00j1 : hj hjþ1 ðhj þ hjþ1 Þ
ð4:8Þ
Using Eqs. (2.4), (4.4), (4.7) and (4.8) in (4.3) and simplifying, we obtain the tridiagonal system Am ¼ F
ð4:9Þ
where A is the matrix of the system and m and F are corresponding vectors. The various entries of this matrix and the components of the r.h.s. vector are given by 9 ðdiagðAÞÞj ¼ rcj ; j ¼ 1; 2; . . . ; n 1; > > > > > > ðsubdiagðAÞÞj ¼ rj ; j ¼ 2; 3; . . . ; n 1; > > > > > þ > ðsupdiagðAÞÞj ¼ rj ; j ¼ 1; 2; . . . ; n 2; > > = þ c F 1 ¼ q1 Z 0 þ q1 Z 1 þ q1 Z 2 r 1 m 0 ; > > > c þ > F j ¼ q Z þ q Z þ q Z ; j ¼ 2; 3; . . . ; n 2; > j1 j jþ1 j j j > > > > þ þ c > > F n1 ¼ q Z þ q Z þ q Z r m ; > n1 n n1 n n1 n2 n1 n1 > > ; m 0 ¼ a0 ; m n ¼ a1 ; where
; 12hjþ1
9 > > > 6 > > > > 3 3 > > 4 4 > h h h þh h j j j jþ1 jþ1 þ > 12hj ; > qj ¼ hjþ1 ðhj þhjþ1 Þ > 6 > > > > > > 3 3 4 4 2 2 2ðh h Þðh h Þðh þh Þ h þh = j jþ1 j j j jþ1 jþ1 jþ1 c ; qj ¼ 2 þ 12hj hjþ1 > n o > > hjþ1 ðhj hjþ1 Þ > ¼ e 1 W ; r þ q j1 > j j > hj ðhj þhjþ1 Þ > > > n o > > hj ðhj hjþ1 Þ þ þ > rj ¼ e 1 þ hjþ1 ðhj þhjþ1 Þ þ qj W jþ1 ; > > > > > n o > 2 > ðhj hjþ1 Þ ; c c rj ¼ e 2 þ hj hjþ1 þ qj W j : h
jþ1 q j ¼ hj ðhj þhjþ1 Þ
h3j h3jþ1
þ
ð4:10Þ
h4j þh4jþ1
ð4:11Þ
It is to be noted that when hj = hj+1 = h (i.e., uniform mesh throughout the region) then (4.11) reduces to ) 2 2 r rþ rcj ¼ 2e þ 5h2 W j =6; j ¼ e þ h W j1 =12; j ¼ e þ h W jþ1 =12; ð4:12Þ 2 h2 q qcj ¼ 5h6 : j ¼ 12 ;
726
K.C. Patidar / Applied Mathematics and Computation 188 (2007) 720–733
Using Eqs. (4.9)–(4.11) or (4.12) we get the approximate solution of V(x) at the grid points xj. Since U(x) is known, therefore the solution to the original problem (1.1) at these grid points will be obtained using (2.2). We call the method consisting of Eqs. (4.9)–(4.11) as the fitted mesh finite difference method (FMFDM) whereas the method consisting of Eqs. (4.9),(4.10) and (4.12) is the standard numerov’s finite difference method (SNFDM). 5. Convergence analysis of FMFDM Throughout the paper M will denote positive constants which may take different values in different equations (inequalities) but are always independent of h and e. Denoting by e L h the discrete operator in the FMFDM, we see that it satisfies e h /i P 0, Lemma 5.1 (Discrete maximum principle). For any mesh function /i satisfying /0 P 0, /n P 0 and L " 0 < i < n, we have /i P 0, " 0 6 i 6 n. and Lemma 5.2 (Uniformly stability estimate). If fi is any mesh function such that f0 = fn = 0, then jfi j 6
1 max j e L h fj j for b 16j6n1
0 6 i 6 n:
To facilitate the analysis properly, we decompose the solution V(x) into its regular (smooth) and singular components and provide the bounds on them and on their derivatives using Lemma [17]: Lemma 5.3. The solution Ve of the problem (2.4) can be decomposed into the form V e :¼ V e;r þ V e;s where for all k 2 [0, 6] and x 2 [0, 1], the regular (smooth) component Ve,r satisfies ðk2Þ=2 Eðx; bÞ; jV ðkÞ e;r ðxÞj 6 M½1 þ e
and the singular component Ve,s satisfies k=2 Eðx; bÞ; jV ðkÞ e;s ðxÞj 6 Me
where 0 < b 6 W ðxÞ
and
n pffiffiffiffiffiffiffi pffiffiffiffiffiffiffio Eðx; bÞ ¼ exp x b=e þ exp ð1 xÞ b=e :
Very often during the analysis of spps, one encounters with the exponential terms involving e and then divided by the terms involving various powers of e, which are the main cause of worry. For their careful consideration while proving the e-uniform convergence, we prove Lemma: Lemma 5.4. For a fixed mesh and for all integers k, we have pffiffi pffiffi exp Mxj = e exp Mð1 xj Þ= e lim max ¼ 0 and lim max ¼ 0; e!0 16j6n1 e!0 16j6n1 ek=2 ek=2 where xj = jh, h = 1/n, "j = 1(1)n 1. Proof. Consider the partition ½0; 1 :¼ f0 ¼ x0 < x1 < x2 < < xn1 < xn ¼ 1g: Clearly, for the interior grid points, we have pffiffi pffiffi pffiffi exp Mxj = e exp ðMx1 = eÞ exp ðMh= eÞ 6 ¼ max 16j6n1 ek=2 ek=2 ek=2
K.C. Patidar / Applied Mathematics and Computation 188 (2007) 720–733
727
and
pffiffi pffiffi pffiffi exp Mð1 xj Þ= e exp ðMð1 xn1 Þ= eÞ exp ðMh= eÞ 6 ¼M : max 16j6n1 ek=2 ek=2 ek=2 (as x1 = h, 1 xn1 = 1 (n 1)h = h). An application of L’Hospital’s rule then gives pffiffi expðMh= eÞ pk k! lim lim ¼ lim ¼ 0; pffi k=2 p!1 ðMhÞk expðMhpÞ e!0 e pð¼1= eÞ!1 expðMhpÞ
which completes the proof.
h
Now, the local truncation error of the scheme Eqs. (4.9)–(4.11) is c c þ 00 c 00 þ 00 ðe L h ðV mÞÞj ¼ ðr j qj W j1 ÞV j1 þ ðr j qj W j ÞV j þ ðr j qj W jþ1 ÞV jþ1 þ eðqj V j1 þ qj V j þ qj V jþ1 Þ:
Expanding Vj1, Vj+1 and their derivatives in terms of Vj and its derivatives, we obtain ð1Þ
ð2Þ
ð3Þ
ð4Þ
ð5Þ
ðe L h ðV mÞÞj ¼ T 0 V j þ T 1 V j þ T 2 V j þ T 3 V j þ T 4 V j þ T 5 V j þ T 6;1 V ð6Þ ðn1j Þ þ T 6;2 V ð6Þ ðn2j Þ þ T 6;3 V ð6Þ ðn3j Þ þ T 6;4 V ð6Þ ðn4j Þ;
ð5:1Þ
where n1j ; n3j 2 ðxj1 ; xj Þ, n2j ; n4j 2 ðxj ; xjþ1 Þ and c þ c þ T 0 ¼ ðr j þ rj þ rj Þ ðqj W j1 þ qj W j þ qj W jþ1 Þ; þ T 1 ¼ hjþ1 ðrþ j qj W jþ1 Þ hj ðrj qj W j1 Þ;
T2 ¼
h2j h2jþ1 þ c þ ðrj q ðrj qþ j W j1 Þ þ j W jþ1 Þ þ eðqj þ qj þ qj Þ; 2 2
T3 ¼ T4 ¼
h3j h3jþ1 þ þ ðrj q ðrj qþ W Þ þ j1 j j W jþ1 Þ þ eðhjþ1 qj hj qj Þ; 6 6
h4j h4jþ1 þ e 2 þ 2 ðrj q ðrj qþ W Þ þ j1 j j W jþ1 Þ þ ðhjþ1 qj þ hj qj Þ; 2 24 24
T5 ¼
ð5:2Þ
h5j h5jþ1 þ e 3 þ 3 ðrj q ðr qþ j W j1 Þ þ j W jþ1 Þ þ ðhjþ1 qj hj qj Þ; 6 120 120 j
T 6;1 ¼
h6j ðr q j W j1 Þ; 720 j
T 6;3 ¼
eh4j q j ; 24
T 6;4 ¼
T 6;2 ¼
h6jþ1 þ ðr qþ j W jþ1 Þ; 720 j
eh4jþ1 qþ j : 24
We have two cases now: either the grid point xj is inside the fine mesh, i.e., j 2 f1; . . . ; j0 1g [ fn j0 þ 1; . . . ; n 1g
ð5:3Þ
or it is inside the coarse mesh, i.e., j 2 fj0 ; . . . ; n j0 g:
ð5:4Þ
We denote by me,r and me,s, respectively, the regular and the singular components of the numerical solution. Then using (4.11) along with Eqs. (3.2), (3.3) into (5.2) and Lemma 5.3, we see that 8 6 6 n in case of ð5:3Þ; > < n ln 16 ; h p ffiffiffi ffi p ffiffiffi ffi e ð5:5Þ jð L ðV e;r me;r ÞÞj j 6 M xj b=e ð1xnj Þ b=e 0 > : n6 e 0 þe e ; in case of ð5:4Þ:
728
K.C. Patidar / Applied Mathematics and Computation 188 (2007) 720–733
(In the above, we have used the fact that the relation 2(1 2d)n1 6 2n1 is true for all step sizes and the quantities 1 + M/e are always less than the quantities M/e for different M.) Similarly, we see that 8 6 6 n in case of ð5:3Þ; > < n ln 16 ; h p ffiffiffi ffi p ffiffiffi ffi e jð L ðV e;s me;s ÞÞj j 6 M ð5:6Þ xj b=e ð1xnj Þ b=e 0 > : n6 e 0 þee2 ; in case of ð5:4Þ: Combining Eqs. (5.5), (5.1)–(5.6) and using Lemma 5.4 followed by the use of uniform stability estimate (Lemma 5.2), we obtain our main result stated in the following Theorem. Theorem 5.5. Let W(x), Z(x) be sufficiently smooth so that V(x) 2 C6[0, 1] and W(x) P b > 0. Let mj, j = 0(1)n, be the approximate solution of (2.4), obtained using Eqs. (4.9)–(4.11) with m0 = V(0), mn = V(1). Then, there is a constant M independent of e and h such that sup max jV j mj j 6 Mn6 ln6 ðn=16Þ 6 Mn4 ;
0
where we have used the fact that ln3(n/16) 6 M n, "n. Remark 5.6. On one hand, the case when n1 6 e then the second parts in the estimates Eqs. (5.5), (5.1)–(5.6) is bounded above, respectively, by Mn5 and Mn4. However, the case when n1 6 e in the region [d, 1 d] is practically not interesting. One the other hand, since xj0 ; ð1 xnj0 Þ 2 ð0; 1Þ, one can see by using Lemma 5.4 that the second parts inpffiffiffi(5.5) (andpffiffiffisimilarly in (5.6)) tends to zero (as e ! 0). However, in general the ffi ffi xj
b=e
ð1xnj Þ
b=e
0 above quantities like e 0 þe e are dominating and deteriorate the results in the pre-asymptotic range of e’s. The readers can see such influence in most of the articles devoted to singularly perturbed problems though not mentioned explicitly in any of the articles so far. But when we say that the method is uniformly convergent (or robust w.r.t. the parameter) then for any value of this parameter (whether it is 101 of 102 or 1010) the errors must remain same for a fix number of partitions of the interval [0, 1]. Theoretically, many authors provide the error bounds as M(hp) or M(hpln (hp)) for the pth order methods with M independent of h and e which appeals everyone very nice but in view of the numerics which they provide, it is really not always right to say so.
6. Schemes of the order higher than four In previous section, we have derived a fourth-order method. The main idea behind deriving the method of order p originated from Eqs. (4.1) and (4.2) in which one can take the expansions up to pth-order derivative terms. Because of variable step sizes, higher the value of p is taken, more the complications arise in deriving such schemes. However, we give here a fifth-order scheme in which we use the relations ðvÞ
ðivÞ
mj ¼ ðmj Þ0 ;
ðivÞ
mj
¼ ðm00j Þ00
and
00 0 m000 j ¼ ðmj Þ :
ð6:1Þ
This along with the procedure used earlier leads to a pentadiagonal system 1 c þ1 þ2 2 1 c þ1 þ2 r2 j mj2 þ rj mj1 þ r j mj þ r j mjþ1 þ rj mjþ2 ¼ qj Z j2 þ qj Z j1 þ qj Z j þ qj Z jþ1 þ qj Z jþ2 ;
ð6:2Þ
where q2 j ¼
hjþ1 h2j ðh5jþ1 h5j Þ ; 60hj hj1
qþ2 j ¼
h2jþ1 hj ðh5jþ1 h5j Þ ; 60hj hjþ1
q1 j
h2jþ1 ðh3jþ1 h3j Þ hjþ1 ðh4jþ1 þ h4j Þ hjþ1 ðh2j h2j1 Þðh5jþ1 h5j Þ h2jþ1 ðhjþ1 þ hj Þðh5jþ1 h5j Þ ¼ þ þ þ ; 6hj 12hj 60hj hj1 60h2j
qþ1 j
h2j ðh3jþ1 h3j Þ hj ðh4jþ1 þ h4j Þ hj ðh2jþ2 h2jþ1 Þðh5jþ1 h5j Þ h2j ðhjþ1 þ hj Þðh5jþ1 h5j Þ ¼ þ þ ; 6hj 12hj 60hj hjþ1 60h2j
K.C. Patidar / Applied Mathematics and Computation 188 (2007) 720–733
729
ðh2jþ1 þ h2j Þ ðh2jþ1 h2j Þðh3jþ1 h3j Þ ðhjþ1 þ hj Þðh4jþ1 þ h4j Þ þ 2 6hj 12hj ( ) ðh5jþ1 h5j Þ h2jþ1 hj ðh2jþ1 h2j Þðhjþ1 þ hj Þ hjþ1 h2j1 þ ; 60 hj hjþ1 hj hj1 h2j ! h2jþ1 ðhjþ1 hj Þ 2 þ2 þ2 1 þ q1 ¼ qj W j2 ; rj ¼ qj W jþ2 ; rj ¼ e 1 þ j W j1 ; hj ! h2j ðhjþ1 hj Þ ðh2jþ1 h2j Þðhjþ1 hj Þ c ¼e 1 W ; r ¼ eð2 þ Þ þ qcj W j þ qþ1 jþ1 j j hj hj
qcj ¼
r2 j rþ1 j
with hj = hjhj+1(hj + hj+1). One can use the scheme (6.2) at the grid points x2 to xn2, whereas at x1 and xn1 one can use the FMFDM derived in the previous section. In a similar manner, one can derive a sixth-order scheme and so on. 7. Test examples and numerical results To illustrate the theory, we give some numerical results corresponding to problem (1.1) for various choices of coefficient functions. Example 7.1. Consider problem (1.1) with aðxÞ ¼ 1 þ x2 ;
bðxÞ ¼ 1 þ xð1 xÞ; pffiffi
pffiffi pffiffi f ðxÞ ¼ 1 þ xð1 xÞ exp x= e xð2x2 3x þ 1Þ 2 e 2x2 x 1 þ e þ 1 pffiffi pffiffi pffiffi
þ exp ð1 xÞ= e x2 ð2x 1Þ þ 2 e 2x2 þ x e þ 1 ; yð0Þ ¼ 0; yð1Þ ¼ 0:
Its exact solution is given by
pffiffi pffiffi yðxÞ ¼ 1 þ ðx 1Þ exp x= e x exp ð1 xÞ= e : Example 7.2. Consider problem (1.1) with aðxÞ ¼ 1; bðxÞ ¼ 1; f ðxÞ ¼ ðcos2 px þ 2ep2 cos 2pxÞ;
yð0Þ ¼ 0;
yð1Þ ¼ 0:
Its exact solution is given by
pffiffi pffiffi pffiffi yðxÞ ¼ exp ð1 xÞ= e þ exp x= e = 1 þ exp 1= e cos2 px: Maximum errors at all the mesh points are evaluated using the formula: En;e :¼ max jyðxj Þ ~mj j; 06j6n
for different values of n and e, where ~mj ~mnj is the approximate solution of (1.1) obtained via Eqs. (2.4) and (2.2). The numerical rates of convergence are computed using the formula [5]: rk rk;e :¼ log2 ðEnk ;e =E2nk ;e Þ;
k ¼ 1; 2; . . . :
Further, we compute En ¼ max En;e 0
whereas the numerical rate of uniform convergence is computed as Rn :¼ log2 ðEn =E2n Þ:
730
K.C. Patidar / Applied Mathematics and Computation 188 (2007) 720–733
For these examples we present the numerical results obtained by the SNFDM (see Tables 1 and 4) followed by those obtained by FMFDM (see Tables 2 and 5). Tables 3 and 6 present the numerical rates of convergence.
Table 1 Results for Example 7.1 (Max. errors using SNFDM) e
n = 64
n = 128
n = 256
n = 512
n = 1024
0.1E01 0.1E02 0.1E03 0.1E04 0.1E05 0.1E08 0.1E12
0.11E05 0.75E04 0.59E02 0.98E01 0.18E+00 0.20E+00 0.20E+00
0.69E07 0.48E05 0.39E03 0.24E01 0.15E+00 0.20E+00 0.20E+00
0.43E08 0.30E06 0.26E04 0.24E02 0.68E01 0.20E+00 0.20E+00
0.27E09 0.19E07 0.16E05 0.15E03 0.12E01 0.20E+00 0.20E+00
0.17E10 0.12E08 0.10E06 0.98E05 0.92E03 0.20E+00 0.20E+00
Table 2 Results for Example 7.1 (Max. errors using FMFDM) e
n = 64
n = 128
n = 256
n = 512
n = 1024
0.1E01 0.1E02 0.1E03 0.1E04 0.1E05 0.1E06 0.1E07 0.1E08 0.1E09 0.1E10 0.1E11 0.1E12
0.11E05 0.46E04 0.63E03 0.45E02 0.10E01 0.13E01 0.14E01 0.15E01 0.15E01 0.15E01 0.15E01 0.15E01
0.69E07 0.48E05 0.84E05 0.91E04 0.40E03 0.67E03 0.80E03 0.84E03 0.86E03 0.86E03 0.86E03 0.86E03
0.43E08 0.30E06 0.10E05 0.11E05 0.11E04 0.31E04 0.43E04 0.48E04 0.49E04 0.50E04 0.50E04 0.50E04
0.27E09 0.19E07 0.15E06 0.15E06 0.19E06 0.11E05 0.22E05 0.27E05 0.29E05 0.29E05 0.30E05 0.30E05
0.17E10 0.12E08 0.20E07 0.19E07 0.19E07 0.28E07 0.97E07 0.15E06 0.17E06 0.17E06 0.18E06 0.18E06
En
0.15E01
0.86E03
0.50E04
0.30E05
0.18E06
Table 3 Results for Example 7.1 (rates of convergence using FMFDM) nk = 64, 128, 256, 512, 1024 e
r1
r2
r3
r4
r5
0.1E01 0.1E02 0.1E03 0.1E04 0.1E05 0.1E06 0.1E07 0.1E08 0.1E09 0.1E10 0.1E11 0.1E12
0.40E+01 0.34E+01 0.62E+01 0.56E+01 0.46E+01 0.43E+01 0.42E+01 0.41E+01 0.41E+01 0.41E+01 0.41E+01 0.41E+01
0.40E+01 0.40E+01 0.33E+01 0.63E+01 0.51E+01 0.44E+01 0.42E+01 0.41E+01 0.41E+01 0.41E+01 0.41E+01 0.41E+01
0.40E+01 0.40E+01 0.27E+01 0.31E+01 0.59E+01 0.47E+01 0.43E+01 0.41E+01 0.41E+01 0.41E+01 0.41E+01 0.41E+01
0.40E+01 0.40E+01 0.29E+01 0.29E+01 0.35E+01 0.53E+01 0.45E+01 0.42E+01 0.41E+01 0.41E+01 0.41E+01 0.41E+01
0.40E+01 0.40E+01 0.31E+01 0.31E+01 0.31E+01 0.38E+01 0.49E+01 0.43E+01 0.41E+01 0.41E+01 0.40E+01 0.39E+01
Rn
0.41E+01
0.41E+01
0.41E+01
0.41E+01
0.39E+01
K.C. Patidar / Applied Mathematics and Computation 188 (2007) 720–733
731
Table 4 Results for Example 7.2 (max. errors using SNFDM) e
n = 64
n = 128
n = 256
n = 512
n = 1024
0.1E01 0.1E02 0.1E03 0.1E04 0.1E05 0.1E08 0.1E12
0.48E06 0.45E04 0.38E02 0.54E01 0.95E01 0.10E+00 0.10E+00
0.30E07 0.28E05 0.27E03 0.15E01 0.78E01 0.10E+00 0.10E+00
0.19E08 0.18E06 0.18E04 0.17E02 0.39E01 0.10E+00 0.10E+00
0.12E09 0.11E07 0.11E05 0.11E03 0.78E02 0.10E+00 0.10E+00
0.74E11 0.70E09 0.70E07 0.69E05 0.67E03 0.99E01 0.10E+00
Table 5 Results for Example 7.2 (max. errors using FMFDM) e
n = 64
n = 128
n = 256
n = 512
n = 1024
0.1E01 0.1E02 0.1E03 0.1E04 0.1E05 0.1E06 0.1E07 0.1E08 0.1E09 0.1E10 0.1E11 0.1E12
0.48E06 0.74E04 0.41E03 0.21E02 0.40E02 0.50E02 0.54E02 0.55E02 0.55E02 0.55E02 0.55E02 0.55E02
0.30E07 0.28E05 0.79E05 0.50E04 0.17E03 0.26E03 0.30E03 0.31E03 0.31E03 0.32E03 0.32E03 0.32E03
0.19E08 0.18E06 0.69E06 0.78E06 0.54E05 0.12E04 0.16E04 0.17E04 0.18E04 0.18E04 0.18E04 0.18E04
0.12E09 0.11E07 0.18E06 0.11E06 0.11E06 0.49E06 0.83E06 0.98E06 0.10E05 0.11E05 0.11E05 0.12E05
0.74E11 0.70E09 0.46E07 0.14E07 0.14E07 0.14E07 0.39E07 0.54E07 0.60E07 0.63E07 0.63E07 0.64E07
En
0.55E02
0.32E03
0.18E04
0.12E05
0.64E07
Table 6 Results for Example 7.2 (rates of convergence using FMFDM) nk = 64, 128, 256, 512, 1024 e
r1
r2
r3
r4
r5
0.1E01 0.1E02 0.1E03 0.1E04 0.1E05 0.1E06 0.1E07 0.1E08 0.1E09 0.1E10 0.1E11 0.1E12
0.40E+01 0.48E+01 0.58E+01 0.54E+01 0.46E+01 0.43E+01 0.42E+01 0.41E+01 0.41E+01 0.41E+01 0.41E+01 0.41E+01
0.40E+01 0.40E+01 0.36E+01 0.60E+01 0.50E+01 0.44E+01 0.42E+01 0.42E+01 0.41E+01 0.41E+01 0.41E+01 0.41E+01
0.40E+01 0.40E+01 0.21E+01 0.30E+01 0.56E+01 0.46E+01 0.43E+01 0.41E+01 0.41E+01 0.41E+01 0.41E+01 0.39E+01
0.40E+01 0.40E+01 0.19E+01 0.29E+01 0.31E+01 0.51E+01 0.44E+01 0.42E+01 0.41E+01 0.41E+01 0.41E+01 0.43E+01
0.40E+01 0.40E+01 0.22E+01 0.31E+01 0.31E+01 0.32E+01 0.47E+01 0.43E+01 0.41E+01 0.41E+01 0.37E+01 0.24E+01
Rn
0.41E+01
0.41E+01
0.39E+01
0.43E+01
0.24E+01
8. Conclusions and future directions Using a fitted mesh finite difference method on a mesh of Shishkin type, we have described a numerical method for solving self-adjoint singular perturbation problems in conservation form. The method has been analyzed for e-uniform convergence. Two examples have been solved to demonstrate the predicted theory. A way to derive the methods of order higher than four is also explained in details.
732
K.C. Patidar / Applied Mathematics and Computation 188 (2007) 720–733 0.16
Diff. of Exact and SNFDM Solutions Diff. of Exact and FMFDM Solutions 0.14
0.12
0.1
0.08
0.06
0.04
0.02
0 0
0.1
0.2
0.3
0.4
0.5
0.6
0.7
0.8
0.9
1
Fig. 1. Errors of SNFDM and FMFDM solutions for Example 7.1 for n = 128 and e = 106.
As is expected the method is advantageous over the conventional methods (e.g., standard finite difference methods, etc). Comparative results obtained by our method are presented along with those obtained by the Standard Numerov’s method. Importance of the fitted mesh method can be seen from the respective tabular results. The solution of such singular perturbation problems possesses boundary (interior) layers in which the solution changes rapidly as e ! 0. On a uniform mesh, the SNFDM fails to capture this singularly perturbed nature of the solution in these layer regions. However, the fitted method takes care of it and the resulting behavior can be seen from Fig. 1. In this figure, we plot the difference between the exact solution and numerical solution(s) obtained via FMFDM and SNFDM. The method presented in this paper takes the advantage of the special nature of the differential Eq. (1.1) that allows one to reduce it to the normal form (2.4). Hence, any suitable method for Eq. (2.4) can be applied for (1.1) after some transformations. However, it would be interesting if such higher-order methods can be designed directly for (1.1). The same would be interesting for the nonlinear versions of it. Note that the higher-order methods for the nonlinear versions of (2.4) (some semi-linear spps) do exist in literature as pointed out earlier. The method proposed in this paper can also be extended for such problems by applying it to the sequence of quasi-linearized equations. However, in view of the current interest of the researchers, instead of this quasilinearization approach, we will prefer to have some direct methods for the general nonlinear problems. Some possible methods in this directions are currently under the investigation of the author. Acknowledgement The author is indebted to Prof. R. Vulanovic´ for his many stimulating suggestions. References [1] D. Allen, R.V. Southwell, Relaxation methods applied to determine the motion in 2-D of a viscous fluid past a fixed cylinder, Q. J. Mech. Appl. Math. VIII (2) (1955) 129–145. [2] N.S. Bakhvalov, Towards optimization of methods for solving boundary value problems in the presence of a boundary layer, Zh. Vychisl. Mat. Mat. Fiz. 9 (1969) 841–859 (in Russian). [3] I.P. Boglaev, A variational difference scheme for a boundary value problem with a small parameter in the highest derivative, USSR Comput. Math. Math. Phys. 21 (4) (1981) 71–81. [4] K.W. Chang, F.A. Howes, Nonlinear Singular Perturbation Phenomena: Theory and Application, Springer-Verlag, New York, 1984. [5] E.P. Doolan, J.J.H. Miller, W.H.A. Schilders, Uniform Numerical Methods for Problems with Initial and Boundary Layers, Boole Press, Dublin, 1980. [6] P.A. Farrell, A. Hegarty, J.J.H. Miller, E. O’Riordan, G.I. Shishkin, Robust Computational Techniques for Boundary Layers, Chapman & Hall/CRC, New York, 2000.
K.C. Patidar / Applied Mathematics and Computation 188 (2007) 720–733
733
[7] D. Herceg, Uniform fourth order difference scheme for a singular perturbation problem, Numer. Math. 56 (1990) 675–693. [8] M.K. Kadalbajoo, K.C. Patidar, A survey of numerical techniques for solving singulalry perturbed ordinary differential equations, Appl. Math. Comp. 130 (2-3) (2002) 457–510. [9] M.K. Kadalbajoo, K.C. Patidar, Singularly perturbed problems in partial differential equations: a survey, Appl. Math. Comp. 134 (23) (2003) 371–429. [10] J. Kevorkian, J.D. Cole, Perturbation Methods in Applied Mathematics, Springer-Verlag, New York, 1981. [11] J. Kevorkian, J.D. Cole, Multiple Scale and Singular Perturbation Methods, Springer-Verlag, New York, 1996. [12] J.D. Lambert, Computational Methods in Ordinary Differential Equations, John Wiley & Sons, New York, 1973. [13] W. Liniger, R.A. Willoughby, Efficient numerical integration methods for stiff system of differential equations, IBM Research Report RC-1970, 1967. [14] T. Linss, H.-G. Roos, R. Vulanovic´, Uniform pointwise convergence on Shishkin type meshes for quasiliner convection diffusion problems, SIAM J. Numer. Anal. 38 (2000) 897–912. [15] J.J.H. Miller, Construction of a FEM for a singularly perturbed problem in 2 dimensions, in Numerische Behandlung von Differentialgleichungen, Band 2 (Tagung, Math. Forschungsinst., Oberwolfach, 1975), Int. Ser. Numer. Math., Birkhau¨ser, Basel, 31 (1976), 165–169. [16] J.J.H. Miller, On the convergence uniformly in e, of difference schemes for a two-point boundary singular perturbation problem, Numerical Analysis of Singular Perturbation Problems, Academic Press, New York, 1979. [17] J.J.H. Miller, E. O’Riordan, G.I. Shishkin, Fitted Numerical Methods for Singular Perturbation Problems, World Scientific, Singapore, 1996. [18] K.W. Morton, Numerical Solution of Convection-Diffusion Problems, Chapman & Hall, London, 1996. [19] A.H. Nayfeh, Perturbation Methods, Wiley, New York, 1973. [20] A.H. Nayfeh, Introduction to Perturbation Techniques, John Wiley & Sons, New York, 1981. [21] K. Niijima, On a three-point difference scheme for a singular perturbation problem without a first derivative term I, Mem. Numer. Math. 7 (1980) 1–10. [22] R.E. O’Malley, Introductions to Singular Perturbations, Academic Press, New York, 1974. [23] R.E. O’Malley, Singular Perturbation Methods for Ordinary Differential Equations, Springer-Verlag, New York, 1991. [24] E. O’Riordan, M. Stynes, A uniformly accurate finite-element method for a singularly perturbed one-dimensional reaction-diffusion problem, Math. Comp. 47 (1986) 555–570. [25] K.C. Patidar, High order fitted operator numerical method for self-adjoint singular perturbation problems, Appl. Math. Comp. 171 (2005) 547–566. [26] H.-G. Roos, Global uniformly convergent schemes for a singularly perturbed boundary-value problem using patched base splinefunctions, J. Comput. Appl. Math. 29 (1) (1990) 69–77. [27] H.-G. Roos, M. Stynes, L. Tobiska, Numerical Methods for Singularly Perturbed Differential Equations, Springer-Verlag, Berlin, 1996. [28] H.-G. Roos, Layer adapted grids for singular perturbation problems, Z. Angew. Math. Mech. 78 (1998) 291–309. [29] H.-G. Roos, T. Linss, Sufficient conditions for uniform convergence on layer adapted grids, Computing 64 (1999) 27–45. [30] A.H. Schatz, L.B. Wahlbin, On the finite element method for singularly perturbed reaction diffusion problems in two and one dimension, Math. Comp. 40 (1983) 47–89. [31] G.I. Shishkin, A difference scheme for a singularly perturbed parabolic equation with a discontinuous boundary condition, Zh. Vychisl. Mat. Mat. Fiz. 28 (1988) 1679–1692 (in Russian). [32] M. Stojanovic´, A uniformly accurate spline collocation method for a singular perturbation problem, Calcolo 27 (1-2) (1990) 81–88. [33] M. Stojanovic´, Exponential splines difference scheme for singular perturbation problem with mixed boundary conditions, Comm. Appl. Numer. Meth. 7 (8) (1991) 625–632. [34] M. Stynes, H.-G. Roos, The midpoint upwind scheme, Appl. Numer. Math. 23 (1997) 361–374. [35] M. Stynes, A jejune heuristic mesh theorem, Comput. Meth. Appl. Math. 3 (2003) 488–492. [36] G. Sun, M. Stynes, An almost fourth order uniformly convergent difference scheme for a semilinear singlularly perturbed reaction diffusion problem, Numer. Math. 17 (1995) 487–500. [37] R. Vulanovic´, On a numerical solution of a type of singularly perturbed boundary value problem by using a special discretization mesh, Univ. Novom. Sadu Zb. Rad. Prirod. Mat. Fak. Ser. Mat 13 (1983) 187–201. [38] R. Vulanovic´, A uniform method for quasilinear singular perturbation problems without turning points, Computing 41 (1989) 97– 106. [39] R. Vulanovic´, A second order numerical method for nonlinear singular perturbation problems without turning points, J. Comput. Math. Math. Phys. 31 (1991) 522–532. [40] R. Vulanovic´, A uniform numerical method for a class of quasilinear turning point problems, in: J.J.H. Miller, R. Vichnevetsky (Eds.), Proceedings of the 13th IMACS World Congress on Computation and Applied Mathematics, Dublin, IMACS, 1991, pp. 493. [41] R. Vulanovic´, P.A. Farrell, Continuous and numerical analysis of a multiple boundary turning point problems, SIAM J. Numer. Anal. 30 (1993) 1400–1418. [42] R. Vulanovic´, Fourth order algorithms for a semilinear singular perturbation problem, Numer. Algor. 16 (1997) 117–128. [43] R. Vulanovic´, A priori meshes for singularly perturbed quasilinear two-point boundary value problems, IMA J. Numer. Anal. 21 (2001) 349–366. [44] R. Vulanovic´, A higher order scheme for quasilinear boundary value problems with two small parameters, Computing 67 (2001) 287– 303.
High order parameter uniform numerical method for singular perturbation problems Kailash C. Patidar Department of Mathematics and Applied Mathematics, University of the Western Cape, Private Bag X17 Bellville 7535, South Africa
Abstract In this paper, we systematically describe how to derive a method of higher order for the numerical solution of singularly perturbed ordinary differential equations. First we apply this idea to derive a fourth-order method for a self-adjoint singularly perturbed two point boundary value problem. This method is uniformly convergent on a piecewise uniform mesh of Shishkin type. After we have developed and analyzed a fourth-order method, we explain with appropriate details, how can one obtain the methods of order higher than four which looks straightforward but has not been seen in the literature so far. Besides these, the fourth-order e-uniformity in the theoretical estimate has been justified by some numerical experiments. 2006 Elsevier Inc. All rights reserved. Keywords: Singular perturbation problems; Ordinary differential equations; Boundary value problems; Finite difference methods; Shishkin Mesh
1. Introduction One question always strikes in every numerical analyst’s mind: ‘‘Why we need higher-order methods?’’ Not only because they converge faster than the lower order methods but also that they provide highly accurate results in fewer number of steps. Whenever the high order accuracy is/was expected, many researchers try/ tried in various different ways to design the higher-order methods. Out of the main two approaches to obtain higher-order methods; the first one is to use the methods based on the idea of acceleration of convergence whereas the other, considered in this paper, is to use the higher-order methods directly. We explore the idea of deriving some higher order schemes (specially those which uses variable step sizes because the one with the uniform step size can be derived from these trivially) by considering the following self-adjoint singularly perturbed two point boundary value problem 0
Ly eðaðxÞy 0 Þ þ bðxÞy ¼ f ðxÞ;
x 2 ½0; 1;
yð0Þ ¼ g0 ;
E-mail address: [email protected] 0096-3003/$ - see front matter 2006 Elsevier Inc. All rights reserved. doi:10.1016/j.amc.2006.10.040
yð1Þ ¼ g1 ;
ð1:1Þ
K.C. Patidar / Applied Mathematics and Computation 188 (2007) 720–733
721
where g0, g1 are given constants and 0 < e 6 1. Further, f(x), a(x) and b(x) are sufficiently smooth functions satisfying the conditions aðxÞ P a > 0;
bðxÞ P b > 0:
Under these conditions, the operator L admits the following maximum principle [5]: Lemma 1.1. Let L be the operator as in (1.1) such that B0y(0) y(0) = g0, B1y(1) y(1) = g1. Suppose /(x) be any smooth function satisfying B0/(0) P 0, B1/(1) P 0 and let L/(x) P 0, " 0 < x < 1 then /(x) P 0, " 0 6 x 6 1. Proof. The proof is by contradiction. Let x* be such that /(x*) = minx2[0,1]/(x) and assume that /(x*) < 0. Clearly, x* 62 {0, 1} and therefore / 0 (x*) = 0 and /00 (x*) P 0. Further, 0
L/ðx Þ ¼ eðaðx Þ/0 ðx ÞÞ þ bðx Þ/ðx Þ < 0; which is a contradiction. It follows that /(x*) P 0 and thus /(x) P 0 "x 2 [0, 1]. This maximum principle implies that the solution is unique and it also exists (as for linear problems, the existence of the solution is implied by its uniqueness). Further, we show below using this principle that the solution of (1.1) is bounded. h Lemma 1.2. Let y(x) be the solution of the problem (1.1), then we have kyk 6 b1 kf k þ maxðjg0 j; jg1 jÞ: Proof. We construct two barrier functions P± defined by P ðxÞ ¼ b1 kf k þ maxðjg0 j; jg1 jÞ yðxÞ: Then we have P ð0Þ ¼ b1 kf k þ maxðjg0 j; jg1 jÞ yð0Þ ¼ b1 kf k þ maxðjg0 j; jg1 jÞ g0 P 0; P ð1Þ ¼ b1 kf k þ maxðjg0 j; jg1 jÞ yð1Þ
¼ b1 kf k þ maxðjg0 j; jg1 jÞ g1 P 0; and we have 0 0
LP ðxÞ ¼ eðaðxÞðP ðxÞÞ Þ þ bðxÞP ðxÞ ¼ bðxÞðb1 kf k þ maxðjg0 j; jg1 jÞÞ LyðxÞ ¼ bðxÞ½b1 kf k þ maxðjg0 j; jg1 jÞ f ðxÞ P ðkf k f ðxÞÞ þ bðxÞ maxðjg0 j; jg1 jÞ P 0;
since kf k P f ðxÞ:
Therefore, by the maximum principle (Lemma 1.1), we obtain P±(x) P 0, for all x 2 [0, 1], which gives the required estimate. h Parameter dependent differential equations occurs naturally in various fields of science and engineering. Singular perturbation problems (Spps) belong to the class of such problems in which a very small positive parameter is multiplied to the highest order derivative term in the differential equation. Most of the conventional methods fail when this small parameter approaches to zero. Numerous books have appeared and are still appearing continuously in this direction (see, e.g., [4,6,10,11,17–20,22,23,27]) which hold the fascination for researchers in this direction. The list of research works on spps is quite long and the interested readers may refer to author’s recently appeared two large survey articles [8,9].
722
K.C. Patidar / Applied Mathematics and Computation 188 (2007) 720–733
Since the solution of the singular perturbation problems possesses boundary and/or interior layer(s), problems have been experienced with the classical methods when applied to the problems of the type (1.1). Except some special cases, most of the classical methods require severe restrictions on the step size to preserve the stability when e is very small. To avoid such restrictions, there are two ways to design the schemes which have small truncation errors inside these layers. The first one is to choose a difference formula (on a uniform mesh) which reflects the behavior of the solution in these layers. This strategy falls under the class of fitted operator methods. The second one is to choose a fine mesh there which forms the class of fitted mesh methods. Both the strategies have their own merits and demerits. On one hand, fitted operator methods work nicely on uniform mesh but are difficult to extend for higher dimensional problems. On the other hand, the fitted mesh methods require the knowledge of the ‘location and exact width’ of the layer(s) but are easy to extend for higher dimensional and nonlinear problems. Sometimes (but not always) the combination of the two methods is useful. A very good discussion using one or both of the above strategies can be found in Miller et al. [17]. Other notable articles for the problems of the type (1.1) are those of Boglaev [3], Miller [16], Niijima [21], O’Riordan and Stynes [24], Roos [26], Schatz and Wahlbin [30], Stojanovic´ [32,33], etc. The fitted operator methods were first suggested by Allen and Southwell [1] for solving the problem of viscous fluid flow past a cylinder. One type of exponentially fitted methods considered in Liniger and Willoughby [13] which is a special class of h-method given in Lambert [12] has been considered by Doolan et al. [5]. The details about the construction of the fitting factor in the above methods has been discussed in [5]. However, it is to be noted that the highest possible order of convergence in most of the above methods was only one (and two in some special cases). Recently, in [25] we have derived some e-uniform second-order and fourth-order (not e-uniform) methods just to fill the gap mentioned above. However, since the method presented in [25] is still not fourth-order e-uniform, we went through the first strategy and obtained, as we expected, the e-uniform results using the variable mesh version of the method (SNFDM) considered in [25]. Towards the use of the fitted mesh methods, tremendous works have been done so far. Whether one need to use a special scheme on a given mesh or a special mesh is to be used for a given scheme, can be answered by taking into account an excellent quote (by Stynes [35]): Miller has moved from [15] the question ‘‘what scheme should one use on a given mesh?’’ to [6] ‘‘what mesh should one use with a given scheme?’’ Hitherto, various mesh selection strategies (either graded or piecewise uniform meshes) have been provided in the literature, however, the most successful and widely used ones are those of Bakhvalov-type and Shishkin-type meshes. Bakhvalov [2] was the first to use an a priori mesh to solve a singular perturbation problem. His mesh was generalized by Vulanovic´ [37]. A mesh of this kind is generated by a suitable smooth mesh generating function which in the layer corresponds to the inverse of the function describing the layer behavior of the continuous solution. Bakhvalov meshes have been successfully applied to many different types of singular perturbation problem, see, e.g., [38–42]. A new impetus to the a priori mesh approach was given by Shishkin [31], who discovered that simple piecewise equidistant meshes can also produce e-uniform results. A mesh of this kind uses a fine equidistant mesh inside each layer, and a coarse, again equidistant, mesh on the outside. The Shishkin meshes are discussed in detail in [27,34,36] as well as in [28], where Roos surveyed results on layer-adapted meshes. It is shown in [42] that though Shishkin meshes are simpler than Bakhvalov meshes, the latter produce much better numerical results. The reason is that Bakhvalov meshes are better adopted to the layer structure. Comparisons of these two meshes made in [14,29,42,43] have shown that the standard three-point difference schemes produce much better results on Bakhvalov meshes than on Shishkin meshes. However, it is also well known that a higher-order scheme is much easier to construct on an equidistant mesh than on a fully nonequidistant one. The paper [44] shows that the real advantage of Shishkin meshes over Bakhvalov meshes is born out when applying more complicated higher order schemes. Since we are also working on higher order schemes, all of these comparisons suggest us therefore to use a Shishkin type mesh for the problem under consideration. To the best of our knowledge, the results presented in this paper are the first, fourth-order results for a self-adjoint spp (which is in conservation form). The systematic way of constructing the methods of any order higher than four (for such problems) is also a plus to this paper. We would like to mention at this stage that higher order uniform methods do exist for some other type of problems, e.g., semi-linear spps, see, [7,36,42].
K.C. Patidar / Applied Mathematics and Computation 188 (2007) 720–733
723
The rest of the paper is organized as follows. In Section 2, we reduce Eq. (1.1) to the normal form. Shishkin mesh (suitable for problem in this paper) has been described in Section 3. Section 4 deals with the design of a fourth-order fitted mesh finite difference scheme. Convergence analysis of this scheme has been given in Section 5. In Section 6, we explain, how can one get the schemes of order higher than four. Two numerical examples demonstrating the e-uniform convergence of the proposed method along with some comparisons with the results obtained via standard method have been given in Section 7. Finally, Section 8 deals with the conclusion and future directions. Notations and terminology Most of the notations and symbols we employ are fairly standard. However, to avoid notational complexities, we omit subscripted e from ye, Le and so on. 2. Reduction to the normal form Re-writing Eq. (1.1) as y 00 þ P ðxÞy 0 þ QðxÞy ¼ RðxÞ;
ð2:1Þ
where P ðxÞ ¼ Let
a0 ðxÞ ; aðxÞ
QðxÞ ¼
bðxÞ eaðxÞ
and
RðxÞ ¼
f ðxÞ : eaðxÞ
Z 1 x U ðxÞ :¼ exp P ðfÞ df : 2 0
Using yðxÞ ¼ U ðxÞV ðxÞ;
ð2:2Þ
we transform (2.1) into the normal form, that is, e ðxÞV ¼ RðxÞ; e V 00 þ W
ð2:3Þ
where e ðxÞ ¼ QðxÞ 1 P 0 ðxÞ 1 ðP ðxÞÞ2 ; W 2 4 Z x 1 e RðxÞ ¼ RðxÞ exp P ðfÞ df ; 2 0 with V ð0Þ ¼
yð0Þ ¼: a0 ; U ð0Þ
V ð1Þ ¼
yð1Þ ¼: a1 ; a0 ; a1 2 R: U ð1Þ
Multiplying (2.3) throughout by e, we get e LV eV 00 þ W ðxÞV ¼ ZðxÞ; V ð0Þ ¼ a0 ;
V ð1Þ ¼ a1 ;
ð2:4Þ
where e ðxÞ and W ðxÞ ¼ e W
e ZðxÞ ¼ e RðxÞ:
It is to be noted that the operator e L also satisfies the continuous maximum principle: Lemma 2.1. Let w(x) be any sufficiently smooth function such that w(0) P 0 and w(1) P 0. Then e LwðxÞ P 0, " x 2 (0, 1) implies that w(x) P 0, " x 2 [0, 1].
724
K.C. Patidar / Applied Mathematics and Computation 188 (2007) 720–733
Proof. Analogous to the proof of Lemma 1.1. To solve problem (2.4), we design the fitted mesh fourth-order method in Section 4. h 3. Description of the mesh Let n be a positive integer and denote by d the width of the boundary layer. Let the interval [0, 1] be divided into three sub-intervals: ½0; 1 :¼ ½0; d [ ½d; 1 d [ ½1 d; 1: The piecewise uniform mesh (of Shishkin type) in these three sub-intervals is designed as follows. Assuming that n = 2m with m P 5, the intervals (0, d) and (1 d, 1) are each divided into n/4 equal mesh elements, while the interval (d, 1 d) is divided into n/2 equal mesh elements. This insures that there are 1 + n/4 equidistant grid points into the intervals [0, d] and [1 d, 1] whereas one has n/2 1 equidistant grid points in (d, 1 d). We define the parameter d by n pffiffiffiffiffiffiffi o d ¼ min 1=4; 4 e=b ln ðn=16Þ : ð3:1Þ where 0 < b 6 W(x), " x 2 [0, 1]. Let j0 ¼ n=4;
xj0 ¼ d;
xnj0 ¼ 1 d
and ½0; 1 :¼ 0 ¼ x0 < x1 < < xj0 < < xnj0 < < xn ¼ 1; with xj = hj hj1, where the mesh spacing is given by 4dn1 ; j ¼ 1; . . . ; j0 ; n j0 þ 1; . . . ; n; ð3:2Þ hj ¼ 2ð1 2dÞn1 ; j ¼ j0 þ 1; . . . ; n j0 : pffiffiffiffiffiffiffi If d = 1/4, i.e., 1=4 < 4 e=b lnðn=16Þ then n1 is very small relative to e. This is very unlikely in practice and in such a case the method can be analyzed using the standard techniques. We therefore assume that pffiffiffiffiffiffiffi d¼4 e=b ln ðn=16Þ: ð3:3Þ 4. Fitted mesh finite difference scheme We use the notations Vj = V(xj), Wj = W(xj) and Zj = Z(xj) and denote the approximations of the Vj at the grid points xj by the unknowns mj. Taking the Taylor series expansion for m around xj, and neglecting the terms containing fifth and higher order terms, we get the following approximations for mj+1 and mj1: mjþ1 mj þ hjþ1 m0j þ
h2jþ1 00 h3jþ1 000 h4jþ1 ðivÞ m þ m þ m 2! j 3! j 4! j
ð4:1Þ
and mj1 mj hj m0j þ
h2j 00 h3j 000 h4j ðivÞ m m þ mj : 2! j 3! j 4!
ð4:2Þ
Adding Eqs. (4.1) and (4.2), we obtain mj1 2mj þ mjþ1 ðhjþ1 hj Þm0j þ
h2jþ1 þ h2j 00 h3jþ1 h3j 000 h4jþ1 þ h4j ðivÞ mj þ mj þ mj : 2! 3! 4!
ð4:3Þ
Neglecting third and fourth-order term in Eqs. (4.1) and (4.2), and solving the resulting system for m0j and m00j , we get
K.C. Patidar / Applied Mathematics and Computation 188 (2007) 720–733
725
m0j
h i 1 h2j mjþ1 ðh2j h2jþ1 Þmj h2jþ1 mj1 hj hjþ1 ðhj þ hjþ1 Þ
ð4:4Þ
m00j
2 ½hj mjþ1 ðhj þ hjþ1 Þmj þ hjþ1 mj1 hj hjþ1 ðhj þ hjþ1 Þ
ð4:5Þ
and
The relations 0
00 m000 j ¼ ðmj Þ
and
ðivÞ
mj
¼ ðm00j Þ
00
then gives the approximations for h i 1 m000 h2j m00jþ1 ðh2j h2jþ1 Þm00j h2jþ1 m00j1 j hj hjþ1 ðhj þ hjþ1 Þ
ð4:6Þ
ð4:7Þ
and ðivÞ
mj
2 ½hj m00jþ1 ðhj þ hjþ1 Þm00j þ hjþ1 m00j1 : hj hjþ1 ðhj þ hjþ1 Þ
ð4:8Þ
Using Eqs. (2.4), (4.4), (4.7) and (4.8) in (4.3) and simplifying, we obtain the tridiagonal system Am ¼ F
ð4:9Þ
where A is the matrix of the system and m and F are corresponding vectors. The various entries of this matrix and the components of the r.h.s. vector are given by 9 ðdiagðAÞÞj ¼ rcj ; j ¼ 1; 2; . . . ; n 1; > > > > > > ðsubdiagðAÞÞj ¼ rj ; j ¼ 2; 3; . . . ; n 1; > > > > > þ > ðsupdiagðAÞÞj ¼ rj ; j ¼ 1; 2; . . . ; n 2; > > = þ c F 1 ¼ q1 Z 0 þ q1 Z 1 þ q1 Z 2 r 1 m 0 ; > > > c þ > F j ¼ q Z þ q Z þ q Z ; j ¼ 2; 3; . . . ; n 2; > j1 j jþ1 j j j > > > > þ þ c > > F n1 ¼ q Z þ q Z þ q Z r m ; > n1 n n1 n n1 n2 n1 n1 > > ; m 0 ¼ a0 ; m n ¼ a1 ; where
; 12hjþ1
9 > > > 6 > > > > 3 3 > > 4 4 > h h h þh h j j j jþ1 jþ1 þ > 12hj ; > qj ¼ hjþ1 ðhj þhjþ1 Þ > 6 > > > > > > 3 3 4 4 2 2 2ðh h Þðh h Þðh þh Þ h þh = j jþ1 j j j jþ1 jþ1 jþ1 c ; qj ¼ 2 þ 12hj hjþ1 > n o > > hjþ1 ðhj hjþ1 Þ > ¼ e 1 W ; r þ q j1 > j j > hj ðhj þhjþ1 Þ > > > n o > > hj ðhj hjþ1 Þ þ þ > rj ¼ e 1 þ hjþ1 ðhj þhjþ1 Þ þ qj W jþ1 ; > > > > > n o > 2 > ðhj hjþ1 Þ ; c c rj ¼ e 2 þ hj hjþ1 þ qj W j : h
jþ1 q j ¼ hj ðhj þhjþ1 Þ
h3j h3jþ1
þ
ð4:10Þ
h4j þh4jþ1
ð4:11Þ
It is to be noted that when hj = hj+1 = h (i.e., uniform mesh throughout the region) then (4.11) reduces to ) 2 2 r rþ rcj ¼ 2e þ 5h2 W j =6; j ¼ e þ h W j1 =12; j ¼ e þ h W jþ1 =12; ð4:12Þ 2 h2 q qcj ¼ 5h6 : j ¼ 12 ;
726
K.C. Patidar / Applied Mathematics and Computation 188 (2007) 720–733
Using Eqs. (4.9)–(4.11) or (4.12) we get the approximate solution of V(x) at the grid points xj. Since U(x) is known, therefore the solution to the original problem (1.1) at these grid points will be obtained using (2.2). We call the method consisting of Eqs. (4.9)–(4.11) as the fitted mesh finite difference method (FMFDM) whereas the method consisting of Eqs. (4.9),(4.10) and (4.12) is the standard numerov’s finite difference method (SNFDM). 5. Convergence analysis of FMFDM Throughout the paper M will denote positive constants which may take different values in different equations (inequalities) but are always independent of h and e. Denoting by e L h the discrete operator in the FMFDM, we see that it satisfies e h /i P 0, Lemma 5.1 (Discrete maximum principle). For any mesh function /i satisfying /0 P 0, /n P 0 and L " 0 < i < n, we have /i P 0, " 0 6 i 6 n. and Lemma 5.2 (Uniformly stability estimate). If fi is any mesh function such that f0 = fn = 0, then jfi j 6
1 max j e L h fj j for b 16j6n1
0 6 i 6 n:
To facilitate the analysis properly, we decompose the solution V(x) into its regular (smooth) and singular components and provide the bounds on them and on their derivatives using Lemma [17]: Lemma 5.3. The solution Ve of the problem (2.4) can be decomposed into the form V e :¼ V e;r þ V e;s where for all k 2 [0, 6] and x 2 [0, 1], the regular (smooth) component Ve,r satisfies ðk2Þ=2 Eðx; bÞ; jV ðkÞ e;r ðxÞj 6 M½1 þ e
and the singular component Ve,s satisfies k=2 Eðx; bÞ; jV ðkÞ e;s ðxÞj 6 Me
where 0 < b 6 W ðxÞ
and
n pffiffiffiffiffiffiffi pffiffiffiffiffiffiffio Eðx; bÞ ¼ exp x b=e þ exp ð1 xÞ b=e :
Very often during the analysis of spps, one encounters with the exponential terms involving e and then divided by the terms involving various powers of e, which are the main cause of worry. For their careful consideration while proving the e-uniform convergence, we prove Lemma: Lemma 5.4. For a fixed mesh and for all integers k, we have pffiffi pffiffi exp Mxj = e exp Mð1 xj Þ= e lim max ¼ 0 and lim max ¼ 0; e!0 16j6n1 e!0 16j6n1 ek=2 ek=2 where xj = jh, h = 1/n, "j = 1(1)n 1. Proof. Consider the partition ½0; 1 :¼ f0 ¼ x0 < x1 < x2 < < xn1 < xn ¼ 1g: Clearly, for the interior grid points, we have pffiffi pffiffi pffiffi exp Mxj = e exp ðMx1 = eÞ exp ðMh= eÞ 6 ¼ max 16j6n1 ek=2 ek=2 ek=2
K.C. Patidar / Applied Mathematics and Computation 188 (2007) 720–733
727
and
pffiffi pffiffi pffiffi exp Mð1 xj Þ= e exp ðMð1 xn1 Þ= eÞ exp ðMh= eÞ 6 ¼M : max 16j6n1 ek=2 ek=2 ek=2 (as x1 = h, 1 xn1 = 1 (n 1)h = h). An application of L’Hospital’s rule then gives pffiffi expðMh= eÞ pk k! lim lim ¼ lim ¼ 0; pffi k=2 p!1 ðMhÞk expðMhpÞ e!0 e pð¼1= eÞ!1 expðMhpÞ
which completes the proof.
h
Now, the local truncation error of the scheme Eqs. (4.9)–(4.11) is c c þ 00 c 00 þ 00 ðe L h ðV mÞÞj ¼ ðr j qj W j1 ÞV j1 þ ðr j qj W j ÞV j þ ðr j qj W jþ1 ÞV jþ1 þ eðqj V j1 þ qj V j þ qj V jþ1 Þ:
Expanding Vj1, Vj+1 and their derivatives in terms of Vj and its derivatives, we obtain ð1Þ
ð2Þ
ð3Þ
ð4Þ
ð5Þ
ðe L h ðV mÞÞj ¼ T 0 V j þ T 1 V j þ T 2 V j þ T 3 V j þ T 4 V j þ T 5 V j þ T 6;1 V ð6Þ ðn1j Þ þ T 6;2 V ð6Þ ðn2j Þ þ T 6;3 V ð6Þ ðn3j Þ þ T 6;4 V ð6Þ ðn4j Þ;
ð5:1Þ
where n1j ; n3j 2 ðxj1 ; xj Þ, n2j ; n4j 2 ðxj ; xjþ1 Þ and c þ c þ T 0 ¼ ðr j þ rj þ rj Þ ðqj W j1 þ qj W j þ qj W jþ1 Þ; þ T 1 ¼ hjþ1 ðrþ j qj W jþ1 Þ hj ðrj qj W j1 Þ;
T2 ¼
h2j h2jþ1 þ c þ ðrj q ðrj qþ j W j1 Þ þ j W jþ1 Þ þ eðqj þ qj þ qj Þ; 2 2
T3 ¼ T4 ¼
h3j h3jþ1 þ þ ðrj q ðrj qþ W Þ þ j1 j j W jþ1 Þ þ eðhjþ1 qj hj qj Þ; 6 6
h4j h4jþ1 þ e 2 þ 2 ðrj q ðrj qþ W Þ þ j1 j j W jþ1 Þ þ ðhjþ1 qj þ hj qj Þ; 2 24 24
T5 ¼
ð5:2Þ
h5j h5jþ1 þ e 3 þ 3 ðrj q ðr qþ j W j1 Þ þ j W jþ1 Þ þ ðhjþ1 qj hj qj Þ; 6 120 120 j
T 6;1 ¼
h6j ðr q j W j1 Þ; 720 j
T 6;3 ¼
eh4j q j ; 24
T 6;4 ¼
T 6;2 ¼
h6jþ1 þ ðr qþ j W jþ1 Þ; 720 j
eh4jþ1 qþ j : 24
We have two cases now: either the grid point xj is inside the fine mesh, i.e., j 2 f1; . . . ; j0 1g [ fn j0 þ 1; . . . ; n 1g
ð5:3Þ
or it is inside the coarse mesh, i.e., j 2 fj0 ; . . . ; n j0 g:
ð5:4Þ
We denote by me,r and me,s, respectively, the regular and the singular components of the numerical solution. Then using (4.11) along with Eqs. (3.2), (3.3) into (5.2) and Lemma 5.3, we see that 8 6 6 n in case of ð5:3Þ; > < n ln 16 ; h p ffiffiffi ffi p ffiffiffi ffi e ð5:5Þ jð L ðV e;r me;r ÞÞj j 6 M xj b=e ð1xnj Þ b=e 0 > : n6 e 0 þe e ; in case of ð5:4Þ:
728
K.C. Patidar / Applied Mathematics and Computation 188 (2007) 720–733
(In the above, we have used the fact that the relation 2(1 2d)n1 6 2n1 is true for all step sizes and the quantities 1 + M/e are always less than the quantities M/e for different M.) Similarly, we see that 8 6 6 n in case of ð5:3Þ; > < n ln 16 ; h p ffiffiffi ffi p ffiffiffi ffi e jð L ðV e;s me;s ÞÞj j 6 M ð5:6Þ xj b=e ð1xnj Þ b=e 0 > : n6 e 0 þee2 ; in case of ð5:4Þ: Combining Eqs. (5.5), (5.1)–(5.6) and using Lemma 5.4 followed by the use of uniform stability estimate (Lemma 5.2), we obtain our main result stated in the following Theorem. Theorem 5.5. Let W(x), Z(x) be sufficiently smooth so that V(x) 2 C6[0, 1] and W(x) P b > 0. Let mj, j = 0(1)n, be the approximate solution of (2.4), obtained using Eqs. (4.9)–(4.11) with m0 = V(0), mn = V(1). Then, there is a constant M independent of e and h such that sup max jV j mj j 6 Mn6 ln6 ðn=16Þ 6 Mn4 ;
0
where we have used the fact that ln3(n/16) 6 M n, "n. Remark 5.6. On one hand, the case when n1 6 e then the second parts in the estimates Eqs. (5.5), (5.1)–(5.6) is bounded above, respectively, by Mn5 and Mn4. However, the case when n1 6 e in the region [d, 1 d] is practically not interesting. One the other hand, since xj0 ; ð1 xnj0 Þ 2 ð0; 1Þ, one can see by using Lemma 5.4 that the second parts inpffiffiffi(5.5) (andpffiffiffisimilarly in (5.6)) tends to zero (as e ! 0). However, in general the ffi ffi xj
b=e
ð1xnj Þ
b=e
0 above quantities like e 0 þe e are dominating and deteriorate the results in the pre-asymptotic range of e’s. The readers can see such influence in most of the articles devoted to singularly perturbed problems though not mentioned explicitly in any of the articles so far. But when we say that the method is uniformly convergent (or robust w.r.t. the parameter) then for any value of this parameter (whether it is 101 of 102 or 1010) the errors must remain same for a fix number of partitions of the interval [0, 1]. Theoretically, many authors provide the error bounds as M(hp) or M(hpln (hp)) for the pth order methods with M independent of h and e which appeals everyone very nice but in view of the numerics which they provide, it is really not always right to say so.
6. Schemes of the order higher than four In previous section, we have derived a fourth-order method. The main idea behind deriving the method of order p originated from Eqs. (4.1) and (4.2) in which one can take the expansions up to pth-order derivative terms. Because of variable step sizes, higher the value of p is taken, more the complications arise in deriving such schemes. However, we give here a fifth-order scheme in which we use the relations ðvÞ
ðivÞ
mj ¼ ðmj Þ0 ;
ðivÞ
mj
¼ ðm00j Þ00
and
00 0 m000 j ¼ ðmj Þ :
ð6:1Þ
This along with the procedure used earlier leads to a pentadiagonal system 1 c þ1 þ2 2 1 c þ1 þ2 r2 j mj2 þ rj mj1 þ r j mj þ r j mjþ1 þ rj mjþ2 ¼ qj Z j2 þ qj Z j1 þ qj Z j þ qj Z jþ1 þ qj Z jþ2 ;
ð6:2Þ
where q2 j ¼
hjþ1 h2j ðh5jþ1 h5j Þ ; 60hj hj1
qþ2 j ¼
h2jþ1 hj ðh5jþ1 h5j Þ ; 60hj hjþ1
q1 j
h2jþ1 ðh3jþ1 h3j Þ hjþ1 ðh4jþ1 þ h4j Þ hjþ1 ðh2j h2j1 Þðh5jþ1 h5j Þ h2jþ1 ðhjþ1 þ hj Þðh5jþ1 h5j Þ ¼ þ þ þ ; 6hj 12hj 60hj hj1 60h2j
qþ1 j
h2j ðh3jþ1 h3j Þ hj ðh4jþ1 þ h4j Þ hj ðh2jþ2 h2jþ1 Þðh5jþ1 h5j Þ h2j ðhjþ1 þ hj Þðh5jþ1 h5j Þ ¼ þ þ ; 6hj 12hj 60hj hjþ1 60h2j
K.C. Patidar / Applied Mathematics and Computation 188 (2007) 720–733
729
ðh2jþ1 þ h2j Þ ðh2jþ1 h2j Þðh3jþ1 h3j Þ ðhjþ1 þ hj Þðh4jþ1 þ h4j Þ þ 2 6hj 12hj ( ) ðh5jþ1 h5j Þ h2jþ1 hj ðh2jþ1 h2j Þðhjþ1 þ hj Þ hjþ1 h2j1 þ ; 60 hj hjþ1 hj hj1 h2j ! h2jþ1 ðhjþ1 hj Þ 2 þ2 þ2 1 þ q1 ¼ qj W j2 ; rj ¼ qj W jþ2 ; rj ¼ e 1 þ j W j1 ; hj ! h2j ðhjþ1 hj Þ ðh2jþ1 h2j Þðhjþ1 hj Þ c ¼e 1 W ; r ¼ eð2 þ Þ þ qcj W j þ qþ1 jþ1 j j hj hj
qcj ¼
r2 j rþ1 j
with hj = hjhj+1(hj + hj+1). One can use the scheme (6.2) at the grid points x2 to xn2, whereas at x1 and xn1 one can use the FMFDM derived in the previous section. In a similar manner, one can derive a sixth-order scheme and so on. 7. Test examples and numerical results To illustrate the theory, we give some numerical results corresponding to problem (1.1) for various choices of coefficient functions. Example 7.1. Consider problem (1.1) with aðxÞ ¼ 1 þ x2 ;
bðxÞ ¼ 1 þ xð1 xÞ; pffiffi
pffiffi pffiffi f ðxÞ ¼ 1 þ xð1 xÞ exp x= e xð2x2 3x þ 1Þ 2 e 2x2 x 1 þ e þ 1 pffiffi pffiffi pffiffi
þ exp ð1 xÞ= e x2 ð2x 1Þ þ 2 e 2x2 þ x e þ 1 ; yð0Þ ¼ 0; yð1Þ ¼ 0:
Its exact solution is given by
pffiffi pffiffi yðxÞ ¼ 1 þ ðx 1Þ exp x= e x exp ð1 xÞ= e : Example 7.2. Consider problem (1.1) with aðxÞ ¼ 1; bðxÞ ¼ 1; f ðxÞ ¼ ðcos2 px þ 2ep2 cos 2pxÞ;
yð0Þ ¼ 0;
yð1Þ ¼ 0:
Its exact solution is given by
pffiffi pffiffi pffiffi yðxÞ ¼ exp ð1 xÞ= e þ exp x= e = 1 þ exp 1= e cos2 px: Maximum errors at all the mesh points are evaluated using the formula: En;e :¼ max jyðxj Þ ~mj j; 06j6n
for different values of n and e, where ~mj ~mnj is the approximate solution of (1.1) obtained via Eqs. (2.4) and (2.2). The numerical rates of convergence are computed using the formula [5]: rk rk;e :¼ log2 ðEnk ;e =E2nk ;e Þ;
k ¼ 1; 2; . . . :
Further, we compute En ¼ max En;e 0
whereas the numerical rate of uniform convergence is computed as Rn :¼ log2 ðEn =E2n Þ:
730
K.C. Patidar / Applied Mathematics and Computation 188 (2007) 720–733
For these examples we present the numerical results obtained by the SNFDM (see Tables 1 and 4) followed by those obtained by FMFDM (see Tables 2 and 5). Tables 3 and 6 present the numerical rates of convergence.
Table 1 Results for Example 7.1 (Max. errors using SNFDM) e
n = 64
n = 128
n = 256
n = 512
n = 1024
0.1E01 0.1E02 0.1E03 0.1E04 0.1E05 0.1E08 0.1E12
0.11E05 0.75E04 0.59E02 0.98E01 0.18E+00 0.20E+00 0.20E+00
0.69E07 0.48E05 0.39E03 0.24E01 0.15E+00 0.20E+00 0.20E+00
0.43E08 0.30E06 0.26E04 0.24E02 0.68E01 0.20E+00 0.20E+00
0.27E09 0.19E07 0.16E05 0.15E03 0.12E01 0.20E+00 0.20E+00
0.17E10 0.12E08 0.10E06 0.98E05 0.92E03 0.20E+00 0.20E+00
Table 2 Results for Example 7.1 (Max. errors using FMFDM) e
n = 64
n = 128
n = 256
n = 512
n = 1024
0.1E01 0.1E02 0.1E03 0.1E04 0.1E05 0.1E06 0.1E07 0.1E08 0.1E09 0.1E10 0.1E11 0.1E12
0.11E05 0.46E04 0.63E03 0.45E02 0.10E01 0.13E01 0.14E01 0.15E01 0.15E01 0.15E01 0.15E01 0.15E01
0.69E07 0.48E05 0.84E05 0.91E04 0.40E03 0.67E03 0.80E03 0.84E03 0.86E03 0.86E03 0.86E03 0.86E03
0.43E08 0.30E06 0.10E05 0.11E05 0.11E04 0.31E04 0.43E04 0.48E04 0.49E04 0.50E04 0.50E04 0.50E04
0.27E09 0.19E07 0.15E06 0.15E06 0.19E06 0.11E05 0.22E05 0.27E05 0.29E05 0.29E05 0.30E05 0.30E05
0.17E10 0.12E08 0.20E07 0.19E07 0.19E07 0.28E07 0.97E07 0.15E06 0.17E06 0.17E06 0.18E06 0.18E06
En
0.15E01
0.86E03
0.50E04
0.30E05
0.18E06
Table 3 Results for Example 7.1 (rates of convergence using FMFDM) nk = 64, 128, 256, 512, 1024 e
r1
r2
r3
r4
r5
0.1E01 0.1E02 0.1E03 0.1E04 0.1E05 0.1E06 0.1E07 0.1E08 0.1E09 0.1E10 0.1E11 0.1E12
0.40E+01 0.34E+01 0.62E+01 0.56E+01 0.46E+01 0.43E+01 0.42E+01 0.41E+01 0.41E+01 0.41E+01 0.41E+01 0.41E+01
0.40E+01 0.40E+01 0.33E+01 0.63E+01 0.51E+01 0.44E+01 0.42E+01 0.41E+01 0.41E+01 0.41E+01 0.41E+01 0.41E+01
0.40E+01 0.40E+01 0.27E+01 0.31E+01 0.59E+01 0.47E+01 0.43E+01 0.41E+01 0.41E+01 0.41E+01 0.41E+01 0.41E+01
0.40E+01 0.40E+01 0.29E+01 0.29E+01 0.35E+01 0.53E+01 0.45E+01 0.42E+01 0.41E+01 0.41E+01 0.41E+01 0.41E+01
0.40E+01 0.40E+01 0.31E+01 0.31E+01 0.31E+01 0.38E+01 0.49E+01 0.43E+01 0.41E+01 0.41E+01 0.40E+01 0.39E+01
Rn
0.41E+01
0.41E+01
0.41E+01
0.41E+01
0.39E+01
K.C. Patidar / Applied Mathematics and Computation 188 (2007) 720–733
731
Table 4 Results for Example 7.2 (max. errors using SNFDM) e
n = 64
n = 128
n = 256
n = 512
n = 1024
0.1E01 0.1E02 0.1E03 0.1E04 0.1E05 0.1E08 0.1E12
0.48E06 0.45E04 0.38E02 0.54E01 0.95E01 0.10E+00 0.10E+00
0.30E07 0.28E05 0.27E03 0.15E01 0.78E01 0.10E+00 0.10E+00
0.19E08 0.18E06 0.18E04 0.17E02 0.39E01 0.10E+00 0.10E+00
0.12E09 0.11E07 0.11E05 0.11E03 0.78E02 0.10E+00 0.10E+00
0.74E11 0.70E09 0.70E07 0.69E05 0.67E03 0.99E01 0.10E+00
Table 5 Results for Example 7.2 (max. errors using FMFDM) e
n = 64
n = 128
n = 256
n = 512
n = 1024
0.1E01 0.1E02 0.1E03 0.1E04 0.1E05 0.1E06 0.1E07 0.1E08 0.1E09 0.1E10 0.1E11 0.1E12
0.48E06 0.74E04 0.41E03 0.21E02 0.40E02 0.50E02 0.54E02 0.55E02 0.55E02 0.55E02 0.55E02 0.55E02
0.30E07 0.28E05 0.79E05 0.50E04 0.17E03 0.26E03 0.30E03 0.31E03 0.31E03 0.32E03 0.32E03 0.32E03
0.19E08 0.18E06 0.69E06 0.78E06 0.54E05 0.12E04 0.16E04 0.17E04 0.18E04 0.18E04 0.18E04 0.18E04
0.12E09 0.11E07 0.18E06 0.11E06 0.11E06 0.49E06 0.83E06 0.98E06 0.10E05 0.11E05 0.11E05 0.12E05
0.74E11 0.70E09 0.46E07 0.14E07 0.14E07 0.14E07 0.39E07 0.54E07 0.60E07 0.63E07 0.63E07 0.64E07
En
0.55E02
0.32E03
0.18E04
0.12E05
0.64E07
Table 6 Results for Example 7.2 (rates of convergence using FMFDM) nk = 64, 128, 256, 512, 1024 e
r1
r2
r3
r4
r5
0.1E01 0.1E02 0.1E03 0.1E04 0.1E05 0.1E06 0.1E07 0.1E08 0.1E09 0.1E10 0.1E11 0.1E12
0.40E+01 0.48E+01 0.58E+01 0.54E+01 0.46E+01 0.43E+01 0.42E+01 0.41E+01 0.41E+01 0.41E+01 0.41E+01 0.41E+01
0.40E+01 0.40E+01 0.36E+01 0.60E+01 0.50E+01 0.44E+01 0.42E+01 0.42E+01 0.41E+01 0.41E+01 0.41E+01 0.41E+01
0.40E+01 0.40E+01 0.21E+01 0.30E+01 0.56E+01 0.46E+01 0.43E+01 0.41E+01 0.41E+01 0.41E+01 0.41E+01 0.39E+01
0.40E+01 0.40E+01 0.19E+01 0.29E+01 0.31E+01 0.51E+01 0.44E+01 0.42E+01 0.41E+01 0.41E+01 0.41E+01 0.43E+01
0.40E+01 0.40E+01 0.22E+01 0.31E+01 0.31E+01 0.32E+01 0.47E+01 0.43E+01 0.41E+01 0.41E+01 0.37E+01 0.24E+01
Rn
0.41E+01
0.41E+01
0.39E+01
0.43E+01
0.24E+01
8. Conclusions and future directions Using a fitted mesh finite difference method on a mesh of Shishkin type, we have described a numerical method for solving self-adjoint singular perturbation problems in conservation form. The method has been analyzed for e-uniform convergence. Two examples have been solved to demonstrate the predicted theory. A way to derive the methods of order higher than four is also explained in details.
732
K.C. Patidar / Applied Mathematics and Computation 188 (2007) 720–733 0.16
Diff. of Exact and SNFDM Solutions Diff. of Exact and FMFDM Solutions 0.14
0.12
0.1
0.08
0.06
0.04
0.02
0 0
0.1
0.2
0.3
0.4
0.5
0.6
0.7
0.8
0.9
1
Fig. 1. Errors of SNFDM and FMFDM solutions for Example 7.1 for n = 128 and e = 106.
As is expected the method is advantageous over the conventional methods (e.g., standard finite difference methods, etc). Comparative results obtained by our method are presented along with those obtained by the Standard Numerov’s method. Importance of the fitted mesh method can be seen from the respective tabular results. The solution of such singular perturbation problems possesses boundary (interior) layers in which the solution changes rapidly as e ! 0. On a uniform mesh, the SNFDM fails to capture this singularly perturbed nature of the solution in these layer regions. However, the fitted method takes care of it and the resulting behavior can be seen from Fig. 1. In this figure, we plot the difference between the exact solution and numerical solution(s) obtained via FMFDM and SNFDM. The method presented in this paper takes the advantage of the special nature of the differential Eq. (1.1) that allows one to reduce it to the normal form (2.4). Hence, any suitable method for Eq. (2.4) can be applied for (1.1) after some transformations. However, it would be interesting if such higher-order methods can be designed directly for (1.1). The same would be interesting for the nonlinear versions of it. Note that the higher-order methods for the nonlinear versions of (2.4) (some semi-linear spps) do exist in literature as pointed out earlier. The method proposed in this paper can also be extended for such problems by applying it to the sequence of quasi-linearized equations. However, in view of the current interest of the researchers, instead of this quasilinearization approach, we will prefer to have some direct methods for the general nonlinear problems. Some possible methods in this directions are currently under the investigation of the author. Acknowledgement The author is indebted to Prof. R. Vulanovic´ for his many stimulating suggestions. References [1] D. Allen, R.V. Southwell, Relaxation methods applied to determine the motion in 2-D of a viscous fluid past a fixed cylinder, Q. J. Mech. Appl. Math. VIII (2) (1955) 129–145. [2] N.S. Bakhvalov, Towards optimization of methods for solving boundary value problems in the presence of a boundary layer, Zh. Vychisl. Mat. Mat. Fiz. 9 (1969) 841–859 (in Russian). [3] I.P. Boglaev, A variational difference scheme for a boundary value problem with a small parameter in the highest derivative, USSR Comput. Math. Math. Phys. 21 (4) (1981) 71–81. [4] K.W. Chang, F.A. Howes, Nonlinear Singular Perturbation Phenomena: Theory and Application, Springer-Verlag, New York, 1984. [5] E.P. Doolan, J.J.H. Miller, W.H.A. Schilders, Uniform Numerical Methods for Problems with Initial and Boundary Layers, Boole Press, Dublin, 1980. [6] P.A. Farrell, A. Hegarty, J.J.H. Miller, E. O’Riordan, G.I. Shishkin, Robust Computational Techniques for Boundary Layers, Chapman & Hall/CRC, New York, 2000.
K.C. Patidar / Applied Mathematics and Computation 188 (2007) 720–733
733
[7] D. Herceg, Uniform fourth order difference scheme for a singular perturbation problem, Numer. Math. 56 (1990) 675–693. [8] M.K. Kadalbajoo, K.C. Patidar, A survey of numerical techniques for solving singulalry perturbed ordinary differential equations, Appl. Math. Comp. 130 (2-3) (2002) 457–510. [9] M.K. Kadalbajoo, K.C. Patidar, Singularly perturbed problems in partial differential equations: a survey, Appl. Math. Comp. 134 (23) (2003) 371–429. [10] J. Kevorkian, J.D. Cole, Perturbation Methods in Applied Mathematics, Springer-Verlag, New York, 1981. [11] J. Kevorkian, J.D. Cole, Multiple Scale and Singular Perturbation Methods, Springer-Verlag, New York, 1996. [12] J.D. Lambert, Computational Methods in Ordinary Differential Equations, John Wiley & Sons, New York, 1973. [13] W. Liniger, R.A. Willoughby, Efficient numerical integration methods for stiff system of differential equations, IBM Research Report RC-1970, 1967. [14] T. Linss, H.-G. Roos, R. Vulanovic´, Uniform pointwise convergence on Shishkin type meshes for quasiliner convection diffusion problems, SIAM J. Numer. Anal. 38 (2000) 897–912. [15] J.J.H. Miller, Construction of a FEM for a singularly perturbed problem in 2 dimensions, in Numerische Behandlung von Differentialgleichungen, Band 2 (Tagung, Math. Forschungsinst., Oberwolfach, 1975), Int. Ser. Numer. Math., Birkhau¨ser, Basel, 31 (1976), 165–169. [16] J.J.H. Miller, On the convergence uniformly in e, of difference schemes for a two-point boundary singular perturbation problem, Numerical Analysis of Singular Perturbation Problems, Academic Press, New York, 1979. [17] J.J.H. Miller, E. O’Riordan, G.I. Shishkin, Fitted Numerical Methods for Singular Perturbation Problems, World Scientific, Singapore, 1996. [18] K.W. Morton, Numerical Solution of Convection-Diffusion Problems, Chapman & Hall, London, 1996. [19] A.H. Nayfeh, Perturbation Methods, Wiley, New York, 1973. [20] A.H. Nayfeh, Introduction to Perturbation Techniques, John Wiley & Sons, New York, 1981. [21] K. Niijima, On a three-point difference scheme for a singular perturbation problem without a first derivative term I, Mem. Numer. Math. 7 (1980) 1–10. [22] R.E. O’Malley, Introductions to Singular Perturbations, Academic Press, New York, 1974. [23] R.E. O’Malley, Singular Perturbation Methods for Ordinary Differential Equations, Springer-Verlag, New York, 1991. [24] E. O’Riordan, M. Stynes, A uniformly accurate finite-element method for a singularly perturbed one-dimensional reaction-diffusion problem, Math. Comp. 47 (1986) 555–570. [25] K.C. Patidar, High order fitted operator numerical method for self-adjoint singular perturbation problems, Appl. Math. Comp. 171 (2005) 547–566. [26] H.-G. Roos, Global uniformly convergent schemes for a singularly perturbed boundary-value problem using patched base splinefunctions, J. Comput. Appl. Math. 29 (1) (1990) 69–77. [27] H.-G. Roos, M. Stynes, L. Tobiska, Numerical Methods for Singularly Perturbed Differential Equations, Springer-Verlag, Berlin, 1996. [28] H.-G. Roos, Layer adapted grids for singular perturbation problems, Z. Angew. Math. Mech. 78 (1998) 291–309. [29] H.-G. Roos, T. Linss, Sufficient conditions for uniform convergence on layer adapted grids, Computing 64 (1999) 27–45. [30] A.H. Schatz, L.B. Wahlbin, On the finite element method for singularly perturbed reaction diffusion problems in two and one dimension, Math. Comp. 40 (1983) 47–89. [31] G.I. Shishkin, A difference scheme for a singularly perturbed parabolic equation with a discontinuous boundary condition, Zh. Vychisl. Mat. Mat. Fiz. 28 (1988) 1679–1692 (in Russian). [32] M. Stojanovic´, A uniformly accurate spline collocation method for a singular perturbation problem, Calcolo 27 (1-2) (1990) 81–88. [33] M. Stojanovic´, Exponential splines difference scheme for singular perturbation problem with mixed boundary conditions, Comm. Appl. Numer. Meth. 7 (8) (1991) 625–632. [34] M. Stynes, H.-G. Roos, The midpoint upwind scheme, Appl. Numer. Math. 23 (1997) 361–374. [35] M. Stynes, A jejune heuristic mesh theorem, Comput. Meth. Appl. Math. 3 (2003) 488–492. [36] G. Sun, M. Stynes, An almost fourth order uniformly convergent difference scheme for a semilinear singlularly perturbed reaction diffusion problem, Numer. Math. 17 (1995) 487–500. [37] R. Vulanovic´, On a numerical solution of a type of singularly perturbed boundary value problem by using a special discretization mesh, Univ. Novom. Sadu Zb. Rad. Prirod. Mat. Fak. Ser. Mat 13 (1983) 187–201. [38] R. Vulanovic´, A uniform method for quasilinear singular perturbation problems without turning points, Computing 41 (1989) 97– 106. [39] R. Vulanovic´, A second order numerical method for nonlinear singular perturbation problems without turning points, J. Comput. Math. Math. Phys. 31 (1991) 522–532. [40] R. Vulanovic´, A uniform numerical method for a class of quasilinear turning point problems, in: J.J.H. Miller, R. Vichnevetsky (Eds.), Proceedings of the 13th IMACS World Congress on Computation and Applied Mathematics, Dublin, IMACS, 1991, pp. 493. [41] R. Vulanovic´, P.A. Farrell, Continuous and numerical analysis of a multiple boundary turning point problems, SIAM J. Numer. Anal. 30 (1993) 1400–1418. [42] R. Vulanovic´, Fourth order algorithms for a semilinear singular perturbation problem, Numer. Algor. 16 (1997) 117–128. [43] R. Vulanovic´, A priori meshes for singularly perturbed quasilinear two-point boundary value problems, IMA J. Numer. Anal. 21 (2001) 349–366. [44] R. Vulanovic´, A higher order scheme for quasilinear boundary value problems with two small parameters, Computing 67 (2001) 287– 303.