# Interpolating Stabilized Element Free Galerkin Method for Neutral Delay Fractional Damped Diffusion-Wave Equation

Authors: Mostafa Abbaszadeh; Mehdi Dehghan; Mahmoud A. Zaky; Ahmed S. Hendy
Journal: Journal of Function Spaces
 (2021)

---

## Abstract

A numerical solution for neutral delay fractional order partial differential equations involving the Caputo fractional derivative is constructed. In line with this goal, the drift term and the time Caputo fractional derivative are discretized by a finite difference approximation. The energy method is used to investigate the rate of convergence and unconditional stability of the temporal discretization. The interpolation of moving Kriging technique is then used to approximate the space derivative, yielding a meshless numerical formulation. We conclude with some numerical experiments that validate the theoretical findings.


---

## Body

## 1. Introduction

Partial differential equations (PDEs) with time delay play an important role in the mathematical modeling of complex phenomena and processes whose states depend not only on a given moment in time but also on one or more previous moments. We can mention a simple scenario involving the hemodynamic behavior of a person suffering from low or high glucose decompensation. This person can then be given intravenous insulin to compensate for the low level. Because the drug must be introduced into the bloodstream for it to take effect, the preceding scenario can be interpreted as a delay problem. As a result, there is a growing interest in studying biological and physical models with delay. The solutions of delay PDEs may represent voltage, concentrations, temperature, or various particle densities such as bacteria, cells, animals, and chemicals [1–3].Delay PDEs with fractional derivatives have recently been studied using various numerical and analytical techniques such as [4–8]. It was pointed out in [9] that the derivatives of the dependent variable in the neutral type delay differential equations are both with and without time delays. Delay differential equations of neutral type appear in a variety of new phenomena, and its theory is even more complicated than the theory of nonneutral delay differential equations. From both a theoretical and practical standpoint, the oscillatory behavior of neutral system solutions is important. For some applications, such as the population growth, motion of radiating electrons, and the spread of epidemics in networks with lossless transmission lines, we refer the interested reader to [9–14].A consideration of the following fractional PDE with a constant delay is the goal of this paper. For that end, we introduce(1)D0CtνΨt,x+Ψtt,x=ΔΨtt,x+ΔΨt,x+ΔΨt−s,x+ft,x,t∈0,T,x∈Λ,with initial and boundary conditions
(2)Ψ0,x=ψx,∂Ψt,x∂tt=0=φx,x∈Λ,Ψt,x=gt,x,x∈∂Λ,t∈0,T,Ψt,x=φ∗t,x,t,x∈−s,0×Λ,where D0Ctν is the Caputo fractional derivative which is defined by
(3)D0CtνΨt,x=∂νΨt,x∂tν=1Γ2−ν∫0t∂2Ψs,x∂s2t−s1−νds,ν∈1,2.A novel interpolating element-free Galerkin approach to approximate the solution of the two-dimensional elastoplasticity problems was constructed in [15] using the interpolating moving least squares scheme for obtaining the shape function. Moreover, an improved element-free Galerkin scheme to solve nonlinear elastic large deformation problems was considered in [16]. The interpolating moving least squares approach using a nonsingular weight function is employed in [17] to approximate the solution of the problem of inhomogeneous swelling of polymer gels, and also the penalty scheme is used to enforce the displacement boundary condition; thus, an improved element-free Galerkin approach was constructed.The interpolating element free Galerkin method has been developed to solve a variety of problems, including two-dimensional elastoplasticity problems [15, 18], two-dimensional potential problems [19], two- and three-dimensional Stokes flow problems [20], two-dimensional large deformation problems [21], incompressible Navier-Stokes equation [22], steady heat conduction problems [23], two-dimensional transient heat conduction problems [24], three-dimensional wave equations [25], two-dimensional Schrödinger equation [26], two-dimensional large deformation of inhomogeneous swelling of gels [27], biological populations [28], two-dimensional elastodynamics problems [29], and two-dimensional unsteady state heat conduction problems [30]. The theoretical analysis for the complex moving least squares approximation, the properties of its shape function, and its stability was analyzed in [31]. In [32], a variational multiscale interpolating element-free Galerkin scheme was established for solving the Darcy flow. For the numerical solution of generalised Oseen problems, a novel variational multiscale interpolating element-free Galerkin scheme was developed in [33] based on moving Kriging interpolation for obtaining shape functions using the Kronecker delta function. Zaky and Hendy [34] constructed a finite difference/Galerkin spectral approach for solving the Higgs boson equation in the de Sitter spacetime universe, which can inherit the discrete energy dissipation property. A high-order efficient difference/Galerkin spectral approach was proposed in [35] for solving the time-space fractional Ginzburg-Landau equation. Hendy and Zaky [36] proposed a finite difference/spectral method based on the L1 formula on nonuniform meshes for time stepping and the Legendre-Galerkin spectral approach for solving a coupled system of nonlinear multiterm time-space fractional diffusion equations.This paper is built up as follows. In Section2, the temporal discretization is discussed. The analysis of the temporal discretization scheme is constructed in Section 3. The moving Kriging technique and its implementation are demonstrated in Section 4. Finally, numerical experiments are presented in Section 5 to illustrate the analysis of the obtained scheme.


