Robust numerical method for singularly perturbed differential equations having both large and small delay

Habtamu Garoma Debela

Arab Journal of Mathematical Sciences

ISSN: 1319-5166

Open Access. Article publication date: 30 March 2021

Issue publication date: 11 January 2022

1311

Abstract

Purpose

The purpose of this study is to develop stable, convergent and accurate numerical method for solving singularly perturbed differential equations having both small and large delay.

Design/methodology/approach

This study introduces a fitted nonpolynomial spline method for singularly perturbed differential equations having both small and large delay. The numerical scheme is developed on uniform mesh using fitted operator in the given differential equation.

Findings

The stability of the developed numerical method is established and its uniform convergence is proved. To validate the applicability of the method, one model problem is considered for numerical experimentation for different values of the perturbation parameter and mesh points.

Originality/value

In this paper, the authors consider a new governing problem having both small delay on convection term and large delay. As far as the researchers' knowledge is considered numerical solution of singularly perturbed boundary value problem containing both small delay and large delay is first being considered.

Keywords

Citation

Garoma Debela, H. (2022), "Robust numerical method for singularly perturbed differential equations having both large and small delay", Arab Journal of Mathematical Sciences, Vol. 28 No. 1, pp. 87-99. https://doi.org/10.1108/AJMS-09-2020-0058

Publisher

:

Emerald Publishing Limited

Copyright © 2020, Habtamu Garoma Debela

License

Published in Arab Journal of Mathematical Sciences. Published by Emerald Publishing Limited. This article is published under the Creative Commons Attribution (CC BY 4.0) licence. Anyone may reproduce, distribute, translate and create derivative works of this article (for both commercial and non-commercial purposes), subject to full attribution to the original publication and authors. The full terms of this licence may be seen at http://creativecommons.org/licences/by/4.0/legalcode


1. Introduction

A differential equation is said to be singularly perturbed delay differential equation, if it includes at least one delay term, involving unknown functions occurring with different arguments, and also, the highest derivative term is multiplied by a small parameter. Such type of delay, differential equations play a very important role in the mathematical models of science and engineering, such as the human pupil light reflex with mixed delay type [], variational problems in control theory with small state problem [], models of HIV infection [] and signal transition [].

Any system involving a feedback control almost involves time delay. The delay occurs because a finite time is required to sense the information and then react to it. Finding the solution of singularly perturbed delay differential equations, whose application mentioned above, is a challenging problem. In response to these, in recent years, there has been a growing interest in numerical methods on singularly perturbed delay differential equations. The authors of [] have developed various numerical schemes on uniform meshes for singularly perturbed second-order differential equations having small delay on convection term. The authors of [] have presented second-order differential equations having large delay.

In this paper, we consider a new governing problem having both small delay on convection term and large delay. Additionally, in recent years the correlative physical phenomena in-depth of the problem under consideration have been done by the authors []. As far as the researchers' knowledge is considered numerical solution of singularly perturbed boundary value problem containing both small delay and large delay is first being considered. Thus, the purpose of this study is to develop stable, convergent and accurate numerical method for solving singularly perturbed differential-difference equations having both small and large delay.

Throughout our analysis C is generic positive constant that are independent of the parameter ε and number of mesh points 2N. We assume that Ω¯=[0,2], Ω=(0,2), Ω1=(0,1), Ω2=(1,2) and Ω*=Ω1Ω2. L1 and L2 are the linear operator associated to the domain Ω1 and Ω2, respectively.

2. Statement of the problem

Consider the following singularly perturbed problem

(1)Ly(x)=εy(x)+a(x)y(x)+b(x)y(x)+c(x)y(x1)+d(x)y(xδ)=f(x),
(2)y(x)=φ(x),x[1,0],y(2)=l,
where δ is small, that is δ=O(ε), 0<ε1, φ(x) is sufficiently smooth on [1,0]. For all xΩ, it is assumed that the sufficient smooth functions a(x),b(x), c(x) and d(x) satisfy at a(x)a1>a>0,b(x)>b0,c(x)γ<0,d(x)ζ0, and 2(a+ζ)+5b+5γ>η>0,a(a1a)+2γ>0. The above assumptions ensure that yX=C0(Ω¯)C1(Ω)C2(Ω1Ω2).

The boundary value problem 1–2 exhibits strong boundary layer at x=2 and interior layer at x=1. Expand y(xδ) about x using the Taylor's expansion and discard higher order terms. Then, the above problem can be approximated by

