# Parameter-uniform fitted mesh higher order finite difference scheme for singularly perturbed problem with an interior turning point

In this paper, a parameter-uniform fitted mesh finite difference scheme is constructed and analyzed for a class of singularly perturbed interior turning point problems. The solution of this class of turning point problem possess two outflow exponential boundary layers. Parameter-explicit theoretical bounds on the derivatives of the analytical solution are given, which are used in the error analysis of the proposed scheme. The problem is discretized by a hybrid finite difference scheme comprises of midpoint-upwind and central difference operator on an appropriate piecewise-uniform fitted mesh. An error analysis has been carried out for the proposed scheme by splitting the solution into regular and singular components and the method has been shown second order uniform convergent except for a logarithmic factor with respect to the singular perturbation parameter. Some relevant numerical examples are also illustrated to verify computationally the theoretical aspects. Numerical experiments show that the proposed method gives competitive results in comparison to those of other methods exist in the literature.

## Authors

• 4 publications
• 1 publication
• 1 publication
06/04/2019

### Parameter uniform essentially first order convergence of a fitted mesh method for a class of parabolic singularly perturbed Robin problem for a system of reaction-diffusion equ

In this paper, a class of linear parabolic systems of singularly perturb...
06/17/2022

### Supercloseness of the local discontinuous Galerkin method for a singularly perturbed convection-diffusion problem

A singularly perturbed convection-diffusion problem posed on the unit sq...
08/18/2021

### Hybrid high-order method for singularly perturbed fourth-order problems on curved domains

We propose a novel hybrid high-order method (HHO) to approximate singula...
09/14/2020

### On construction of a global numerical solution for a semilinear singularly–perturbed reaction diffusion boundary value problem

A class of different schemes for the numerical solving of semilinear sin...
10/19/2020

### Numerical approximations to a singularly perturbed convection-diffusion problem with a discontinuous initial condition

A singularly perturbed parabolic problem of convection-diffusion type wi...
07/27/2021

### Parameter-uniform numerical methods for singularly perturbed linear transport problems

Pointwise accurate numerical methods are constructed and analysed for th...
03/11/2020

### A Non-Standard Finite Difference Scheme for MHD Boundary Layer Fluid Flow

This paper deals with a non-standard finite difference scheme defined on...
##### This week in AI

Get the week's most popular data science and artificial intelligence research sent straight to your inbox every Saturday.

