1 Introduction
Singularly perturbed problems arise often in the modeling of various modern complicated processes, such as viscous flow problems with large Reynolds numbers [9], convective heat transport problems with large Péclet numbers [10], drift diffusion equation of semiconductor device modelling [21], electromagnetic field problems in moving media [8], financial modelling [5] and turbulence models [15] etc. Most of the singularly perturbed problems cannot be completely solved by analytical techniques. Consequently, numerical techniques are getting much attention to get some useful insights on the solutions of singularly perturbed problems. In general, two classes of methods, namely, fitted operator methods and fitted mesh methods have been used to solve such problems.
Those singularly perturbed convection-diffusion problems, in which the convection coefficient vanishes at some points of the domain of the problem, are called singularly perturbed turning point problems (SPTPPs), and zeros of the convection coefficient are said to be turning points. Here, we consider the following class of singularly perturbed two-point boundary value problems with an interior turning point at [12]:
(1.1) |
where is a small perturbation parameter satisfying , and are given constants, and are sufficiently smooth functions. We impose the following restriction to ensure that the solution of Eq. (1.1) exhibits twin boundary layers
(1.2) |
Moreover, for some constant there exists a positive constant , such that
(1.3) |
Also is required to be bounded below by some positive constant , ,
(1.4) |
to guarantee that the operator is inverse monotone on and to exclude the so-called resonance phenomena [2]. We also impose the following restriction to ensure that there are no other turning points in the interval :
(1.5) |
This class of singularly perturbed turning point problem (SPTPP) (1.1) has a unique solution possess twin outflow boundary layers of exponential type at both end points , under the assumptions (1.2)- (1.5).
It is very difficult to deal singularly perturbed turning point problems analytically. The study of these problems received much attention in the literature due to the complexity involved in finding uniformly valid asymptotic expansions unlike non-turning problems. Some authors, such as, Jingde [7], O’Malley [16, 17], Wasow [24] studied qualitative aspects of these problems, namely, existence, uniqueness and asymptotic behavior of the solution.
In general, since the convection coefficient has zero inside the domain therefore numerical treatment of singularly perturbed turning point problem becomes more difficult than the singularly perturbed non-turning point problems. Abrahamsson [1], Berger et al. [3] and Farrell [6] establish a priori bounds for interior turning point problems; in particular it is shown that a bound is independent of singular perturbation parameter if and only if reaction coefficient is greater than zero at the turning point. It is also shown there how the ratio of reaction coefficient and first derivative of convection coefficient , , at the turning point plays a key role in determining the behavior of the solution [3]. It is shown that for , the solution is smooth near turning point and two outflow boundary layers of exponential type exhibits at both the endpoints of the domain. In this case the turning point is sometimes called a diverging flow or expansion turning point. On the other hand, if , there is in general no boundary layers exhibited and an interior layer appears at the turning point, the nature of which depends in a fundamental way on . For , the interior layer is called cusp layer because it can be approximately modelled by a cusp-like function. Interior layer turning point is sometimes called a converging flow or compression turning point. In inviscid fluid dynamics, the diverging flow turning point corresponds to a sonic point and converging flow turning point to a shock point. For the case at the turning point, the solution exhibits a very interesting phenomenon called Ackerberg-O’Malley resonance phenomenon [2].
Berger et al. [3] also show that the modified version of El Mistikawy Werle scheme is uniformly convergent of in the norm using the analytic bounds obtained in [3]. Farrell [6] obtained a set of sufficient conditions for uniform convergence in the discrete norm on uniform mesh, not only for exponentially fitted schemes, but also for a large class of schemes of upwinded type. Kadalbajoo and Patidar [11] gave a numerical scheme based on cubic spline approximation with nonuniform mesh for SPTPP (1.1)-(1.5) and established second order -uniform convergence. Natesan et al. [20] proposed a numerical method based on the classical upwind finite difference scheme on a Shishkin mesh and proved that the proposed scheme is uniformly convergent of almost order one. In [12], Kadalbajoo and Gupta derived asymptotic bounds for the derivatives of the analytical solution of SPTPP (1.1)-(1.5) and proposed a computational method comprises B-spline collocation scheme on a non-uniform Shiskin mesh. They shown that this scheme is second order accurate in the maximum norm. Kadalbajoo et al. [13] also suggested B-spline collocation with artificial viscosity on uniform mesh for the same class of SPTPP (1.1)-(1.5). In [18], Munyakazi and Patidar conclude that convergence acceleration Richardson extrapolation technique on existing numerical schemes for the above class of turning point problem does not improve the rate of convergence. However, Becher and Roos [4] show that Richardson extrapolation on upwind scheme with piecewise-uniform Shishkin mesh works fine and improves the accuracy to under the assumption . Recently, Munyakazi et al. [19] proposed a fitted operator finite difference scheme for singularly perturbed turning point problem having an interior layer and also shown that with Richardson extrapolation technique, accuracy and order of convergence of the scheme can be improved upto two. For a general review of existing literature on asymptotic and numerical analysis of turning point problems, one can see [22].
In this paper, we focus to devise a second order uniformly
convergent finite difference scheme for SPTPP (1.1) on
piecewise uniform mesh of Shishkin type without using any convergence acceleration technique like Richardson extrapolation. The proposed method
combines the midpoint upwind difference scheme and classical
central finite difference scheme on piecewise uniform mesh. The
requirements of higher order truncation error and monotonicity
play a vital role in the construction of this scheme. One can
observe the fact that the classical central difference scheme is
monotone if is relatively large than the convection coefficient , if
, where is the mesh width and has second order truncation error
on uniform mesh. On the other hand, midpoint upwind difference
operator
is monotone for all value of and for relatively large convection coefficient than the reaction coefficient such
that . Moreover, midpoint upwind operator possess second order truncation error away from the boundary layer region.
Also, Shiskin mesh equally distribute the number of mesh points inside and outside the boundary layers, therefore one can gets a coarse mesh region outside the boundary layer
and fine mesh region inside the boundary layer. Utilizing these facts, we employ midpoint upwind difference scheme in coarse
mesh region and central difference operator in fine mesh region of
Shishkin mesh. Since, central difference operator yields first
order truncation error at transition points, we use midpoint
upwind operator on transition points. Such type of higher order
scheme for singularly perturbed non-turning convection-diffusion
problem was introduced by Stynes and Roos [23]. To
analyze the proposed scheme theoretically, we split the numerical
solution into regular and singular components and analyze them
separately by using tools such as
truncation error bounds, discrete
minimum principle and appropriate choices of barrier functions.
Notation. Throughout the paper we use as a generic
positive constant independent of and mesh
parameters. For any given function
( a non-negative integer), is a global maximum norm
over the domain defined by
2 A-priori Estimates for Continuous Problem
In this section some bounds of the exact solution and its derivatives are discussed. These bounds will be needed for error analysis of proposed numerical scheme in later sections. Derivation of these bounds are well known and can be found in [12]. Systematically, we use minimum principle to derive these bounds.
Lemma 2.1.
([12].) (Minimum Principle) Let and . Then implies that
Since the concerned SPTPP (1.1)-(1.5) is linear, minimum principle ensure the existence and uniqueness of the classical solution. Using the above minimum principle, one can easily prove the following uniform stability estimate for the differential operator .
Lemma 2.2.
To exclude the turning point and to obtain the bounds for the solution and its derivatives in the non-turning point region of the domain, we divide the domain into three subdomains as , and such that , where . Further, following theorem gives bounds for the derivatives of in the subintervals and individually.
Theorem 2.1.
Next, we state a theorem, which gives the bounds for the derivatives of the solution in the turning point region and deduce that the solution is smooth in subdomain .
Theorem 2.2.
It turns out that the bounds for continuous solution given in Theorem 2.1 and Theorem 2.2 are not adequate to obtain -uniform error estimate for the proposed scheme. Therefore, to analyze the proposed scheme correctly, we need to derive more precise bounds on these derivatives by decomposing the solution into regular component and singular component as
where the smooth component satisfies homogeneous problem and singular component satisfies homogeneous problem with appropriate boundary conditions. Using the technique given in [12], we get the following bounds for smooth and singular components in the region :
In the same manner, we can obtain analogous estimates for subinterval , while the solution and its derivatives are smooth in the subinterval . Hence, on the whole domain , the bounds on and , and their derivatives are given in the following theorem:
Theorem 2.3.
3 Fitted Mesh Higher-Order Scheme
In this section, first we construct fitted piecewise-uniform mesh of Shishkin type to discretize the domain and then employ a specially designed finite difference scheme on this mesh to discretize the SPTPP (1.1)-(1.2). The fitted mesh is constructed by dividing into three subintervals and such that . For be an integer, divides each of the subintervals and into mesh intervals and with mesh intervals such that . Here, the transition parameter is obtained by taking
The constant is independent of the parameter and the number of mesh points and will be chosen later on during the analysis of proposed scheme. This mesh is coarse on and fine on and on . If and are fine and coarse mesh width respectively, then mesh width is defined as
One can easily observe that
Since, convection coefficient changes its sign at the turning point , therefore, we construct a finite difference scheme to discretize the SPTPP (1.1) in the following manner
(3.1) |
where,
Here, we used the following definition to construct above scheme
It is clear that proposed finite difference operator in scheme (3.1) is a combination of central difference operator and midpoint upwind difference operator , which is constructed by using knowledge judiciously about the sign of the convection term, location of the turning point and truncation error behavior of these operators. After simplifying the terms in (3.1), the difference scheme takes the form , where the coefficients are given by
4 Uniform Convergence
Here, in this section first we shall establish the consistency and stability estimate through discrete minimum principle and then analyze proposed numerical method (3.1) for -uniform convergence by analogous decomposition of discrete solution into smooth and singular components as of continuous solution.
Lemma 4.1.
( Discrete Minimum Principle) Let us suppose that , where
(4.1) |
Then the operator defined by (3.1) satisfies a discrete minimum principle, if is a mesh function that satisfies and for , then for .
Proof. In order to establish the discrete minimum principle, We simply check that the associated system matrix is -matrix with the choice of the midpoint upwind and central difference operator used in the definition of the difference scheme (3.1). It allow us to establish the following inequalities on the coefficients of the difference operator :
(4.2) |
In the case of central difference operator , conditions in (4.2) are satisfied if if , then one can check and for and . For the case of midpoint upwind operator , the conditions in (4.2) are satisfied if if From these sign patterns on the coefficients of associated system matrix, one can deduce that operator is of negative type and therefore satisfies a discrete minimum principle. Moreover, it ensures that the operator is uniformly stable in the maximum norm.
Lemma 4.2.
Let be any mesh function such that Then for all we have
Proof. Let us introduce two comparison functions defined by
Clearly one can notice that since Furthermore, for we have
as Therefore, discrete minimum principle (4.1) implies that , which gives desired result.
Further, using the valid Taylor’s series expansion, we obtained the following truncation error estimates for different finite difference operator employed in the operator : On a uniform mesh with step size , we have
On an arbitrary non-uniform mesh, we have
Here, one can notice that order of truncation error is reduced to one only if the central difference operator is employed on arbitrary non-uniform mesh instead of uniform mesh. Moreover, We have the following truncation error bounds corresponding to the midpoint upwind difference operator, which are valid for both uniform and non-uniform mesh:
Note that the order of truncation error is higher by one in the convection term for midpoint upwind operator than the centered difference operator on a non-uniform mesh. This is the reason to apply midpoint upwind scheme at the transition points and of proposed mesh.
Further the solution of the discrete problem can be decomposed in an analogous manner as that of the continuous solution into the following sum
(4.3a) | |||
where, | |||
(4.3b) | |||
(4.3c) |
Therefore, the error can be written in the form
so the errors in the smooth and singular components of the solution can be estimated separately.
Lemma 4.3.
(Error in smooth component) Assume that satsifies the assumption (4.1). Then the regular component of the error satisfies the following error bound
Proof. Using the usual truncation error estimates given above and bounds for the smooth component given in Theorem (2.3), we have
and applying Lemma 4.2, we obtain the required result.
Since and , we consider both the region and individually to get the error estimates for the layer component . Therefore, we consider the following barrier functions for a positive constant :
(4.4) |
First we prove the following technical result.
Lemma 4.4.
If , the barrier functions satsisfy the inequalities
Proof. We begin with the left hand barrier function and analyze each of the different discretizations used in the definition of the operator . First, in the case of midpoint upwind operator with , we have . Using the properties, and , and with the condition , one can easily observe that
In the case of central difference operator with , we have
Similarly, applying the midpoint upwind operator for the case , we have
In the same manner if we use central difference operator with , we also get It completes the proof.
Lemma 4.5.
The barrier functions and and layer component satisfy
Moreover, following bounds are valid for the layer component in no layer region
Proof. Construct the barrier functions . By Lemma 4.4, we have . Now using the discrete minimum principle we obtain the requred bound. Furthermore, to obtain the bound for in no layer region , we have for :
for the choice of . Here, we have used the inequality with to prove . Using similar argument for barrier function , we obtain desired bounds for in the domain .
Lemma 4.6.
(Error in singular component) Assume that satsifies the assumption (4.1) and . Then the singular component of the error satisfies the following error estimates
Proof. We split our discussion into the two cases of boundary layer region and no boundary layer region to analyze the singular component of the error. Since , it is sufficient to consider only the subinterval and using same argument one can get similar estimate for the subinterval . Both and are small in , therefore we will use triangle inequlaity, Theorem 2.3, Lemma 4.5 instead of the usual truncaton error argument, to get the required error bounds on layer component in . For , using triangle inequality, we have
(4.5) |
Proceeding in a similar manner in subinterval , one can prove
(4.6) |
We now consider the boundary layer region to estimate the singular component of the error. In this case, we obtain the following singular component of the local truncation error estimates for :
(4.7) |
From the Eq. (4), , also we have . Therefore, if we choose