# Local RBF-FD-Based Mesh-free Scheme for Singularly Perturbed Convection-Diffusion-Reaction Models with Variable Coefficients

Authors: Ram Jiwari; Sukhveer Singh; Paramjeet Singh
Journal: Journal of Mathematics
 (2022)

---

## Abstract

This work analyze singularly perturbed convection-diffusion-reaction (CDR) models with two parameters and variable coefficients by developing a mesh-free scheme based on local radial basis function-finite difference (LRBF-FD) approximation. In the evolvement of the scheme, time derivative is discretized by forward finite difference. After that, LRBF-FD approximation is used for spatial discretization, and we obtained a system of linear equations. Then, the obtained linear system is solved by LU decomposition method in MATLAB. For numerical simulation, four singularly perturbed models are pondered to check the efficiency and chastity of the proposed scheme.


---

## Body

## 1. Introduction

Singularly perturbed models (SPMs) can be seen in different areas of science, medicine, and engineering. Herein, we consider the following class of parabolic SPMs with two small parameters and variable coefficients:(1)Lε,μwx,t≡εwxxx,t+μax,twxx,t−bx,twtx,t−cx,tw=fx,t,x,t∈Ω=0,1×0,T,with initial and boundary conditions (BCs)(2)wx,0=w0x,x∈0,1,w0,t=g1t,w1,t=g2t,t∈0,T,where ε and μ are small positive parameters and the variable coefficients ax,t,bx,t,cx,t, and fx,t are sufficiently smooth functions on Ω¯ and satisfies the following conditions:(3)ax,t≥α>0,bx,t≥β>0,cx,t≥γ>0,∀x,t∈Ω¯.We also suppose that the initial and BCs have enough smoothness on∂Ω and sufficient compatibility conditions on the edges of the domain Ω. These conditions on the initial and BCs ensure that there exists a unique classical solution wx,t∈C4,2Ω¯ [1].In the SP case, the solutions of the abovementioned models may possess a multiscale character, presenting rapid variations in some narrow regions adjacent to both the left and right lateral boundary depending on the magnitude of the parametersε and μ. For μ = 1, ε ⟶ 0, model (1) becomes convection dominated model (CDM) and for μ = 0, ε ⟶ 0, model (1) was converted into the reaction-diffusion model (RDM). There are various numerical techniques available in literature for the simulation of CDM and RDM, e.g., adaptive B-spline collocation method [2], ε-uniform schemes [3], parameter uniform difference scheme [4], uniformly convergent B-spline collocation method [5], finite difference fitted schemes [6], uniformly convergent difference schemes [7, 8], SDFEM [9], parameter-uniform hybrid finite difference scheme [10], finite difference domain decomposition algorithms [11], high order methods [12], Layer-adapted meshes and FEM [13], parameter uniform approximations [14], uniformly convergent scheme [15], and finite difference scheme [16]. Except these, there are many schemes for these types of problems [17–22].The singularly perturbed CDR models occurred in mathematical modelling of physical problems, such as fluid dynamics, elasticity, quantum mechanics, heat and mass transfer in mechanical engineering, theory of plates and shells, oil and gas reservoir simulation, and magneto-hydrodynamic flow. Two parameter problems arise in chemical flow reactor theory [23] as well as in the case of boundary layers controlled by suction (or blowing) of some fluids [17, 18, 24]. It is important to study these CDR models with two parameters and variable coefficients so that other information can be clarify by the interaction between convection and diffusion.Recently, singularly perturbed convention diffusion problems are solved by Ahmad and Kalarestaghi [25] using LRBFs methods. But, these problems are steady state and with constant coefficients. Secondly, upwind LRBFs and central LRBFs methods are used for solving these problems. This motivates the authors to develop LRBF-FD scheme for unsteady-state singularly perturbed CDR models. To achieve this aim, we present a mesh-free scheme based on LRBF-FD approximation for the unsteady-state singularly perturbed CDR models with two parameters and variable coefficients. The novelty of the work is that we solved unsteady-state singularly perturbed problems with two parameters and variable coefficients. Second, we used LRBF-FD approximation in place of upwind LRBFs and central LRBFs methods. In the evolvement of the scheme, time derivative is discretized by forward finite difference. After that, RBF-FD approximation is used for spatial discretization, and we obtained a system of linear equations. Finally, the obtained linear system is solved by the LU decomposition method in MATLAB. For numerical simulation, four problems are pondered to check the efficiency and chastity of the proposed scheme.Rest of the paper is organized as follows: in Section2, the details of RBFs and Local RBF-FD approximation has been described. Mesh-free numerical scheme based on LRBF-FD approximation will be introduced in Section 3. In Section 4, MATLAB implementation steps are given in form of algorithm of the developed scheme. Numerical simulation and discussion are carried out in Section 5. Concluding remarks of the work are given in the final Section 6.


