# Fitted Operator Method Using Multiple Fitting Factors for Two Parameters Singularly Perturbed Parabolic Problems

Authors: Dagnachew Mengstie Tefera; Awoke Andargie Tiruneh; Getachew Adamu Derese
Journal: Mathematical Problems in Engineering
 (2022)

---

## Abstract

In this paper, we produceϵ,μ− uniform numerical method for a singularly perturbed parabolic differential equation with two parameters. To approximate the solution, we consider the implicit Euler method for time direction, the finite difference method for spatial direction on a uniform mesh, and the fitted operator method with multiple fitting factors. To accelerate the convergence of the method, the Richardson extrapolation method is applied. The results show that the proposed method is second-order convergent in both temporal and spatial directions. The convergence of the scheme is insensitive to the two perturbation parameters. Two model examples are considered to validate the applicability of the proposed method and produced more accurate results compared to some methods that appear in the literature. Matlab software is used to manipulate the results.


---

## Body

## 1. Introduction

The parabolic singularly perturbed convection-reaction-diffusion equation with two parameters is a type of partial differential equation in which the highest order and the convection term are multiplied by a small parameter. The small parameters change the solution rapidly in a small region (the boundary layer region). Due to the presence of the boundary layer, the numerical methods on a uniform mesh for the singular perturbed problems are insufficient and fail to provide better accuracy unless an unacceptably large number of mesh points are utilized to make the step size at least as small as the perturbation parameters (which is not easy to implement practically). This was the driving force for the development of theϵ-uniform numerical scheme.A variety of situations in which one small parameter is multiplied by the highest derivative term (singularly perturbed equation) can be found in the survey article by Kadalbajoo and Gupta [1] and many researchers. However, two-parameter singularly perturbed problems have been studied by a few researchers. Two-parameter singularly perturbed problems occur in chemical reactor theory by Chen and O’Malley [2], lubrication theory by DiPrima [3], and in some other works such as fluid dynamics, quantum mechanics, chemical reactor theory, meteorology, oceanography, reaction-diffusion process, Navier–Stokes equation of fluid flow at high Reynolds number, Riccati and logistic equations Khader and Adel [4]. Boundary layer occurs in different fields such as in temperature distribution Khader and Babatin [5], in heat generation or absorption Khader [6].During the last two decades, many scholars have presented a variety ofϵ− uniformly convergent numerical methods for one-parameter singularly perturbed parabolic problems [7–13]. Abd-Elhameed et al. [14] analyzed numerical spectral solutions of linear and nonlinear second-order singular, singularly perturbed, and Bratu-type equations using Petrov–Galerkin and collocation methods. Abdelhakem Youssri [15] developed two methodical spectral Legendre’s derivative algorithms to numerically attack the Lane-Emden, Bratu’s, and singularly perturbed type equations.A few scholars have proposedϵ,μ− uniformly convergent numerical methods for two-parameter singularly perturbed problems such as Tesfaye et al. The authors of [16] presented a robust finite difference method based on a nonstandard finite difference operator for solving singularly perturbed parabolic initial boundary value problems with two parameters. Tariku and Gemechis [17] designed the numerical method for singularly perturbed parabolic convection-diffusion problems with two small positive parameters using the Crank-Nicolson method over the nonstandard finite difference method. Das and Mehrmann [18] introduced singularly perturbed parabolic convection-reaction-diffusion problems with two small parameters using a moving mesh-adaptive algorithm in which the meshes are generated by equi-distribution of a special positive monitor function. Brdar and Zarin [19] considered the numerical solution for singularly perturbed problem with two small parameters using a Bakhvalov type mesh via a Galerkin finite element method with piecewise linear functions. Kadalbajoo [20] developed parameter-uniform numerical methods for a class of singularly perturbed parabolic reaction-diffusion problems with two small parameters by means of Rothe’s method and the finite element method in spatial direction on a piecewise uniform mesh of the Shishkin type. Gupta et al. [21] considered the implicit Euler method for time stepping on a uniform mesh and the hybrid monotone difference operator for spatial discretization on a Shishkin mesh for solving singularly perturbed problems with two parameters.In this paper, we present a fitted operator method using multiple fitting factors for two parameters, the parabolic convection-reaction-diffusion initial-boundary value problems. The novelty of the proposed method is that it produces two fitting factors for the two parameters, thus getting a more accurate numerical solution than methods that appear in the literature. The order of convergence of the present method is shown to be two in both the spatial and temporal variables.