## 2. Temporal Discretization

Assume thatτ=s/m such that m is a positive integer. Take N=T/τ and tn=nτ,∀n∈ℕ+∪0. Also, to make t=s,2s,⋯ being grid points, the time-variable step size should be surrounded by s=mτ instead of τ=T/N1 for N1∈ℤ+. Thus, tn=nτ for n=−m,−m+1,⋯,0. Here, we present a time-discrete scheme for Equation (1). For any function ξn=ξx,y,tn, we set
(4)ξn−1/2=12ξn+ξn−1,δtξn−1/2=1τξn−ξn−1.Lemma 1 (see [37]).
Assumeϕt∈C20,tn and ν∈1,2. Then
(5)∫0tnϕ′ttn−tν−1dt−BCνϕtn,ϕt0≤Cνmaxϕ′′tτ3−ν0≤t≤tn,in which
(6)BCνϕtn,ϕt0=1τλ0ϕtn−∑k=1n−1λn−k−1−λn−kϕtk−λn−1ϕt0,λk=τ2−ν2−νk+12−ν−k2−ν.LetΨ be the exact solution of (1) and
(7)wt,x,y,z=1Γ2−ν∫0tt−s1−ν∂vs,x,y,z∂sds,where vt,x,y,z=∂Ψt,x,y,z/∂t. Thus, Equation (1) at tn,x,y,z can be rewritten as
(8)wtn−1/2,x+vtn−1/2,x=Δvtn−1/2,x+ΔΨtn−1/2,x+ΔΨtn−1/2−m,x+ftn−1/2,x,n≥0.Making use of Taylor expansion yields(9)vn−12=δtΨn−12+e1n−12,(10)wn−1/2+vn−1/2=Δvn−1/2+ΔΨn−1/2+ΔΨn−1/2−s+fn−1/2+e2n−1/2.Employing Lemma1 and putting v0=vx,0=φx=φ give
(11)wn−12=1Γ2−ντBCνvn−12,φ+e3n−12.Furthermore, there is a constantc>0 that
(12)e1n−12≤cτ2,e2n−12≤cτ2,e3n−12≤cτ3−ν.Substituting the above result into (10) arrives at
(13)1Γ2−ντBCνδtΨn−1/2,φ+δtΨn−1/2=ΔδtΨn−1/2+ΔΨn−1/2+ΔΨn−1/2−m+fn−1/2+Eνn−1/2,n≥1,in which there exists C∈ℝ+ such that
(14)Eνn−1/2≤Cτ3−ν.RemovingEνn−1/2 yields
(15)1Γ2−ντBCνδtUn−1/2,φ+δtUn−1/2=ΔδtUn−1/2+ΔUn−1/2+ΔUn−1/2−m+fn−1/2,n≥1.In the current paper,Un is an approximation of exact solution Ψn.


## 3. Analysis of the Temporal Discretization