## 2. RBFs and Local RBF-FD Approximation

Herein, we presented the brief introduction of RBFs and LRBF-FD approximation inℜn space as follows.


### 2.1. RBFs

For scattered pointsΧ∈Ξ⊂ℜn, the approximation of a function χX can be written as a linear combination of RBFs [26, 27].(4)χΧ≅∑j=1MλjφΧ−Χj2+ψΧ,where M is the number of centers or sometimes called knotsΧ, Χ=x1,x2,…,xn,n is the dimension of the problem, λ′s are coefficients to be determined and φ is a RBF. (4) can be written without the additional polynomial ψΧ. If Ψqn denotes the space of n-variate polynomials of order not exceeding q, and letting the polynomials P1,P2,…,Pm be the basis of Ψqn in ℜn, then the polynomial ψΧ, in (4), is usually written in the following form:(5)ψX=∑i=1mζiPiX,where m=q−1+n!/n!q−1!.To determine the coefficientsλ1,λ2,…,λN and ζ1,ζ2,…,ζm extra m equations are required in addition to the M equations resulting from the collocating equation (4) at M knots. This is insured by the m conditions for (4), viz.,(6)∑j=1MλjPiXj=0,i=1,2,…,m.There are many RBFs (expression ofφ) available in literature such as Multiquadric (MQ), Inverse Multiquadric (IMQ), Gaussian, Thin Plate spline [27–34] etc.


### 2.2. LRBF-FD Method Formulation

In this work, we used LRBF-FD approximation for the discretization of spatial derivatives. LRBF-FD approximation is more general form of classical finite difference method. To explain the concept of LRBF-FD approximation, assume that there areM scattered points in the bounded domain of the problem i.e. Ξ⊂ℜn. In this approximation, the differential operator ℓ. is approximated at a particular point xc∈Ξ as a weighted linear combination of function values at closest ns neighbor stencil of the point xc as follows:(7)ℓwxc≅∑j∈Ixcαjwj,where Ixc=j∈1,2,…,M:xj−xc2≤δ and δ is the radius of support domain.Now, the main task is to calculate the unknown coefficientsαj occurred in (7). For this purpose, assume the approximation of the unknown solution w. is as follows [35, 36].(8)wx=∑k∈Ixβkφx−xk2,x∈Ξ,where φ. is MQ-RBF discussed in above subsection.Apply the differential operatorℓ. over (8) and take x=xc, we get(9)ℓwxc=∑k∈Ixβkℓφxc−xk2,xc∈Ξ.Now, using the (7)–(9), we can write as follows:(10)∑k∈Ixβkℓφxc−xk2=∑j∈Ixcαjwj=∑j∈Ixcαj∑k∈Ixβkφxj−xk2.To compute the unknown weightsαj, the system algebraic (10) can be written as(11)φx1−x1φx2−x1…φxIxc−x1φx1−x2φx2−x2…φxIxc−x2........φx1−xIxcφx2−xIxc…φxIxc−xIxcα1α2..αIxc=ℓφxc−x1ℓφxc−x1..ℓφxc−Ixc.Solve the (11) and substitute in the (7) to the approximation of the differential operator ℓ. at the point xc.To simulate the unsteady SPMs model (1), we need to approximate the 1st and 2nd order derivatives. For this aim, we will divide the M scattered points Χ=x1,x2,…,xM into MB boundary points and MI interior points and M=MB+MI. Then, the weights associated to differential operators ∂/∂x,∂2/∂x2 in node xi based on stencil Ixi for internal points, will be determined by using (11) as follows:(12)wxxi=αi11wi1+αi21wi2+⋯+αiM1wiM,wxxxi=αi12wi1+αi22wi2+⋯+αiM2wiM.The Sparsity pattern of the matricesα1 and α2 is plotted in Figure 1.Figure 1
Sparsity pattern of matricesα1 (a) and α2 (b) for N = 100 and ns=15.
(a)(b)


## 2.1. RBFs