## 2. Description of the Problem

Consider the singularly perturbed parabolic initial-boundary value problem with two parameters in a domainD=Dx×Dt with Dx=0,1 and Dt=0,T of the form:(1)ϵuxx+μpx,tux−qx,tux,t−ut=gx,t,x,t∈D.Subject to initial condition and boundary conditions:(2)ux,0=sx,x∈D¯x,u0,t=z1t,u1,t=z2t,t∈D¯t,where 0<ϵ≪1 and 0<μ≪1 are the two parameters. The functions px,t, qx,t, gx,t, sx, z1t, and z2t are assumed to be sufficiently smooth, bounded on D¯, and satisfy px,t≥α>0 and qx,t≥β>0; α and β are real constants. This problem includes both the convection-diffusion problem when μ=1 and the reaction-diffusion problem when μ=0.The widths of the boundary layer depends on the relationship between the two parametersϵ and μ. That is, if μ2<ϵ as ϵ⟶0, then two boundary layers of each width Oϵ appear in the neighborhood of both left and right of the solution domain. When ϵ<μ2 as μ⟶0, a layer of width Oϵ/μ occurs on the left side and a layer of width Oμ appears on the right side.Under the assumption of sufficiently smooth and suitable compatibility conditions on the data, the problem (1) has a unique solution ux,t, which exhibits exponential (or parabolic) boundary layers on both the left and right side of the lateral surface of the domain, depending on the size of the parameters [18]. We set the compatibility condition to s0=z10 and s1=z20, and the data match at the two corners 0,0 and 1,0. This condition ensures that there is a constant C1 independent of the perturbation parameter, such that ∀x,t∈D¯x×D¯t.(3)ux,t−ux,0≤C1t,ux,t−z1t=ux,t≤C11−x,for the proof of equation (3) (see Roos et al. [22]).Lemma 1.
Continuous Minimum Principle: assume thatπx,t∈C2,1D¯ such that πx,t≥0, ∀x,t∈∂D. Then, Lϵ,μπx,t≤0, ∀x,t∈D implies that πx,t≥0, ∀x,t∈D¯, where Lϵ,μπx,t≡ϵπxxx,t+μpx,tπxx,t−qx,tπx,t−πtx,t.Proof.
Assumex∗,t∗∈D¯−∂D such that πx∗,t∗=minx,t∈D¯πx,t and suppose πx∗,t∗<0. This gives that πtx∗,t∗=0, πxx∗,t∗=0, and πxxx∗,t∗≥0, this implies(4)Lϵ,μπx∗,t∗=ϵπxxx∗,t∗+μpx∗,t∗πxx∗,t∗−qx∗,t∗πx∗,t∗−πtx∗,t∗=ϵπxxx∗,t∗−qx∗,t∗πx∗,t∗≥0.
We haveLϵ,μπx,t≥0 which contradicts our assumption.
Thus,πx∗,t∗≥0 which leads to πx,t≥0, ∀x,t∈D¯.Lemma 2.
The bound on the solutionux,t of equation (1) is given by(5)ux,t≤C,x,t∈D¯.Proof.
Using equation (3), we have(6)ux,t≤C1t,x,t∈D¯.
Sincet∈0,T, we got(7)ux,t≤C,x,t∈D¯.Lemma 3.
The bound on the derivative ofu with respect to t is given by(8)utx,t≤C,x,t∈D¯.Proof.
Under the assumption of Lemma1 and Lemma 2, we can prove Lemma 3 using the mean value theorem [20].Lemma 4.
Stability Estimation: letux,t be the solution of the continuous problem for equations (1) and (2). Then, we have(9)ux,t≤β−1g+maxz1t,z2t,where qx,t≥β>0,∀x,t∈D.Proof.
Let the two barrier functionsϕ± and defined as(10)ϕ±=β−1g+maxz1t,z2t±ux,t.
Then,(11)ϕ±0,0=β−1g+maxz1t,z2t±0,0≥maxz10,z20±z10≥0ϕ±1,0=β−1g+maxz1t,z2t±u1,0≥maxz10,z20±z20≥0.
Furthermore, we have(12)Lϵ,μϕ±x,t=ϵϕxx±x,t+μpx,tϕx±x,t−qx,tϕ±x,t−ϕt±x,t=−qx,tβ−1g+maxzt,z2t±Lϵ,μux,t=−qx,tβ−1g+maxz1t,z2t±gx,t≤−g∓gx,t−β−1maxz1t,z2t≤−g∓gx,t≤0.
Therefore, using Lemma1, we obtain ϕ±x,t≥0, ∀x,t∈D¯.Lemma 5.
.For all nonnegative integersi,j such that 0≤i+3j≤4, the solution ux,t satisfies(13)∂i+ju∂xi∂tjD¯≤C1ϵi,ifμ2≤Cϵμϵi,ifμ2≥Cϵ,where the constant C is independent of ϵ and μ and is dependent only on ∂vp/∂xvD¯,∂vq/∂xvD¯, ∂vg/∂xvD¯, for v=0,1,2.Proof.
The reader can refer the proof in Gupta et al. [21].