In the current section, we check the stability of the numerical procedure.Lemma 2 (see [38]).
Letωs be a nonnegative sequence, and the sequence χs satisfies
(16)χ0≤a0,χs≤a0+∑r=0s−1br+∑r=0s−1ωrχr,s≥1,
Then, fora0≥0 and b0≥0, we have
(17)χs≤a0+∑r=0s−1brexp∑r=0s−1ωr,s≥1.Lemma 3 (see [37]).
For anyK=K1,K2,⋯ and Q, we obtain
(18)∑r=1Nλ0Kr−∑s=1r−1λr−s−1−λr−sKs−λr−1QKr≥tN1−ν2τ∑r=1NKr2−tN2−ν22−νQ2.Theorem 4.
LetΨs∈H01Λ; then scheme (15) is unconditionally stable.Proof.
We defineWs=Ψs−Us. Now, we have
(19)1τΓ2−νλ0δtWs−1/2−∑k=1s−1λs−k−1−λs−kδtWk−1/2+δtWs−1/2=ΔWs−1/2+ΔWs−1/2−m.
Multiplying relation (19) by τδtWs−1/2, integrating over Λ and then summing from s=1 to M give
(20)1Γ2−ν∑s=1Mλ0δtWs−1/2,δtWs−1/2−∑k=1s−1λs−k−1−λs−kδtWk−1/2,δtWs−1/2+τ∑s=1MδtWs−1/2,δtWs−1/2=τ∑s=1MΔδtWs−1/2,δtWs−1/2+τ∑s=1MΔWs−1/2,δtWs−1/2+τ∑s=1MΔWs−1/2−m,δtWs−1/2.
Recalling the left hand side of the above relation, invoking Schwartz inequality and Lemma3 yields
(21)1Γ2−ν∑s=1Mλ0δtWs−1/2L2Λ2−∑k=1s−1λs−k−1−λs−kδtWk−1/2,δtWs−1/2≥1Γ2−ν∑s=1Mλ0δtWs−1/2L2Λ2−∑k=1s−1λs−k−1−λs−kδtWk−1/2L2ΛδtWs−1/2L2Λ=1Γ2−ν∑s=1MBCνδtWs−1/2L2Λ,0δtWs−1/2L2Λ≥τtM1−ν2Γ2−ν∑s=1MδtWs−1/2L2Λ2.
Moreover, for the first term in the right hand side of Equation (20), we have
(22)τ∑s=1MΔδtWs−12,δtWs−12=−τ∑s=1M∇δtWs−12,∇δtWs−12=−τ∑s=1M∇δtWs−12L2Λ2.
On the other hand, according to some simple mathematical actions, we have(23)τ∑s=1MΔWs−1/2,δtWs−1/2=−τ∑s=1M∇Ws−1/2,∇δtWs−1/2=−τ∑s=1M∫Λ∇Ws+∇Ws−12∇Ws−∇Ws−1τdΛ=−12∑s=1M∫Λ∇Ws2−∇Ws−12dΛ=12∑s=1M∫Λ∇Ws−12dΛ−∫Λ∇Ws2dΛ=12∑s=1M∇Ws−1L2Λ2−∇WsL2Λ2=12∇W0L2Λ2−∇WML2Λ2.
Also, for the delay term, we arrive at(24)τ∑s=1MΔWs−m−12,δtWs−12=−τ∑s=1M∇Ws−m−12,∇δtWs−12≤τ∑s=1M∇Ws−m−12L2Λ∇δtWs−12L2Λ≤τ2∑s=1M∇Ws−m−1/2L2Λ2+τ2∑s=1M∇δtWs−1/2L2Λ2.
Replacing the above relations in Equation (20) yields
(25)τtM1−ν2Γ2−ν∑s=1MδtWs−12L2Λ2+τ∑s=1MδtWs−12L2Λ2≤−τ∑s=1M∇δtWs−12L2Λ2+12∇W0L2Λ2−∇WML2Λ2+τ2∑s=1M∇Ws−m−12L2Λ2+τ2∑s=1M∇δtWs−12L2Λ2,or
(26)τtM1−νΓ2−ν∑s=1MδtWs−12L2Λ2+2τ∑s=1MδtWs−12L2Λ2≤−2τ∑s=1M∇δtWs−12L2Λ2+∇W0L2Λ2−∇WML2Λ2+τ∑s=1M∇Ws−m−12L2Λ2+τ∑s=1M∇δtWs−12L2Λ2.
Now, Equation (26) can be simplified as
(27)∇WML2Λ2≤∇W0L2Λ2+τ∑s=1M∇Ws−m−12L2Λ2.
Changing index fromM to s arrives at
(28)∇WsL2Λ≤∇W0L2Λ+2τ∑k=1s∇Wk−mL2Λ=∇W0L2Λ+2τ∑k=1s∇φ∗k−mL2Λ.
The use of Equation (29) and Lemma 2 yields
(29)∇WsL2Λ≤∇W0L2Λ+2τ∑k=1s∇φ∗k−mL2Λexp2sτ≤∇W0L2Λ+2τ∑k=1s∇φ∗k−mL2Λexp2T.
Thus, there existsC∈ℝ+ that
(30)∇WsL2Λ≤C∇W0L2Λ+2τ∑k=1s∇φ∗k−mL2Λ.


## 4. Moving Kriging Interpolation and Its Implementation