For scattered pointsΧ∈Ξ⊂ℜn, the approximation of a function χX can be written as a linear combination of RBFs [26, 27].(4)χΧ≅∑j=1MλjφΧ−Χj2+ψΧ,where M is the number of centers or sometimes called knotsΧ, Χ=x1,x2,…,xn,n is the dimension of the problem, λ′s are coefficients to be determined and φ is a RBF. (4) can be written without the additional polynomial ψΧ. If Ψqn denotes the space of n-variate polynomials of order not exceeding q, and letting the polynomials P1,P2,…,Pm be the basis of Ψqn in ℜn, then the polynomial ψΧ, in (4), is usually written in the following form:(5)ψX=∑i=1mζiPiX,where m=q−1+n!/n!q−1!.To determine the coefficientsλ1,λ2,…,λN and ζ1,ζ2,…,ζm extra m equations are required in addition to the M equations resulting from the collocating equation (4) at M knots. This is insured by the m conditions for (4), viz.,(6)∑j=1MλjPiXj=0,i=1,2,…,m.There are many RBFs (expression ofφ) available in literature such as Multiquadric (MQ), Inverse Multiquadric (IMQ), Gaussian, Thin Plate spline [27–34] etc.


## 2.2. LRBF-FD Method Formulation

In this work, we used LRBF-FD approximation for the discretization of spatial derivatives. LRBF-FD approximation is more general form of classical finite difference method. To explain the concept of LRBF-FD approximation, assume that there areM scattered points in the bounded domain of the problem i.e. Ξ⊂ℜn. In this approximation, the differential operator ℓ. is approximated at a particular point xc∈Ξ as a weighted linear combination of function values at closest ns neighbor stencil of the point xc as follows:(7)ℓwxc≅∑j∈Ixcαjwj,where Ixc=j∈1,2,…,M:xj−xc2≤δ and δ is the radius of support domain.Now, the main task is to calculate the unknown coefficientsαj occurred in (7). For this purpose, assume the approximation of the unknown solution w. is as follows [35, 36].(8)wx=∑k∈Ixβkφx−xk2,x∈Ξ,where φ. is MQ-RBF discussed in above subsection.Apply the differential operatorℓ. over (8) and take x=xc, we get(9)ℓwxc=∑k∈Ixβkℓφxc−xk2,xc∈Ξ.Now, using the (7)–(9), we can write as follows:(10)∑k∈Ixβkℓφxc−xk2=∑j∈Ixcαjwj=∑j∈Ixcαj∑k∈Ixβkφxj−xk2.To compute the unknown weightsαj, the system algebraic (10) can be written as(11)φx1−x1φx2−x1…φxIxc−x1φx1−x2φx2−x2…φxIxc−x2........φx1−xIxcφx2−xIxc…φxIxc−xIxcα1α2..αIxc=ℓφxc−x1ℓφxc−x1..ℓφxc−Ixc.Solve the (11) and substitute in the (7) to the approximation of the differential operator ℓ. at the point xc.To simulate the unsteady SPMs model (1), we need to approximate the 1st and 2nd order derivatives. For this aim, we will divide the M scattered points Χ=x1,x2,…,xM into MB boundary points and MI interior points and M=MB+MI. Then, the weights associated to differential operators ∂/∂x,∂2/∂x2 in node xi based on stencil Ixi for internal points, will be determined by using (11) as follows:(12)wxxi=αi11wi1+αi21wi2+⋯+αiM1wiM,wxxxi=αi12wi1+αi22wi2+⋯+αiM2wiM.The Sparsity pattern of the matricesα1 and α2 is plotted in Figure 1.Figure 1
Sparsity pattern of matricesα1 (a) and α2 (b) for N = 100 and ns=15.
(a)(b)


## 3. Mesh-free Numerical Scheme Based on LRBF-FD Approximation for the Parabolic Singularly Perturbed Models

Herein, we will develop a mesh-free scheme based on forward finite difference (FFD) and LRBF-FD for the simulation of parabolic SPMs (1). The steps of the scheme are as follows


### 3.1. Semidiscretization