## 3. Derivation of the Numerical Method

In this section, we develop a fitted operator finite difference method with multiple fitting factors to solve equation (1) together with equation (2).


### 3.1. Time Discretization

We begin by semidiscretizing the problem in the time direction by employing the implicit Euler method with uniform time steps ofk=T/n in the time interval 0,T. This results in a time-mesh Dtn≡tj=jk,j=0,1,2,…,n, and the method results in a system of boundary value problems.(14)ϵuxxj+1x+μpj+1xuxj+1x−qj+1xuj+1x−uj+1x−ujxk=gj+1x,where ux,tj+1=uj+1x, px,tj+1=pj+1x, qx,tj+1=qj+1x, gx,tj+1=gj+1x. We can rewrite equation (14) as a differential equation in the spatial variable,(15)ϵuxxj+1x+μpj+1xuxj+1x−q¯j+1xuj+1x=g¯j+1x,uj+10=z1tj+1,uj+10=z1tj+1,where q¯j+1x=qj+1x+1/k and g¯j+1x=gj+1x−ujx/k.Lemma 6.
Local Error Estimate: the local error estimate in the time discretization is given by(16)Tj+1∞≤Ck2.Proof.
The reader can see the proof in Clavero et al. [23].Lemma 7.
Global Error Estimate: under the assumption of Lemma6, we have(17)Ej∞≤C3k.Proof.
Using Lemma6, the global error estimate at the jth time level is given by(18)Ej∞=∑l=1jTl∞,≤T1∞+T2∞+⋯+Tj∞,≤Cjkk,using Lemma6,≤CTk,usingj≤Tk≤C3k,where C,C3 are positive real constants independent of perturbation parameters ϵ and μ.
Equation (15) at j+1th time level can be written in terms of spatial variables as(19)L∗y≡ϵy″x+μp∗xy′x−q∗xyx=g∗xy0=z1tj+1,y1=z2tj+1,,where yx=uj+1x,p∗x=pj+1x,q∗x=q¯j+1x,andg∗x=g¯j+1x.
Equation (19) satisfies the following continuous minimum principle and stability estimate for semidiscretization.Lemma 8.
Continuous Minimum Principle for Semidiscretization: assume thatφx∈C2Dx¯ such that φ0≥0, φ1≥0. Then, L∗φx≤0, ∀x∈Dx implies that φx≥0, ∀x∈Dx¯.Proof.
This lemma can be proved using the same approach as Lemma1.Lemma 9.
Stability Estimate for Semidiscretization: Letyx be the solution of the continuous problem for equation (19). Then, we have(20)yx≤β−1g∗+maxz1tj+1,z2tj+1,where q∗x≥β>0,∀x∈Dx.Proof.
This lemma can be proved using the same approach as Lemma4.


### 3.2. Spatial Discretization

LetDxm be the partition of the spatial domain 0,1 into m subinterval of a uniform mesh as(21)xi=ih,fori=0,1,2,…,m,where h=1/m is step length in the spatial direction such that h>ϵ. The resulting boundary value problems are treated using the multiple exponentially fitted operator finite difference method. For equation (15), multiple exponential fitting parameters are first developed, followed by the spatial discretization.


### 3.3. Computing Multiple Exponential Fitting Parameters

