- Email: [email protected]
A uniformly convergent numerical method on non-uniform mesh for singularly perturbed unsteady Burger–Huxley equation
Available online at www.sciencedirect.com
Applied Mathematics and Computation 195 (2008) 688–706 www.elsevier.com/locate/amc
A uniformly convergent numerical method on non-uniform mesh for singularly perturbed unsteady Burger–Huxley equation q Aditya Kaushik a
a,*
, M.D. Sharma
b
INRIA, Mathe´matiques Applique´es de Bordeaux, University Bordeaux 1, 351 Cours de Libration F-33405, Cedex Talence, France b Department of Mathematics, Kurukshetra University, Kurukshetra 136 119, Haryana, India
Abstract In this paper, we present a numerical method which deals with the solvability of initial boundary value problem and related general properties of singularly perturbed Burger–Huxley equation. Singular perturbation problems are the problems in which a very small positive parameter is multiplied to the highest order derivative term. When this parameter approaches the limiting value of interest, i.e. ! 0, the problem exhibits boundary layers and most of the conventional methods fails to capture this effect. Thus, the quest for some new numerical techniques that may handle the difficulties occurring due to the presence of perturbation parameter and nonlinearity in the problem earns relevance. In this paper, one such parameter robust numerical scheme is constructed. The first step in this direction is the discretization of the time variable using the Euler implicit method with a constant time step. This produces a set of nonlinear stationary singularly perturbed semidiscrete problem class which is linearized further using the quasi-linearization process followed by spatial discretization using upwind finite difference operator on a piecewise uniform mesh. An extensive amount of analysis is carried out that uses a suitable decomposition of the error into smooth and singular component and a comparison principle combined with appropriate barrier functions. The error estimates are obtained, which ensures uniform convergence of the method. A set of numerical experiments is carried out in support of the predicted theory which validates the theoretical results computationally. Ó 2007 Elsevier Inc. All rights reserved. Keywords: Burger–Huxley equation; Piecewise uniform mesh; Uniform convergence; Singular perturbation; Quasi-linearization
1. Introduction In this paper, we construct a numerical method for solving time dependent singularly perturbed Burger– Huxley equation, L ½u ut þ auux uxx bð1 uÞðu cÞu ¼ 0; q
ð1:1aÞ
This work was completed with financial support by the UGC jointly with CSIR and HRDG Unit under Junior Research Fellowship vide No. F.17-39/98(SA-I). * Corresponding author. E-mail address: [email protected] (A. Kaushik). 0096-3003/$ - see front matter Ó 2007 Elsevier Inc. All rights reserved. doi:10.1016/j.amc.2007.05.067
A. Kaushik, M.D. Sharma / Applied Mathematics and Computation 195 (2008) 688–706
689
D fðx; tÞ 2 X ð0; T ð0; 1Þ ð0; T g; uðx; 0Þ ¼ UðxÞ; x 2 X ð0; 1Þ;
ð1:1bÞ ð1:1cÞ
uð0; tÞ ¼ uð1; tÞ ¼ 0;
ð1:1dÞ
t 2 ð0; T ;
where a; b P 0 and c 2 ð0; 1Þ and 0 < 1 are the parameters. The equation describes the interaction between convection, diffusion, and reaction. When a = 0 and = 1, Eq. (1.1) reduces to Huxley equation [1]. This equation describes nerve pulse propagation in nerve fibers and wall motion in liquid crystals [2–6]: ut uxx ðx; tÞ ¼ bð1 uÞðu cÞu: ð1:2Þ In the case when the epsilon is very small, i.e. ! 0, qualitative analysis of the model, which closely mimics the ionic processes at a real nerve membrane, is performed by the authors, by means of a singular perturbation theory in [7]. On the other hand, when b = 0, Eq. (1.1) reduces to the Burgers equation at high Reynolds number that establishes a balance between time evolution, nonlinearity and diffusion: ut þ auux ¼ uxx ðx; tÞ:
ð1:3Þ
This equation was first introduced by Burger [8] primarily to throw light on turbulence described by the interaction of two opposite effects of convection and diffusion. Burger equation arises in many physical problems including one dimensional turbulence, sound waves in viscous medium [9], shock waves in viscous medium, waves in a medium with fluid filled viscoelastic tubes and magnetohydrodynamic waves in a medium with finite electrical conductivity. This paper deals with the numerical study of generalized time dependent Burger–Huxley equation which is singularly perturbed from mathematical perspective. Parameter dependent differential equations such as (1.1) occur in diverse area of science and engineering. The concept of singular perturbation is not new, indeed, it has been a formidable tool in the solution of some important applied mathematical problems. In accordance with the informal principle, the behavior of solutions is governed primarily by the highest order terms, a solution u of the perturbed problem will often behave analytically quite differently from a solution of the original equation. When the perturbation parameter specifying the problem tends to zero, a breakdown occurs in the solution of the problem, and gives rise to boundary or the transition layers in the outflow boundary regions [10]. For higher values of a number of solution methodologies existing in the literature, but for sufficiently small , the existing solution methodology fails and a discrepancy occurs in the literature. In order to approximate the solution of such type of problems, numerical and asymptotic are two principal approaches. There are a wide variety of asymptotic expansion methods available for solving the problems of above type. But there can be difficulties in applying these asymptotic expansion methods, such as finding the appropriate asymptotic expansions in the inner and the outer regions, which are not routine exercises but require skill, insight and experimentation. Then one comes down ultimately to the numerical methods. It is a well established fact that nonlinear diffusion equations (1.1)–(1.3) with a small parameter plays an important role in nonlinear physics. Also it is of great practical interest to study the nonlinear phenomena. In order to solve Eq. (1.1), various numerical approaches were adopted by the researchers. Ismail et al. [11] studied the adomian decomposition method for Burger–Huxley and Burger–Fisher equation. Wang et al. [1] studied the solitary wave solutions of the generalized Burger–Huxley equation and Estevez [12] presented the non-classical symmetries and the singular modified solutions of the Burger and Burger–Huxley equation. More recently, Javidi [13] proposed a numerical solution of generalized Burger–Huxley equation by pseudospectral method and Darvishi’s preconditioning. However, it is well known for singular perturbation problem, that the numerical methods consisting of standard finite difference operators lead to severe restriction on the step size in order to attain the stability for small values of . To avoid these restrictions on the step size, a parameter uniform numerical scheme is constructed. The method can easily be extended to higher dimensions and for the nonlinear problems, but the only thing required in the process is the location and the exact width of the layer. 2. Bounds on the solution and time derivatives Since we shall use such a priori estimates later on in several different situations, we present in this section the a priori information needed in sequel. Moreover, the estimates theorems will be seen to be useful as a guide in the effective construction of the convergence and stability of the time discretization.
690
A. Kaushik, M.D. Sharma / Applied Mathematics and Computation 195 (2008) 688–706
Lemma 2.1. Let Uðx; tÞ; Wðx; tÞ 2 C 2;1 , i.e., twice continuously differential in space and once in time, and satisfy the relation jL½Uj < L½W;
ð2:4Þ
where L is the operator defined in Eq. (1.1). If jUj 6 W
ð2:5Þ
holds on the boundary, then jUj 6 W holds throughout the domain. Proof. Assume that ½U W attains the positive maximum at an interior point say (x, t) of the domain. Then according to the maximum principle ðL ½U WÞðx;tÞ 6 0, a contradiction to the assumption (2.4). Thus our assumption that ½U W attains positive maximum in the interior of the domain is wrong and hence ½U W does not have positive maximum within the domain, also from (2.5), along the boundary of the domain U W 6 0, it follows that U W 6 0 in D:
ð2:6Þ
Now suppose that U W attains a positive maximum in an interior point of the domain, then similarly it follows that U W 6 0
ð2:7Þ
in D:
Now from (2.7), (2.6) we have, jUj 6 W holding throughout the domain. Hence the result follows.
h
For obvious reasons the function W is called a barrier function for U. We impose the compatibility conditions u0 ð0; 0Þ ¼ u0 ð1; 0Þ ¼ 0, so that the data match at the corners (0, 0) and (1, 0) of the domain. As an application of this lemma we immediately obtain Theorem 2.2. There exists a number C independent of the perturbation parameter such that for all sufficiently small positive values of ; juðx; tÞ UðxÞj 6 Ct; t 2 ð0; T . Furthermore, juðx; tÞj 6 C;
ðx; tÞ 2 X ð0; T :
ð2:8Þ
Proof. We set gðx; tÞ ¼ uðx; tÞ UðxÞ, then g satisfies the differential equation L gðx; tÞ ¼ Uxx ðxÞ aUUx þ bð1 UÞðU cÞ ¼ H ðx; tÞ;
ð2:9Þ
where gðx; 0Þ ¼ 0; 0 < x < 1, gð0; tÞ ¼ uð0; tÞ Uð0Þ ¼ Uð0Þ;
gð1; tÞ ¼ uð1; tÞ Uð1Þ ¼ Uð1Þ;
0 6 t 6 1:
Now, clearly L ðCtÞ 6 C:
ð2:10Þ
Now, it can be immediately obtained from Lemma 2.1 and Eq. (2.10) that Ct is the barrier function for g for all small provided C is chosen so large that C P supðx;tÞ2D jH ðx; tÞj and so that Ct > jgj on the boundary. Thus jgðx; tÞj ¼ juðx; tÞ UðxÞj 6 Ct
ð2:11Þ
is uniform D and small . Using the properties of modulus, Eq. (2.11) is written as juðx; tÞj jUðxÞj 6 juðx; tÞ UðxÞj 6 Ct;
t 2 ð0; T :
Now since UðxÞ is sufficiently smooth and x and t also lie in the bounded interval, it is clear from the above relations that the solution u is bounded. h
Lemma 2.3. By keeping x fixed along the line fðx; tÞ : 0 6 t 6 T g, the bound ut is given as jut ðx; tÞj 6 C:
ð2:12Þ
A. Kaushik, M.D. Sharma / Applied Mathematics and Computation 195 (2008) 688–706
691
Proof. We assume that the solution uðx; tÞ is sufficiently smooth in the domain X ½0; T and by mean value theorem, there exits a t* in the interval ðt; t þ kÞ along the line fðx; tÞ : 0 6 t 6 T g such that uðx; t þ kÞ uðx; tÞ ; k 2juðx; tÞj jut ðx; t Þj 6 : k
ut ðx; t Þ ¼
ð2:13Þ
Using the Eq. (2.8), we get jut ðx; tÞj 6 C:
ð2:14Þ
Similarly we get the bounds of the utt ðx; tÞ and uttt ðx; tÞ along the line fðx; tÞ : 0 6 t 6 T g, i.e. we have i o uðx; tÞ ð2:15Þ oti 6 C for i ¼ 0; 1; 2; 3: 3. Time semidiscretization and quasi-linearization The first step in the numerical solution consists of discretizing the time variable with the Euler’s implicit method with constant time step Dt. This produces a set of stationary singularly perturbed problems of type Dt ujþ1 xx þ aDt½ujþ1 ½ujþ1 x þ ½ujþ1 bDtð1 ½ujþ1 Þð½ujþ1 cÞ½ujþ1 ¼ ½uj ; ð3:16aÞ ujþ1 ð0Þ ¼ 0
ujþ1 ð1Þ ¼ 0;
and
0 6 j 6 M 1:
ð3:16bÞ
Now for each j, (3.16) is a singularly perturbed nonlinear ordinary differential equation. To develop a numerical scheme for the boundary value problem Eq. (1.1), we first linearize the nonlinear problem (3.16) by using the quasi-linearization process. The nonlinear differential equation is linearized around a nominal solution of the nonlinear differential equation which satisfies the boundary conditions. Suppose uðkÞ ðxÞ is the nominal solution of the problem (3.16). Their application of the quasi-linearization process to the nonlinear problem (3.16) yields a sequence huðkÞ i of linear equations determined by the following recurrence relation: ðu00jþ1 Þ
ðkþ1Þ
þ ukþ1 ¼ f ðkÞ þ ððu0jþ1 Þ j
ðkþ1Þ
ðkÞ
ðkÞ
ðu0jþ1 Þ Þfu0
jþ1
þ ððujþ1 Þ
ðkþ1Þ
ðkÞ
ðujþ1 Þ ÞfuðkÞ ; jþ1
k ¼ 0; 1; . . .
Thus (3.16) followed by quasi-linearization leads to Lx; ujþ1 ¼ hðxÞ
ð3:17aÞ
with boundary condition ujþ1 ð0Þ ¼ 0
and
ujþ1 ð1Þ ¼ 0;
0 6 j 6 M 1:
ð3:17bÞ
where Lx; ¼
o2 o þ aðxÞ þ bðxÞ; ox ox2 h i
aðxÞ ¼ ak ¼ a ukjþ1 P b > 0; h h i i h i2 h i 1 P h > 0; bðxÞ ¼ bk ¼ þ a ukjþ1 b 2 ukjþ1 3 ukjþ1 c þ 2c ukjþ1 x Dt ½ukþ1 j þ a½ukjþ1 ½ukjþ1 x þ bð2½ukjþ1 3 ð1 þ cÞ½ukjþ1 2 Þ: hðxÞ ¼ Dt Remark 3.1. At first sight, the parameter Dt seems to behave like another perturbation parameter into the time semidiscretized problem, but under enough smoothness and compatibility requirements on the data of the continuous problem we have seen that this parameter does not have any severe effect as that of the perturbation parameter , to the multiscale character of the semidiscrete solutions.
692
A. Kaushik, M.D. Sharma / Applied Mathematics and Computation 195 (2008) 688–706 ðkþ1Þ
Now the boundary value problem (3.17) is linear in ujþ1 and now instead of solving the original nonlinear problem (3.16), we will solve the sequence of boundary value problems for singularly perturbed second order linear differential equations given by (3.17), for k ¼ 0; 1; 2; . . . and 0 6 j 6 M 1. Analytically, we requires as ðkÞ k tends to 1; ujþ1 ðxÞ to converge to the solution ujþ1 ðxÞ of the original nonlinear differential equation while numerically, we require that ðkþ1Þ ðkÞ ujþ1 ðxÞ ujþ1 ðxÞ < Tol:; 0 6 x 6 1; where Tol. is a small tolerance prescribed by us. If the above condition is satisfied, we terminate the iteration ðkþ1Þ and the profile ujþi ðxÞ is the numerical solution of the nonlinear boundary value problem (3.16). ðkÞ Now we prove the convergence of the sequence hujþ1 i. For the sake of convenience, consider eu00jþ1 ¼ f ðujþ1 Þ;
ð3:18aÞ
ujþ1 ð0Þ ¼ ujþ1 ðbÞ ¼ 0:
ð3:18bÞ ðkÞ
After quasi-linearization, we obtain a sequence hujþ1 i of linear equations determined by the recurrence eðu00jþ1 Þ
ðkþ1Þ
ðkþ1Þ
ðkÞ
ðkþ1Þ
ðkÞ
ðkÞ
ðxÞ ¼ f ðujþ1 Þ þ ðujþ1 ujþ1 Þf 0 ðujþ1 Þ;
ð3:19Þ
ðkþ1Þ
ujþ1 ð0Þ ¼ ujþ1 ðbÞ ¼ 0;
ð3:20Þ
where f 0 ðujþi Þ ¼ df ðujþi Þ=dujþi . Suppose ðk1Þ
ðkÞ
ð0Þ ujþ1 ðxÞ
ðk1Þ
is an initial approximation. Now from Eq. (3.19), we have ðk1Þ
eðu00jþ1 ÞðkÞ ðxÞ ¼ f ðujþ1 Þ þ ðujþ1 ujþ1 Þf 0 ðujþ1 Þ:
ð3:21Þ
Subtracting Eq. (3.21) from Eq. (3.19), we obtain ðkþ1Þ
ðkÞ
ðkÞ
00
ðk1Þ
ðkÞ
ðk1Þ
ðk1Þ
ðkþ1Þ
ðkÞ
ðkÞ
eðujþ1 ujþ1 Þ ¼ f ðujþ1 Þ f ðujþ1 Þ ðujþ1 ujþ1 Þf 0 ðujþ1 Þ þ ðujþ1 ujþ1 Þf 0 ðujþ1 Þ: ðkþ1Þ
ð3:22Þ
ðkÞ
Eq. (3.22) is a second order differential equation in ðujþ1 ujþ1 Þ. Converting this into an integral equation using Green’s function, we get Z b h i ðkþ1Þ ðkÞ ðkÞ ðk1Þ ðkÞ ðk1Þ ðk1Þ ðkþ1Þ ðkÞ ðkÞ eðujþ1 ujþ1 Þ ¼ Gðx; sÞ f ðujþ1 Þ f ðujþ1 Þ ðujþ1 ujþ1 Þf 0 ðujþ1 Þ þ ðujþ1 ujþ1 Þf 0 ðujþ1 Þ ds; 0
ð3:23Þ where Gðx; sÞ is the Green’s function and is defined by ( ðxbÞs ; 0 6 s 6 x 6 b; b Gðx; sÞ ¼ xðsbÞ ; 0 6 x 6 s 6 b; b
ð3:24Þ
and max jGðx; sÞj ¼ b=4:
ð3:25Þ
x;s
By the mean value theorem, we have ðkÞ
ðk1Þ
ðujþ1 ujþ1 Þ2 00 f ðnÞ; ð3:26Þ 2 ðk1Þ ðkÞ ðkÞ ðk1Þ where ujþ1 6 n 6 ujþ1 . Using Eq. (3.26) for f ðujþ1 Þ f ðujþ1 Þ in Eq. (3.23) followed by a simplification yields Z b ðkþ1Þ ðkÞ ðkÞ ðk1Þ 2 ðkþ1Þ ðkÞ ðkÞ eðujþ1 ujþ1 Þ ¼ Gðx; sÞ½ðujþ1 ujþ1 Þ f 00 ðnÞ=2 þ ðujþ1 ujþ1 Þf 0 ðujþ1 Þ ds: ð3:27Þ ðkÞ
ðk1Þ
ðkÞ
ðk1Þ
ðk1Þ
f ðujþ1 Þ f ðujþ1 Þ ¼ ðujþ1 ujþ1 Þf 0 ðujþ1 Þ þ
0
Let
max f 0 ðujþ1 Þ ¼ p;
kujþ1 k61
max f 00 ðujþ1 Þ ¼ q:
kujþ1 k61
ð3:28Þ
A. Kaushik, M.D. Sharma / Applied Mathematics and Computation 195 (2008) 688–706
693
Now taking the maximum of the moduli over the domain of consideration on both the sides of Eq. (3.27), we get ðkþ1Þ
ðkÞ
e max jujþ1 ujþ1 j 6
Z
06x6b
b 0
ðkÞ ðk1Þ 2 ðkþ1Þ ðkÞ ðkÞ max jGðx; sÞj maxðujþ1 ujþ1 Þ max jf 00 ðnÞj=2 þ max jujþ1 ujþ1 j max jf 0 ðujþ1 Þj ds: x;s
x
x
ð3:29Þ
Using Eqs. (3.25) and (3.28) in Eq. (3.29), we obtain b2 q ðkþ1Þ ðkÞ ðkÞ ðk1Þ 2 ðkþ1Þ ðkÞ maxðujþ1 ujþ1 Þ þ p max jujþ1 ujþ1 j : max jujþ1 ujþ1 j 6 x 06x6b 4e 2 x A rearrangement of terms yields b2 q ðkÞ ðk1Þ 2 maxðujþ1 ujþ1 Þ : x 8eð1 b2 p=4eÞ x This proves that the sequence hukjþ1 i of linear equations converges quadratically provided ðkþ1Þ
ðkÞ
max jujþ1 ujþ1 j 6
ð3:30Þ
b2 q < 1: 8eð1 b2 p=4eÞ
ð3:31Þ
The inequality (3.31) holds for sufficiently small b. If the interval ½0; b is too large, we can still keep ð1Þ ð0Þ ð0Þ jujþ1 ðxÞ ujþ1 ðxÞj sufficiently small by choosing a judicious initial approximation ujþ1 ðxÞ, which satisfies the ðkþ1Þ ðkÞ sufficient condition for the convergence, i.e., maxx jujþ1 ujþ1 j is small enough for some k. Remark 3.2. Thus we have shown that the sequence given by (3.17) of linear boundary value problems obtained after linearizing the nonlinear problem (3.16) converges to the original nonlinear boundary value problem. Therefore to obtain the approximate solution of the nonlinear problem (3.16), it is sufficient to approximate the linearized boundary value problem given by (3.17). 3.1. Convergence and stability i
Lemma 3.3. If j oto i uðx; tÞj 6 C, for ðx; tÞ 2 X ð0; T and 0 6 i 6 2. Then the local truncation error satisfies the estimate 2
kejþ1 k1 6 CðDtÞ : Proof. Now, since the solution of (3.17) is smooth enough, it holds Z tjþ1 ou o2 u uðtj Þ ¼ uðtjþ1 Þ Dt ðtjþ1 Þ þ ðtj nÞ 2 ðnÞ dn ot ot tj ¼ uðtjþ1 Þ Dtðuxx aðxÞux bðxÞuðx; tÞ þ f ðx; tÞÞðtjþ1 Þ þ OðDt2 Þ:
ð3:32Þ
Subtracting (3.16) from (3.32), and note that the local truncation error ejþi uðtjþ1 Þ ^ujþ1 at ðj þ 1Þth time step is the solution of boundary value problem of type Lx; ejþ1 ¼ OðDt2 Þ;
ð3:33aÞ
ejþ1 ð0Þ ¼ ejþ1 ð1Þ ¼ 0;
ð3:33bÞ
where ^ ujþ1 is the solution of the boundary value problem (3.17). Now, the operator Lx; clearly satisfies the maximum principle: consequently kL1 x; k1 6 C which ensures the stability of the semidiscretized scheme, where C is the positive constant independent of Dt. Thus we have kejþ1 k1 6 CðDtÞ2 :
ð3:34Þ
694
A. Kaushik, M.D. Sharma / Applied Mathematics and Computation 195 (2008) 688–706
Now, using the local error estimates up to the ðj þ 1Þth time step, followed by classical combination of stability and consistency property of the scheme, leads us to j X ek 6 CDt; ð3:35Þ kEj k1 ¼ k¼1 1
since jDt 6 T , Here C is a positive constant independent of and Dt. Thus we have proved the following: Theorem 3.4. Under the hypothesis of the preceding lemma, the global truncation error satisfies the estimate kEj k1 6 CDt:
ð3:36Þ
4. A priori estimates for semidiscretized problem In this section of the paper, we meet the important and recurring concept of an a priori estimate. Also we discuss qualitative properties of the solution of the semidiscretized problem. We now give to the maximum principle a mathematical basis, that can be a useful tool for deriving a priori bounds on the solutions of the differential equations and their derivatives. Rewriting the semidiscretized equation k k kþ1 kþ1 k Lx; ½ukþ1 jþ1 xx þ a ðxÞ½ujþ1 x þ b ðxÞ½ujþ1 ¼ f ðxÞ;
ujþ1 ð0Þ ¼ 0
and
ujþ1 ð1Þ ¼ 0;
06j6M 1
and the differential operator satisfies the following continuous maximum principal: Lemma 4.1 (Maximum Principle). If a function U, twice continuously differentiable in space and once in time, i.e. U 2 C 2;1 ðXÞ be such that Uðx; tÞ P 0, for all ðx; tÞ 2 oD, then Lx; Uðx; tÞ P 0 for all ðx; tÞ 2 D implies that Uðx; tÞ P 0 for all ðx; tÞ 2 D. Proof. Let ðx ; t Þ be such that Uðx ; t Þ ¼ minx2D Uðx; tÞ, and suppose that Uðx ; t Þ < 0. It is then clear from the hypothesis that ðx ; t Þ 62 oD, from this it follows ðx ; t Þ 2 D, hence Ux ðx ; t Þ ¼ 0, and Uxx ðx ; t Þ P 0, further Lx; Uðx ; t Þ ¼ Uxx ðx ; t Þ þ ak ðxÞUx ðx ; t Þ þ bk ðxÞUðx ; t Þ 6 0;
ð4:37Þ
which is a contradiction to our assumption, hence our assumption Uðx ; t Þ < 0 is wrong. Hence, Uðx; tÞ P 0 for all x 2 D. Hence the required result holds. h An immediate and important consequence of the maximum principle is the uniqueness and continuous dependence of the solution on the boundary values. For let hvi;jþ1 i and hwi;jþ1 i be the two sets of solutions of the semidiscretized difference equations satisfying the boundary conditions. Suppose zi;jþ1 ¼ hvi;jþ1 i hwi;jþ1 i, then since the semidiscrete operator is linear, we have zi;jþ1 ¼ hvi;jþ1 i hwi;jþ1 i satisfies Lx; zi;jþ1 ¼ Lx; vi;jþ1 Lx; wi;jþ1 ¼ 0 and obviously, z0;jþ1 ¼ 0; zN ;jþ1 ¼ 0. Thus the mesh function zi;jþ1 satisfies the hypothesis of discrete maximum principle that leads to zi;jþ1 P 0
for all i ¼ 1; 2; . . . ; N :
ð4:38Þ
On the other hand, if we set zi;jþ1 ¼ ðhvi;jþ1 i hwi;jþ1 iÞ, then a similar treatment to hzi;jþ1 i with an application of maximum principle leads us to zi;jþ1 P 0:
ð4:39Þ
Combining Eqs. (4.38) and (4.39), we have zi;jþ1 ¼ 0 and uniqueness follows immediately. Further since the problem under consideration is linear, the existence follows from uniqueness. The maximum principle can also be used to derive estimates on the derivatives of the solution provided some additional conditions are placed on the equation. In the next two lemmas we establish a priori bound on the exact solution and its derivatives, of time semidiscretized problem class (3.17). It is convenient to use the terminology suggested by the comparison principle combined with appropriate barrier functions.
A. Kaushik, M.D. Sharma / Applied Mathematics and Computation 195 (2008) 688–706
695
Lemma 4.2. Let u be the solution of any problem from semidiscretized problem class (3.17), then ku k1 6 C:
ð4:40Þ
Proof. Consider the barrier functions, f ¼ 1h khk1 u ðx; tÞ. Note that the functions f are non-negative at x ¼ 0; 1 for all t 2 ð0; T , and that for all ðx; tÞ Lx; f ðx; tÞ ¼ bðxÞh1 khk1 hðxÞ P 0;
ð4:41Þ
since bðxÞ P h > 0 and khk1 P hðxÞ. Thus maximum principle gives f ðx; tÞ P 0 and so, 1 ku k1 6 khk1 ¼ C h
for all ðx; tÞ 2 D:
Lemma 4.3. The derivatives uðkÞ of the solution u of semidiscretized problem satisfy the following bounds: pffiffi pffiffi kuk kX 6 Cð Þk ð1 þ ð Þk Þ maxfkuk; khkg; k ¼ 1; 2; pffiffi 3 pffiffi 3 ku3 kX 6 Cð Þ ð1 þ ð Þ Þ maxfkuk; khk; kh0 kg:
ð4:42Þ ð4:43Þ
Proof. Given any x 2 ð0; 1Þ, we can construct the neighborhood N x ¼ ðp; p þ rÞ such that x 2 N x and N x ð0; 1Þ. The mean value theorem implies that there exists x 2 N x such that u0 ðx Þ ¼
uðp þ rÞ uðpÞ r
and so ju0 ðx Þj 6
2kuk : r
Now we have 0
0
u ðxÞ ¼ u ðx Þ þ
ð4:44Þ
Z
x
u00 ðsÞ ds x
and therefore from the time semidiscretized differential Eq. (3.17), we have Z x Z x u0 ðxÞ ¼ u0 ðx Þ 1 hðsÞ ds þ 1 ½aðsÞu0 ðsÞ þ bðsÞuðsÞ ds: x
Now integration by parts leads us to Z Z x aðsÞu0 ðsÞ ds ¼ aðsÞujx a0 ðsÞuðsÞ : x
ð4:45Þ
x
ð4:46Þ
x
Now using the fact that the maximum norm of a function is always greater than the value of the over the domain of consideration, followed by calculation, yields Z x 0 0 ð4:47Þ aðsÞu ðsÞ ds 6 ð2kak þ ka kÞkuk: x
Substituting (4.48) in (4.45), we have Z x Z x Z x x u0 ðxÞ ¼ u0 ðx Þ 1 hðsÞ ds þ 1 bðsÞuðsÞ ds þ 1 aðsÞujx a0 ðsÞuðsÞ ds : x
x
x
Now since jx x j 6 r, and from (4.44) followed by (4.47) on simplification leads to 1 r 1 r 1 r 1 0 0 ju ðxÞj 6 C þ þ kuk þ kf k 6 ju ðxÞj 6 C þ þ maxfkuk; kf kg: r r
ð4:48Þ
696
A. Kaushik, M.D. Sharma / Applied Mathematics and Computation 195 (2008) 688–706
pffiffi If we choose r ¼ , then the right-hand side of the above expression is minimized with respect to r, and we obtain the required result for k = 1. Using the differential Eq. (3.17) for u, we can obtain the required bounds for k = 2 and on differentiating again (3.17) the result for k = 3 follows. h In order to have parameter uniform convergence results one needs to prove the stronger pointwise error bounds. To find these sharper pointwise bounds the solution of the semidiscretized problem has to be decomposed into smooth and singular components. The solution u of the semidiscrete problem is decomposed into a smooth component v and a singular component w as follows: u ¼ v þ w ; where v can be written in the form v ¼ v0 þ v1 þ 2 v2 ; and v0 ; v1 and v2 are defined, respectively, to be the solutions of the problems aðxÞv00 ðxÞ þ bðxÞv0 ðxÞ ¼ hðxÞ;
x 2 X;
aðxÞv01 ðxÞ
v1 ð1Þ ¼ 0;
v002 ðxÞ
þ bðxÞv1 ðxÞ ¼
þ
aðxÞv02 ðxÞ
v000 ;
þ bðxÞv2 ðxÞ ¼
v0 ð1Þ ¼ u ð1Þ ¼ 0;
v001 ;
ð4:49aÞ ð4:49bÞ
v2 ð0Þ ¼ v2 ð1Þ ¼ 0:
ð4:49cÞ
Thus the smooth component is the solution of Lx; v ¼ hðxÞ;
v ð0Þ ¼ v0 ð0Þ þ v1 ð0Þ;
v ð1Þ ¼ u ð1Þ ¼ 0
ð4:50Þ
and consequently singular component satisfies Lx; w ¼ 0;
w ð0Þ ¼ u ð0Þ v ð0Þ;
w ð1Þ ¼ 0:
ð4:51Þ
Lemma 4.4. For 0 < k 6 3; v0 satisfies the following estimates: ðkÞ
kv0 k 6 C: Proof. The problem (4.49a) can be written as v00 ðtÞ þ pðtÞv0 ðtÞ ¼ f ðtÞpðtÞ;
v0 ð1Þ ¼ 0;
ð4:52Þ
1
where pðtÞ ¼ bðxÞaðxÞ . The problem (4.52) isR a first order linear differential equation in v0, in order to obtain the solution we multiply Eq. (4.52) by exp pðtÞ dt , which followed by integration from x to 1, for some x 2 ð0; 1Þ, yields Z 1 Z 1 Z v0 ðtÞ exp pðtÞ dt ¼ f ðtÞ exp pðtÞ dt pðtÞ dt x
x
which on simplification reduces to Z Z 1 1 R v0 ðxÞ ¼ f ðtÞ exp pðtÞ dt pðtÞ dt: exp pðtÞ dt t¼x x
ð4:53Þ
Now since aðtÞ is sufficiently smooth and we have aðtÞ > h > 0 for all t 2 ½0; 1 also f ðtÞ is bounded for all t 2 ½0; 1. Therefore this boundedness of v0 follows immediately. Again from Eq. (4.49a), we have v00 ðxÞ ¼ f ðxÞpðxÞ pðxÞv0 ðxÞ and the boundedness of v00 ðxÞ follows directly from the boundedness of v0. Bounds on v000 ðxÞ and v000 0 ðxÞ can be easily obtained by successively differentiating Eq. (4.49a). Thus for 0 < k 6 3; v0 satðkÞ isfies the estimates kv0 k 6 C. h Lemma 4.5. For 0 < k 6 3; v1 satisfies the following estimates: ðkÞ
kv1 k 6 C:
A. Kaushik, M.D. Sharma / Applied Mathematics and Computation 195 (2008) 688–706
697
Proof. Since v1 is the solution of the first order linear differential Eq. (4.49b) and since v0 is bounded, using the classical argument used in Lemma 4.4 it is straightforward to prove that for 0 < k 6 3; v1 satisfies the required estimates. h Theorem 4.6. For 0 < k 6 3; v2 satisfies the following estimates: hð1 xÞ ðkÞ k kv2 k 6 C 1 þ exp for all x 2 X:
ð4:54Þ
Now we need to obtain the bound on the singular component of the solution. For if, introducing the two barrier functions hð1 xÞ w :¼ C exp w ;
then we have w ð0Þ ¼ C exp hð1Þ 0 P 0, and we choose C so that it holds w ð1Þ ¼ C w P 0. Further h hð1 xÞ hð1 xÞ ðnÞ ðnÞ L W ðxÞ ¼ h þ a ðxÞ exp þ b ðxÞ exp > 0; ð4:55Þ since aðnÞ ðxÞ P h > 0 and bn ðxÞ P b > 0. Therefore from Lemma 4.1 we obtain w P 0 for all x 2 X which gives us hð1 xÞ kw ðxÞk 6 C exp for all x 2 X; ð4:56Þ Z x Z x ðnÞ ðaðw ÞxÞðsÞ ds ¼ aðnÞ w ðxÞ ; 6 C exp hð1 xÞ : a ðw Þ ð4:57Þ ðsÞ ds x x x 0 0 By mean value theorem in the interval ð0; Þ, there exist z 2 ð0; Þ such that w ðÞ w ð0Þ 6 C1 : jw ðzÞj ¼ Integrating Eq. (4.51) w.r.t. x from 0 to x, we have Z x Z x ðnÞ ðw Þx ðxÞ þ ðw Þx ð0Þ þ a ðw Þx ðsÞ ds þ bðnÞ w ðsÞ ds ¼ 0: 0
ð4:59Þ
0
The above equation gives Z z Z z ðnÞ ðnÞ jðw Þx ð0Þj ¼ ðw Þx ðzÞ a ðw Þx ðsÞ ds b w ðsÞ ds: 0
ð4:58Þ
ð4:60Þ
0
Making use of Eqs. (4.57) and (4.58), we have hð1 Þ jðw Þx ð0Þj 6 C exp 6 C:
ð4:61Þ
Now Eq. (4.59) gives the following estimate with the help of the estimates given by (4.57) and (4.61) jðw ðxÞÞx j 6 C1 ;
x 2 X:
ð4:62Þ
Using the above estimate in Eq. (4.51), we get the estimate for ðw ðxÞÞxx , i.e. jðw Þxx j 6 C2
ð4:63Þ
Similarly by differentiating Eq. (4.51) w.r.t. x and using the estimates given by Eqs. (4.62) and (4.63), we get the estimate for ðw Þxxx ðkÞ
jðw Þ j 6 Ck
for k ¼ 1; 2; 3:
Now, the above results proven by us are summarized in the following theorem:
ð4:64Þ
698
A. Kaushik, M.D. Sharma / Applied Mathematics and Computation 195 (2008) 688–706
Theorem 4.7. The solution u of the semidiscretized problem admits the decomposition u :¼ v þ w ; where for all k; 0 6 k 6 3 and for all x 2 X, the smooth component of the solution v satisfies ðkÞ
jðv Þ ðxÞj 6 Cð1 þ 2k Þ
ð4:65Þ
and the singular component of the error w satisfies ðkÞ
jðw Þ j 6 Ck
ð4:66Þ
for some constant C independent of the perturbation parameter. 5. The spatial discretization 5.1. Classical difference scheme In order to have a better insight into the problem, we first consider the classical finite difference method that constitutes the base of the construction of the parameter uniform scheme later on. So let us first introduce a non-uniform mesh of the nodal points xi in the spatial direction. Let K 1 be the number of the intervals. Define hi ¼ xiþ1 xi ¼ max hi : i
For the continuous problem we use the following difference scheme [14]: Lj ¼ f ðxÞ:
ð5:67Þ
Here, Lj ¼ Dþ D ui;jþ1 þ aðxi ÞDþ ui;jþ1 þ bðxi Þui;jþ1 ; þ
þ
ð5:68aÞ
D D ui;jþ1 ¼ 2ðD ui;jþ1 D ui;jþ1 Þ=ðhi þ hiþ1 Þ;
ð5:68bÞ
D ui;jþ1 ¼ ðui;jþ1 ui1;jþ1 Þ=hi ; Dþ ui;jþ1 ¼ ðuiþ1;jþ1 ui;jþ1 Þ=hi þ 1:
ð5:68cÞ ð5:68dÞ
The difference scheme so defined is monotone. By means of the maximum principle and taking into account the a priori estimates on the derivative and the solution of the semidiscretized scheme, we find that the solution of the difference scheme converges for a fixed value of kuðx; tÞ jðx; tÞk 6 Cð2 N 1 þ M 1 Þ:
ð5:69Þ
This error bound for the classical scheme is clearly not -uniform. The proof of (5.69) follows the lines of the classical convergence proof for the difference schemes [14]. 5.2. The -uniformly convergent scheme In this paper, we are considering the upwind discretization on piecewise uniform fitted meshes. Since in practice we often want to construct a numerical method whose system matrix is a M-matrix. In this section, we discretize the boundary value problem (3.17) using the fitted mesh finite difference method composed of a standard upwind finite difference operator on a piecewise uniform mesh condensing at the boundary points. N The fitted piecewise uniform mesh XNx ¼ fxi gi¼0 on the interval [0, 1] is constructed by partitioning the interval into two subintervals, ð0; 1 dÞ; ð1 d; 1Þ. The resulting piecewise uniform mesh depends on one parameter which is called the transition parameter d and is chosen such that d min½0:5; C ln N ;
ð5:70Þ k
where C is the constant whose value depends upon the method applied. Also it is assumed that N ¼ 2 , where k P 2 is an integer, which guarantees that there is at least one point in the boundary regions. On each of these subintervals a uniform mesh is then placed.
A. Kaushik, M.D. Sharma / Applied Mathematics and Computation 195 (2008) 688–706
699
Define xi ¼
ihi
for 0 6 i 6 N =2;
ð1 dÞ þ ði N =2Þhi
for N =2 þ 1 6 i 6 N ;
where N is the number of mesh points. 2ð1 dÞ=N for 0 6 i 6 N =2; hi ¼ 2d=N for N =2 þ 1 6 i 6 N : The fitted mesh finite difference method for the time semidiscretized problem (4) at ðj þ 1Þth time step on the piecewise uniform mesh XNx is defined by LNx; ui;jþ1 ¼ hðxi Þ;
ð5:71Þ
u0;jþ1 ¼ uN ;jþ1 ¼ 0;
ð5:72Þ
where the discrete operator is defined as LNx; ui;jþ1 ¼ Dþ D ui;jþ1 þ aðxi ÞDþ ui;jþ1 þ bðxi Þui;jþ1 :
ð5:73Þ
6. Convergence and stability analysis 6.1. Stability analysis A discrete maximum principle can also be established directly, without appealing to the properties of the elements in the system matrix. Analogously to the continuous case, this is usually proved by the method of contradiction as we now show. Theorem 6.1 (Discrete Maximum Principle). Assume that mesh function Wi;jþ1 satisfies W0;jþ1 P 0 and WN ;jþ1 P 0. Then ðLNx; ÞWi;jþ1 P 0 for all xi 2 XNx implies that Wi;jþ1 P 0 for all xi 2 XNx . Proof. Assume that there exists a positive integer k such that Wk;jþ1 ¼ min06i6N Wi;jþ1 and assume Wk;jþ1 < 0. Then we have W0;jþ1 P 0 and WN ;jþ1 P 0, therefore it follows from the hypothesis that k 62 f0; N g. Also Wk;jþ1 Wk1;jþ1 6 0 and Wkþ1;jþ1 Wk;jþ1 P 0 and note that ðLNx; ÞWk;jþ1 < 0 since Wk;jþ1 < 0 by assumption and bðxi Þ > h > 0, which is a contradiction to our assumption. Therefore the assumption Wk;jþ1 < 0 is wrong, hence Wk;jþ1 P 0. We have chosen k to be fixed and arbitrary, so Wi;jþ1 P 0 for all xi 2 XNx . h As an immediate application of the maximum and minimum principle, we have the following uniform stability estimates: Lemma 6.2. Let Z i;jþ1 be any mesh function such that. Z 0;jþ1 ¼ Z N ;jþ1 ¼ 0, then for all i; 0 6 i 6 N jZ i;jþ1 j 6 h1 max jLNx; Z k;jþ1 j: 16k6N 1
Proof. To prove the estimate we construct two mesh functions W i;jþ1 ¼ Mxi;jþ1 Z i;jþ1 ; where M ¼ h1 max16k6N 1 jðLNx ÞZ k;jþ1 j. Then we clearly have 1 W max jLNx; Z k;jþ1 jxi;jþ1 Z 0;jþ1 ¼ 0; 0;jþ1 ¼ h
ð6:74aÞ
16k6N 1
1 max jLNx; Z k;jþ1 jxN ;jþ1 Z N ;jþ1 ¼ h1 max jLNx; Z k;jþ1 jxN ;jþ1 P 0 W N ;jþ1 ¼ h 16k6N 1
16k6N 1
ð6:74bÞ
700
A. Kaushik, M.D. Sharma / Applied Mathematics and Computation 195 (2008) 688–706
and for 1 6 i 6 N 1, ðLNx; ÞW i;jþ1 ¼ ðaðxi Þ þ bðxi Þxi;jþ1 Þ
h1
max jðLNx; ÞZ k;jþ1 j LNx; Z i;jþ1 P 0;
16k6N 1
Thus an application of the discrete maximum principle leads to the required estimates.
h
6.2. Convergence analysis The discrete problem is decomposed in an analogous manner as the second decomposition of the solution u of the continuous problem, in order to estimate the error in the smooth and singular component separately. Thus uN ðxi;jþ1 Þ ¼ vN ðxi;jþ1 Þ þ wN ðxi;jþ1 Þ; where vN is the solution of the inhomogeneous problem LNx; vN ðxi;jþ1 Þ ¼ f ðxi;jþ1 Þ; vN ðx0;jþ1 Þ ¼ 0 and
vN ðx1;jþ1 Þ ¼ 0
and wN ðxi;jþ1 Þ is the solution of the homogeneous problem LNx; wN ðxi;jþ1 Þ ¼ 0; wN ðx0;jþ1 Þ ¼ 0
and
wN ðx1;jþ1 Þ ¼ 0:
Then the error at each time step (j + 1) can be written in the form ðuN u Þðxi;jþ1 Þ ¼ ðvN v Þðxi;jþ1 Þ þ ðwN w Þðxi;jþ1 Þ: Now in the next two results, we estimate the error in the smooth and singular component, respectively. Lemma 6.3. At each time step (j + 1) and mesh point xi 2 XNx , the smooth component of error satisfies jðvN v Þðxi;jþ1 Þj 6 CN 1 :
ð6:75Þ
Proof. This is obtained using the following standard stability and consistency argument, we consider the local truncation error 2 o d þ LNx; ðvN v Þðxi;jþ1 Þ ¼ D D D ðx Þ þ aðx Þ v v ðxi;jþ1 Þ: i;jþ1 i dx ox2 At each time step, to reduce the order of differentiation in the integral, integration by parts leads us to Z xi;jþ1 o 1 D ðxi1;jþ1 sÞv00 ðsÞ ds: v ðxi;jþ1 Þ ¼ ox xi;jþ1 xi1;jþ1 xi1;jþ1 Taking the modulus on both the sides, we have Z xi;jþ1 00 o 1 6 jv ðxi;jþ1 Þj D ðx Þ ðs xi1;jþ1 Þ ds ¼ jv00 ðxi;jþ1 Þjðxi;jþ1 xi1;jþ1 Þ: v i;jþ1 ox x 2 x i;jþ1
i1;jþ1
ð6:76Þ
xi1;jþ1
Also, integrating by parts twice we have 2 3 R xi;jþ1 2 000 1 2 ðx sÞ v ðsÞ ds iþ1;jþ1 x o 1 x x 4 iþ1;jþ1 i;jþ1 Riþ1;jþ1 5: Dþ D 2 v ðxi;jþ1 Þ ¼ xi1;jþ1 1 xiþ1;jþ1 xi1;jþ1 ox ðs xi1;jþ1 Þ2 v000 ðsÞ ds xi;jþ1 xi1;jþ1
xi;jþ1
From this it follows that 3 R xiþ1;jþ1 000 2 2 2 1 ðxiþ1;jþ1 sÞ ds v ðxi;jþ1 Þ o xiþ1;jþ1 xi;jþ1 xi;jþ1 þ 4 5 R xi;jþ1 ox2 D D v ðxi;jþ1 Þ 6 x 2 1 iþ1;jþ1 xi1;jþ1 þ ðs xi1;jþ1 Þ ds xi;jþ1 xi1;jþ1
1 6 ðxiþ1;jþ1 xi1;jþ1 Þv000 ðxi;jþ1 Þ : 3
xi1;jþ1
ð6:77Þ
A. Kaushik, M.D. Sharma / Applied Mathematics and Computation 195 (2008) 688–706
Substituting (6.76) and (6.77) into (6.76) yields 00 jLNx; ðvN v Þðxi;jþ1 Þj 6 Cðxiþ1;jþ1 xi1;jþ1 Þðjv000 ðxi;jþ1 Þj þ jv ðxi;jþ1 ÞjÞ:
701
ð6:78Þ
1
Now, since xiþ1;jþ1 xi1;jþ1 6 2N holds for all choice of piecewise uniform meshes, and the bounds on the smooth components lead to jLNx; ðvN v Þðxi;jþi Þj 6 C N 1 by an application of the preceding Lemma 6.2, we have jðvN v Þðxi;jþ1 Þj 6 jLNx; ðvN v Þðxi;jþ1 Þj 6 CN 1 . h Lemma 6.4. For all N P 4, at each time step (j þ 1) and mesh point xi 2 XNx the singular component of the error satisfies the estimates jðwN w Þðxi;jþ1 Þj 6 CN 1 ln N : ð6:79Þ For total discrete uniform convergence, we have the following parameter uniform convergence result. Theorem 6.5. The fitted mesh finite difference method with the standard upwind finite difference operator and the piecewise uniform fitted mesh XNx , condensing at the boundary points, is -uniform for the continuous problem Eq. (1.1) provided that d is chosen to satisfy condition (5.70). Moreover, the solution u of the continuous problem Eq. (1.1) and the solution uN of the discrete problem satisfy the following -uniform error estimates: sup kuN u k1;XNx 6 CðN 1 ðln N Þ2 þ M 1 Þ;
0<61
where C is a constant independent of . 7. Computational results and discussion To examine the performance of the method, a set of numerical experiments is carried out. Choice of the range of the parameter reflects our interests in the perturbed case. In order to tabulate the maximum error and the order of convergence of the method, we use the double mesh principal [15]. These are defined as EN ¼ max juN ðxi ; tj Þ u2N ðxi ; tj Þj; EN ¼ max EN ; 06i;j6N ;M
where superscript indicates the number of mesh points used in the spatial direction, and tj ¼ jDt and Dt is the time step followed by the tabulation of computed order of convergence pN , and the computed parameter uniform order of convergence pN on the same lines, defined as logðEN =E2N logðEN =E2N Þ Þ ; pN ¼ : pN ¼ log 2 log 2 Example 1. ut þ uux uxx ¼ ð1 uÞðu 0:5Þu, under the homogeneous boundary conditions uð0; tÞ ¼ 0; uð1; tÞ ¼ 0; 0 6 t 6 T , and initial condition u0 ðx; 0Þ ¼ sinðpxÞ 0 6 x 6 1. Example 2. ut þ uux ¼ uxx , under the homogeneous boundary conditions uð0; tÞ ¼ 0; uð1; tÞ ¼ 0; 0 6 t 6 T , and initial condition u0 ðx; 0Þ ¼ xð1 x2 Þ; 0 6 x 6 1. 8. Conclusion In this paper, the initial boundary value problem for generalized singularly perturbed time dependent Burger– Huxley equation is considered. There are two factors in the problem that lead to severe difficulties in approximating the solution of the problem. The first one is due to the presence of the perturbation parameter and the second one is due to the nonlinearity in the problem. Further, the current interest in the problem is due to its great practical value. The equation admits two particular cases each for a = 0 and b = 0, respectively. In the first case, the equation reduces to an initial boundary value problem of reaction diffusion type, and in the second case, the resultant is a convection diffusion type problem. Now in the limiting case both the cases will have different behaviors. For the problem of reaction diffusion type, the solution admits the layer at both ends of the outflow boundary region, while for the convection diffusion problem the solution exibits only one layer on either side of the interval. The first case requires a different analysis to deal with, so it is differently treated in [7].
702
A. Kaushik, M.D. Sharma / Applied Mathematics and Computation 195 (2008) 688–706
In order to overcome the difficulties, on a specially generated piecewise uniform mesh, a uniformly convergent numerical scheme using the standard finite difference upwind operator is proposed to analyze the layer behavior. An extensive amount of theoretical analysis is carried out to obtain the stability and the error estimates. It was found that the method proposed is unconditionally stable and the convergence obtained is
Table 1 The maximum absolute error for Example 1, when Dt ¼ 0:01 e
N 8
16
32
64
128
256
512
22 24 26 28 210 212 214 216 218 220
0.046959 0.064927 0.069603 0.071201 0.071695 0.071818 0.071849 0.071857 0.071859 0.071859
0.027434 0.039037 0.042693 0.043571 0.043786 0.043839 0.043852 0.043856 0.043856 0.043857
0.014867 0.021442 0.023719 0.024304 0.024446 0.024481 0.024490 0.024492 0.024493 0.024493
0.007745 0.011273 0.012526 0.012873 0.012957 0.012977 0.012983 0.012984 0.012984 0.012984
0.003955 0.005787 0.006450 0.006632 0.006678 0.006689 0.006692 0.006693 0.006693 0.006693
0.001999 0.002932 0.003274 0.003368 0.003391 0.003397 0.003399 0.003399 0.003399 0.003399
0.001005 0.001476 0.001650 0.001697 0.001709 0.001712 0.001713 0.001713 0.001713 0.001713
EN
0.071859
0.043857
0.024493
0.012984
0.006693
0.003399
0.001713
Table 2 The maximum absolute error for Example 1, when Dt ¼ 0:001 e
N 8
16
32
64
128
256
512
22 24 26 28 210 212 214 216 218 220
0.045165 0.061563 0.066019 0.067127 0.067403 0.067472 0.067489 0.067494 0.067495 0.067495
0.026306 0.036563 0.039974 0.040840 0.041057 0.041111 0.041124 0.041128 0.041129 0.041129
0.014235 0.020058 0.021987 0.022530 0.022675 0.022711 0.022720 0.022722 0.022723 0.022723
0.007411 0.010533 0.011611 0.011904 0.011982 0.012002 0.012006 0.012008 0.012008 0.012008
0.003783 0.005402 0.005969 0.006126 0.006166 0.006176 0.006179 0.006179 0.006179 0.006179
0.001911 0.002736 0.003028 0.003108 0.003129 0.003134 0.003135 0.003136 0.003136 0.003136
0.000961 0.001377 0.001525 0.001566 0.001576 0.001579 0.001580 0.001580 0.001580 0.001580
EN
0.067495
0.041129
0.022723
0.012008
0.006179
0.003136
0.001580
Table 3 The maximum absolute error for Example 1, when Dt ¼ 0:0001 e
N 8
16
32
64
128
256
512
22 24 26 28 210 212 214 216 218 220
0.044977 0.061213 0.065642 0.066747 0.067022 0.067091 0.067108 0.067112 0.067113 0.067114
0.026188 0.036307 0.039685 0.040547 0.040763 0.040817 0.040831 0.040834 0.040835 0.040835
0.014169 0.019914 0.021829 0.022338 0.022482 0.022518 0.022527 0.022529 0.022530 0.022530
0.007377 0.010457 0.011520 0.011803 0.011881 0.011900 0.011905 0.011906 0.011907 0.011907
0.003765 0.005362 0.005921 0.006074 0.006114 0.006125 0.006127 0.006128 0.006128 0.006128
0.001902 0.002716 0.003003 0.003082 0.003103 0.003108 0.003109 0.003109 0.003110 0.003110
0.000956 0.001367 0.001512 0.001553 0.001563 0.001566 0.001566 0.001566 0.001566 0.001566
EN
0.067114
0.040835
0.022530
0.011907
0.006128
0.003110
0.001566
A. Kaushik, M.D. Sharma / Applied Mathematics and Computation 195 (2008) 688–706
703
parameter uniform. Further, in Section 7, a set of numerical experiments is carried out in support of the predicted theory, which were found in good accordance with the predicted theory. The maximum absolute error is tabulated for both the examples for different time steps and can be seen from Tables 1–3 for example 1 and Tables 4–6 for example 2. It was observed from the tables that the order Table 4 The maximum absolute error for Example 2, when Dt ¼ 0:01 N
e
8
16
32
64
128
256
512
2 24 26 28 210 212 214 216 218 220
0.007691 0.012118 0.014024 0.014374 0.014449 0.014467 0.014471 0.014472 0.014472 0.014472
0.004348 0.006814 0.007701 0.007865 0.007919 0.007931 0.007935 0.007935 0.007936 0.007936
0.002311 0.003621 0.004086 0.004174 0.004194 0.004199 0.004200 0.004201 0.004201 0.004201
0.001195 0.001867 0.002099 0.002148 0.002160 0.002163 0.002163 0.002163 0.002163 0.002163
0.000607 0.000948 0.001065 0.001089 0.001095 0.001097 0.001097 0.001097 0.001097 0.001097
0.000306 0.000477 0.000536 0.000549 0.000552 0.000552 0.000553 0.000553 0.000553 0.000553
0.000154 0.000240 0.000269 0.000275 0.000277 0.000277 0.000277 0.000277 0.000277 0.000277
EN
0.014472
0.007936
0.004201
0.002163
0.001097
0.000553
0.000277
2
Table 5 The maximum absolute error for Example 2, when At Dt ¼ 0:001 e
N 8
16
32
64
128
256
512
22 24 26 28 210 212 214 216 218 220
0.007368 0.011327 0.012945 0.013244 0.013309 0.013325 0.013328 0.013329 0.013330 0.013330
0.004169 0.006351 0.007084 0.007205 0.007232 0.007239 0.007241 0.007241 0.007241 0.007241
0.002217 0.003369 0.003736 0.003810 0.003828 0.003832 0.003833 0.003833 0.003833 0.003833
0.001147 0.001735 0.001917 0.001953 0.001963 0.001965 0.001966 0.001966 0.001966 0.001966
0.000583 0.000880 0.000971 0.000990 0.000995 0.000996 0.000996 0.000996 0.000996 0.000996
0.000294 0.000443 0.000489 0.000498 0.000501 0.000501 0.000502 0.000502 0.000502 0.000502
0.000147 0.000222 0.000245 0.000250 0.000251 0.000252 0.000252 0.000252 0.000252 0.000252
EN
0.013330
0.007241
0.003833
0.001966
0.000996
0.000502
0.000252
Table 6 The maximum absolute error for Example 2, when At Dt ¼ 0:0001 e
N 8
16
32
64
128
256
512
22 24 26 28 210 212 214 216 218 220
0.007334 0.011247 0.012836 0.013130 0.013194 0.013209 0.013213 0.013214 0.013214 0.013214
0.004150 0.006304 0.007022 0.007141 0.007168 0.007175 0.007176 0.007177 0.007177 0.007177
0.002207 0.003343 0.003701 0.003774 0.003791 0.003795 0.003796 0.003796 0.003797 0.003797
0.001142 0.001722 0.001898 0.001935 0.001943 0.001946 0.001946 0.001946 0.001946 0.001946
0.000580 0.000873 0.000962 0.000980 0.000985 0.000986 0.000986 0.000987 0.000987 0.000987
0.000292 0.000440 0.000484 0.000493 0.000496 0.000496 0.000497 0.000497 0.000497 0.000497
0.000147 0.000221 0.000243 0.000248 0.000249 0.000249 0.000249 0.000249 0.000249 0.000249
EN
0.013214
0.007177
0.003797
0.001946
0.000987
0.000497
0.000249
704
A. Kaushik, M.D. Sharma / Applied Mathematics and Computation 195 (2008) 688–706
of convergence is very good. But the convergence of the time discretization is a little low, the reason lying in the fact that the choice of the temporal step size is governed by the accuracy and stability considerations. In the ideal case, accuracy and stability restrictions on the time step should be about the same. If the time step is much more restricted by stability than by accuracy, efficiency may be gained by switching to a possibly more
1
0.8
ε=2 ε=2
0.9
0.7
ε=2
ε=2
0.8
0.6 ε=2 Computed Solutions
Computed Solutions
0.7 ε=2
0.6 0.5 0.4
0.5 0.4
0.3
0.3 0.2 ε=2 0.2 0.1
0.1 0
0
0.5 T=0.1
0
1
ε=2 0
0.5 T=0.9
1
Fig. 1. Computed solutions of Example 1 for different values of , with N = 128 and M = 10.
0.6
1 T=0.1
T=0.3
0.5
0.8 0.7
0.4
Computed Solutions
Computed Solutions
T=0.1
0.9
0.3
T=0.3
0.2
T=0.6
0.6 T=0.9 0.5 0.4 0.3 0.2
0.1 0.1
T=0.6 0
T=0.9 0
0.5 =2
1
0
0
0.5 =2
1
Fig. 2. Computed solution of Example 1 for different values of T, with N = 128 and M = 10.
A. Kaushik, M.D. Sharma / Applied Mathematics and Computation 195 (2008) 688–706
705
expansive but more stable time discretization. Further, the effect of perturbation parameter and the final time on the layer behavior can be seen from Figs. 1 and 2, and Figs. 3 and 4, for example 1 and example 2, respectively. It can be seen from the figures that this scheme faithfully mimics the dynamics of the corresponding differential equation.
0.4
0.35
ε=2
ε=2 0.35
ε=2
ε=2
0.3
ε=2
0.3 0.25
Computed Solutions
Computed Solutions
0.25 ε=2
0.2 0.15
0.15
0.1
0.1
ε=2
0.05
0.05 0
0.2
0
0.5 T=0.1
ε=2
0
1
0
0.5 T=0.9
1
Fig. 3. Computed solutions of Example 2 for different values of , with N = 128 and M = 10.
0.35
0.4 T=0.1 0.35
0.3
0.3
T=0.3 T=0.6 T=0.9
T=0.3
Computed Solutions
Computed Solutions
0.25
T=0.1
0.2
0.15 T=0.6
0.25
0.2
0.15
0.1 0.1
T=0.9 0.05
0
0
0.05
0.5 ε=2
1
0
0
0.5 ε=2
1
Fig. 4. Computed solution of Example 2 for different values of T, with N = 128 and M = 10.
706
A. Kaushik, M.D. Sharma / Applied Mathematics and Computation 195 (2008) 688–706
References [1] X.Y. Wang, Z.S. Zhu, Y.K. Lu, Solitary wave solutions of the generalized Burger’s–Huxley equation, J. Phys. A: Math. Gen. 23 (1990) 271–274. [2] A.C. Scott, Neurophysics, John Wiley, New York, 1977. [3] J. Satsuma, Explicit solutions of nonlinear equations with density dependent diffusion, J. Phys. Soc. Jpn. 56 (1987) 1947–1950. [4] X.Y. Wang, Nerve propagation and wall in liquid crystals, Phy. Lett. 112A (1985) 402–406. [5] X.Y. Wang, Brochard–Lager wall in liquid crystals, Phys. Rev. A 34 (1986) 5179–5182. [6] X.Y. Wang, Z.S. Zhu, Y.K. Lu, Solitary Waves solutions of the generalized Burger–Huxley equations, J. Phys.: Math. Gen. 23 (1990) 271–274. [7] A. Kaushik, M.D. Sharma, Numerical analysis of a mathematical model describes propagation of electrical pulses in a neuron, Numer. Meth. Partial Differ. Equation, in press. [8] J.M. Burgers, A mathematical model illustrating the theory of turbulence, Adv. Appl. Mech. 1 (1948) 171–199. [9] M.J. Lighthill, Viscosity effects in sound waves of finite amplitude, in: G.K. Batchelor, R.M. Davies (Eds.), Surveys in Mechanics, Cambridge University Press, Cambridge, 1956, pp. 250–351. [10] H.G. Roos, M. Stynes, L. Tobiska, Numerical Methods for Singularly Perturbed Differential Equations. Convection–Diffusion and Flow Problems, Springer-Verlag, Berlin, 1996. [11] H.N.A. Ismail, K. Raslan, A.A.A. Rabboh, Adomian decomposition method for Burgers’–Huxley and Burger’s–Fisher equations, Appl. Math. Comput. 159 (2004) 291–301. [12] P.G. Estevez, Non-classical symmetries and the singular modified the Burger’s and Burgers’–Huxley equation, J. Phys. A 27 (1994) 2113–2127. [13] M. Javidi, A numerical solution of the generalized Burger’s–Huxley equation by pseudospectral method and Darvishi’s preconditioning, Appl. Math. Comput. 175 (2006) 1619–1628. [14] A.A. Samarskii, Theory of Difference Schemes, third ed., Nauka, Moscow, 1989. [15] 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.
Applied Mathematics and Computation 195 (2008) 688–706 www.elsevier.com/locate/amc
A uniformly convergent numerical method on non-uniform mesh for singularly perturbed unsteady Burger–Huxley equation q Aditya Kaushik a
a,*
, M.D. Sharma
b
INRIA, Mathe´matiques Applique´es de Bordeaux, University Bordeaux 1, 351 Cours de Libration F-33405, Cedex Talence, France b Department of Mathematics, Kurukshetra University, Kurukshetra 136 119, Haryana, India
Abstract In this paper, we present a numerical method which deals with the solvability of initial boundary value problem and related general properties of singularly perturbed Burger–Huxley equation. Singular perturbation problems are the problems in which a very small positive parameter is multiplied to the highest order derivative term. When this parameter approaches the limiting value of interest, i.e. ! 0, the problem exhibits boundary layers and most of the conventional methods fails to capture this effect. Thus, the quest for some new numerical techniques that may handle the difficulties occurring due to the presence of perturbation parameter and nonlinearity in the problem earns relevance. In this paper, one such parameter robust numerical scheme is constructed. The first step in this direction is the discretization of the time variable using the Euler implicit method with a constant time step. This produces a set of nonlinear stationary singularly perturbed semidiscrete problem class which is linearized further using the quasi-linearization process followed by spatial discretization using upwind finite difference operator on a piecewise uniform mesh. An extensive amount of analysis is carried out that uses a suitable decomposition of the error into smooth and singular component and a comparison principle combined with appropriate barrier functions. The error estimates are obtained, which ensures uniform convergence of the method. A set of numerical experiments is carried out in support of the predicted theory which validates the theoretical results computationally. Ó 2007 Elsevier Inc. All rights reserved. Keywords: Burger–Huxley equation; Piecewise uniform mesh; Uniform convergence; Singular perturbation; Quasi-linearization
1. Introduction In this paper, we construct a numerical method for solving time dependent singularly perturbed Burger– Huxley equation, L ½u ut þ auux uxx bð1 uÞðu cÞu ¼ 0; q
ð1:1aÞ
This work was completed with financial support by the UGC jointly with CSIR and HRDG Unit under Junior Research Fellowship vide No. F.17-39/98(SA-I). * Corresponding author. E-mail address: [email protected] (A. Kaushik). 0096-3003/$ - see front matter Ó 2007 Elsevier Inc. All rights reserved. doi:10.1016/j.amc.2007.05.067
A. Kaushik, M.D. Sharma / Applied Mathematics and Computation 195 (2008) 688–706
689
D fðx; tÞ 2 X ð0; T ð0; 1Þ ð0; T g; uðx; 0Þ ¼ UðxÞ; x 2 X ð0; 1Þ;
ð1:1bÞ ð1:1cÞ
uð0; tÞ ¼ uð1; tÞ ¼ 0;
ð1:1dÞ
t 2 ð0; T ;
where a; b P 0 and c 2 ð0; 1Þ and 0 < 1 are the parameters. The equation describes the interaction between convection, diffusion, and reaction. When a = 0 and = 1, Eq. (1.1) reduces to Huxley equation [1]. This equation describes nerve pulse propagation in nerve fibers and wall motion in liquid crystals [2–6]: ut uxx ðx; tÞ ¼ bð1 uÞðu cÞu: ð1:2Þ In the case when the epsilon is very small, i.e. ! 0, qualitative analysis of the model, which closely mimics the ionic processes at a real nerve membrane, is performed by the authors, by means of a singular perturbation theory in [7]. On the other hand, when b = 0, Eq. (1.1) reduces to the Burgers equation at high Reynolds number that establishes a balance between time evolution, nonlinearity and diffusion: ut þ auux ¼ uxx ðx; tÞ:
ð1:3Þ
This equation was first introduced by Burger [8] primarily to throw light on turbulence described by the interaction of two opposite effects of convection and diffusion. Burger equation arises in many physical problems including one dimensional turbulence, sound waves in viscous medium [9], shock waves in viscous medium, waves in a medium with fluid filled viscoelastic tubes and magnetohydrodynamic waves in a medium with finite electrical conductivity. This paper deals with the numerical study of generalized time dependent Burger–Huxley equation which is singularly perturbed from mathematical perspective. Parameter dependent differential equations such as (1.1) occur in diverse area of science and engineering. The concept of singular perturbation is not new, indeed, it has been a formidable tool in the solution of some important applied mathematical problems. In accordance with the informal principle, the behavior of solutions is governed primarily by the highest order terms, a solution u of the perturbed problem will often behave analytically quite differently from a solution of the original equation. When the perturbation parameter specifying the problem tends to zero, a breakdown occurs in the solution of the problem, and gives rise to boundary or the transition layers in the outflow boundary regions [10]. For higher values of a number of solution methodologies existing in the literature, but for sufficiently small , the existing solution methodology fails and a discrepancy occurs in the literature. In order to approximate the solution of such type of problems, numerical and asymptotic are two principal approaches. There are a wide variety of asymptotic expansion methods available for solving the problems of above type. But there can be difficulties in applying these asymptotic expansion methods, such as finding the appropriate asymptotic expansions in the inner and the outer regions, which are not routine exercises but require skill, insight and experimentation. Then one comes down ultimately to the numerical methods. It is a well established fact that nonlinear diffusion equations (1.1)–(1.3) with a small parameter plays an important role in nonlinear physics. Also it is of great practical interest to study the nonlinear phenomena. In order to solve Eq. (1.1), various numerical approaches were adopted by the researchers. Ismail et al. [11] studied the adomian decomposition method for Burger–Huxley and Burger–Fisher equation. Wang et al. [1] studied the solitary wave solutions of the generalized Burger–Huxley equation and Estevez [12] presented the non-classical symmetries and the singular modified solutions of the Burger and Burger–Huxley equation. More recently, Javidi [13] proposed a numerical solution of generalized Burger–Huxley equation by pseudospectral method and Darvishi’s preconditioning. However, it is well known for singular perturbation problem, that the numerical methods consisting of standard finite difference operators lead to severe restriction on the step size in order to attain the stability for small values of . To avoid these restrictions on the step size, a parameter uniform numerical scheme is constructed. The method can easily be extended to higher dimensions and for the nonlinear problems, but the only thing required in the process is the location and the exact width of the layer. 2. Bounds on the solution and time derivatives Since we shall use such a priori estimates later on in several different situations, we present in this section the a priori information needed in sequel. Moreover, the estimates theorems will be seen to be useful as a guide in the effective construction of the convergence and stability of the time discretization.
690
A. Kaushik, M.D. Sharma / Applied Mathematics and Computation 195 (2008) 688–706
Lemma 2.1. Let Uðx; tÞ; Wðx; tÞ 2 C 2;1 , i.e., twice continuously differential in space and once in time, and satisfy the relation jL½Uj < L½W;
ð2:4Þ
where L is the operator defined in Eq. (1.1). If jUj 6 W
ð2:5Þ
holds on the boundary, then jUj 6 W holds throughout the domain. Proof. Assume that ½U W attains the positive maximum at an interior point say (x, t) of the domain. Then according to the maximum principle ðL ½U WÞðx;tÞ 6 0, a contradiction to the assumption (2.4). Thus our assumption that ½U W attains positive maximum in the interior of the domain is wrong and hence ½U W does not have positive maximum within the domain, also from (2.5), along the boundary of the domain U W 6 0, it follows that U W 6 0 in D:
ð2:6Þ
Now suppose that U W attains a positive maximum in an interior point of the domain, then similarly it follows that U W 6 0
ð2:7Þ
in D:
Now from (2.7), (2.6) we have, jUj 6 W holding throughout the domain. Hence the result follows.
h
For obvious reasons the function W is called a barrier function for U. We impose the compatibility conditions u0 ð0; 0Þ ¼ u0 ð1; 0Þ ¼ 0, so that the data match at the corners (0, 0) and (1, 0) of the domain. As an application of this lemma we immediately obtain Theorem 2.2. There exists a number C independent of the perturbation parameter such that for all sufficiently small positive values of ; juðx; tÞ UðxÞj 6 Ct; t 2 ð0; T . Furthermore, juðx; tÞj 6 C;
ðx; tÞ 2 X ð0; T :
ð2:8Þ
Proof. We set gðx; tÞ ¼ uðx; tÞ UðxÞ, then g satisfies the differential equation L gðx; tÞ ¼ Uxx ðxÞ aUUx þ bð1 UÞðU cÞ ¼ H ðx; tÞ;
ð2:9Þ
where gðx; 0Þ ¼ 0; 0 < x < 1, gð0; tÞ ¼ uð0; tÞ Uð0Þ ¼ Uð0Þ;
gð1; tÞ ¼ uð1; tÞ Uð1Þ ¼ Uð1Þ;
0 6 t 6 1:
Now, clearly L ðCtÞ 6 C:
ð2:10Þ
Now, it can be immediately obtained from Lemma 2.1 and Eq. (2.10) that Ct is the barrier function for g for all small provided C is chosen so large that C P supðx;tÞ2D jH ðx; tÞj and so that Ct > jgj on the boundary. Thus jgðx; tÞj ¼ juðx; tÞ UðxÞj 6 Ct
ð2:11Þ
is uniform D and small . Using the properties of modulus, Eq. (2.11) is written as juðx; tÞj jUðxÞj 6 juðx; tÞ UðxÞj 6 Ct;
t 2 ð0; T :
Now since UðxÞ is sufficiently smooth and x and t also lie in the bounded interval, it is clear from the above relations that the solution u is bounded. h
Lemma 2.3. By keeping x fixed along the line fðx; tÞ : 0 6 t 6 T g, the bound ut is given as jut ðx; tÞj 6 C:
ð2:12Þ
A. Kaushik, M.D. Sharma / Applied Mathematics and Computation 195 (2008) 688–706
691
Proof. We assume that the solution uðx; tÞ is sufficiently smooth in the domain X ½0; T and by mean value theorem, there exits a t* in the interval ðt; t þ kÞ along the line fðx; tÞ : 0 6 t 6 T g such that uðx; t þ kÞ uðx; tÞ ; k 2juðx; tÞj jut ðx; t Þj 6 : k
ut ðx; t Þ ¼
ð2:13Þ
Using the Eq. (2.8), we get jut ðx; tÞj 6 C:
ð2:14Þ
Similarly we get the bounds of the utt ðx; tÞ and uttt ðx; tÞ along the line fðx; tÞ : 0 6 t 6 T g, i.e. we have i o uðx; tÞ ð2:15Þ oti 6 C for i ¼ 0; 1; 2; 3: 3. Time semidiscretization and quasi-linearization The first step in the numerical solution consists of discretizing the time variable with the Euler’s implicit method with constant time step Dt. This produces a set of stationary singularly perturbed problems of type Dt ujþ1 xx þ aDt½ujþ1 ½ujþ1 x þ ½ujþ1 bDtð1 ½ujþ1 Þð½ujþ1 cÞ½ujþ1 ¼ ½uj ; ð3:16aÞ ujþ1 ð0Þ ¼ 0
ujþ1 ð1Þ ¼ 0;
and
0 6 j 6 M 1:
ð3:16bÞ
Now for each j, (3.16) is a singularly perturbed nonlinear ordinary differential equation. To develop a numerical scheme for the boundary value problem Eq. (1.1), we first linearize the nonlinear problem (3.16) by using the quasi-linearization process. The nonlinear differential equation is linearized around a nominal solution of the nonlinear differential equation which satisfies the boundary conditions. Suppose uðkÞ ðxÞ is the nominal solution of the problem (3.16). Their application of the quasi-linearization process to the nonlinear problem (3.16) yields a sequence huðkÞ i of linear equations determined by the following recurrence relation: ðu00jþ1 Þ
ðkþ1Þ
þ ukþ1 ¼ f ðkÞ þ ððu0jþ1 Þ j
ðkþ1Þ
ðkÞ
ðkÞ
ðu0jþ1 Þ Þfu0
jþ1
þ ððujþ1 Þ
ðkþ1Þ
ðkÞ
ðujþ1 Þ ÞfuðkÞ ; jþ1
k ¼ 0; 1; . . .
Thus (3.16) followed by quasi-linearization leads to Lx; ujþ1 ¼ hðxÞ
ð3:17aÞ
with boundary condition ujþ1 ð0Þ ¼ 0
and
ujþ1 ð1Þ ¼ 0;
0 6 j 6 M 1:
ð3:17bÞ
where Lx; ¼
o2 o þ aðxÞ þ bðxÞ; ox ox2 h i
aðxÞ ¼ ak ¼ a ukjþ1 P b > 0; h h i i h i2 h i 1 P h > 0; bðxÞ ¼ bk ¼ þ a ukjþ1 b 2 ukjþ1 3 ukjþ1 c þ 2c ukjþ1 x Dt ½ukþ1 j þ a½ukjþ1 ½ukjþ1 x þ bð2½ukjþ1 3 ð1 þ cÞ½ukjþ1 2 Þ: hðxÞ ¼ Dt Remark 3.1. At first sight, the parameter Dt seems to behave like another perturbation parameter into the time semidiscretized problem, but under enough smoothness and compatibility requirements on the data of the continuous problem we have seen that this parameter does not have any severe effect as that of the perturbation parameter , to the multiscale character of the semidiscrete solutions.
692
A. Kaushik, M.D. Sharma / Applied Mathematics and Computation 195 (2008) 688–706 ðkþ1Þ
Now the boundary value problem (3.17) is linear in ujþ1 and now instead of solving the original nonlinear problem (3.16), we will solve the sequence of boundary value problems for singularly perturbed second order linear differential equations given by (3.17), for k ¼ 0; 1; 2; . . . and 0 6 j 6 M 1. Analytically, we requires as ðkÞ k tends to 1; ujþ1 ðxÞ to converge to the solution ujþ1 ðxÞ of the original nonlinear differential equation while numerically, we require that ðkþ1Þ ðkÞ ujþ1 ðxÞ ujþ1 ðxÞ < Tol:; 0 6 x 6 1; where Tol. is a small tolerance prescribed by us. If the above condition is satisfied, we terminate the iteration ðkþ1Þ and the profile ujþi ðxÞ is the numerical solution of the nonlinear boundary value problem (3.16). ðkÞ Now we prove the convergence of the sequence hujþ1 i. For the sake of convenience, consider eu00jþ1 ¼ f ðujþ1 Þ;
ð3:18aÞ
ujþ1 ð0Þ ¼ ujþ1 ðbÞ ¼ 0:
ð3:18bÞ ðkÞ
After quasi-linearization, we obtain a sequence hujþ1 i of linear equations determined by the recurrence eðu00jþ1 Þ
ðkþ1Þ
ðkþ1Þ
ðkÞ
ðkþ1Þ
ðkÞ
ðkÞ
ðxÞ ¼ f ðujþ1 Þ þ ðujþ1 ujþ1 Þf 0 ðujþ1 Þ;
ð3:19Þ
ðkþ1Þ
ujþ1 ð0Þ ¼ ujþ1 ðbÞ ¼ 0;
ð3:20Þ
where f 0 ðujþi Þ ¼ df ðujþi Þ=dujþi . Suppose ðk1Þ
ðkÞ
ð0Þ ujþ1 ðxÞ
ðk1Þ
is an initial approximation. Now from Eq. (3.19), we have ðk1Þ
eðu00jþ1 ÞðkÞ ðxÞ ¼ f ðujþ1 Þ þ ðujþ1 ujþ1 Þf 0 ðujþ1 Þ:
ð3:21Þ
Subtracting Eq. (3.21) from Eq. (3.19), we obtain ðkþ1Þ
ðkÞ
ðkÞ
00
ðk1Þ
ðkÞ
ðk1Þ
ðk1Þ
ðkþ1Þ
ðkÞ
ðkÞ
eðujþ1 ujþ1 Þ ¼ f ðujþ1 Þ f ðujþ1 Þ ðujþ1 ujþ1 Þf 0 ðujþ1 Þ þ ðujþ1 ujþ1 Þf 0 ðujþ1 Þ: ðkþ1Þ
ð3:22Þ
ðkÞ
Eq. (3.22) is a second order differential equation in ðujþ1 ujþ1 Þ. Converting this into an integral equation using Green’s function, we get Z b h i ðkþ1Þ ðkÞ ðkÞ ðk1Þ ðkÞ ðk1Þ ðk1Þ ðkþ1Þ ðkÞ ðkÞ eðujþ1 ujþ1 Þ ¼ Gðx; sÞ f ðujþ1 Þ f ðujþ1 Þ ðujþ1 ujþ1 Þf 0 ðujþ1 Þ þ ðujþ1 ujþ1 Þf 0 ðujþ1 Þ ds; 0
ð3:23Þ where Gðx; sÞ is the Green’s function and is defined by ( ðxbÞs ; 0 6 s 6 x 6 b; b Gðx; sÞ ¼ xðsbÞ ; 0 6 x 6 s 6 b; b
ð3:24Þ
and max jGðx; sÞj ¼ b=4:
ð3:25Þ
x;s
By the mean value theorem, we have ðkÞ
ðk1Þ
ðujþ1 ujþ1 Þ2 00 f ðnÞ; ð3:26Þ 2 ðk1Þ ðkÞ ðkÞ ðk1Þ where ujþ1 6 n 6 ujþ1 . Using Eq. (3.26) for f ðujþ1 Þ f ðujþ1 Þ in Eq. (3.23) followed by a simplification yields Z b ðkþ1Þ ðkÞ ðkÞ ðk1Þ 2 ðkþ1Þ ðkÞ ðkÞ eðujþ1 ujþ1 Þ ¼ Gðx; sÞ½ðujþ1 ujþ1 Þ f 00 ðnÞ=2 þ ðujþ1 ujþ1 Þf 0 ðujþ1 Þ ds: ð3:27Þ ðkÞ
ðk1Þ
ðkÞ
ðk1Þ
ðk1Þ
f ðujþ1 Þ f ðujþ1 Þ ¼ ðujþ1 ujþ1 Þf 0 ðujþ1 Þ þ
0
Let
max f 0 ðujþ1 Þ ¼ p;
kujþ1 k61
max f 00 ðujþ1 Þ ¼ q:
kujþ1 k61
ð3:28Þ
A. Kaushik, M.D. Sharma / Applied Mathematics and Computation 195 (2008) 688–706
693
Now taking the maximum of the moduli over the domain of consideration on both the sides of Eq. (3.27), we get ðkþ1Þ
ðkÞ
e max jujþ1 ujþ1 j 6
Z
06x6b
b 0
ðkÞ ðk1Þ 2 ðkþ1Þ ðkÞ ðkÞ max jGðx; sÞj maxðujþ1 ujþ1 Þ max jf 00 ðnÞj=2 þ max jujþ1 ujþ1 j max jf 0 ðujþ1 Þj ds: x;s
x
x
ð3:29Þ
Using Eqs. (3.25) and (3.28) in Eq. (3.29), we obtain b2 q ðkþ1Þ ðkÞ ðkÞ ðk1Þ 2 ðkþ1Þ ðkÞ maxðujþ1 ujþ1 Þ þ p max jujþ1 ujþ1 j : max jujþ1 ujþ1 j 6 x 06x6b 4e 2 x A rearrangement of terms yields b2 q ðkÞ ðk1Þ 2 maxðujþ1 ujþ1 Þ : x 8eð1 b2 p=4eÞ x This proves that the sequence hukjþ1 i of linear equations converges quadratically provided ðkþ1Þ
ðkÞ
max jujþ1 ujþ1 j 6
ð3:30Þ
b2 q < 1: 8eð1 b2 p=4eÞ
ð3:31Þ
The inequality (3.31) holds for sufficiently small b. If the interval ½0; b is too large, we can still keep ð1Þ ð0Þ ð0Þ jujþ1 ðxÞ ujþ1 ðxÞj sufficiently small by choosing a judicious initial approximation ujþ1 ðxÞ, which satisfies the ðkþ1Þ ðkÞ sufficient condition for the convergence, i.e., maxx jujþ1 ujþ1 j is small enough for some k. Remark 3.2. Thus we have shown that the sequence given by (3.17) of linear boundary value problems obtained after linearizing the nonlinear problem (3.16) converges to the original nonlinear boundary value problem. Therefore to obtain the approximate solution of the nonlinear problem (3.16), it is sufficient to approximate the linearized boundary value problem given by (3.17). 3.1. Convergence and stability i
Lemma 3.3. If j oto i uðx; tÞj 6 C, for ðx; tÞ 2 X ð0; T and 0 6 i 6 2. Then the local truncation error satisfies the estimate 2
kejþ1 k1 6 CðDtÞ : Proof. Now, since the solution of (3.17) is smooth enough, it holds Z tjþ1 ou o2 u uðtj Þ ¼ uðtjþ1 Þ Dt ðtjþ1 Þ þ ðtj nÞ 2 ðnÞ dn ot ot tj ¼ uðtjþ1 Þ Dtðuxx aðxÞux bðxÞuðx; tÞ þ f ðx; tÞÞðtjþ1 Þ þ OðDt2 Þ:
ð3:32Þ
Subtracting (3.16) from (3.32), and note that the local truncation error ejþi uðtjþ1 Þ ^ujþ1 at ðj þ 1Þth time step is the solution of boundary value problem of type Lx; ejþ1 ¼ OðDt2 Þ;
ð3:33aÞ
ejþ1 ð0Þ ¼ ejþ1 ð1Þ ¼ 0;
ð3:33bÞ
where ^ ujþ1 is the solution of the boundary value problem (3.17). Now, the operator Lx; clearly satisfies the maximum principle: consequently kL1 x; k1 6 C which ensures the stability of the semidiscretized scheme, where C is the positive constant independent of Dt. Thus we have kejþ1 k1 6 CðDtÞ2 :
ð3:34Þ
694
A. Kaushik, M.D. Sharma / Applied Mathematics and Computation 195 (2008) 688–706
Now, using the local error estimates up to the ðj þ 1Þth time step, followed by classical combination of stability and consistency property of the scheme, leads us to j X ek 6 CDt; ð3:35Þ kEj k1 ¼ k¼1 1
since jDt 6 T , Here C is a positive constant independent of and Dt. Thus we have proved the following: Theorem 3.4. Under the hypothesis of the preceding lemma, the global truncation error satisfies the estimate kEj k1 6 CDt:
ð3:36Þ
4. A priori estimates for semidiscretized problem In this section of the paper, we meet the important and recurring concept of an a priori estimate. Also we discuss qualitative properties of the solution of the semidiscretized problem. We now give to the maximum principle a mathematical basis, that can be a useful tool for deriving a priori bounds on the solutions of the differential equations and their derivatives. Rewriting the semidiscretized equation k k kþ1 kþ1 k Lx; ½ukþ1 jþ1 xx þ a ðxÞ½ujþ1 x þ b ðxÞ½ujþ1 ¼ f ðxÞ;
ujþ1 ð0Þ ¼ 0
and
ujþ1 ð1Þ ¼ 0;
06j6M 1
and the differential operator satisfies the following continuous maximum principal: Lemma 4.1 (Maximum Principle). If a function U, twice continuously differentiable in space and once in time, i.e. U 2 C 2;1 ðXÞ be such that Uðx; tÞ P 0, for all ðx; tÞ 2 oD, then Lx; Uðx; tÞ P 0 for all ðx; tÞ 2 D implies that Uðx; tÞ P 0 for all ðx; tÞ 2 D. Proof. Let ðx ; t Þ be such that Uðx ; t Þ ¼ minx2D Uðx; tÞ, and suppose that Uðx ; t Þ < 0. It is then clear from the hypothesis that ðx ; t Þ 62 oD, from this it follows ðx ; t Þ 2 D, hence Ux ðx ; t Þ ¼ 0, and Uxx ðx ; t Þ P 0, further Lx; Uðx ; t Þ ¼ Uxx ðx ; t Þ þ ak ðxÞUx ðx ; t Þ þ bk ðxÞUðx ; t Þ 6 0;
ð4:37Þ
which is a contradiction to our assumption, hence our assumption Uðx ; t Þ < 0 is wrong. Hence, Uðx; tÞ P 0 for all x 2 D. Hence the required result holds. h An immediate and important consequence of the maximum principle is the uniqueness and continuous dependence of the solution on the boundary values. For let hvi;jþ1 i and hwi;jþ1 i be the two sets of solutions of the semidiscretized difference equations satisfying the boundary conditions. Suppose zi;jþ1 ¼ hvi;jþ1 i hwi;jþ1 i, then since the semidiscrete operator is linear, we have zi;jþ1 ¼ hvi;jþ1 i hwi;jþ1 i satisfies Lx; zi;jþ1 ¼ Lx; vi;jþ1 Lx; wi;jþ1 ¼ 0 and obviously, z0;jþ1 ¼ 0; zN ;jþ1 ¼ 0. Thus the mesh function zi;jþ1 satisfies the hypothesis of discrete maximum principle that leads to zi;jþ1 P 0
for all i ¼ 1; 2; . . . ; N :
ð4:38Þ
On the other hand, if we set zi;jþ1 ¼ ðhvi;jþ1 i hwi;jþ1 iÞ, then a similar treatment to hzi;jþ1 i with an application of maximum principle leads us to zi;jþ1 P 0:
ð4:39Þ
Combining Eqs. (4.38) and (4.39), we have zi;jþ1 ¼ 0 and uniqueness follows immediately. Further since the problem under consideration is linear, the existence follows from uniqueness. The maximum principle can also be used to derive estimates on the derivatives of the solution provided some additional conditions are placed on the equation. In the next two lemmas we establish a priori bound on the exact solution and its derivatives, of time semidiscretized problem class (3.17). It is convenient to use the terminology suggested by the comparison principle combined with appropriate barrier functions.
A. Kaushik, M.D. Sharma / Applied Mathematics and Computation 195 (2008) 688–706
695
Lemma 4.2. Let u be the solution of any problem from semidiscretized problem class (3.17), then ku k1 6 C:
ð4:40Þ
Proof. Consider the barrier functions, f ¼ 1h khk1 u ðx; tÞ. Note that the functions f are non-negative at x ¼ 0; 1 for all t 2 ð0; T , and that for all ðx; tÞ Lx; f ðx; tÞ ¼ bðxÞh1 khk1 hðxÞ P 0;
ð4:41Þ
since bðxÞ P h > 0 and khk1 P hðxÞ. Thus maximum principle gives f ðx; tÞ P 0 and so, 1 ku k1 6 khk1 ¼ C h
for all ðx; tÞ 2 D:
Lemma 4.3. The derivatives uðkÞ of the solution u of semidiscretized problem satisfy the following bounds: pffiffi pffiffi kuk kX 6 Cð Þk ð1 þ ð Þk Þ maxfkuk; khkg; k ¼ 1; 2; pffiffi 3 pffiffi 3 ku3 kX 6 Cð Þ ð1 þ ð Þ Þ maxfkuk; khk; kh0 kg:
ð4:42Þ ð4:43Þ
Proof. Given any x 2 ð0; 1Þ, we can construct the neighborhood N x ¼ ðp; p þ rÞ such that x 2 N x and N x ð0; 1Þ. The mean value theorem implies that there exists x 2 N x such that u0 ðx Þ ¼
uðp þ rÞ uðpÞ r
and so ju0 ðx Þj 6
2kuk : r
Now we have 0
0
u ðxÞ ¼ u ðx Þ þ
ð4:44Þ
Z
x
u00 ðsÞ ds x
and therefore from the time semidiscretized differential Eq. (3.17), we have Z x Z x u0 ðxÞ ¼ u0 ðx Þ 1 hðsÞ ds þ 1 ½aðsÞu0 ðsÞ þ bðsÞuðsÞ ds: x
Now integration by parts leads us to Z Z x aðsÞu0 ðsÞ ds ¼ aðsÞujx a0 ðsÞuðsÞ : x
ð4:45Þ
x
ð4:46Þ
x
Now using the fact that the maximum norm of a function is always greater than the value of the over the domain of consideration, followed by calculation, yields Z x 0 0 ð4:47Þ aðsÞu ðsÞ ds 6 ð2kak þ ka kÞkuk: x
Substituting (4.48) in (4.45), we have Z x Z x Z x x u0 ðxÞ ¼ u0 ðx Þ 1 hðsÞ ds þ 1 bðsÞuðsÞ ds þ 1 aðsÞujx a0 ðsÞuðsÞ ds : x
x
x
Now since jx x j 6 r, and from (4.44) followed by (4.47) on simplification leads to 1 r 1 r 1 r 1 0 0 ju ðxÞj 6 C þ þ kuk þ kf k 6 ju ðxÞj 6 C þ þ maxfkuk; kf kg: r r
ð4:48Þ
696
A. Kaushik, M.D. Sharma / Applied Mathematics and Computation 195 (2008) 688–706
pffiffi If we choose r ¼ , then the right-hand side of the above expression is minimized with respect to r, and we obtain the required result for k = 1. Using the differential Eq. (3.17) for u, we can obtain the required bounds for k = 2 and on differentiating again (3.17) the result for k = 3 follows. h In order to have parameter uniform convergence results one needs to prove the stronger pointwise error bounds. To find these sharper pointwise bounds the solution of the semidiscretized problem has to be decomposed into smooth and singular components. The solution u of the semidiscrete problem is decomposed into a smooth component v and a singular component w as follows: u ¼ v þ w ; where v can be written in the form v ¼ v0 þ v1 þ 2 v2 ; and v0 ; v1 and v2 are defined, respectively, to be the solutions of the problems aðxÞv00 ðxÞ þ bðxÞv0 ðxÞ ¼ hðxÞ;
x 2 X;
aðxÞv01 ðxÞ
v1 ð1Þ ¼ 0;
v002 ðxÞ
þ bðxÞv1 ðxÞ ¼
þ
aðxÞv02 ðxÞ
v000 ;
þ bðxÞv2 ðxÞ ¼
v0 ð1Þ ¼ u ð1Þ ¼ 0;
v001 ;
ð4:49aÞ ð4:49bÞ
v2 ð0Þ ¼ v2 ð1Þ ¼ 0:
ð4:49cÞ
Thus the smooth component is the solution of Lx; v ¼ hðxÞ;
v ð0Þ ¼ v0 ð0Þ þ v1 ð0Þ;
v ð1Þ ¼ u ð1Þ ¼ 0
ð4:50Þ
and consequently singular component satisfies Lx; w ¼ 0;
w ð0Þ ¼ u ð0Þ v ð0Þ;
w ð1Þ ¼ 0:
ð4:51Þ
Lemma 4.4. For 0 < k 6 3; v0 satisfies the following estimates: ðkÞ
kv0 k 6 C: Proof. The problem (4.49a) can be written as v00 ðtÞ þ pðtÞv0 ðtÞ ¼ f ðtÞpðtÞ;
v0 ð1Þ ¼ 0;
ð4:52Þ
1
where pðtÞ ¼ bðxÞaðxÞ . The problem (4.52) isR a first order linear differential equation in v0, in order to obtain the solution we multiply Eq. (4.52) by exp pðtÞ dt , which followed by integration from x to 1, for some x 2 ð0; 1Þ, yields Z 1 Z 1 Z v0 ðtÞ exp pðtÞ dt ¼ f ðtÞ exp pðtÞ dt pðtÞ dt x
x
which on simplification reduces to Z Z 1 1 R v0 ðxÞ ¼ f ðtÞ exp pðtÞ dt pðtÞ dt: exp pðtÞ dt t¼x x
ð4:53Þ
Now since aðtÞ is sufficiently smooth and we have aðtÞ > h > 0 for all t 2 ½0; 1 also f ðtÞ is bounded for all t 2 ½0; 1. Therefore this boundedness of v0 follows immediately. Again from Eq. (4.49a), we have v00 ðxÞ ¼ f ðxÞpðxÞ pðxÞv0 ðxÞ and the boundedness of v00 ðxÞ follows directly from the boundedness of v0. Bounds on v000 ðxÞ and v000 0 ðxÞ can be easily obtained by successively differentiating Eq. (4.49a). Thus for 0 < k 6 3; v0 satðkÞ isfies the estimates kv0 k 6 C. h Lemma 4.5. For 0 < k 6 3; v1 satisfies the following estimates: ðkÞ
kv1 k 6 C:
A. Kaushik, M.D. Sharma / Applied Mathematics and Computation 195 (2008) 688–706
697
Proof. Since v1 is the solution of the first order linear differential Eq. (4.49b) and since v0 is bounded, using the classical argument used in Lemma 4.4 it is straightforward to prove that for 0 < k 6 3; v1 satisfies the required estimates. h Theorem 4.6. For 0 < k 6 3; v2 satisfies the following estimates: hð1 xÞ ðkÞ k kv2 k 6 C 1 þ exp for all x 2 X:
ð4:54Þ
Now we need to obtain the bound on the singular component of the solution. For if, introducing the two barrier functions hð1 xÞ w :¼ C exp w ;
then we have w ð0Þ ¼ C exp hð1Þ 0 P 0, and we choose C so that it holds w ð1Þ ¼ C w P 0. Further h hð1 xÞ hð1 xÞ ðnÞ ðnÞ L W ðxÞ ¼ h þ a ðxÞ exp þ b ðxÞ exp > 0; ð4:55Þ since aðnÞ ðxÞ P h > 0 and bn ðxÞ P b > 0. Therefore from Lemma 4.1 we obtain w P 0 for all x 2 X which gives us hð1 xÞ kw ðxÞk 6 C exp for all x 2 X; ð4:56Þ Z x Z x ðnÞ ðaðw ÞxÞðsÞ ds ¼ aðnÞ w ðxÞ ; 6 C exp hð1 xÞ : a ðw Þ ð4:57Þ ðsÞ ds x x x 0 0 By mean value theorem in the interval ð0; Þ, there exist z 2 ð0; Þ such that w ðÞ w ð0Þ 6 C1 : jw ðzÞj ¼ Integrating Eq. (4.51) w.r.t. x from 0 to x, we have Z x Z x ðnÞ ðw Þx ðxÞ þ ðw Þx ð0Þ þ a ðw Þx ðsÞ ds þ bðnÞ w ðsÞ ds ¼ 0: 0
ð4:59Þ
0
The above equation gives Z z Z z ðnÞ ðnÞ jðw Þx ð0Þj ¼ ðw Þx ðzÞ a ðw Þx ðsÞ ds b w ðsÞ ds: 0
ð4:58Þ
ð4:60Þ
0
Making use of Eqs. (4.57) and (4.58), we have hð1 Þ jðw Þx ð0Þj 6 C exp 6 C:
ð4:61Þ
Now Eq. (4.59) gives the following estimate with the help of the estimates given by (4.57) and (4.61) jðw ðxÞÞx j 6 C1 ;
x 2 X:
ð4:62Þ
Using the above estimate in Eq. (4.51), we get the estimate for ðw ðxÞÞxx , i.e. jðw Þxx j 6 C2
ð4:63Þ
Similarly by differentiating Eq. (4.51) w.r.t. x and using the estimates given by Eqs. (4.62) and (4.63), we get the estimate for ðw Þxxx ðkÞ
jðw Þ j 6 Ck
for k ¼ 1; 2; 3:
Now, the above results proven by us are summarized in the following theorem:
ð4:64Þ
698
A. Kaushik, M.D. Sharma / Applied Mathematics and Computation 195 (2008) 688–706
Theorem 4.7. The solution u of the semidiscretized problem admits the decomposition u :¼ v þ w ; where for all k; 0 6 k 6 3 and for all x 2 X, the smooth component of the solution v satisfies ðkÞ
jðv Þ ðxÞj 6 Cð1 þ 2k Þ
ð4:65Þ
and the singular component of the error w satisfies ðkÞ
jðw Þ j 6 Ck
ð4:66Þ
for some constant C independent of the perturbation parameter. 5. The spatial discretization 5.1. Classical difference scheme In order to have a better insight into the problem, we first consider the classical finite difference method that constitutes the base of the construction of the parameter uniform scheme later on. So let us first introduce a non-uniform mesh of the nodal points xi in the spatial direction. Let K 1 be the number of the intervals. Define hi ¼ xiþ1 xi ¼ max hi : i
For the continuous problem we use the following difference scheme [14]: Lj ¼ f ðxÞ:
ð5:67Þ
Here, Lj ¼ Dþ D ui;jþ1 þ aðxi ÞDþ ui;jþ1 þ bðxi Þui;jþ1 ; þ
þ
ð5:68aÞ
D D ui;jþ1 ¼ 2ðD ui;jþ1 D ui;jþ1 Þ=ðhi þ hiþ1 Þ;
ð5:68bÞ
D ui;jþ1 ¼ ðui;jþ1 ui1;jþ1 Þ=hi ; Dþ ui;jþ1 ¼ ðuiþ1;jþ1 ui;jþ1 Þ=hi þ 1:
ð5:68cÞ ð5:68dÞ
The difference scheme so defined is monotone. By means of the maximum principle and taking into account the a priori estimates on the derivative and the solution of the semidiscretized scheme, we find that the solution of the difference scheme converges for a fixed value of kuðx; tÞ jðx; tÞk 6 Cð2 N 1 þ M 1 Þ:
ð5:69Þ
This error bound for the classical scheme is clearly not -uniform. The proof of (5.69) follows the lines of the classical convergence proof for the difference schemes [14]. 5.2. The -uniformly convergent scheme In this paper, we are considering the upwind discretization on piecewise uniform fitted meshes. Since in practice we often want to construct a numerical method whose system matrix is a M-matrix. In this section, we discretize the boundary value problem (3.17) using the fitted mesh finite difference method composed of a standard upwind finite difference operator on a piecewise uniform mesh condensing at the boundary points. N The fitted piecewise uniform mesh XNx ¼ fxi gi¼0 on the interval [0, 1] is constructed by partitioning the interval into two subintervals, ð0; 1 dÞ; ð1 d; 1Þ. The resulting piecewise uniform mesh depends on one parameter which is called the transition parameter d and is chosen such that d min½0:5; C ln N ;
ð5:70Þ k
where C is the constant whose value depends upon the method applied. Also it is assumed that N ¼ 2 , where k P 2 is an integer, which guarantees that there is at least one point in the boundary regions. On each of these subintervals a uniform mesh is then placed.
A. Kaushik, M.D. Sharma / Applied Mathematics and Computation 195 (2008) 688–706
699
Define xi ¼
ihi
for 0 6 i 6 N =2;
ð1 dÞ þ ði N =2Þhi
for N =2 þ 1 6 i 6 N ;
where N is the number of mesh points. 2ð1 dÞ=N for 0 6 i 6 N =2; hi ¼ 2d=N for N =2 þ 1 6 i 6 N : The fitted mesh finite difference method for the time semidiscretized problem (4) at ðj þ 1Þth time step on the piecewise uniform mesh XNx is defined by LNx; ui;jþ1 ¼ hðxi Þ;
ð5:71Þ
u0;jþ1 ¼ uN ;jþ1 ¼ 0;
ð5:72Þ
where the discrete operator is defined as LNx; ui;jþ1 ¼ Dþ D ui;jþ1 þ aðxi ÞDþ ui;jþ1 þ bðxi Þui;jþ1 :
ð5:73Þ
6. Convergence and stability analysis 6.1. Stability analysis A discrete maximum principle can also be established directly, without appealing to the properties of the elements in the system matrix. Analogously to the continuous case, this is usually proved by the method of contradiction as we now show. Theorem 6.1 (Discrete Maximum Principle). Assume that mesh function Wi;jþ1 satisfies W0;jþ1 P 0 and WN ;jþ1 P 0. Then ðLNx; ÞWi;jþ1 P 0 for all xi 2 XNx implies that Wi;jþ1 P 0 for all xi 2 XNx . Proof. Assume that there exists a positive integer k such that Wk;jþ1 ¼ min06i6N Wi;jþ1 and assume Wk;jþ1 < 0. Then we have W0;jþ1 P 0 and WN ;jþ1 P 0, therefore it follows from the hypothesis that k 62 f0; N g. Also Wk;jþ1 Wk1;jþ1 6 0 and Wkþ1;jþ1 Wk;jþ1 P 0 and note that ðLNx; ÞWk;jþ1 < 0 since Wk;jþ1 < 0 by assumption and bðxi Þ > h > 0, which is a contradiction to our assumption. Therefore the assumption Wk;jþ1 < 0 is wrong, hence Wk;jþ1 P 0. We have chosen k to be fixed and arbitrary, so Wi;jþ1 P 0 for all xi 2 XNx . h As an immediate application of the maximum and minimum principle, we have the following uniform stability estimates: Lemma 6.2. Let Z i;jþ1 be any mesh function such that. Z 0;jþ1 ¼ Z N ;jþ1 ¼ 0, then for all i; 0 6 i 6 N jZ i;jþ1 j 6 h1 max jLNx; Z k;jþ1 j: 16k6N 1
Proof. To prove the estimate we construct two mesh functions W i;jþ1 ¼ Mxi;jþ1 Z i;jþ1 ; where M ¼ h1 max16k6N 1 jðLNx ÞZ k;jþ1 j. Then we clearly have 1 W max jLNx; Z k;jþ1 jxi;jþ1 Z 0;jþ1 ¼ 0; 0;jþ1 ¼ h
ð6:74aÞ
16k6N 1
1 max jLNx; Z k;jþ1 jxN ;jþ1 Z N ;jþ1 ¼ h1 max jLNx; Z k;jþ1 jxN ;jþ1 P 0 W N ;jþ1 ¼ h 16k6N 1
16k6N 1
ð6:74bÞ
700
A. Kaushik, M.D. Sharma / Applied Mathematics and Computation 195 (2008) 688–706
and for 1 6 i 6 N 1, ðLNx; ÞW i;jþ1 ¼ ðaðxi Þ þ bðxi Þxi;jþ1 Þ
h1
max jðLNx; ÞZ k;jþ1 j LNx; Z i;jþ1 P 0;
16k6N 1
Thus an application of the discrete maximum principle leads to the required estimates.
h
6.2. Convergence analysis The discrete problem is decomposed in an analogous manner as the second decomposition of the solution u of the continuous problem, in order to estimate the error in the smooth and singular component separately. Thus uN ðxi;jþ1 Þ ¼ vN ðxi;jþ1 Þ þ wN ðxi;jþ1 Þ; where vN is the solution of the inhomogeneous problem LNx; vN ðxi;jþ1 Þ ¼ f ðxi;jþ1 Þ; vN ðx0;jþ1 Þ ¼ 0 and
vN ðx1;jþ1 Þ ¼ 0
and wN ðxi;jþ1 Þ is the solution of the homogeneous problem LNx; wN ðxi;jþ1 Þ ¼ 0; wN ðx0;jþ1 Þ ¼ 0
and
wN ðx1;jþ1 Þ ¼ 0:
Then the error at each time step (j + 1) can be written in the form ðuN u Þðxi;jþ1 Þ ¼ ðvN v Þðxi;jþ1 Þ þ ðwN w Þðxi;jþ1 Þ: Now in the next two results, we estimate the error in the smooth and singular component, respectively. Lemma 6.3. At each time step (j + 1) and mesh point xi 2 XNx , the smooth component of error satisfies jðvN v Þðxi;jþ1 Þj 6 CN 1 :
ð6:75Þ
Proof. This is obtained using the following standard stability and consistency argument, we consider the local truncation error 2 o d þ LNx; ðvN v Þðxi;jþ1 Þ ¼ D D D ðx Þ þ aðx Þ v v ðxi;jþ1 Þ: i;jþ1 i dx ox2 At each time step, to reduce the order of differentiation in the integral, integration by parts leads us to Z xi;jþ1 o 1 D ðxi1;jþ1 sÞv00 ðsÞ ds: v ðxi;jþ1 Þ ¼ ox xi;jþ1 xi1;jþ1 xi1;jþ1 Taking the modulus on both the sides, we have Z xi;jþ1 00 o 1 6 jv ðxi;jþ1 Þj D ðx Þ ðs xi1;jþ1 Þ ds ¼ jv00 ðxi;jþ1 Þjðxi;jþ1 xi1;jþ1 Þ: v i;jþ1 ox x 2 x i;jþ1
i1;jþ1
ð6:76Þ
xi1;jþ1
Also, integrating by parts twice we have 2 3 R xi;jþ1 2 000 1 2 ðx sÞ v ðsÞ ds iþ1;jþ1 x o 1 x x 4 iþ1;jþ1 i;jþ1 Riþ1;jþ1 5: Dþ D 2 v ðxi;jþ1 Þ ¼ xi1;jþ1 1 xiþ1;jþ1 xi1;jþ1 ox ðs xi1;jþ1 Þ2 v000 ðsÞ ds xi;jþ1 xi1;jþ1
xi;jþ1
From this it follows that 3 R xiþ1;jþ1 000 2 2 2 1 ðxiþ1;jþ1 sÞ ds v ðxi;jþ1 Þ o xiþ1;jþ1 xi;jþ1 xi;jþ1 þ 4 5 R xi;jþ1 ox2 D D v ðxi;jþ1 Þ 6 x 2 1 iþ1;jþ1 xi1;jþ1 þ ðs xi1;jþ1 Þ ds xi;jþ1 xi1;jþ1
1 6 ðxiþ1;jþ1 xi1;jþ1 Þv000 ðxi;jþ1 Þ : 3
xi1;jþ1
ð6:77Þ
A. Kaushik, M.D. Sharma / Applied Mathematics and Computation 195 (2008) 688–706
Substituting (6.76) and (6.77) into (6.76) yields 00 jLNx; ðvN v Þðxi;jþ1 Þj 6 Cðxiþ1;jþ1 xi1;jþ1 Þðjv000 ðxi;jþ1 Þj þ jv ðxi;jþ1 ÞjÞ:
701
ð6:78Þ
1
Now, since xiþ1;jþ1 xi1;jþ1 6 2N holds for all choice of piecewise uniform meshes, and the bounds on the smooth components lead to jLNx; ðvN v Þðxi;jþi Þj 6 C N 1 by an application of the preceding Lemma 6.2, we have jðvN v Þðxi;jþ1 Þj 6 jLNx; ðvN v Þðxi;jþ1 Þj 6 CN 1 . h Lemma 6.4. For all N P 4, at each time step (j þ 1) and mesh point xi 2 XNx the singular component of the error satisfies the estimates jðwN w Þðxi;jþ1 Þj 6 CN 1 ln N : ð6:79Þ For total discrete uniform convergence, we have the following parameter uniform convergence result. Theorem 6.5. The fitted mesh finite difference method with the standard upwind finite difference operator and the piecewise uniform fitted mesh XNx , condensing at the boundary points, is -uniform for the continuous problem Eq. (1.1) provided that d is chosen to satisfy condition (5.70). Moreover, the solution u of the continuous problem Eq. (1.1) and the solution uN of the discrete problem satisfy the following -uniform error estimates: sup kuN u k1;XNx 6 CðN 1 ðln N Þ2 þ M 1 Þ;
0<61
where C is a constant independent of . 7. Computational results and discussion To examine the performance of the method, a set of numerical experiments is carried out. Choice of the range of the parameter reflects our interests in the perturbed case. In order to tabulate the maximum error and the order of convergence of the method, we use the double mesh principal [15]. These are defined as EN ¼ max juN ðxi ; tj Þ u2N ðxi ; tj Þj; EN ¼ max EN ; 06i;j6N ;M
where superscript indicates the number of mesh points used in the spatial direction, and tj ¼ jDt and Dt is the time step followed by the tabulation of computed order of convergence pN , and the computed parameter uniform order of convergence pN on the same lines, defined as logðEN =E2N logðEN =E2N Þ Þ ; pN ¼ : pN ¼ log 2 log 2 Example 1. ut þ uux uxx ¼ ð1 uÞðu 0:5Þu, under the homogeneous boundary conditions uð0; tÞ ¼ 0; uð1; tÞ ¼ 0; 0 6 t 6 T , and initial condition u0 ðx; 0Þ ¼ sinðpxÞ 0 6 x 6 1. Example 2. ut þ uux ¼ uxx , under the homogeneous boundary conditions uð0; tÞ ¼ 0; uð1; tÞ ¼ 0; 0 6 t 6 T , and initial condition u0 ðx; 0Þ ¼ xð1 x2 Þ; 0 6 x 6 1. 8. Conclusion In this paper, the initial boundary value problem for generalized singularly perturbed time dependent Burger– Huxley equation is considered. There are two factors in the problem that lead to severe difficulties in approximating the solution of the problem. The first one is due to the presence of the perturbation parameter and the second one is due to the nonlinearity in the problem. Further, the current interest in the problem is due to its great practical value. The equation admits two particular cases each for a = 0 and b = 0, respectively. In the first case, the equation reduces to an initial boundary value problem of reaction diffusion type, and in the second case, the resultant is a convection diffusion type problem. Now in the limiting case both the cases will have different behaviors. For the problem of reaction diffusion type, the solution admits the layer at both ends of the outflow boundary region, while for the convection diffusion problem the solution exibits only one layer on either side of the interval. The first case requires a different analysis to deal with, so it is differently treated in [7].
702
A. Kaushik, M.D. Sharma / Applied Mathematics and Computation 195 (2008) 688–706
In order to overcome the difficulties, on a specially generated piecewise uniform mesh, a uniformly convergent numerical scheme using the standard finite difference upwind operator is proposed to analyze the layer behavior. An extensive amount of theoretical analysis is carried out to obtain the stability and the error estimates. It was found that the method proposed is unconditionally stable and the convergence obtained is
Table 1 The maximum absolute error for Example 1, when Dt ¼ 0:01 e
N 8
16
32
64
128
256
512
22 24 26 28 210 212 214 216 218 220
0.046959 0.064927 0.069603 0.071201 0.071695 0.071818 0.071849 0.071857 0.071859 0.071859
0.027434 0.039037 0.042693 0.043571 0.043786 0.043839 0.043852 0.043856 0.043856 0.043857
0.014867 0.021442 0.023719 0.024304 0.024446 0.024481 0.024490 0.024492 0.024493 0.024493
0.007745 0.011273 0.012526 0.012873 0.012957 0.012977 0.012983 0.012984 0.012984 0.012984
0.003955 0.005787 0.006450 0.006632 0.006678 0.006689 0.006692 0.006693 0.006693 0.006693
0.001999 0.002932 0.003274 0.003368 0.003391 0.003397 0.003399 0.003399 0.003399 0.003399
0.001005 0.001476 0.001650 0.001697 0.001709 0.001712 0.001713 0.001713 0.001713 0.001713
EN
0.071859
0.043857
0.024493
0.012984
0.006693
0.003399
0.001713
Table 2 The maximum absolute error for Example 1, when Dt ¼ 0:001 e
N 8
16
32
64
128
256
512
22 24 26 28 210 212 214 216 218 220
0.045165 0.061563 0.066019 0.067127 0.067403 0.067472 0.067489 0.067494 0.067495 0.067495
0.026306 0.036563 0.039974 0.040840 0.041057 0.041111 0.041124 0.041128 0.041129 0.041129
0.014235 0.020058 0.021987 0.022530 0.022675 0.022711 0.022720 0.022722 0.022723 0.022723
0.007411 0.010533 0.011611 0.011904 0.011982 0.012002 0.012006 0.012008 0.012008 0.012008
0.003783 0.005402 0.005969 0.006126 0.006166 0.006176 0.006179 0.006179 0.006179 0.006179
0.001911 0.002736 0.003028 0.003108 0.003129 0.003134 0.003135 0.003136 0.003136 0.003136
0.000961 0.001377 0.001525 0.001566 0.001576 0.001579 0.001580 0.001580 0.001580 0.001580
EN
0.067495
0.041129
0.022723
0.012008
0.006179
0.003136
0.001580
Table 3 The maximum absolute error for Example 1, when Dt ¼ 0:0001 e
N 8
16
32
64
128
256
512
22 24 26 28 210 212 214 216 218 220
0.044977 0.061213 0.065642 0.066747 0.067022 0.067091 0.067108 0.067112 0.067113 0.067114
0.026188 0.036307 0.039685 0.040547 0.040763 0.040817 0.040831 0.040834 0.040835 0.040835
0.014169 0.019914 0.021829 0.022338 0.022482 0.022518 0.022527 0.022529 0.022530 0.022530
0.007377 0.010457 0.011520 0.011803 0.011881 0.011900 0.011905 0.011906 0.011907 0.011907
0.003765 0.005362 0.005921 0.006074 0.006114 0.006125 0.006127 0.006128 0.006128 0.006128
0.001902 0.002716 0.003003 0.003082 0.003103 0.003108 0.003109 0.003109 0.003110 0.003110
0.000956 0.001367 0.001512 0.001553 0.001563 0.001566 0.001566 0.001566 0.001566 0.001566
EN
0.067114
0.040835
0.022530
0.011907
0.006128
0.003110
0.001566
A. Kaushik, M.D. Sharma / Applied Mathematics and Computation 195 (2008) 688–706
703
parameter uniform. Further, in Section 7, a set of numerical experiments is carried out in support of the predicted theory, which were found in good accordance with the predicted theory. The maximum absolute error is tabulated for both the examples for different time steps and can be seen from Tables 1–3 for example 1 and Tables 4–6 for example 2. It was observed from the tables that the order Table 4 The maximum absolute error for Example 2, when Dt ¼ 0:01 N
e
8
16
32
64
128
256
512
2 24 26 28 210 212 214 216 218 220
0.007691 0.012118 0.014024 0.014374 0.014449 0.014467 0.014471 0.014472 0.014472 0.014472
0.004348 0.006814 0.007701 0.007865 0.007919 0.007931 0.007935 0.007935 0.007936 0.007936
0.002311 0.003621 0.004086 0.004174 0.004194 0.004199 0.004200 0.004201 0.004201 0.004201
0.001195 0.001867 0.002099 0.002148 0.002160 0.002163 0.002163 0.002163 0.002163 0.002163
0.000607 0.000948 0.001065 0.001089 0.001095 0.001097 0.001097 0.001097 0.001097 0.001097
0.000306 0.000477 0.000536 0.000549 0.000552 0.000552 0.000553 0.000553 0.000553 0.000553
0.000154 0.000240 0.000269 0.000275 0.000277 0.000277 0.000277 0.000277 0.000277 0.000277
EN
0.014472
0.007936
0.004201
0.002163
0.001097
0.000553
0.000277
2
Table 5 The maximum absolute error for Example 2, when At Dt ¼ 0:001 e
N 8
16
32
64
128
256
512
22 24 26 28 210 212 214 216 218 220
0.007368 0.011327 0.012945 0.013244 0.013309 0.013325 0.013328 0.013329 0.013330 0.013330
0.004169 0.006351 0.007084 0.007205 0.007232 0.007239 0.007241 0.007241 0.007241 0.007241
0.002217 0.003369 0.003736 0.003810 0.003828 0.003832 0.003833 0.003833 0.003833 0.003833
0.001147 0.001735 0.001917 0.001953 0.001963 0.001965 0.001966 0.001966 0.001966 0.001966
0.000583 0.000880 0.000971 0.000990 0.000995 0.000996 0.000996 0.000996 0.000996 0.000996
0.000294 0.000443 0.000489 0.000498 0.000501 0.000501 0.000502 0.000502 0.000502 0.000502
0.000147 0.000222 0.000245 0.000250 0.000251 0.000252 0.000252 0.000252 0.000252 0.000252
EN
0.013330
0.007241
0.003833
0.001966
0.000996
0.000502
0.000252
Table 6 The maximum absolute error for Example 2, when At Dt ¼ 0:0001 e
N 8
16
32
64
128
256
512
22 24 26 28 210 212 214 216 218 220
0.007334 0.011247 0.012836 0.013130 0.013194 0.013209 0.013213 0.013214 0.013214 0.013214
0.004150 0.006304 0.007022 0.007141 0.007168 0.007175 0.007176 0.007177 0.007177 0.007177
0.002207 0.003343 0.003701 0.003774 0.003791 0.003795 0.003796 0.003796 0.003797 0.003797
0.001142 0.001722 0.001898 0.001935 0.001943 0.001946 0.001946 0.001946 0.001946 0.001946
0.000580 0.000873 0.000962 0.000980 0.000985 0.000986 0.000986 0.000987 0.000987 0.000987
0.000292 0.000440 0.000484 0.000493 0.000496 0.000496 0.000497 0.000497 0.000497 0.000497
0.000147 0.000221 0.000243 0.000248 0.000249 0.000249 0.000249 0.000249 0.000249 0.000249
EN
0.013214
0.007177
0.003797
0.001946
0.000987
0.000497
0.000249
704
A. Kaushik, M.D. Sharma / Applied Mathematics and Computation 195 (2008) 688–706
of convergence is very good. But the convergence of the time discretization is a little low, the reason lying in the fact that the choice of the temporal step size is governed by the accuracy and stability considerations. In the ideal case, accuracy and stability restrictions on the time step should be about the same. If the time step is much more restricted by stability than by accuracy, efficiency may be gained by switching to a possibly more
1
0.8
ε=2 ε=2
0.9
0.7
ε=2
ε=2
0.8
0.6 ε=2 Computed Solutions
Computed Solutions
0.7 ε=2
0.6 0.5 0.4
0.5 0.4
0.3
0.3 0.2 ε=2 0.2 0.1
0.1 0
0
0.5 T=0.1
0
1
ε=2 0
0.5 T=0.9
1
Fig. 1. Computed solutions of Example 1 for different values of , with N = 128 and M = 10.
0.6
1 T=0.1
T=0.3
0.5
0.8 0.7
0.4
Computed Solutions
Computed Solutions
T=0.1
0.9
0.3
T=0.3
0.2
T=0.6
0.6 T=0.9 0.5 0.4 0.3 0.2
0.1 0.1
T=0.6 0
T=0.9 0
0.5 =2
1
0
0
0.5 =2
1
Fig. 2. Computed solution of Example 1 for different values of T, with N = 128 and M = 10.
A. Kaushik, M.D. Sharma / Applied Mathematics and Computation 195 (2008) 688–706
705
expansive but more stable time discretization. Further, the effect of perturbation parameter and the final time on the layer behavior can be seen from Figs. 1 and 2, and Figs. 3 and 4, for example 1 and example 2, respectively. It can be seen from the figures that this scheme faithfully mimics the dynamics of the corresponding differential equation.
0.4
0.35
ε=2
ε=2 0.35
ε=2
ε=2
0.3
ε=2
0.3 0.25
Computed Solutions
Computed Solutions
0.25 ε=2
0.2 0.15
0.15
0.1
0.1
ε=2
0.05
0.05 0
0.2
0
0.5 T=0.1
ε=2
0
1
0
0.5 T=0.9
1
Fig. 3. Computed solutions of Example 2 for different values of , with N = 128 and M = 10.
0.35
0.4 T=0.1 0.35
0.3
0.3
T=0.3 T=0.6 T=0.9
T=0.3
Computed Solutions
Computed Solutions
0.25
T=0.1
0.2
0.15 T=0.6
0.25
0.2
0.15
0.1 0.1
T=0.9 0.05
0
0
0.05
0.5 ε=2
1
0
0
0.5 ε=2
1
Fig. 4. Computed solution of Example 2 for different values of T, with N = 128 and M = 10.
706
A. Kaushik, M.D. Sharma / Applied Mathematics and Computation 195 (2008) 688–706
References [1] X.Y. Wang, Z.S. Zhu, Y.K. Lu, Solitary wave solutions of the generalized Burger’s–Huxley equation, J. Phys. A: Math. Gen. 23 (1990) 271–274. [2] A.C. Scott, Neurophysics, John Wiley, New York, 1977. [3] J. Satsuma, Explicit solutions of nonlinear equations with density dependent diffusion, J. Phys. Soc. Jpn. 56 (1987) 1947–1950. [4] X.Y. Wang, Nerve propagation and wall in liquid crystals, Phy. Lett. 112A (1985) 402–406. [5] X.Y. Wang, Brochard–Lager wall in liquid crystals, Phys. Rev. A 34 (1986) 5179–5182. [6] X.Y. Wang, Z.S. Zhu, Y.K. Lu, Solitary Waves solutions of the generalized Burger–Huxley equations, J. Phys.: Math. Gen. 23 (1990) 271–274. [7] A. Kaushik, M.D. Sharma, Numerical analysis of a mathematical model describes propagation of electrical pulses in a neuron, Numer. Meth. Partial Differ. Equation, in press. [8] J.M. Burgers, A mathematical model illustrating the theory of turbulence, Adv. Appl. Mech. 1 (1948) 171–199. [9] M.J. Lighthill, Viscosity effects in sound waves of finite amplitude, in: G.K. Batchelor, R.M. Davies (Eds.), Surveys in Mechanics, Cambridge University Press, Cambridge, 1956, pp. 250–351. [10] H.G. Roos, M. Stynes, L. Tobiska, Numerical Methods for Singularly Perturbed Differential Equations. Convection–Diffusion and Flow Problems, Springer-Verlag, Berlin, 1996. [11] H.N.A. Ismail, K. Raslan, A.A.A. Rabboh, Adomian decomposition method for Burgers’–Huxley and Burger’s–Fisher equations, Appl. Math. Comput. 159 (2004) 291–301. [12] P.G. Estevez, Non-classical symmetries and the singular modified the Burger’s and Burgers’–Huxley equation, J. Phys. A 27 (1994) 2113–2127. [13] M. Javidi, A numerical solution of the generalized Burger’s–Huxley equation by pseudospectral method and Darvishi’s preconditioning, Appl. Math. Comput. 175 (2006) 1619–1628. [14] A.A. Samarskii, Theory of Difference Schemes, third ed., Nauka, Moscow, 1989. [15] 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.