In this scheme, the first step is semidiscretization of model (1) in time space. For this purpose, define tn=nτ,n=0,1,…,K, where τ=T/K is the time step. Now, approximate time derivative by FFD and apply Crank–Nicolson (CN) on other part of the model, then we have(13)εwxxn+1+wxxn2−cn+1wn+1+wn2+μan+1wxn+1+wxn2−bn+1wn+1−wnτ=fx,tn+1+sn,0≤n≤K−1,where wn+1=wx,tn+1,an+1=ax,tn+1,bn+1=bx,tn+1,cn+1=cx,tn+1, and sn is the truncation error (T. E.) in semidiscrete problem.With initial and BCs(14)w0=w0x,wn+1=0.The system (13) can written as follows:(15)ετwxxn+1+τμan+1wxn+1−2bn+1+τcn+1wn+1=−ετwxxn−τμan+1wxn+τcn+1−2bn+1wn+2τfx,tn+1+2τsn.Equation (15) is the system of linear ordinary differential equations (ODEs) for each value of n and if we drop out the small T. E. sn, we have(16)ετwxxn+1+τμan+1wxn+1−2bn+1+τcn+1wn+1=−ετwxxn−τμan+1wxn+τcn+1−2bn+1wn+2τfx,tn+1.


### 3.2. Fully Discretization

For fully discretization, we used LRBF-FD approximation. First of all discretize the given domain by the uniform nodexj,j=1,2,…,M. Then, discretize the equation (16) at x=xi and the spatial derivatives are approximated as follows:(17)wxn+1xi≅αi11wi1+αi21wi2+⋯+αiM1wiM,wxxn+1xi≅αi12wi1+αi22wi2+⋯+αiM2wiM.After that use the above approximation in (16), we have the following system of linear equations:(18)ετ∑j=1Mαij2wjn+1+τμan+1∑j=1Mαij1wjn+1−2bn+1+τcn+1win+1=−ετ∑j=1Mαij2wjn−τμan+1∑j=1Mαij1wjn+τcn+1−2bn+1win+2τfx,tn+1,j=1,2,…,M.The above system can be written in matrix form as(19)AM×MWM×1=BM×1,where the elements of the matrices are(20)Aij=εταij2+τμan+1αij1,i≠j,εταij2+τμan+1αij1−2bn+1+τcn+1,i=j,Bi=−ετ∑j=1Mαij2wjn−τμan+1∑j=1Mαij1wjn+τcn+1−2bn+1win+2τfx,tn+1.Now, apply the BCs and solved the obtained linear system of orderM−2×M−2 by the LU decomposition method.For simulation purpose, we used the uniform nodes instead of nonuniform nodes since uniform nodes gives better results than nonuniform nodes for the proposed algorithm.


## 3.1. Semidiscretization

In this scheme, the first step is semidiscretization of model (1) in time space. For this purpose, define tn=nτ,n=0,1,…,K, where τ=T/K is the time step. Now, approximate time derivative by FFD and apply Crank–Nicolson (CN) on other part of the model, then we have(13)εwxxn+1+wxxn2−cn+1wn+1+wn2+μan+1wxn+1+wxn2−bn+1wn+1−wnτ=fx,tn+1+sn,0≤n≤K−1,where wn+1=wx,tn+1,an+1=ax,tn+1,bn+1=bx,tn+1,cn+1=cx,tn+1, and sn is the truncation error (T. E.) in semidiscrete problem.With initial and BCs(14)w0=w0x,wn+1=0.The system (13) can written as follows:(15)ετwxxn+1+τμan+1wxn+1−2bn+1+τcn+1wn+1=−ετwxxn−τμan+1wxn+τcn+1−2bn+1wn+2τfx,tn+1+2τsn.Equation (15) is the system of linear ordinary differential equations (ODEs) for each value of n and if we drop out the small T. E. sn, we have(16)ετwxxn+1+τμan+1wxn+1−2bn+1+τcn+1wn+1=−ετwxxn−τμan+1wxn+τcn+1−2bn+1wn+2τfx,tn+1.


## 3.2. Fully Discretization

For fully discretization, we used LRBF-FD approximation. First of all discretize the given domain by the uniform nodexj,j=1,2,…,M. Then, discretize the equation (16) at x=xi and the spatial derivatives are approximated as follows:(17)wxn+1xi≅αi11wi1+αi21wi2+⋯+αiM1wiM,wxxn+1xi≅αi12wi1+αi22wi2+⋯+αiM2wiM.After that use the above approximation in (16), we have the following system of linear equations:(18)ετ∑j=1Mαij2wjn+1+τμan+1∑j=1Mαij1wjn+1−2bn+1+τcn+1win+1=−ετ∑j=1Mαij2wjn−τμan+1∑j=1Mαij1wjn+τcn+1−2bn+1win+2τfx,tn+1,j=1,2,…,M.The above system can be written in matrix form as(19)AM×MWM×1=BM×1,where the elements of the matrices are(20)Aij=εταij2+τμan+1αij1,i≠j,εταij2+τμan+1αij1−2bn+1+τcn+1,i=j,Bi=−ετ∑j=1Mαij2wjn−τμan+1∑j=1Mαij1wjn+τcn+1−2bn+1win+2τfx,tn+1.Now, apply the BCs and solved the obtained linear system of orderM−2×M−2 by the LU decomposition method.For simulation purpose, we used the uniform nodes instead of nonuniform nodes since uniform nodes gives better results than nonuniform nodes for the proposed algorithm.


