# Analysis of Two-Level Mesh Method for Partial Integro-Differential Equation

Authors: Quan Tang; Ziyang Luo; Xindong Zhang; Juan Liu
Journal: Journal of Function Spaces
 (2022)

---

## Abstract

In this paper, we present two-level mesh scheme to solve partial integro-differential equation. The proposed method is based on a finite difference method. For the first step, we use finite difference method in time and global radial basis function (RBF) finite difference (FD) in space. For the second step, we use the finite difference method to solve the proposed problem. This two-level mesh scheme is obtained by combining the radial basis function with finite difference. We prove the stability and convergence of scheme and show that the convergence order isOτ2+h2, where τ and h are the time step size and space step size, respectively. The results of numerical examples are compared with analytical solutions to show the efficiency of proposed scheme. The numerical results are in good agreement with theoretical ones.


---

## Body

## 1. Introduction

For partial integro-differential equation,(1)∂ux,t∂t=Aux,t+∫0tκt,sBux,sds+fx,t,x,t∈Ω×0,T,where A and B are linear or nonlinear operators, the kernel κt,s is the weakly singular function, and fx,t is given function. T is the terminal time, and Ω denotes a bounded domain with boundary Γ, where Ω⊂ℝd with d=1,2,3.This kind partial integro-differential equation is indicated by Volterra [1]; the new results of Volterra partial integro-differential equation can be found in [2, 3]. The problem Equation (1) can be found in the model of heat flow in a conducting materials with memory [4–7] and in viscoelastic mechanics [8–10].In this paper, our particular concern is to obtain stability and convergence results for finite difference scheme and then use two-level mesh method to solve the following partial integro-differential equation. The following equation can be obtained by Equation (1), as Aux,t=Δux,t=∂2ux,t/∂x2, Bux,s=∂2ux,s/∂x2, and κt,s=t−sβ−1.
(2)∂ux,t∂t=∂2ux,t∂x2+∫0tt−sβ−1∂2ux,s∂x2ds+fx,t,x,t∈0,L×0,T,ux,0=φx,x∈0,L,u0,t=uL,t=0,t∈0,T,where 0<β<1, fx,t, and φx are known functions. Throughout this paper, we assume that Equation (2) has a unique solution in 0,L×0,T and suppose that ∂4ux,t/∂x4 and ∂3ux,t/∂t3 are continuous in 0,L×0,T.Some attention has been paid to numerical methods of Equation (2), such as spectral collocation approximation by Kim and Choi [11], orthogonal spline collocation by Yan and Fairweather [12], alternating direction implicit difference scheme by Qiao and Xu [13], finite difference scheme by Tang [14], and compact difference scheme by Luo et al. [15]. If A is nonlinear operator, the numerical investigations for Equation (1) have been studied by several authors; one can refer to [16–18]. If we replace ∂ux,t/∂t with a fractional differential operator, then the fractional integro-differential equations can be obtained. The numerical solutions of fractional integro-differential equations have been studied by several authors; one can refer to [19–23].Radial basis function, as a meshless method, has advantage of the ability to handle scattered node layouts and has the characteristics of simple implementation and high precision. The mainly purpose of this paper is to construct a numerical scheme using RBF and FD to solve partial integro-differential equation. In recent years, there are many researchers that studied the meshless method for solving partial differential equations. For example, Oruç developed a meshless approach based on radial basis function-finite difference (RBF-FD) method for solving one-dimensional, two-dimensional, and three-dimensional Schrödinger system [24, 25], two-dimensional fractional cable equation [26], and wave equations [27, 28]. Nikan and Avazzadeh [29] constructed an improved localized radial basis-pseudospectral numerical method for solving fractional reaction-subdiffusion equation.This paper is outlined as follows. In Section2, we introduce some notations and lemma and propose a two-level mesh discrete scheme based on finite difference. In Section 3, we discuss the stability and error estimation of the constructed discrete scheme. In Section 4, we consider the two-level algorithm which will be used in practical calculation. For more details on a two-level method, see [30, 31]. In Section 5, the results of numerical experiments are compared with exact solution to confirm the efficiency of the proposed discrete scheme. Finally, the conclusion and future studies are given for this paper in Section 6.