We consider the problem in which the boundary layers are on the right and left sides of the domain. The solution of equation (15) can be described by the roots of the characteristic equation at the j+1th time level:(22)ϵλj+1x2+μpj+1xλj+1x−q¯j+1x=0.At thej+1th time level, equation (22) yields two continuous functions.(23)λ1j+1x=−μpj+1x2ϵ−μpj+1x2ϵ2+q¯j+1xϵ,(24)λ2j+1x=−μpj+1x2ϵ+μpj+1x2ϵ2+q¯j+1xϵ.To handle the effect of the singular perturbation parameters, multiply the diffusion term and the convection term by exponentially fitting parametersσρ and τρ, respectively, and apply a finite difference scheme for spatial variable, and then, equation (15) becomes(25)ϵσiρD+D−uij+1+μτiρpj+1xiD+uij+1−q¯j+1xiuij+1=g¯j+1xi,uj+10=z1tj+1,uj+11=z2tj+1,fori=1,2,…,m−1,where D+D−uij+1=ui−1j+1−2uij+1+ui+1j+1/h2,D+uij+1=ui+1j+1−uij+1/h,andρ=h/ϵ.The fitting factors atj+1th time level are obtained by substituting equation (23) and equation (24) in the corresponding homogeneous difference equation (25) as shown in the following equation:(26)σiρ=−q¯j+1xiρh4epj+1xih/2ϵsinhλ1j+1xih/2sinhλ2j+1xih/2,and(27)τiρ=q¯j+1xih2pj+1xicothλ1j+1xih2+cothλ2j+1xih2,for j=0,1,2,…,n and i=1,2,…,m−1.The scheme (25) can be written as three recurrence relation terms of the form:(28)Eij+1ui−1j+1+Fij+1uij+1+Gij+1ui+1j+1=Hij+1,where(29)Eij+1=ϵσiρh2,Fij+1=−2ϵσiρh2−μτiρhpj+1xi−q¯j+1xi,Gij+1=ϵσiρh2+μτiρhpj+1xi,Hij+1=g¯j+1xi,for j=0,1,2,…,nandi=1,2,…,m−1.The tridiagonal system of equation (28) can be solved by the Thomas algorithm [11].Lemma 10.
Discrete Minimum Principle: let the discrete operator of equation (27) is denoted by Lϵ,μm,n, and φij+1 is any mesh that satisfies φj+1x0≥0 and φj+1xm≥0. If Lϵ,μφxi≤0 for 1≤i≤m−1, then φj+1xi>0,for alli,0≤i≤m.Proof.
Assume thatφz+1xy=min1≤i≤m−1φj+1xi and suppose that φz+1xy<0, then(30)Lϵ,μφz+1xy=ϵσyρφz+1xy−1−2φz+1xy+φz+1xy+1h2+μτyρpz+1xyφz+1xy+1−φz+1xyh−q¯z+1xyφz+1xy>0,which contradicts our assumption. Thus, φj+1xi>0,for alli,0≤i≤m.Lemma 11.
Discrete Stability Estimate: the solutionUij+1 of the discrete method satisfies the bound(31)Uij+1≤β−1maxLϵ,μUij+1+maxz1tj+1,z2tj+1,where q¯j+1xi≥β>0,∀i,j.Proof.
Define the barrier functionψi,j+1± as(32)ψi,j+1±=β−1maxLϵ,μUij+1+maxz1tj+1,z2tj+1,∞.
Using the barrier function at the boundary conditions, we obtain(33)ψj+1±0=β−1maxLϵ,μUij+1+maxz1tj+1,z2tj+1±z1tj+1≥maxz1tj+1,z2tj+1±z1tj+1≥0,ψj+1±1=β−1maxLϵ,μUij+1+maxz1tj+1,z2tj+1±z2tj+1≥maxz1tj+1,z2tj+1±z2tj+1≥0.
Furthermore, we have(34)Lϵ,μψi,j+1±=Lϵ,μβ−1maxLϵ,μUij+1+maxz1tj+1,z2tj+1±Lϵ,μUij+1=−q¯j+1xiβ−1maxLϵ,μUij+1+maxz1tj+1,z2tj+1±Lϵ,μUij+1≤−maxLϵ,μUij+1±Lϵ,μUij+1≤0.
Using Lemma10, we have that ψi,j+1±≥0 for all xi,tj+1∈D¯. Thus, the required bound in Lemma 11 is satisfied.


## 3.1. Time Discretization