## 4. MATLAB Implementation Steps of the Developed Scheme

The implementation steps of the scheme in MATLAB are as follows:A1. Partition the spatial domain (0, 1) into grid points with uniform step length and the grid points arexj=j−1/N−1,j=1,2,...,N and enter the constant shape parameterA2. Compute the local first and second order approximation named as sparse differentiation matricesα1 and α2A3. Enter the linear system (19) to get the solution at grid pointA4. Solve the entered system in MATLAB via LU decomposition method and iterates uptoT (final time)


## 5. Numerical Simulation and Discussion

This section ponders some particular SPMs of (1) for simulation purpose. To check the chastity and competence of the developed numerical scheme, maximum absolute error L∞ and relative error L2 are computed by the formulas:(21)L∞=uanalytic−uappx.∞=maxjujanalytic−ujappx.,L2=Δx∑k=1Nukanalytic−ukappx.2,where uanalytic and uappx. denote analytic and approximated solutions respectively.The simulation work has been done in the MATLAB environment on desktop computer (Intel(R) Core(TM) i3-6100T CPU, 4 GB RAM).The shape parameter used in RBFs is calculated by the following algorithm [37]. (Algorithm 1)Algorithm 1:Local RBF-FD scheme algorithm.
Input kmin, kmax, incrementk = 1;while (k < kmin or k > kmax). produce interpolation matrix say A; use [U, S,V] = svd(A); K = max(S)/min(S)ifK < minKshape = shape − increment;elseifK > maxKshape = shape + increment;endend


### 5.1. SP Model 1 (Reaction-Diffusion, Forμ=0)

Consider the one parameter SP reaction-diffusion model aswt−εwxx+1+xe−tw=fx,t,wx,0=0,w0,t=0,w1,t=0, over the domain x,t∈0,t×0,1 which is particular case of the model (1). The function fx,t is taken from exact solution given by [38].(22)wx,t=te−x/ε+e−1−x/ε1+e−1/ε−cos2πx.The results of the model are presented in Tables1 and 2 and Figure 2. The Table 1 submits the maximum absolute L∞ and L2 errors at different time and parameter up to ε=2−35. The Table concludes that as we decrease the parameter the errors decrease. For ε=2−30 at time t=1.0, the present L∞ error with M=150, ε=2−30 is 6.660E − 06 while the same error with N=256, 6.72E − 06 (hybrid method) and 2.80E − 04 (backward Euler method in time and central difference approximation in space) [38]. This concludes that the proposed scheme produces good results with small number of grid points. Table 2 reports the CPU time at different time and perturbation parameter ε. Figure 2 depicts the physical behavior of the model at different times with parameter ε=2−10. The figure shows that, as time increases, the solution decreases and comes closer and closer near x = 0.5.Table 1
Maximum absoluteL∞ and L2 errors at different times and parameters with M=150,ns=3.
TimeErrorsε=2−10ε=2−15ε=2−20ε=2−25ε=2−30ε=2−350.1L∞2.809E − 041.869E − 049.654E − 052.411E − 061.466E − 061.469E − 06L21.573E − 051.919E − 069.011E − 076.007E − 085.797E − 085.814E − 080.3L∞2.300E − 031.500E − 038.173E − 042.411E − 053.545E − 063.547E − 06L21.281E − 041.525E − 057.483E − 062.625E − 071.446E − 071.458E − 070.5L∞5.800E − 033.700E − 032.100E − 032.491E − 054.891E − 064.897E − 06L23.264E − 043.788E − 051.926E − 052.625E − 072.0502E − 072.078E − 070.7L∞1.060E − 036.400E − 033.900E − 031.261E − 045.795E − 065.805E − 06L25.927E − 046.725E − 053.525E − 056.330E − 072.487E − 072.530E − 071.0L∞1.980E − 021.130E − 027.300E − 032.386E − 046.660E − 066.678E − 06L21.100E − 031.205E − 046.556E − 052.146E − 062.954E − 073.006E − 07Table 2
CPU time at different time andε.
Corresponding CPU timeTimeε=2−10ε=2−15ε=2−20ε=2−25ε=2−30ε=2−350.10.4760.4350.4590.4100.3710.4200.51.571.341.431.331.341.541.02.672.782.732.792.562.91Figure 2
Solution behavior of the mode at different time with parameterε=2−10 and M=150.