## 2. Construction of Numerical Scheme

In this section, we mainly consider the construction of discrete scheme. The following notations will be used in this section and others. LetM and N be positive integers. Let xj≔jh with j=0,1,⋯,M and tk≔kτ with k=0,1,⋯,N, where h≔L/M and τ≔T/N are the spatial and temporal step sizes, respectively. Let Vh=VV=V0,V1,⋯,VM be the grid function space. For any functions V,W∈Vh, we define the following notations:
(3)δxVj=Vj−Vj−1h,δx2Vj=Vj+1−2Vj+Vj−1h2,j=1,2,⋯,M−1,V,W=h∑j=1M−1VjWj,V2=V,V,V∞=max1≤j≤M−1Vj.Throughout this paper, the values ofC in different places are allowed to be different. Lemma 1 will be used in the construction of numerical scheme.Lemma 1 (see [32]).
Ifgt∈C20,T, then,
(4)∫0tktk−sβ−1gsds−∑i=0kλi,kgti≤Cmax0≤t≤Tg′′ttkβτ2,1≤k≤N,where C is a constant depending only on β0<β<1 and
(5)λi,k=τβββ+1∗k−1β+1−k−1−βkβ,i=0,k−i+1β+1−2k−iβ+1+k−i−1β+1,1≤i≤k−1,1,i=k.Now, we will consider the terms of Equation (2) one by one, and construct the discrete scheme of Equation (2).For the term∂ux,t/∂t in the left hand side of Equation (2), by the Taylor expansion, we obtain that
(6)∂uxj,tk∂t=3uxj,tk−4uxj,tk−1+uxj,tk−22τ+Oτ2,1≤j≤M,1≤k≤N.For the first term∂2ux,t/∂x2 in the right hand side of Equation (2), by the Taylor expansion, we obtain that
(7)∂2uxj,tk∂x2=uxj+1,tk−2uxj,tk+uxj−1,tkh2+Oh2,1≤j≤M−1,1≤k≤N.For the second term∫0tt−sβ−1∂2ux,s/∂x2ds in the right hand side of Equation (2), by Lemma 1 and the Taylor expansion, we obtain that
(8)∫0tktk−sβ−1∂2uxj,s∂x2ds=∑i=0kλi,k∂2uxj,ti∂x2+Oτ2=∑i=0kλi,kuxj+1,ti−2uxj,ti+uxj−1,tih2+Oτ2+h2,1≤j≤M−1,1≤k≤N.Forj=0,1,⋯,M and k=0,1,⋯,N, we use ujk and fjk to denote the exact solutions of uxj,tk and fxj,tk, respectively. If substituting Equations (6), (7), and (8) into Equation (2), then, we have
(9)3ujk−4ujk−1+ujk−22τ=δx2ujk+∑i=0kλi,kδx2uji+fjk+Rjk,1≤j≤M−1,1≤k≤N,uj0=φ1xj,0≤j≤M,u0k=uMk=0,1≤k≤N,where Rjk≤Cτ2+h2 and C is a positive constant with independent τ and h.LetUj−1 and Ujk be the numerical approximation solutions of uxj,t−1 and uxj,tk, respectively. Omitting small term Rjk in Equation (9), the discrete scheme of Equation (2) can be obtained as follows:
(10)3Ujk−4Ujk−1+Ujk−22τ=δx2Ujk+∑i=0kλi,kδx2Uji+fjk,1≤j≤M−1,1≤k≤N,Uj0=φ1xj,0≤j≤M,U0k=UMk=0,1≤k≤N.


## 3. Stability and Convergence