## 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]:

 {Lu(x)≡εu′′(x)+a(x)u′(x)−b(x)u(x)=f(x),x∈Ω=(−1,1),u(−1)=A,u(1)=B, (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

 a(0)=0,a′(0)<0. (1.2)

Moreover, for some constant there exists a positive constant , such that

 |a(x)|≥α>0,δ≤|x|≤1. (1.3)

Also is required to be bounded below by some positive constant , ,

 b(x)≥β>0,x∈¯Ω=[−1,1], (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 :

 |a′(x)|≥∣∣∣a′(0)2∣∣∣,x∈¯Ω=[−1,1]. (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

 ||g||=max¯Ω|g(x)|.

## 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.

([12].) (Uniform Stability Estimate) , solution of the SPTPP (1.1)-(1.5), satisfies the following stability estimate:

 ||u(x)||≤||f||β+max(|A|,|B|),∀x∈¯Ω.

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.

([12].) If and , then solution of the SPTPP (1.1)-(1.2) satisfies the following bounds for any :

 |uj(x)|≤C(1+ε−jexp(−α(1+x)ε)),j=1,⋯,m+1,x∈Ω1,
 |uj(x)|≤C(1+ε−jexp(−α(1−x)ε)),j=1,⋯,m+1,x∈Ω3,

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.

([3].) Let be the solution of SPTPP defined from the equations(1.1)-(1.5), and . Then for and sufficiently small , there exists a positive constant such that

 |u(j)(x)|≤C,j=1,2,…,m,∀x∈Ω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

 u(x)=v(x)+w(x),∀x∈¯Ω,

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 :

 |v(j)(x)|≤C(1+ε((m−1)−j)e−α(1+x)/ε),∀x∈Ω1,
 |w(j)(x)|≤Cε−je−α(1+x)/ε.∀x∈Ω1.

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.

([12].) Let and then for all the smooth component satisfies

 |v(j)(x)|≤C(1+ε((m−1)−j)(exp(−α(1+x)ε)+exp(−α(1−x)ε))),∀x∈¯Ω,

and the singular component satisfies

 |w(i)(x)|≤Cε−i(exp(−α(1+x)ε)+exp(−α(1−x)ε)),∀x∈¯Ω.

## 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

 τ=min{14,τ0εlnN}.

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

 hi=⎧⎪⎨⎪⎩h=4τ/N, i=1,2,…,N/4,H=4(1−τ)/N, i=N/4+1,…,3N/4,h=4τ/N, i=3N/4+1,…,N.

One can easily observe that

 N−1≤H≤4N−1,h=4τ0εN−1lnN

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

 LNU(xi)≡⎧⎪ ⎪ ⎪ ⎪ ⎪⎨⎪ ⎪ ⎪ ⎪ ⎪⎩LNcU≡εδ2Ui+aiD0Ui−biUi=fi,i=1,2,…,N/4−1,LNmpU≡εδ2Ui+ai±1/2D±Ui−(bU)i±1/2=fi±1/2,i=N/4,…,3N/4,LNcU≡εδ2Ui+aiD0Ui−biUi=fi,i=3N/4+1,…,N−1,U0=A,UN=B, (3.1)

where,

 LNmpU≡{εδ2Ui+ai+1/2D+Ui−(bU)i+1/2=fi+1/2,if ai>0εδ2Ui+ai−1/2D−Ui−(bU)i−1/2=fi−1/2,% if ai<0.

Here, we used the following definition to construct above scheme

 vi=v(xi),vi+1/2=vi+vi+12,vi−1/2=vi−1+vi2,ˆhi=hi+hi+12,
 D+vi=vi+1−vihi+1, ,D−vi=vi−vi−1hi, D0vi=vi+1−vi−12ˆhi, δ2vi=(D+vi−D−vi)ˆhi.

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

 pli=(εhiˆhi−ai2ˆhi),pci=(−pli−pri−bi), pri=(εhi+1ˆhi+ai2ˆhi), i=1,2,…,N/4−1, 3N/4+1…,N−1, pli=(εhiˆhi), pci=(−pli−pri−bi+1/2), pri=(εhi+1ˆhi+ai+1/2hi+1−bi+12), if ai>0 pli=(εhiˆhi−ai−1/2hi−bi−12), pci=(−pli−pri−bi−1/2), pri=(εhi+1ˆhi), if ai<0, i=N/4,…,3N/4.

## 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

 h||a||2ε<1,  i.e.,  2τ0||a||

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 :

 pli>0,pri>0,pli+pci+pri<0,i=1,2,…,N−1. (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

 |ZNi|≤1βmax1≤j≤N−1|LNZNj|.

Proof. Let us introduce two comparison functions defined by

 Ψ±i=1βmax1≤j≤N−1|LNZNj|±ZNi.

Clearly one can notice that since Furthermore, for we have

 LNΨ±i=−bβmax1≤j≤N−1|LNZNj|±LNZNi≤0,

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

 |LNcui−(Lu)(xi)|≤C(ε~h2|u(iv)|+~h2|u(iii)|).

On an arbitrary non-uniform mesh, we have

 |LNcui−(Lu)(xi)|≤C(ε(hi+hi+1)|u(iii)|+(hi+hi+1)|u(ii)|).

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:

 |LNmpui−(Lu)(xi−1/2)|≤{C(ε(hi+hi+1)|u(iii)|+h2i+1(|u(iii)|+|u(ii)|+|ui|)),if a(x)>0,C(ε(hi+hi+1)|u(iii)|+h2i(|u(iii)|+|u(ii)|+|ui|)),if a(x)<0.

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

 U=V+W, (4.3a) where, LNV=f,V(−1)=v(−1),V(1)=v(1), (4.3b) LNW=0,W(−1)=w(−1),W(1)=w(1). (4.3c)

Therefore, the error can be written in the form

 U−u=(V−v)+(W−w),

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

 |(V−v)(xi)|≤{CN−2,∀i=0,1,…,N/4−1,3N/4+1,…,N,CN−1(ε+N−1)∀i=N/4,N/4+1,…,3N/4.

Proof. Using the usual truncation error estimates given above and bounds for the smooth component given in Theorem (2.3), we have

 |LN(V−v)(xi)| ≤{CN−2(ε|v(iv)|+|v(iii)|),∀i=0,1,…,N/4−1,3N/4+1,…,N,CN−1(ε|v(iii)|+N−1(|v(iii)|+|v(ii)|+|vi|)),∀i=N/4,N/4+1,…,3N/4. ≤{CN−2,∀i=0,1,…,N/4−1,3N/4+1,…,N,CN−1(ε+N−1),∀i=N/4,N/4+1,…,3N/4,

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 :

 ΦLi=⎧⎪⎨⎪⎩∏ij=1(1+γhjε)−1,i=1,…N/2,1,i=0.ΦRi=⎧⎪⎨⎪⎩∏Nj=i+1(1+γhjε)−1,i=N/2,…N−1,1,i=N. (4.4)

First we prove the following technical result.

###### Lemma 4.4.

If , the barrier functions satsisfy the inequalities

 LNΦLi≤0,∀i=1,2,…,N/2,LNΦRi≤0,∀i=N/2,…,N−1.

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

 LNmpΦLi =(γ2εhi+1ˆhi−ai+1/2γε−bi+12)ΦLi+1−bi2ΦLi =(2γ2ε(hi+12ˆhi−1)+(2γ2ε−ai+1/2γε−bi+12)−bi2(1+γhi+1ε))ΦLi+1≤0.

In the case of central difference operator with , we have

 LNcΦLi =(2γ2ε(hi+12ˆhi−1)+(2γ2ε−aiγεhi+12ˆhi))ΦLi+1−(aiγεhi2ˆhi+bi)ΦLi≤0.

Similarly, applying the midpoint upwind operator for the case , we have

 LNmpΦRi=(2γ2ε(hi2ˆhi−1)+(2γ2ε+ai−1/2γε−bi−12)−bi2(1+γhiε))ΦLi−1≤0.

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

 |Wi|≤CΦLi,∀i=0,1,…N/2,|Wi|≤CΦRi,∀i=N/2,…,N.

Moreover, following bounds are valid for the layer component in no layer region

 |Wi|≤CN−2,∀i=N/4,…,3N/4.

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 :

 ΦLi≤ΦLN4 =N/4∏j=1(1+γhjε)−1=(1+γhε)−N/4=(1+4γτεN)−N/4=(1+4γτ0N−1lnN)−N/4 =(1+8N−1lnN)−N/4=((1+8N−1lnN)−N/8)2≤CN−2,

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

 |(W−w)(xi)|≤{CN−2(lnN)2,∀i=0,1,…,N/4−1,3N/4+1,…,N,CN−2,∀i=N/4,N/4+1,…,3N/4.

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

 |(W−w)(xi)| ≤|W(xi)|+|w(xi)| ≤Ci∏j=1(1+γhjε)−1+Cexp(−α(1+xi)ε) ≤Ci∏j=1(1+γhjε)−1(since e−α(1+xi)/ε≤ΦLi) ≤CN−2(Using Lemma~{}???). (4.5)

Proceeding in a similar manner in subinterval , one can prove

 |(W−w)(xi)|≤CN−2,∀i=N/2,…,3N/4. (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 :

 |LN(W−w)(xi)| ≤Ch2(ε|w(iv)|+|w(iii)|) =16CN−2τ2(ε|w(iv)|+|w(iii)|) ≤CN−2ε2(lnN)2(ε−3exp(−α(1+xi)ε)) =C(N−2(lnN)2εexp(−α(1+xi)ε)) ≤C(N−2(lnN)2εΦLi). (4.7)

From the Eq. (4), , also we have . Therefore, if we choose

 Ψ±(xi)=CN−2(1+(lnN)2ΦLi)±(W−w)(xi),∀i=0,1,…