### 5.2. SP Model 2 (Reaction-Diffusion, forμ=0)

Consider the one parameter SP reaction-diffusion model aswt−εwxx+1+xe−tw=fx,t,wx,0=0,w0,t=0,w1,t=0, over the domain x,t∈0,t×0,1 which is same as the above example except the known function fx,t. The function fx,t is taken from exact solution given by [38].(23)wx,t=1−e−te−x/ε+e−1−x/ε1+e−1/ε−cos2πx.Tables3 and 4 and Figure 3 offer the results of the model. Table 3 submits the maximum absolute L∞ and L2 errors at different time and parameter up to ε=2−35. The table concludes that, as we decrease the parameter the errors decrease. For ε=2−30, at time t=1.0, the present L∞ error with M=150, ε=2−30 is 3.854E − 05, while the same error with N=256, 5.67E − 03 (hybrid method) [38]. This concludes that the proposed scheme produces good results with small number of grid points. Table 4 reports the CPU time at different time and perturbation parameter ε. Figure 3 depicts the physical behavior in 2D and 3D of the model at different times with parameter ε=2−20. The figure shows that, as time increases the solution decreases and comes closer and closer between x = 0.4 to x = 0.6. Second, the solution becomes steeper near the right and left boundaries.Table 3
Maximum absoluteL∞ and L2 errors at different times and parameters with M=150,ns=3.
TimeErrorsε=2−10ε=2−15ε=2−20ε=2−25ε=2−30ε=2−350.1L∞9.634E − 051.007E − 047.792E − 051.417E − 051.395E − 051.340E − 05L25.337E − 061.120E − 068.957E − 075.493E − 075.489E − 075.430E − 070.3L∞6.895E − 049.171E − 046.553E − 043.210E − 053.039E − 053.043E − 05L23.894E − 058.625E − 066.103E − 061.232E − 061.248E − 061.249E − 060.5L∞1.600E − 031.800E − 031.600E − 034.166E − 053.763E − 053.756E − 05L29.198E − 051.666E − 051.468E − 051.591E − 061.594E − 061.597E − 060.7L∞2.800E − 032.900E − 032.800E − 036.528E − 053.986E − 053.966E − 05L21.555E − 042.711E − 052.518E − 051.802E − 061.734E − 061.738E − 061.0L∞4.700E − 034.600E − 034.800E − 031.292E − 043.854E − 053.818E − 05L22.587E − 044.282E − 054.248E − 052.055E − 061.729E − 061.736E − 06Table 4
CPU time at different time andε.
Corresponding CPU timeTimeε=2−10ε=2−15ε=2−20ε=2−25ε=2−30ε=2−350.10.4560.4400.4480.4080.3670.4180.51.321.321.341.281.291.521.02.562.672.652.672.472.88Figure 3
Solution behavior of the mode at different time with parameterε=2−20 and M=150.


### 5.3. SP Model 3

Consider the two parameter SP reaction-diffusion model with initial and BCs as [39, 40].(24)εwxxx,t+μ1+xwxx,t−wtx,t−wx,t=16x21−x2,and wx,0=0,w0,t=0,w1,t=0 over the domain x,t∈0,t×0,T. The exact solution of the model is not given in literature but the solution behaviour is given in [39, 40] up to time t=1.0. Here, we present numerical results in 2D and 3D form of the model in Figure 4 with parameter ε=10−4,μ=10−10, M=150. The figure exhibits that as time increases up to t=1.0 the solution decreases and goes upto −0.6316. Second, the solution comes closer and closer between the right and left boundaries.Figure 4
Solution behavior of the mode at different time with parameterε=10−4,μ=10−10 and M=150.


### 5.4. SP Model 4