We begin by semidiscretizing the problem in the time direction by employing the implicit Euler method with uniform time steps ofk=T/n in the time interval 0,T. This results in a time-mesh Dtn≡tj=jk,j=0,1,2,…,n, and the method results in a system of boundary value problems.(14)ϵuxxj+1x+μpj+1xuxj+1x−qj+1xuj+1x−uj+1x−ujxk=gj+1x,where ux,tj+1=uj+1x, px,tj+1=pj+1x, qx,tj+1=qj+1x, gx,tj+1=gj+1x. We can rewrite equation (14) as a differential equation in the spatial variable,(15)ϵuxxj+1x+μpj+1xuxj+1x−q¯j+1xuj+1x=g¯j+1x,uj+10=z1tj+1,uj+10=z1tj+1,where q¯j+1x=qj+1x+1/k and g¯j+1x=gj+1x−ujx/k.Lemma 6.
Local Error Estimate: the local error estimate in the time discretization is given by(16)Tj+1∞≤Ck2.Proof.
The reader can see the proof in Clavero et al. [23].Lemma 7.
Global Error Estimate: under the assumption of Lemma6, we have(17)Ej∞≤C3k.Proof.
Using Lemma6, the global error estimate at the jth time level is given by(18)Ej∞=∑l=1jTl∞,≤T1∞+T2∞+⋯+Tj∞,≤Cjkk,using Lemma6,≤CTk,usingj≤Tk≤C3k,where C,C3 are positive real constants independent of perturbation parameters ϵ and μ.
Equation (15) at j+1th time level can be written in terms of spatial variables as(19)L∗y≡ϵy″x+μp∗xy′x−q∗xyx=g∗xy0=z1tj+1,y1=z2tj+1,,where yx=uj+1x,p∗x=pj+1x,q∗x=q¯j+1x,andg∗x=g¯j+1x.
Equation (19) satisfies the following continuous minimum principle and stability estimate for semidiscretization.Lemma 8.
Continuous Minimum Principle for Semidiscretization: assume thatφx∈C2Dx¯ such that φ0≥0, φ1≥0. Then, L∗φx≤0, ∀x∈Dx implies that φx≥0, ∀x∈Dx¯.Proof.
This lemma can be proved using the same approach as Lemma1.Lemma 9.
Stability Estimate for Semidiscretization: Letyx be the solution of the continuous problem for equation (19). Then, we have(20)yx≤β−1g∗+maxz1tj+1,z2tj+1,where q∗x≥β>0,∀x∈Dx.Proof.
This lemma can be proved using the same approach as Lemma4.


## 3.2. Spatial Discretization

LetDxm be the partition of the spatial domain 0,1 into m subinterval of a uniform mesh as(21)xi=ih,fori=0,1,2,…,m,where h=1/m is step length in the spatial direction such that h>ϵ. The resulting boundary value problems are treated using the multiple exponentially fitted operator finite difference method. For equation (15), multiple exponential fitting parameters are first developed, followed by the spatial discretization.


## 3.3. Computing Multiple Exponential Fitting Parameters