Following [39, 40], we will invoke the technique of moving Kriging. Up to our knowledge and armed by the fact of the advantage of less CPU time consuming needed for the element free Galerkin approach based on the shape functions of moving Kriging than what needed for the element free Galerkin approach based on the shape functions of moving least squares approximation. In the meantime, the shape functions of moving Kriging interpolation can be deduced, which is analogous to moving least squares approximation over subdomain Λ1⊂Λ. Let Ψhx is the approximate solution of Ψx on Λ. The local approximation is formulated for any subdomain as
(31)Ψhx=∑r=1mqrxar+Sx=qTxa+Sx,such that qr and ar are monomial basis functions and monomial coefficients, respectively. Also, Sx be the realization of a stochastic process. The covariance matrix of Sx is given as
(32)covSxi,Sxr=σ2EExi,xr,in which
(i)
EExi,xr is the correlation matrix(ii)
Exi,xr is the correlation function between any pair of nodes located at xi and xrThe correlation function is defined as [39, 40]
(33)Exi,xr=exp−θrir2,rir=xi−xr,such thatθ>0 is a value of the correlation parameter [39, 40]. Using the best linear unbiased (BLUP) [39], we can write Equation (31) as follows [39, 40]
(34)Ψhx=qTxη+rTxE−1u−Qη,in which
(35)η=QTE−1Q−1QTE−1u.We will introduce some notations. The vector of knownm functions can be written as follows [39, 40]
(36)qx=q1xq2x⋮qmx1×m,and the matrix of defined function values at the nodesx1,x2,⋯,xn has the following representation [39, 40]
(37)Q=q1x1q2x1⋯qmx1q1x2q2x2⋯qmx2⋮⋮⋱⋮q1xnq2xn⋯qmxnn×m.The correlation matrix is given as [39, 40]
(38)EExi,xr=1Ex1,x2⋯Ex1,xnEx2,x11⋯Ex2,xn⋮⋮⋱⋮Exn,x1Exn,x2⋯1n×n.The correlation vector at the nodesx1,x2,⋯,xn has the following form
(39)rx=Ex1,xEx2,x⋮Exn,x.The matricesA and B are given as
(40)A=QTE−1Q−1QTE−1,B=E−1I−QA,whereI is the n×n identity matrix. Accordingly, Equation (34) can be written as follows [39, 40]
(41)Ψhx=qTxAu+rTxBu,or
(42)Ψhx=qTxA+rTxBu=∑k=1nϕkxΨk=φxu,where the moving Kriging approach’s shape functions are as follows [39, 40]:
(43)φx=qTxA+rTxBr=ϕ1,ϕ2,⋯,ϕnT.Now, we are ready to implement that kind of interpolation to the problem under consideration. Let the approximation solution of this equation be(44)Unx,y,z=∑j=1Nϖjnϕjx,y,z,in whichϕjx are shape functions of moving Kriging approximation. Substituting Equation (44) in relation (15) gives
(45)1Γ2−ντBCνδt∑j=1Nϖjn−1/2ϕjx,y,z,φ+∑j=1Nδtϖjn−1/2ϕjx,y,z=∑j=1Nδtϖjn−1/2Δϕjx,y,z+∑j=1Nϖjn−1/2Δϕjx,y,z+∑j=1Nδtϖjn−1/2−mΔϕjx,y,z+fn−1/2.By collocating a set of arbitrary distributed nodesxi,yi,zii=1N in the computational domain Λ concludes
(46)1Γ2−ντBCνδt∑j=1Nϖjn−1/2ϕjxi,yi,zi,φ+∑j=1Nδtϖjn−1/2ϕjxi,yi,zi=∑j=1Nδtϖjn−1/2Δϕjxi,yi,zi+∑j=1Nϖjn−1/2Δϕjxi,yi,zi+∑j=1Nδtϖjn−1/2−mΔϕjxi,yi,zi+fn−1/2xi,yi,zi.After doing some simplifications, we have(47)μλ0+1∑j=1Nϖjnϕjxi,yi,zi−1+τ2∑j=1NϖjnΔϕjxi,yi,zi=μλ0+1∑j=1Nϖjn−1ϕjxi,yi,zix,y,z+−1+τ2∑j=1Nϖjn−1Δϕjxi,yi,zi+μτλn−1φ+μτ∑k=1n−1λn−k−1−λn−k∑j=1Nϖjkϕjxi,yi,zi−∑j=1Nϖjk−1ϕjxi,yi,zi+τ2∑j=1Nϖjn−mϕjxi,yi,zi−∑j=1Nϖjn−m−1ϕjxi,yi,zi+τfn−1/2,whereμ=1/Γ2−ντ. Now, the above formulation yields the following system of equations
(48)AΛn=Fn,in which
(49)Aij=μλ0+1ϕjxi,yi,zi−1+τ2Δϕjxi,yi,zi,Fi=μλ0+1∑j=1Nϖjn−1ϕjx,y,z+−1+τ2∑j=1Nϖjn−1Δϕjx,y,z+μτ∑k=1n−1λn−k−1−λn−k∑j=1Nϖjkϕjx,y,z−∑j=1Nϖjk−1ϕjx,y,z+μτλn−1φ+τ2∑j=1Nϖjn−mϕjx,y,z−∑j=1Nϖjn−m−1ϕjx,y,z+τfn−1/2.