Consider the two parameter SP reaction-diffusion model with initial and BCs as [40, 41].(25)εwxxx,t+μ1+x1−x+t2wxx,t−wtx,t−1+5xtwx,t=x1−xet−1,and wx,0=0,w0,t=0,w1,t=0 over the domain x,t∈0,t×0,T. The exact solution of the model is not given in literature but the solution behaviour is given in [40, 41] up to time t=1.0. Here, we present numerical results in 2D and 3D form of the model in Figure 5 with parameter ε=10−4,μ=10−10, M=150. The figure exhibits that, as time increases up to t=2.0, the solution decreases and goes upto −0.274287. Second, the solution comes closer and closer between the right and left boundaries. Since for the above model, we do not know the analytical solution; therefore to explain the performance of the proposed scheme, we depicted the solution behavior of the model.Figure 5
Solution behavior of the mode at different time with parameterε=10−4,μ=10−10 and M=150.


## 5.1. SP Model 1 (Reaction-Diffusion, Forμ=0)

Consider the one parameter SP reaction-diffusion model aswt−εwxx+1+xe−tw=fx,t,wx,0=0,w0,t=0,w1,t=0, over the domain x,t∈0,t×0,1 which is particular case of the model (1). The function fx,t is taken from exact solution given by [38].(22)wx,t=te−x/ε+e−1−x/ε1+e−1/ε−cos2πx.The results of the model are presented in Tables1 and 2 and Figure 2. The Table 1 submits the maximum absolute L∞ and L2 errors at different time and parameter up to ε=2−35. The Table concludes that as we decrease the parameter the errors decrease. For ε=2−30 at time t=1.0, the present L∞ error with M=150, ε=2−30 is 6.660E − 06 while the same error with N=256, 6.72E − 06 (hybrid method) and 2.80E − 04 (backward Euler method in time and central difference approximation in space) [38]. This concludes that the proposed scheme produces good results with small number of grid points. Table 2 reports the CPU time at different time and perturbation parameter ε. Figure 2 depicts the physical behavior of the model at different times with parameter ε=2−10. The figure shows that, as time increases, the solution decreases and comes closer and closer near x = 0.5.Table 1
Maximum absoluteL∞ and L2 errors at different times and parameters with M=150,ns=3.
TimeErrorsε=2−10ε=2−15ε=2−20ε=2−25ε=2−30ε=2−350.1L∞2.809E − 041.869E − 049.654E − 052.411E − 061.466E − 061.469E − 06L21.573E − 051.919E − 069.011E − 076.007E − 085.797E − 085.814E − 080.3L∞2.300E − 031.500E − 038.173E − 042.411E − 053.545E − 063.547E − 06L21.281E − 041.525E − 057.483E − 062.625E − 071.446E − 071.458E − 070.5L∞5.800E − 033.700E − 032.100E − 032.491E − 054.891E − 064.897E − 06L23.264E − 043.788E − 051.926E − 052.625E − 072.0502E − 072.078E − 070.7L∞1.060E − 036.400E − 033.900E − 031.261E − 045.795E − 065.805E − 06L25.927E − 046.725E − 053.525E − 056.330E − 072.487E − 072.530E − 071.0L∞1.980E − 021.130E − 027.300E − 032.386E − 046.660E − 066.678E − 06L21.100E − 031.205E − 046.556E − 052.146E − 062.954E − 073.006E − 07Table 2
CPU time at different time andε.
Corresponding CPU timeTimeε=2−10ε=2−15ε=2−20ε=2−25ε=2−30ε=2−350.10.4760.4350.4590.4100.3710.4200.51.571.341.431.331.341.541.02.672.782.732.792.562.91Figure 2
Solution behavior of the mode at different time with parameterε=2−10 and M=150.


## 5.2. SP Model 2 (Reaction-Diffusion, forμ=0)