In this section, our purpose is to investigate the stability and convergence of the discrete scheme Equation (10). Some lemmas will be given for the proof of the main results. For convenience, we omit the subscript j from this section.Lemma 2 (see [16]).
IfV,W∈Vh, then δx2V,W=−δxV,δxW.Lemma 3 (see [16]).
IfV,W∈Vh, then δx2V,W≤4/h2VW.Lemma 4. (see [33]).
Letλi,k be defined as in (5). For 0<β<1 and 0≤i≤k, each of the following holds.
(i)
λi,k>0, in particular, 0<λ0,k<τβ/β+1,1≤k≤N(ii)
∑i=1k−1λi,k≤τβ/β,2≤k≤NLetut≔∂ux,t/∂t, in particular, denote ux,t−1≈u−1=ux,0−τutx,0, where t−1=−τ. By Taylor expansion, Lemma 5 can be obtained easily. We omit the proof of this lemma.Lemma 5.
Letu−1=u0−τut0 and ε−1=ux,t−1−u−1. Assume that ut∈C20,T; then ε−1≤Cτ2.Theorem6 is about stability of discrete scheme Equation (10).Theorem 6.
LetUk be the numerical solution of Equation (10) with the given initial and boundary conditions; if τ≤C1h2, then,
(11)Uk≤CU0+max1≤s≤kfs,1≤k≤N,where C1 and C are positive constants.Proof.
Sinceu−1=u0−τut, by the assumption that ut is continuous in the given domains, we get that u−1≤Cu0+τut<∞, where C is a positive constant.
By Equation (10), we get that
(12)Uk=4Uk−1−Uk−2+2τ31+λk,kδx2Uk+2τ3∑i=0k−1λi,kδx2Ui+2τ3fk.
Multiplying both sides of Equation (12) by Uk and integrating on 0,L, we obtain that
(13)Uk,Uk=4Uk−1−Uk−2,Uk+2τ31+λk,kδx2Uk,Uk+2τ3∑i=0k−1λi,kδx2Ui,Uk+2τ3fk,Uk.
By Lemma2, we have 2τ/31+λk,kδx2Uk,Uk=−2τ/31+λk,kδxUk,δxUk≤0. Thus, Equation (13) can be rewritten as
(14)Uk,Uk≤4Uk−1−Uk−2,Uk+2τ3∑i=0k−1λi,kδx2Ui,Uk+2τ3fk,Uk.
In the following, we will use mathematical induction to prove our theorem. Fork=1, by Equation (14) and λ0,1=τβ/β+1, we get that
(15)U1,U1≤4U0−U−1,U1+2τ1+β3β+1δx2U0,U1+2τ3f1,U1.
Using Cauchy-Schwarz inequality and Lemma3, we get
(16)U12≤4U0U1+U−1U1+8τ1+β3β+1h2U0U1+2τ3f1U1,
i.e.,(17)U1≤4U0+U−1+C18τβ3β+1U0+2τ3f1≤CU0+f1.
Assuming thatUk≤CU0+max1≤s≤kfs holds for k≤N−1, we want to prove that
(18)UN≤CU0+max1≤s≤Nfs.
Similarly, fork=N, by Equation (14), we get that
(19)UN,UN≤4UN−1−UN−2,UN+2τ3∑i=0N−1λi,Nδx2Ui,UN+2τ3fN,UN.
Using Cauchy-Schwarz inequality, Lemma3, Lemma 4 and τ/h2≤C1, we get
(20)UN2≤4UN−1−UN−2UN+8τ3h2∑i=0N−1λi,NUiUN+2τ3fNUN≤4UN−1UN+UN−2UN+8τ3h2∑i=0N−1λi,NUiUN+2τ3fNUN,
i.e.,(21)UN≤4UN−1+UN−2+8τ3h2∑i=0N−1λi,NUi+2τ3fN≤C2U0+max1≤s≤N−1fs+C3U0+max1≤s≤N−2fs+83C1λ0,1U0+83C1∑i=1N−1λi,NC4U0+max1≤s≤ifs+2τ3fN≤C2+C3+83C1λ0,1+83C1C4∑i=1N−1λi,NU0+C2+C3+83C1λ0,1+83C1C4∑i=1N−1λi,N+2τ3max1≤s≤Nfs≤CU0+max1≤s≤Nfs.
This proves the theorem.In order to establish the convergence of proposed method, the following notations and results are needed. ForV0∈Vh, let Φ0=u0−V0, Ψ0=V0−U0, Φt0=ut0−Vt0, and Ψt0=Vt0−Ut0; then we obtain that
(22)Φ0≤Ch2u0,Φt0≤Ch2ut0,Ψ0=Ψt0=0.The convergence theorem and an error bound will been given for discrete scheme Equation (10) as follows.Theorem 7.
Letuk be the exact solution of Equation (2) and Uk be the numerical solution of Equation (10). Letting ek=uk−Uk, then, for τ≤C1h2 and the given initial and boundary conditions, the following holds:
(23)ek≤Cτ2+h2,1≤k≤N,where C1 and C are positive constants.Proof.
Subtracting Equation(10) from Equation(9), we obtain
(24)ek=4ek−1−ek−2+2τ3δx2ek+2τ3∑i=0kλi,kδx2ei+2τ3Rk,
wheree0=0 and e−1 will be discussed in the following.
Multiplying both sides of Equation (25) by ek and integrating on 0,L, by Lemmas 1, 2, and 4, we obtain
(25)ek,ek=4ek−1−ek−2,ek+2τ3δx2ek,ek+2τ3∑i=0kλi,kδx2ei,ek+2τ3Rk,ek≤4ek−1−ek−2,ek+2τ3∑i=0k−1λi,kδx2ei,ek+2τ3Rk,ek.
Using Cauchy-Schwarz inequality, Lemma3 and τ/h2≤C1, we have
(26)ek2≤4ek−1−ek−2ek+8τ3h2∑i=0k−1λi,keiek+2τ3Rkek.≤4ek−1ek+ek−2ek+83C1∑i=0k−1λi,keiek+2τ3Rkek,
i.e.,(27)ek≤4ek−1+ek−2+83C1∑i=0k−1λi,kei+2τ3Rk.
In order to prove this theorem, we first prove thate−1=ux,t−1−U−1≤Cτ2+h2. By Lemma 5 and (22), we get that
(28)e−1=ux,t−1−U−1=ux,t−1−u−1+u−1−U−1=ε−1+u0−τut0−U0−τUt0=ε−1+u0−U0−τut0−Ut0=ε−1+u0−V0+V0−U0−τut0−Vt0+Vt0−Ut0=ε−1+Φ0+Ψ0−τΦt0+Ψt0≤ε−1+Φ0+Ψ0+τΦt0+Ψt0≤Cτ2+h2.
Next, we use mathematical induction to get the convergence theorem and error bound. Fork=1 in (27), by (28) and Rk=Rjk≤Cτ2+h2, we get that
(29)e1≤4e0+e−1+83C1λ0,1e0+2τ3Rk≤Cτ2+h2.
Assuming thatek≤Cτ2+h2 holds for k≤N−1, we want to prove that eN≤Cτ2+h2. By Lemma 4 and (27), we get that
(30)eN≤4eN−1+eN−2+83C1∑i=1N−1λi,Nei+2τ3RN≤4C2τ2+h2+C3τ2+h2+83C1C4∑i=1N−1λi,Nτ2+h2+2τ3C5τ2+h2≤Cτ2+h2.
This completes the proof.