## 5. Numerical Verification

In the current section, we investigate the convergence, capability, and stability of the developed numerical procedure. Also, the computational rate is calculated by(50)C−order=log2−1×logEh,2τEh,τ.We consider the following problem(51)1Γ2−ν∫0t∂2Ψx,y,s∂s2dst−sν−1+∂Ψx,y,t∂t=∂∂t∂2Ψx,y,t∂x2+∂2Ψx,y,t∂y2+∂2Ψx,y,t∂x2+∂2Ψx,y,t∂y2+∂2Ψx,y,t−s∂x2+∂2Ψx,y,t−s∂y2+fx,y,t,in which
(52)fx,y,t=cosπxcosπy2t+2π2s−t2+4π2t+2π2t2,the initial conditions are
(53)Ψx,y,tt=0=0,∂Ψx,y,t∂tt=0=0,x,y∈Λ,and also
(54)Ψx,y,t=t2cosπxcosπy,x,y,t∈Λ×−s,0,with no-flux boundary condition. The exact solution is Ψx,y,t=t2cosπxcosπy.Table1 shows the results obtained based on the 500 collocation points, ν=1.1, ν=1.2, ν=1.9, and different values of τ. Table 1 confirms that the theoretical order (TO) in temporal direction is near to the computational order, i.e., 3−ν. Figure 1 demonstrates the approximate solutions (a) and its absolute errors (b) on square domains 0,1×0,1 (top figures), 0,4×0,4 (middle figures), and 0,8×0,8 (bottom figures) with τ=0.001, ν=1.5, and also 1000 collocation points. Figure 2 illustrates the approximate solutions (a) and its absolute errors (b) on irregular domains
(55)rθ=310n21+2n+n.2−n+1cosnθ,where n=4 (top figures), n=8 (middle figures), and n=12 (bottom figures) with τ=0.001, ν=1.5, and also 1000 collocation points. Figure 3 presents the approximate solutions (a) and its absolute errors (b) on irregular domains
(56)rθ=310n21+2n+n.2−n+1cosnθ,with τ=0.001, ν=1.5, and also 1000 collocation points. Figure 4 presents the approximate solutions (a) and its absolute errors (b) on irregular domains
(57)rθ=0.3cos2θ+sqrt1.1−sinθ2,with τ=0.001, ν=1.5, and also 1000 collocation points.Table 1
Results obtained with 500 collocation points for Example 1.
τν=1.1ν=1.2ν=1.9CPU time(s)L∞Cτ-orderL∞Cτ-orderL∞Cτ-order0.11.20×10−1—1.04×10−1—7.74×10−1—0.50.053.00×10−22.002.69×10−21.946.02×10−10.360.90.0257.51×10−31.997.20×10−31.903.31×10−10.863.20.01251.88×10−31.991.97×10−31.871.61×10−11.0410.40.006254.70×10−41.995.46×10−41.857.37×10−21.1337.10.0031251.17×10−41.991.53×10−41.833.29×10−21.1687.50.00156252.94×10−51.994.36×10−51.811.45×10−21.18163.40.000781257.38×10−61.991.25×10−51.806.35×10−31.19277.3TO1.91.81.1Figure 1
Approximate solution and its absolute error on square domains.Figure 2
Approximate solution and its absolute error on irregular domains.Figure 3
Approximate solution and its absolute error on irregular domains.Figure 4
Approximate solution and its absolute error on irregular domains.


## 6. Conclusion

The current paper presented a new numerical procedure for solving fractional damped diffusion-wave equations with delay. In this process, the time derivative is discretized by a finite difference scheme, and we constructed a time-discrete scheme. The stability and convergence of the proposed numerical formulation are studied, analytically and numerically. Then, the moving Kriging interpolation technique, as a meshless method, is used to get a fully discrete scheme. The proposed numerical method is flexible to simulate a wide range of PDEs including delay PDEs on irregular computational domains. Finally, an example is provided to demonstrate the stability and convergence of the new technique.



---

*Source: 6665420-2021-11-19.xml*