We consider the problem in which the boundary layers are on the right and left sides of the domain. The solution of equation (15) can be described by the roots of the characteristic equation at the j+1th time level:(22)ϵλj+1x2+μpj+1xλj+1x−q¯j+1x=0.At thej+1th time level, equation (22) yields two continuous functions.(23)λ1j+1x=−μpj+1x2ϵ−μpj+1x2ϵ2+q¯j+1xϵ,(24)λ2j+1x=−μpj+1x2ϵ+μpj+1x2ϵ2+q¯j+1xϵ.To handle the effect of the singular perturbation parameters, multiply the diffusion term and the convection term by exponentially fitting parametersσρ and τρ, respectively, and apply a finite difference scheme for spatial variable, and then, equation (15) becomes(25)ϵσiρD+D−uij+1+μτiρpj+1xiD+uij+1−q¯j+1xiuij+1=g¯j+1xi,uj+10=z1tj+1,uj+11=z2tj+1,fori=1,2,…,m−1,where D+D−uij+1=ui−1j+1−2uij+1+ui+1j+1/h2,D+uij+1=ui+1j+1−uij+1/h,andρ=h/ϵ.The fitting factors atj+1th time level are obtained by substituting equation (23) and equation (24) in the corresponding homogeneous difference equation (25) as shown in the following equation:(26)σiρ=−q¯j+1xiρh4epj+1xih/2ϵsinhλ1j+1xih/2sinhλ2j+1xih/2,and(27)τiρ=q¯j+1xih2pj+1xicothλ1j+1xih2+cothλ2j+1xih2,for j=0,1,2,…,n and i=1,2,…,m−1.The scheme (25) can be written as three recurrence relation terms of the form:(28)Eij+1ui−1j+1+Fij+1uij+1+Gij+1ui+1j+1=Hij+1,where(29)Eij+1=ϵσiρh2,Fij+1=−2ϵσiρh2−μτiρhpj+1xi−q¯j+1xi,Gij+1=ϵσiρh2+μτiρhpj+1xi,Hij+1=g¯j+1xi,for j=0,1,2,…,nandi=1,2,…,m−1.The tridiagonal system of equation (28) can be solved by the Thomas algorithm [11].Lemma 10.
Discrete Minimum Principle: let the discrete operator of equation (27) is denoted by Lϵ,μm,n, and φij+1 is any mesh that satisfies φj+1x0≥0 and φj+1xm≥0. If Lϵ,μφxi≤0 for 1≤i≤m−1, then φj+1xi>0,for alli,0≤i≤m.Proof.
Assume thatφz+1xy=min1≤i≤m−1φj+1xi and suppose that φz+1xy<0, then(30)Lϵ,μφz+1xy=ϵσyρφz+1xy−1−2φz+1xy+φz+1xy+1h2+μτyρpz+1xyφz+1xy+1−φz+1xyh−q¯z+1xyφz+1xy>0,which contradicts our assumption. Thus, φj+1xi>0,for alli,0≤i≤m.Lemma 11.
Discrete Stability Estimate: the solutionUij+1 of the discrete method satisfies the bound(31)Uij+1≤β−1maxLϵ,μUij+1+maxz1tj+1,z2tj+1,where q¯j+1xi≥β>0,∀i,j.Proof.
Define the barrier functionψi,j+1± as(32)ψi,j+1±=β−1maxLϵ,μUij+1+maxz1tj+1,z2tj+1,∞.
Using the barrier function at the boundary conditions, we obtain(33)ψj+1±0=β−1maxLϵ,μUij+1+maxz1tj+1,z2tj+1±z1tj+1≥maxz1tj+1,z2tj+1±z1tj+1≥0,ψj+1±1=β−1maxLϵ,μUij+1+maxz1tj+1,z2tj+1±z2tj+1≥maxz1tj+1,z2tj+1±z2tj+1≥0.
Furthermore, we have(34)Lϵ,μψi,j+1±=Lϵ,μβ−1maxLϵ,μUij+1+maxz1tj+1,z2tj+1±Lϵ,μUij+1=−q¯j+1xiβ−1maxLϵ,μUij+1+maxz1tj+1,z2tj+1±Lϵ,μUij+1≤−maxLϵ,μUij+1±Lϵ,μUij+1≤0.
Using Lemma10, we have that ψi,j+1±≥0 for all xi,tj+1∈D¯. Thus, the required bound in Lemma 11 is satisfied.


## 4. Convergence Analysis

The truncation error of the fitted operator finite difference scheme (28) is given by(35)T.E.=Lϵ,μuxi,tj+1−uij+1.Consider the following Taylor series expansion at the pointxi,tj+1(36)ui+1j+1=uij+1+hux,ij+1+h22uxx,ij+1+h33!uxxx,ij+1+⋯ui−1j+1=uij+1−hux,ij+1+h22uxx,ij+1−h33!uxxx,ij+1+⋯,and(37)uij=uij+1−kut,ij+1+k22utt,ij+1−k33!uttt,ij+1+⋯.By using the above expansions, the truncation error becomes(38)T.E.=ϵh2uij+1+hux,ij+1+h22uxx,ij+1+⋯−2uij+1+uij+1−hux,ij+1+h22uxx,ij+1−⋯+μpj+1xihuij+1+hux,ij+1+h22uxx,ij+1+⋯−uij+1+qj+1xiuij+1−1kuij+1−uij+1−kut,ij+1+k22utt,ij+1−⋯−gxi,tj+1.After simplifying the terms, truncation error becomes(39)T.E.=k2utt,ij+1+ϵh212uxxxx,ij+1+μpj+1xih2uxx,ij+1+Oh2+k2.The major outcome of this work is described in the following theorem.Theorem 1.
Let the exact solution and the numerical solution using the current method of the equations (1) and (2), respectively, are uxi,tj+1 and Uij+1. Then, there exists a constant C independent of ϵ, μ, k, and h such that(40)max0<i<m,0<j<nuxi,tj+1−Uij+1≤Ch+k.This result indicates that the proposed method is first-order convergent in both time and spatial variables.Proof.
The reader can refer the detailed proof in Das and Mehrmann [18].


## 5. Richardson Extrapolation Method