## 4. Two-Level Mesh Method for Integro-Differential Equation

In this section, we will consider the two-level algorithm which will be used in practical calculation. The considered two-level algorithm is given as follows.Algorithm 8.
(Step 1)
On the coarse grid, we use the global RBF method to solve the following problem with Newton iterative:(31)3UHk−4UHn−1+UHn−22τ=ΔUHk+∑i=0kλi,kδx2UHi+fHk,where k∈ℕ+ is the number of iteration. For given n with n≥1, we use UHn−1 to denote the initial value on the coarse mesh, that is, UHn−1=UH0, and use UHn−2 to denote the value of UH−1, by Lemma 5; we have UHn−2=UH0−τUt0=UH−1. We denote UHn=UHk after the iteration stops, and UHk will be used in Step 2.
(Step 2)
On the fine grid, FD method will be used to solve the proposed problem, and we construct the following discrete scheme:(32)3Uhk−4Uhn−1+Uhn−22τ=δx2Uhk+∑i=0kλi,kδx2Uhi+fhk.In Step 1, the global RBF method will be used to solve the problem Equation (31). For the global RBF method, one can refer to [34–37] for details. Let x=x1,x2,⋯,xp,xp+1,⋯,xM′, where x1,x2,⋯,xp within Ω and xp+1,⋯,xM on Γ (Γ is the boundary), and xj=xj1,xj2,⋯,xjM′ represents coordinates of the jth node. Then, we can get that
(33)uHkx=∑j=1Mγjkϕx−xj,where · is the standard Euclidean distance norm, ϕ· is RBF, and γjk is interpolation coefficient. By Equations (31) and (33), we get that
(34)∑j=1Mγjkϕxℓ−xj1−2τ3λk,k−2τ3Δϕxℓ−xj=2τ34UHn−1−UHn−2+∑i=0k−1λi,kδx2UHi+fHk,forxℓ∈Ω\Γ,∑j=1Nγjkϕxℓ−xj=0,forxℓ∈Γ.The matrix form of Equation (34) is
(35)AΛ=G,where A is the coefficient matrix and the element Aℓ,j has the following form:
(36)Aℓ,j=ϕxℓ−xj1−2τ3λk,k−2τ3Δϕxℓ−xj,ℓ=1,⋯,p,j=1,⋯,M,ϕxℓ−xj,ℓ=p+1,⋯,M,j=1,2,⋯,M.The element ofG has the following form:
(37)Gℓ,1=2τ34UHn−1xℓ−UHn−2xℓ+∑i=0k−1λi,kδx2UHixℓ+fHkxℓ,ℓ=1,⋯,p,0,ℓ=p+1,⋯,M,and Λ=γ1k,γ2k,⋯,γMk′. The coefficients γjk will be obtained by solving the system of Equation (35). By Equation (33), we can get the approximate solution of u for any point, which will be used in Step 2. In Step 2, FD method will be used to solve the proposed problem. By Equation (32), we can get the approximate solution on fine grid, that is, U=Ukx1,Ukx2,⋯,UkxM′.