(3)cε,δ(x)y(x)+p(x)y(x)+b(x)y(x)+c(x)y(x1)=f(x),
where cε,δ(x)=εδd(x) and p(x)=a(x)+d(x),
(4)y(x)=φ(x),x[1,0],y(2)=l.

As we observed from , the values of y(x1) is known for the domain Ω1 and unknown for the domain Ω2 due to the large delay at x=1. So, it is impossible to treat the problem throughout the domain (Ω¯). Thus, we have to treat the problem at Ω1 and Ω2 separately. are equivalent to

(5)Ly(x)=R(x),
where
(6)Ly(x)={L1y(x)=cε,δ(x)y(x)+p(x)y(x)+b(x)y(x),xΩ1,L2y(x)=cε,δ(x)y(x)+p(x)y(x)+b(x)y(x)+c(x)y(x1),xΩ2.
(7)R(x)={f(x)c(x)φ(x1),xΩ1,f(x),xΩ2.

with boundary conditions

(8){y(x)=φ(x),x[1,0],y(1)=y(1+),y(1)=y(1+),y(2)=l.

3. Properties of continuous solution

Lemma 3.1.

(Maximum principle) Let ψ(x) be any function in X such that ψ(0)0,ψ(2)0,L1ψ(x)0,xΩ1,L2ψ(x)0,xΩ2 and [ψ](1)0 then ψ(x)0,xΩ¯.

Proof. For the proof refer [] ▪

Lemma 3.2.

(Stability result) The solution y(x) of problem , satisfies the bound

|y(x)|Cmax{|y(0)|,|y(2)|,supxΩ*|Ly(x)|},xΩ¯.

Proof. For the proof refer [] ▪
Lemma 3.3.

Let y(x) be the solution of . Then we have the following bounds

||y(k)(x)||Ω*Cεk,k=1,2,3.

Proof. For the proof refer [] ▪

4. Numerical scheme formulation

We divide the interval [0,2] into 2N equal parts with constant mesh length h. Let 0=x0,x1,,xN=1,xN+1,xN+2,,x2N=2 be the mesh points. Then we have xi=ih,i=0,1,2,,2N.

Consider a uniform mesh with interval [0,1] in which 0=x0<x1<<xN=1 where h=1N and xi=ih,i=0,1,2,,N.

We can rewrite as

(9)cε,δ(x)y(x)+p(x)y(x)+b(x)y(x)=Q(x),xΩ1,
where Q(x)=f(x)c(x)y(x1).

For each segment [xi,xi+1],i=1,2,,N1 the non-polynomial cubic spline SΔ(x) has the following form

(10)SΔ(x)=ai+bi(xxi)+ci(ew(xxi)+ew(xxi))+di(ew(xxi)ew(xxi)),
where ai,bi,ci and di are unknown coefficients, and w0 arbitrary parameter which will be used to increase the accuracy of the method.

To determine the unknown coefficients in in terms of yi,yi+1,Mi and Mi+1 first we define

(11){SΔ(xi)=yi,SΔ(xi+1)=yi+1,SΔ(xi)=Mi,SΔ(xi+1)=Mi+1.

The coefficients in are determined as

(12){ai=yiMiw2,bi=yi+1yih+MiMi+1wθ,ci=Mi+1w2(eθeθ)Mi(eθ+eθ)2w2(eθeθ),di=Mi2w2,
where θ=wh.

Using the continuity condition of the first derivative at xi, SΔ1(xi)=SΔ(xi), we have

(13)bi1+wci1(eθ+eθ)+wdi1(eθeθ)=bi+2wci.

Reducing indices of by one and substituting into , we obtain yiyi1h+MiMi+1wθ+w(2Mi(eθ+eθ)Mi12w2(eθ+eθ))=yi+1yih+MiMi+1wθ+2w(Mi+1w2(eθeθ)Mi(eθ+eθ)2w2(eθeθ)),

(14)yi12yi+yi+1h2=αMi1+2βMi+αMi+1,
where
α=1θ2(12θ(eθeθ)),β=1θ2(θ(eθ+eθ)(eθeθ)1).

If h0, then θ=hk0. Thus, using L'Hopitals rule we have limh0α=16 and limh0β=13.

Using SΔ(xi)=yi=Mi in to (9), we get

(15){cε,δ(x)Mi=Qipiyibiyi,cε,δ(x)Mi1=Qi1pi1yi1bi1yi1,cε,δ(x)Mi+1=Qi+1pi+1yi+1bi+1yi+1.

Using Taylorâ™s series expansions of yi1,yi+1,yi1 and yi+1 simplifying, we have

(16){yi=yi+1yi12h+T1,yi1=yi+1+4yi3yi12h+T2,yi+1=3yi+14yi+yi12h+T2,
where

T1=h26y(ξ) and T2=h212y(ξ), for ξ(xi1,xi).

Substituting in to , we obtain

(17){Mi=1cε,δ(x){Qipi(yi+1yi12h+T1)biyi},Mi1=1cε,δ(x){Qi1pi1(yi+1+4yi3yi12h+T2)bi1yi1},Mi+1=1cε,δ(x){Qi+1pi+1(3yi+14yi+yi12h+T2)bi+1yi+1}.

Substituting into and rearranging, we get

(18)cε,δ(x)h2(yi12yi+yi+1)+αpi12h(yi+14yi3yi1)+2βpi2h(yi+1yi1)+αpi+12h(3yi+13yi+yi1)=α(Qi1bi1yi1+Qi+1bi+1yi+1)+2β(Qibiyi)+T,
where, T=(4βpiαpi1αpi+1)h212y(ξ) is the local truncation error.

From the theory of singular perturbations described in [] and the Taylorâ™s series expansion of y(x) about the point ‘0’ in the asymptotic solution of the problem in , we have

y(xi)y0(xi)+(φ0y0(0))ep(0)ihcε,δ(x),
and letting ρ=hcε,δ(x), we get
limh0y(ih)y0(ih)+(φ0y0(0))ep(0)iρ,
since xi=x0+ih.

Introducing a fitting factor σ(ρ) in to , we get

(19)σ(ρ)cε,δ(x)h2(yi12yi+yi+1)+αpi12h(yi+14yi3yi1)+2βpi2h(yi+1yi1)+αpi+12h(3yi+13yi+yi1)=α(Qi1bi1yi1+Qi+1bi+1yi+1)+2β(Qibiyi)+T.

Multiplying by h and taking a limit as h0 we get

(20)σρlimh0(yi12yi+yi+1)+αp(0)hlimh0(yi+14yi3yi1)+βp(0)limh0(yi+1yi1)+αp(0)2limh0(3yi+13yi+yi1)=0.

Thus, we consider two cases of the boundary layers.

  • Case 1: For p(x)>0 (Left-end boundary layer), we have

(21){limh0(yi12yi+yi+1)=(φ0y0(0))ep(0)iρ(ep(0)ρ+ep(0)ρ2),limh0(yi+14yi3yi1)=(φ0y0(0))ep(0)iρ(3ep(0)ρep(0)ρ+4),limh0(yi+1yi1)=(φ0y0(0))ep(0)iρ(ep(0)ρ+3ep(0)ρ4),limh0(3yi+13yi+yi1)=(φ0y0(0))ep(0)iρ(ep(0)ρep(0)ρ).

Substituting into and simplifying, we get

σ0=ρp(0)(α+β)coth(p(0)ρ2).
  • Case 2: For p(x)<0 (Right-end boundary layer), we have

(22){limh0(yi12yi+yi+1)=(ϕy0(1))ep(1)iρ(ep(1)ρ+ep(1)ρ2),limh0(yi+14yi3yi1)=(ϕy0(1))ep(1)iρ(3ep(1)ρep(1)ρ+4),limh0(yi+1yi1)=(ϕy0(1))ep(1)iρ(ep(1)ρ+3ep(1)ρ4),limh0(3yi+13yi+yi1)=(ϕy0(1))ep(1)iρ(ep(1)ρep(1)ρ).

Substituting into and simplifying, we obtain

(23)σ1=ρp(1)(α+β)coth(p(1)ρ2).

In general, we take a variable fitting parameter as

(24)σi=ρip(xi)(α+β)coth(p(xi)ρi2),
where, ρi=hcε,δ(x).

Thus, can be written as

(25){cε,δ(x)σih23αpi12h+αbi1βpih+αpi+12h}yi1{2cε,δ(x)σih22αpi1h2βbi+2αpi+1h}yi+{cε,δ(x)σih2αpi12h+αbi+1+βpih+3αpi+12h}yi+1=α(Qi1+Qi+1)+2βQi.

Simplifying , for the domain Ω1=(0,1), we get the tri-diagonal system of the equation of the form

(26)LNEiyi1Fiyi+Giyi+1=Hi,i=1,2,,N1,
where
{Ei=cε,δ(x)σih23αpi12h+αbi1βpih+αpi+12h,Fi=2cε,δ(x)σih22αpi1h2βbi+2αpi+1h,Gi=cε,δ(x)σih2αpi12h+αbi+1+βpih+3αpi+12h,Hi=α(Qi1+Qi+1)+2βQi.

Similarly, if we consider the domain Ω2=(1,2), from , we have

(27){cε,δ(x)y(x)+p(x)y(x)+b(x)y(x)+c(x)y(x1)=f(x),xΩ2y(1)=θ,y(2)=l.

Using SΔ(xi)=yi=Mi for , we get

(28){cε,δ(xi)Mi=fipiyibiyiciy(xi1),cε,δ(xi)Mi1=fi1pi1yi1bi1yi1ci1y(xi11),cε,δ(xi)Mi+1=fi+1pi+1yi+1bi+1yi+1ci+1y(xi+11).

Substituting in to , we obtain

(29){Mi=1cε,δ(xi){fipi(yi+1yi12h+T1)biyiciy(xi1)},Mi1=1cε,δ(xi){fi1pi1(yi+1+4yi3yi12h+T2)bi1yi1ci1y(xi11)},Mi+1=1cε,δ(xi){fi+1pi+1(3yi+14yi+yi12h+T2)bi+1yi+1ci+1y(xi+11)}.

Substituting in to , introducing fitting factor and rearranging, we get

(30)LNEiyi1Fiyi+Giyi+1+Ti=Hi,i=N+1,N+2,,2N1,
where
{Ei=cε,δ(xi)σih23αpi12h+αbi1βpih+αpi+12h,Fi=2cε,δ(xi)σih22αpi1h2βbi+2αpi+1h,Gi=cε,δ(xi)σih2αpi12h+αbi+1+βpih+3αpi+12hHi=α(fi1+fi+1)+2βfi,Ti=α{ci1y(xi11)+ci+1y(xi+11)+2βciy(xi1)}.

Therefore, on the whole domain Ω¯=[0,2], the basic schemes to solve are the schemes given in and .

5. Stability and convergence analysis

5.1 Truncation error

Let expand the terms yi±1 and Mi±1 from , using Taylor's series as

(31){yi+1=yi+hyi+h22!yi+h33!yi+h44!yi(4)+h55!yi(5)+h66!yi(6)+O(h7),yi1=yihyi+h22!yih33!yi+h44!yi(4)h55!yi(5)+h66!yi(6)+O(h7),Mi+1=yi+1=yi+hyi+h22!yi(4)+h33!yi(5)+h44!yi(6)+O(h7),Mi1=yi1=yihyi+h22!yi(4)h33!yi(5)+h44!yi(6)+O(h7).

The local truncation error Ti(h) obtained from as

(32)Ti(h)=yi12yi+yi+1h2α(Mi1+Mi+1)2βMi.

Substituting the series of yi±1 and Mi±1 from and collecting like terms gives

(33)Ti(h)=(12(α+β))yi+h2(112α)yi(4)+O(h4).

But from the values of α=16 and β=13, becomes

Ti(h)=h2(112)yi(4)+O(h4),
which implies
(34)Ti(h)Ch2,
where C=112|yi(4)|.

This establishes that the developed method is second order accurate or its order of convergence is O(h2).

5.2 Convergence analysis

Local truncation errors refer to the differences between the original differential equation and its finite difference approximation at a mesh points. Finite difference scheme is called consistent if the limit of truncation error (Ti(h)) is equal to zero as the mesh size h goes to zero. Hence, the proposed method in with local truncation error in satisfies the definition of consistency as

(35)limh0Ti(h)=limh0Ch2=0.

Thus, the proposed scheme is consistent.

5.3 Stability analysis

Consider the developed scheme in ,

(36)Eiyi1Fiyi+Giyi+1=Hi,
where the coefficients Ei, Fi and Gi are as in . If we multiply both sides of by h2 and consider the values of Ei, Fi and Gi for sufficiently small h, we get
(37)Ei=Gi=εσ,Fi=2εσ,

Considering into the one which is multiplied by h2 the developed scheme can be written in a matrix form

(38)AY=B
where the matrices A=(2εσεσ00εσ2εσεσ00mtdmtd0εσ2εσ),Y=(y1y2yN2yN1) and B=(h2H1E1y0h2H2h2HN2h2HN1GN1yN).

Here, the coefficient matrix A is a tri-diagonal matrix with size (N1)×(N1). Matrix A is irreducible if its codiagonals contain nonzero elements only. The codiagonals contain Ei and Gi. It is clearly seen that, for sufficiently small h both Ei0 and Gi0 for i=1,2,,N1. Hence, A is irreducible.

Again we can see that all |Ei|, |F|i, |Gi|>0 for i=1,2,,N1 and in each row of A, the modulus of diagonal element is greater than or equal to the sum of the modulus of the two codiagonal elements (i.e. |Fi||Ei|+|Gi|). This implies that A is diagonally dominant. Under this condition, the Thomas algorithm is stable for sufficiently small h.

As discussed in [] the eigenvalues of a tri-diagonal matrix A are given by

λs=2εσ+2{(εσ)(εσ)}cossπN,s=1(1)N1.

Hence, the eigenvalues of matrix A in are

λs=2εσ+2{(εσ)2}cossπN=2εσ(1cossπN),s=1(1)N1.

But from trigonometric identity, we have 1cossπN=2sin2sπ2N. Thus, the eigenvalues of A

(39)λs=2εσ(2sin2sπ2N)=4εσsin2sπ2N4εσ.

A finite difference method for the boundary value problems is stable if A is nonsingular and A1C, for 0<h<h0 .where, C and h0 are two constants that are independent of h.

Since A is real and symmetric it follows that A1 is also real and symmetric so that, its eigenvalues are real and given by 1λs. Hence, as [] the stability condition of the method will be satisfied when; A1=|1λs|=|14εσ|=14εσC , where, C is independent of h . Thus the developed scheme in is stable. A consistent and stable finite difference method is convergent by []. Hence as we have shown above, the proposed method is satisfying both the criteria of consistency and stability which are equivalent to convergence of the method.

6. Numerical examples and results

In this section, one example is given to illustrate the numerical method discussed above. The exact solutions of the test problem is not known. Therefore, we use the double mesh principle to estimate the error and compute the experiment rate of convergence to the computed solution. For this we put

(40)EεN=max0i2N|YiNY2i2N|,
where YiN and Y2i2N are the ith components of the numerical solutions on meshes of N and 2N, respectively. We compute the uniform error and the rate of convergence as
(41)EN=maxεEεN,andRN=log2(ENE2N).

The numerical results are presented for the values of the perturbation parameter ε{104,108,,1020}.

Example 6.1.

Consider the model singularly perturbed boundary value problem:

εy(x)+10y(x)y(x1)+y(xε)=xx(0,1)(1,2),

Subject to the boundary conditions

y(x)=1,x[1,0],y(2)=2.

7. Discussion and conclusion

This study introduces a fitted nonpolynomial spline method for singularly perturbed differential equations having both small and large delay. The numerical scheme is developed on uniform mesh using fitted operator in the given differential equation. The stability of the developed numerical method is established and its uniform convergence is proved. To validate the applicability of the method, one model problem is considered for numerical experimentation for different values of the perturbation parameter and mesh points. The numerical results are tabulated in terms of maximum absolute errors, numerical rate of convergence and uniform errors (see ). Further, behavior of the numerical solution (), point-wise absolute error () and the ε -uniform convergence of the method is shown by the log-log plot (). The method is shown to be ε -uniformly convergent with order of convergence O(h2). The proposed method gives accurate, stable and ε -uniform numerical result.

Figures

The behavior of the numerical solution for Example 6.1 at ε=10−12 and N=32

Figure 1

The behavior of the numerical solution for at ε=1012 and N=32

Point wise absolute error of Example 6.1 at ε=10−12 with different mesh point N

Figure 2

Point wise absolute error of at ε=1012 with different mesh point N

ε-uniform convergence with fitted operator in log-log scale for Example 6.1

Figure 3

ε-uniform convergence with fitted operator in log-log scale for

Maximum absolute errors and rate of convergence for at different perturbation parameter and number of mesh points

εN = 32N = 64N = 128N = 256N = 512
1041.9799e-041.0004e-045.0281e-052.5206e-051.2619e-05
1081.9799e-041.0004e-045.0281e-052.5206e-051.2619e-05
10121.9799e-041.0004e-045.0281e-052.5206e-051.2619e-05
10161.9799e-041.0004e-045.0281e-052.5206e-051.2619e-05
10201.9799e-041.0004e-045.0281e-052.5206e-051.2619e-05
EN1.9799e-041.0004e-045.0281e-052.5206e-051.2619e-05
RN0.98490.99250.99620.9982

References

[1]Longtin A, Milton J. Complex oscillations in the human pupil light reflex with mixed and delayed feedback. Math Biosci. 1988; 90: 183-99.

[2]Glizer VY. Asymptotic analysis and Solution of a finite-horizon H control problem for singularly perturbed linear systems with small state delay, J Optim Theory Appl. 2003; 117: 295-325.

[3]Culshaw RV, Ruan S. A delay differential equation model of HIV infection of C D4+ T-cells, Math Biosci. 2000; 165: 27-39.

[4]Elsâ EL. ™golâ™ts, Qualitative methods in mathematical analysis. American Mathematical Society. 1964; 12.

[5]Reddy YN, Soujanya GBSL, Phaneendra K. Numerical integration method for singularly perturbed delay differential equations, Int J Appl Sci Eng. 2012; 10(3): 249-61.

[6]Sirisha CL, Reddy YN. Numerical integration of singularly perturbed delay differential equations using exponential integrating factor, Math Commun. 2017; 22(2): 251-64.

[7]Gadissa G, File G. Fitted fourth order scheme for singularly perturbed delay convection-diffusion equations, Ethiopian J Edu Sci. 2019; 14(2): 102-18.

[8]Subburayan V, Ramanujam N. An initial value technique for singularly perturbed convection–diffusion problems with a negative shift, J Optim Theor Appl. 2013; 158(1): 234-50.

[9]Geng FZ, Qian SP. A hybrid method for singularly perturbed delay boundary value problems exhibiting a right boundary layer. Bull Iranian Math Soc. 2015; 41(5): 1235-247.

[10]Kumar NS, Rao RN. A second order stabilized central difference method for singularly perturbed differential equations with a large negative shift. Diff Eq Dyn Sys. 2020: 1-18, doi: 10.1007/s12591-020-00532-w.

[11]Kumar D., Kumari P. Parameter-uniform numerical treatment of singularly perturbed initial-boundary value problems with large delay. Appl Numer Math. 2020; 153: 412-29.

[12]Debela HG, Duressa GF. Accelerated fitted operator finite difference method for singularly perturbed delay differential equations with non-local boundary condition. J Egy Math Soc. 2020; 28: 1-16.

[13]Liu Z, Papageorgiou NS. Positive solutions for resonant (p, q)-equations with convection. Adv Nonlinear Anal. 2021; 10: 217-32.

[14]Perera K, Shivaji R, Sim I. A class of semipositone p-Laplacian problems with a critical growth reaction term. Adv Nonlinear Anal. 2020; 9: 516-25.

[15]Tang X, Chen S. Singularly perturbed Choquard equations with nonlinearity satisfying Berestycki-Lions assumptions. Adv Nonlinear Anal. 2020; 9: 413-37.

[16]Jiang S, Liang L, Sun M, Su F., Uniform high-order convergence of multiscale finite element computation on a graded recursion for singular perturbation. Electron Res Arch. 2020; 28: 935-49.

[17]Varma VD, Chamakuri N, Nadupuri SK. Discontinuous Galerkin solution of the convection-diffusion-reaction equations in fluidized beds. Appl Numer Math. 2020; 153: 188-201.

[18]O'Malley RE. Singular perturbation methods for ordinary differential equations. New York: Springe, 1991; 89: 8-225.

[19]Siraj MK, Duressa GF, Bullo TA. Fourth-order stable central difference with Richardson extrapolation method for second-order self-adjoint singularly perturbed boundary value problems. J Egyp Math Soc. 2019; 27(1): 50.

[20]Smith GD. Numerical solution of partial differential equations: finite difference methods. New York: Oxford university press. 1985.

Acknowledgements

The authors wish to express their thanks to Jimma University, College of Natural Sciences, for technical support and the authors of the literature for the provided scientific aspects and idea for this work.

Corresponding author

Habtamu Garoma Debela can be contacted at: habte200@gmail.com

Related articles