The Richardson extrapolation approach is designed to improve the accuracy of the computed solutions in the basic scheme. From equation (40), we have(41)uxi,tj+1−Uij+1≤Ch+k,where uxi,tj+1 and Uij+1 are the exact and numerical solution for Dmn mesh interval, respectively, and C is a constant independent of ϵ, μ, k and h.LetD2m2n be the mesh found by doubling the mesh interval of Dmn and the numerical solution using D2m2n mesh interval is U¯ij+1. Equation (41) with mesh size h,k≠0, becomes(42)uxi,tj+1−Uij+1=Ch+k+Rmn,xi,tj+1∈Dmn,and equation (42) with mesh size h/2,k/2≠0, becomes(43)uxi,tj+1−U¯ij+1=Ch2+k2+R2m2n,xi,tj+1∈D2m2n,where the remainders Rmn and R2m2n are Oh2+k2. Subtracting equation (42) from (43) to obtain the following extrapolation formula:(44)Uij+1ext=2U¯ij+1−Uij+1,where Uij+1ext is the approximate solution of uxi,tj+1 using the Richardson extrapolation method with truncation error.(45)uxi,tj+1−Uij+1ext≤Ch2+k2.Thus, using the Richardson extrapolation method, the proposed method becomes second-order convergent in both time and spatial variables. Furthermore, the reader can refer to Das [24], and Das and Natesan [25].


## 6. Numerical Examples and Results

To determine the efficiency of the current scheme, we look at model problems that had been addressed in the literature and had approximate solutions that could be compared.We use the double mesh principle to estimate the absolute maximum error of the current approach when the exact solution to the given problem is unknown Doolan et al. [26]. We use the following formula to estimate the absolute maximum error at the selected mesh points:Case 1.
If the exact solution is known,(46)Eϵ,μm,n=maxxi,tj∈Duxi,tj−uijext.Case 2.
If the exact solution is unknown,(47)Eϵ,μm,n=maxxi,tj∈Duij−u2i−12j−1ext,where uijext is the numerical solution obtained by using m number of subintervals of space variable and n number of subintervals of the time variable, and u2i−12j−1ext is the numerical solution obtained by using 2m number of subintervals of the space variable and 2n number of subintervals of the time variable after applying Richardson extrapolation method.
We also evaluate the corresponding rate of convergence using the following formula:(48)Rϵ,μm,n=logEϵ,μm,n−logEϵ,μ2m,2nlog2.Example 1.
Consider equation (1) with px,t=1+x, qx,t=1 and gx,t=16x21−x2∀x,t∈D=0,1×0,1 subject to the initial and boundary conditions(49)ux,0=0,∀x∈0,1andu0,t=0=u1,t,∀t∈0,1.Example 2.
Consider equation (1) with px,t=1+x1−x+t2, qx,t=1+5xt, and gx,t=x1−xet−1∀x,t∈D=0,1×0,1 subject to the initial and boundary conditions(50)ux,0=0,∀x∈0,1andu0,t=0=u1,t,∀t∈0,1.


## 7. Discussion and Conclusions