Consider the one parameter SP reaction-diffusion model aswt−εwxx+1+xe−tw=fx,t,wx,0=0,w0,t=0,w1,t=0, over the domain x,t∈0,t×0,1 which is same as the above example except the known function fx,t. The function fx,t is taken from exact solution given by [38].(23)wx,t=1−e−te−x/ε+e−1−x/ε1+e−1/ε−cos2πx.Tables3 and 4 and Figure 3 offer the results of the model. Table 3 submits the maximum absolute L∞ and L2 errors at different time and parameter up to ε=2−35. The table concludes that, as we decrease the parameter the errors decrease. For ε=2−30, at time t=1.0, the present L∞ error with M=150, ε=2−30 is 3.854E − 05, while the same error with N=256, 5.67E − 03 (hybrid method) [38]. This concludes that the proposed scheme produces good results with small number of grid points. Table 4 reports the CPU time at different time and perturbation parameter ε. Figure 3 depicts the physical behavior in 2D and 3D of the model at different times with parameter ε=2−20. The figure shows that, as time increases the solution decreases and comes closer and closer between x = 0.4 to x = 0.6. Second, the solution becomes steeper near the right and left boundaries.Table 3
Maximum absoluteL∞ and L2 errors at different times and parameters with M=150,ns=3.
TimeErrorsε=2−10ε=2−15ε=2−20ε=2−25ε=2−30ε=2−350.1L∞9.634E − 051.007E − 047.792E − 051.417E − 051.395E − 051.340E − 05L25.337E − 061.120E − 068.957E − 075.493E − 075.489E − 075.430E − 070.3L∞6.895E − 049.171E − 046.553E − 043.210E − 053.039E − 053.043E − 05L23.894E − 058.625E − 066.103E − 061.232E − 061.248E − 061.249E − 060.5L∞1.600E − 031.800E − 031.600E − 034.166E − 053.763E − 053.756E − 05L29.198E − 051.666E − 051.468E − 051.591E − 061.594E − 061.597E − 060.7L∞2.800E − 032.900E − 032.800E − 036.528E − 053.986E − 053.966E − 05L21.555E − 042.711E − 052.518E − 051.802E − 061.734E − 061.738E − 061.0L∞4.700E − 034.600E − 034.800E − 031.292E − 043.854E − 053.818E − 05L22.587E − 044.282E − 054.248E − 052.055E − 061.729E − 061.736E − 06Table 4
CPU time at different time andε.
Corresponding CPU timeTimeε=2−10ε=2−15ε=2−20ε=2−25ε=2−30ε=2−350.10.4560.4400.4480.4080.3670.4180.51.321.321.341.281.291.521.02.562.672.652.672.472.88Figure 3
Solution behavior of the mode at different time with parameterε=2−20 and M=150.


## 5.3. SP Model 3

Consider the two parameter SP reaction-diffusion model with initial and BCs as [39, 40].(24)εwxxx,t+μ1+xwxx,t−wtx,t−wx,t=16x21−x2,and wx,0=0,w0,t=0,w1,t=0 over the domain x,t∈0,t×0,T. The exact solution of the model is not given in literature but the solution behaviour is given in [39, 40] up to time t=1.0. Here, we present numerical results in 2D and 3D form of the model in Figure 4 with parameter ε=10−4,μ=10−10, M=150. The figure exhibits that as time increases up to t=1.0 the solution decreases and goes upto −0.6316. Second, the solution comes closer and closer between the right and left boundaries.Figure 4
Solution behavior of the mode at different time with parameterε=10−4,μ=10−10 and M=150.


## 5.4. SP Model 4

Consider the two parameter SP reaction-diffusion model with initial and BCs as [40, 41].(25)εwxxx,t+μ1+x1−x+t2wxx,t−wtx,t−1+5xtwx,t=x1−xet−1,and wx,0=0,w0,t=0,w1,t=0 over the domain x,t∈0,t×0,T. The exact solution of the model is not given in literature but the solution behaviour is given in [40, 41] up to time t=1.0. Here, we present numerical results in 2D and 3D form of the model in Figure 5 with parameter ε=10−4,μ=10−10, M=150. The figure exhibits that, as time increases up to t=2.0, the solution decreases and goes upto −0.274287. Second, the solution comes closer and closer between the right and left boundaries. Since for the above model, we do not know the analytical solution; therefore to explain the performance of the proposed scheme, we depicted the solution behavior of the model.Figure 5
Solution behavior of the mode at different time with parameterε=10−4,μ=10−10 and M=150.


## 6. Concluding Remarks

In this article, the authors developed a mesh-free numerical scheme based on FFD and LRBF-FD for the simulation of unsteady-state 1D parabolic singularly perturbed CDR models with two parameters and variable coefficients. The novelty of the work is that we solved unsteady-state singularly perturbed problems with two parameters and variable coefficients. The proposed numerical scheme is tested by taking several numerical examples and the results found satisfactory. The proposed mesh-free scheme is new for the 1D parabolic singularly perturbed CDR models and provides similar results as given in [38–41]. The developed scheme easily worked for small perturbed parameter ε=10−35 and variable coefficients and produced comparable results with [38–41]. The main advantage of the scheme is that it is easy to implement and an attractive alternative to the FDM and FEM which are classical solution schemes for the higher dimensional nonlinear PDEs.



---

*Source: 3119482-2022-03-15.xml*