## 5. Numerical Results

In this section, we present numerical results for compact finite difference which supports the theoretical analysis in Section3. All numerical experiments are implemented on an Intel(R) Core(TM) i7-10510U CPU with 1.80 GHz and 16Gb of RAM, using MATLAB(R2020b) software. We will use multiquadric MQ,ϕr=1+cr2 radial basis functions for the numerical experiments, where c=1.Letτ=1/N and h=1/M, where N,M are the numbers of meshes in time and space. Define the L∞-Error and L∞-Rate as the following forms:
(38)L∞‐Error=max1≤n≤Nuxj,tn−ujn∞,L∞‐Rate=logE1/E2/logh1/h2,in which E1 and E2 are errors correspond to grids with mesh size h1 and h2, respectively.Example 9.
In this example, we consider the problem (2) with L=T=1, ux,0=0 and u0,t=u1,t=0. The source term as follows:
(39)fx,t=sinπx4t3+π2t4+24π2ββ+1β+2β+3β+4tβ+4.
Then, its exact solution has the following form:ux,t=t4sinπx.The results of this example will be shown in Tables1 and 2. In Table 1, the convergence rates for space are obtained as Oh2. In Table 2, the convergence rates for time are obtained as Oτ2.Table 1
Errors, convergence rates for space, and CPU times withβ=0.2 and 0.5.
βhτL∞-ErrorL∞-RateCPU time (s)0.21/101/206.7888×10−3—0.0124741/201/401.6680×10−32.030.0145441/401/804.1145×10−42.020.0198291/801/1601.0170×10−42.020.0788131/1601/3202.5176×10−52.010.9658190.51/101/206.4482×10−3—0.0135821/201/401.6016×10−32.010.0143171/401/803.9907×10−42.000.0216611/801/1609.9555×10−52.000.0781901/1601/3202.4852×10−52.000.970974Table 2
Errors, convergence rates for time, and CPU times withβ=0.6 and 0.9.
βhτL∞-ErrorL∞-RateCPU time (s)0.61/201/101.4174×10−3—0.0120511/401/203.4654×10−42.030.0096081/801/408.5053×10−52.030.0159161/1601/802.0928×10−52.020.0922961/3201/1605.1652×10−62.020.7303720.91/201/103.6375×10−3—0.0098431/401/209.3300×10−41.960.0096891/801/402.3642×10−41.980.0145981/1601/805.9514×10−51.990.0938961/3201/1601.4930×10−52.000.712324Table1 shows the maximum error, convergence order in the spatial direction, and CPU times (second) in the case of β=0.2 and 0.5 with M=2N. It can be obviously found that the convergence order in the spatial direction is Oh2, which is in accordance with the convergence analysis.Table2 shows the maximum error, convergence order in the temporal direction, and CPU times (second) of finite difference scheme Equation (10) for β=0.6 and 0.9 with N=2M. It is observed that the convergence order in the temporal direction is Oτ2, which are almost as the theoretical results in Theorem 7.Figure1 shows the exact solution, numerical solution, absolute error, and contour plot of absolute error with β=0.5 at N=1000 and M=500. It demonstrates that the numerical solution of Equation (10) is in excellent agreement with the exact solution.Example 10.
In the second example, we consider the problem (2) with L=1,T=4 and u0,t=u1,t=0, the source term as follows:
(40)fx,t=π2+tββπ2sinπx−4t3+4π2t4+96Γβπ2tβ+4Γβ+5sin2πx.Figure 1
Exact solution, numerical solution, absolute error, and contour plot absolute error for Example9 with β=0.5, M=500, and N=1000.
(a)
Exact solution(b)
Numerical solution(c)
Absolute error(d)
Contour plot absolute errorThen the corresponding initial term isux,0=sinπx, and the exact solution is ux,t=sinπx−t4sin2πx.In Figures2(a) and 2(b), we show the errors in L∞-norm attain second order of accuracy in temporal direction for β=0.2 and β=0.5, respectively. In Figures 2(c) and 2(d), we show the errors in L∞-norm attain second order of accuracy in spatial direction for β=0.2 and β=0.5, respectively. It is observed that the convergence order in the temporal and spatial directions are Oτ2 and Oh2, respectively, which are almost as the theoretical results in Theorem 7.Figure 2
log(error) as a function of 2log(τ) and 2log(h), respectively, for given β.
(a)
β=0.2(b)
β=0.5(c)
β=0.2(d)
β=0.5


## 6. Conclusion

In this paper, we present two-level mesh method based on finite difference scheme to solve partial integro-differential equation. For the first step on coarse grid, we use finite difference method in time and global radial basis function finite difference in space. For the second step on fine grid, the finite difference method is used to solve the proposed problem, and the numerical solutions obtained by first step will be used as the initial conditions. The stability and error are discussed, and the convergence rate is obtained. Numerical examples show the efficiency of proposed scheme. It shows that the numerical results are in good agreement with theoretical analysis.In our future work, we will study the numerical scheme of fractional derivatives, such as Riemann-Liouville derivative, Caputo derivative, and Caputo-Fabrizio derivative. Meanwhile, we would like to investigate fractional derivative in both time and space. In addition, we are interested to combine other numerical methods and RBF to consider the numerical solutions of high dimensional fractional integro-differential equations, irregular region problems, and nonlinear problems.



---

*Source: 4557844-2022-07-13.xml*