We proposed a fitted operator numerical method using the multiple fitting factors for solving the singularly perturbed two-parameter parabolic partial differential equations with initial and boundary conditions. The basic mathematical procedures are defining the model problem, approximating the time variable with implicit Euler’s method, approximating the second- and first-order differential equations with central and forward finite difference methods, multiplying the two parameters with exponential fitting factor, and reducing them into a tridiagonal system of equations, and this can be solved using Thomas algorithm.To validate the applicability of the proposed method, we used two model examples. Comparisons with other methods that appear in the literature are shown in Tables1 and 2. The maximum absolute error and the rate of convergence are shown in Tables 3 and 4 with different values of ϵ,μ,m,andn. The physical behaviors of the solution are shown in Figures 1 and 2.Table 1
Comparison of absolute maximum errors for example1 at m=n=32.
μ↓ϵ⟶10−410−610−810−1010−12Present method10−22.6245e-053.0339e-053.0339e-053.0339e-053.0339e-0510−42.6245e-053.0339e-053.0339e-053.0339e-053.0339e-0510−62.6245e-053.0339e-053.0339e-053.0339e-053.0339e-0510−82.6245e-053.0339e-053.0339e-053.0339e-053.0339e-0510−102.6245e-053.0339e-053.0339e-053.0339e-053.0339e-05Gupta et al. [21]10−21.1635e-022.1259e-022.1397e-022.1398e-022.1398e-0210−42.2124e-022.2039e-022.2101e-022.2100e-022.2100e-0210−62.2120e-022.2033e-022.2057e-022.2060e-022.2060e-0210−82.2120e-022.2033e-022.2057e-022.2059e-022.2060e-0210−102.2120e-022.2033e-022.2057e-022.2059e-022.2060e-02Table 2
Comparison of absolute maximum errors and rate of convergence for example2 at μ=10−7.
ϵ↓m,n⟶m = 64, n = 16m = 128, n = 32m = 256, n = 64m = 512, n = 128Present method10−66.1013e-051.5300e-053.8299e-069.5763e-071.99561.99821.999810−76.1013e-051.5300e-053.8299e-069.5802e-071.99561.99821.999210−86.1013e-051.5300e-053.8299e-069.5802e-071.99561.99821.999210−96.1013e-051.5300e-053.8299e-069.5802e-071.99561.99821.9992Das and Mehrmann [18]10−69.6949e-044.9906e-042.5233e-041.2824e-040.95800.98400.976410−79.8712e-045.0050e-042.5485e-041.2853e-040.95800.97370.987610−89.5128e-045.0026e-042.5237e-041.2781e-040.92720.98710.982610−99.6746e-045.0012e-042.5237e-041.2804e-040.95190.97390.9919Table 3
Maximum absolute point-wise error and rate of convergence for example1 before and after Richardson extrapolation method applied, when ϵ=10−10 where L≥10.
μ↓ M,N ⟶16, 1632, 3264, 64128, 128256, 256After10−21.1762e-043.0339e-057.7088e-061.9427e-064.8766e-0710−41.1762e-043.0339e-057.7088e-061.9427e-064.8766e-0710−61.1762e-043.0339e-057.7088e-061.9427e-064.8766e-0710−81.1762e-043.0339e-057.7088e-061.9427e-064.8766e-0710−101.1762e-043.0339e-057.7088e-061.9427e-064.8766e-07Rate of convergence10−L1.95491.97661.98841.9942Before10−25.5124e-032.8165e-031.4228e-037.1498e-043.5838e-0410−45.5124e-032.8165e-031.4228e-037.1498e-043.5838e-0410−65.5124e-032.8165e-031.4228e-037.1498e-043.5838e-0410−85.5124e-032.8165e-031.4228e-037.1498e-043.5838e-0410−105.5124e-032.8165e-031.4228e-037.1498e-043.5838e-04Rate of convergence10−L0.96880.98520.99280.9964Table 4
Maximum absolute point-wise error and rate of convergence for example2 before and after Richardson extrapolation method applied, when ϵ=10−10 where L≥10.
μ↓m, n⟶16, 1632, 3264, 64128, 128256, 256After10−26.0984e-051.5300e-053.8286e-069.5798e-072.3957e-0710−46.0984e-051.5300e-053.8286e-069.5798e-072.3957e-0710−66.0984e-051.5300e-053.8286e-069.5798e-072.3957e-0710−86.0984e-051.5300e-053.8286e-069.5798e-072.3957e-0710−106.0984e-051.5300e-053.8286e-069.5798e-072.3957e-07Rate of convergence10−L1.99491.99861.99871.9995Before10−21.8034e-038.9261e-044.4396e-042.2137e-041.1053e-0410−41.8034e-038.9261e-044.4396e-042.2137e-041.1053e-0410−61.8034e-038.9261e-044.4396e-042.2137e-041.1053e-0410−81.8034e-038.9261e-044.4396e-042.2137e-041.1053e-0410−101.8034e-038.9261e-044.4396e-042.2137e-041.1053e-04Rate of convergence10−L1.01461.00761.00441.0020Figure 1
Log-log plot of maximum absolute errors before and after Richardson extrapolation method applied for example1 for μ=10−2.Figure 2
The physical behavior of the solutions for example2 at m = n = 32, when ϵ=10−10 and μ=10−6.To improve the accuracy of the method, we used the Richardson extrapolation technique and show that the resulting finite difference scheme isϵ,μ-uniformly convergent (the convergence of the method is insensitive to the two perturbation parameters). And the method has order two convergence in both spatial and time variables. The suggested scheme’s performance is investigated by comparing the results, and it is shown that the numerical results’ accuracy is better than that of existing difference schemes that appear in the literature.



---

*Source: 6267522-2022-07-14.xml*