# Investigations for a Type of Variable Coefficient Fractional Subdiffusion Equation with Multidelay

Authors: Wei Gu
Journal: Journal of Function Spaces
 (2020)

---

## Abstract

A difference scheme is constructed for a type of variable coefficient time fractional subdiffusion equation with multidelay. Stability and convergence results of the scheme are obtained, and theoretical results are proved by two numerical tests.


---

## Body

## 1. Introduction

This paper will consider the following variable coefficient time fractional subdiffusion equation with arbitrary multidelay(1)∂αu∂tα−Dxuxx=fux,t,ux,t−s1,ux,t−sq,x,t,x,t∈0,1×0,T,ux,t=ϕx,t,x∈0,1,t∈−s,0,u0,t=γt,u1,t=βt,t∈0,T,where 0<c0≤Dx≤c1 is the variable diffusion coefficient, sj>0,j=1,2,⋯,q are delay terms, and s=max1≤i≤qsi; time fractional partial derivative ∂αu/∂tα0<α<1 is defined as
(2)∂αu∂tα=1Γ1−α∫0tt−ξ−α∂ux,ξ∂ξdξ,where Γ· is the gamma function.For simplicity, we consider the following equation with arbitrary two different delays:(3)∂αu∂tα−Dxuxx=fux,t,ux,t−s1,ux,t−s2,x,t,x,t∈0,1×0,T,(4)ux,t=ϕx,t,x∈0,1,t∈−s,0,(5)u0,t=γt,u1,t=βt,t∈0,T,where s2>s1>0 are two delay terms and satisfy s1=m1+δ1τ, s2=m2+δ2τ, m2>m1>0 are integers, and δ1,δ2∈0,1, s=m2+1τ.Applications of delay differential equations (DDEs) can be found in many fields, such as electrotechnics, population dynamics, biology, and economics [1–5]. However, one can obtain the analytical solutions of DDEs in few cases. The researchers turn to numerical methods for solving DDEs in general cases [6–9]. Most studies about DDEs consider one delay; in fact, more than one delays should be considered in some systems. Tian and Kuang considered H-method and GP-stability for DDEs with several delay terms [10]; the applications of such DDEs can be found in [11–14].Fractional delay differential equations (FDDEs) are widely used in automatic control, population dynamics, finance, etc. [15–17]. Numerical solutions for FDDEs can be found in [18–25], where [18–20] considered numerical solutions for FDDEs, and [21–25] considered numerical solutions for fractional delay partial differential equations (FDPDEs). However, the researchers seldom considered the numerical methods for FDDEs with multidelays; some studies considered such FDDEs without theoretical analysis.In this paper, there is an effective difference scheme for solving systems (3)–(5) with multidelay. Some stability and convergence results of the scheme are obtained by mathematical proof, and the theoretical results are proved by two numerical tests.We organize the paper as the follows. Section2 constructs a numerical scheme to solve (3)–(5). Section 3 provides the stability and convergence results by proof. Section 4 presents two numerical tests. Section 5 gives a brief discussion.


## 2. The Construction of the Second Finite Difference Method

The following assumptions are assumed to be true in this paper.(H1)
AssumeDx to be sufficiently smooth function which fulfills 0<c0≤Dx≤c1(H2)
Suppose thatfux,t,ux,t−s1,ux,t−s2,x,t is sufficiently smooth and satisfies [26, 27](6)fμ+ε1,ν+ε2,ι+ε3,x,t−fμ,ν,ι,x,t≤c2ε1+c3ε2+c4ε3,where ε1,ε2,ε3 are arbitrary real numbers and c2>0,c3>0, and c4>0 are three constants.TakingM>0 and N>0 to be two integers, we have space and time stepsize h=1/M and τ=T/N, respectively, and discrete points xi=ih and tk=kτ. Assume Ωhτ=Ωh×Ωτ, Ωh=xi∣0≤i≤M, Ωτ=tk∣−m2+1≤k≤N, and Uik=uxi,tk. Define the following grid function space on Ωhτ:
(7)W=vik∣0≤i≤M,−m2+1≤k≤Nrgb1.00,0.00,0.00,We have(8)δxvi+1/2k=vi+1k−vikh,δx2vik=vi+1k−2vik+vi−1kh2,δxDδxvik=Di+1/2δxvi+1/2k−Di−1/2δxvi−1/2k/h,δtαvk=τ−αΓ2−αd0vk−∑j=1k−1dk−j−1−dk−jvj−dk−1v0,where dj=j+11−α−j1−α,j≥0.We introduce the following two Lemmas,Lemma 1 [28].
Suppose0<α<1, y∈C20,tk; it holds that
(9)1Γ1−α∫0tky′sdstk−sα−τ−αΓ2−αd0ytk−∑j=1k−1dk−j−1−dk−jytj−dk−1yt0≤1Γ2−α1−α12+22−α2−α−1+2−αmax0≤t≤tky″tτ2−α.Lemma 2 [29].
Assume0<α<1; then, it holds that
(i)
dj decrease monotonically with j increases, and 0<dj≤1(ii)
d0=1, ∑j=1k−1dk−j−1−dk−j=d0−dk−1Considering (3) at the point xi,tk, we have
(10)∂αu∂tαxi,tk−∂∂xDx∂ux,t∂xxi,tk=fuxi,tk,uxi,tk−s1,uxi,tk−s2,xi,tk,0≤i≤M,0≤k≤N.From Lemma1, we obtain
(11)∂αu∂tαxi,tk=δtαUik+rik,where
(12)rik=1Γ2−α1−α12+22−α2−α−1+2−αmax0≤t≤tk∂2uxi,tk∂t2τ2−α.From Taylor’s expansion, we have(13)Di+1/2δxUi+1/2k=Di+1/2∂u∂xxi+1/2,tk+h224Di+1/2∂3u∂x3xi+1/2,tk+Oh3,(14)Di−1/2δxUi−1/2k=Di−1/2∂u∂xxi−1/2,tk+h224Di−1/2∂3u∂x3xi−1/2,tk+Oh3.Subtracting (13) and (14) and dividing the result by h, we obtain
(15)δxDδxUik=1hDi+1/2δxUi+1/2k−Di−1/2δxUi−1/2k=1hDi+1/2∂u∂xxi+1/2,tk−Di−1/2∂u∂xxi−1/2,tk+Oh2=∂∂xDx∂ux,t∂xxi,tk+Oh2,(16)fuxi,tk,uxi,tk−s1,uxi,tk−s2,xi,tk=f2Uik−1−Uik−2,Uik−m1−δ1,Uik−m2−δ2,xi,tk+τ2∂2u∂t2xi,ηikfμρk,Uik−m1−δ1,Uik−m2−δ2,xi,tk.where Uik−m1−δ1=1−δ1Uik−m1+δ1Uik−m1+1, Uik−m2−δ2=1−δ2Uik−m2+δ2Uik−m2+1, ηik∈tk−s2,tk, and ρk is a real number between uxi,tk and 2Uik−1−Uik−2.Substituting (11), (15), and (16) into (10), we obtain
(17)δtαUik−δxDδxUik=f2Uik−1−Uik−2,Uik−m1−δ1,Uik−m2−δ2,xi,tk+Rik,where
(18)Rik=rik+τ2∂2u∂t2xi,ηikfμρk,Uik−m1−δ1,Uik−m2−δ2,xi,tk+Oh2.From Assumptions (H1) and (H2), we can easily obtain(19)Rik≤CRτ2−α+h2,0≤i≤M,1≤k≤N.We have the discretization of (4) and (5):
(20)Uik=ϕxi,tk,0≤i≤M,−m2+1≤k≤0,(21)U0k=γtk,UMk=βtk,1≤k≤N.ReplacingUik by uik in (17), (20), and (21), and omitting Rik, we have the following numerical scheme:
(22)δtαuik−δxDδxuik=f2uik−1−uik−2,uik−m1−δ1,uik−m2−δ2,xi,tk,(23)uik=ϕxi,tk,0≤i≤M,−m2+1≤k≤0,(24)u0k=γtk,uMk=βtk,1≤k≤N.


## 3. The Solvability, Convergence, and Stability of the Difference Scheme

Assume the following to be defined onΩh(25)Vh,0=u∣u=u0,u1,⋯,uM,u0=uM=0.Ifu,v∈Vh,0, the following inner products and corresponding norms are introduced
(26)u,v=h∑i=1M−1uivi,u=u,u,u∞=max1≤i≤M−1ui,δxu,δxv=h∑i=0M−1δxui+1/2δxvi+1/2,δxu1=δxu,δxu,δxu,δxvD=h∑i=0M−1Dxi+1/2δxui+1/2δxvi+1/2,δxu1D=δxu,δxuD.We can obtain Lemma3 from Assumption (H1):Lemma 3.
For∀u∈Vh,0, we have c0δxu1≤δxu1D≤c1δxu1.Lemma 4 [30].
For∀u∈Vh,0, one can obtain
(27)u∞≤δxu1/2,u≤δxu1/6.We also introduce Lemma5 to be utilized in the following proof.Lemma 5 [30].
Assuming the sequenceFk≥0∣k≥0 satisfies
(28)Fk+1≤A+Bτ∑i=1kFl,k=0,1,⋯,then
(29)Fk+1≤AexpBkτ,k=0,1,2,⋯,where A≥0 and B≥0 are two constants.Theorem 6.
Under the assumptions of (H1) and (H2), the solution of the difference schemes (22)–(24) is unique.Proof.
Difference schemes (22)–(24) can be reformed as follows:
Whenk=1, we have
(30)−Di−1/2h2ui−11+1λ+Di−1/2+Di+1/2h2ui1−Di+1/2h2ui+11=1λui0+fi1,Whenk≥2, we have
(31)−Di−1/2h2ui−1k+1λ+Di−1/2+Di+1/2h2uik−Di+1/2h2ui+1k=1λ∑j=1k−1dk−j−1−dk−juij+dk−1ui0+fik,where λ=ταΓ2−α, fik=f2uik−1−uik−2,uik−m1−δ1,uik−m2−δ2,xi,tk,k=1,2,⋯. From (30) and (31), we can see that the schemes (22)–(24) are a linear tridiagonal system with strictly diagonally dominant coefficient matrix. Thus, schemes (22)–(31) have a unique solution.Denoteeik=Uik−uik,0≤i≤M,−m2+1≤k≤N, subtracting (22)–(24) from (17), (20), and (21), respectively, the following error equations can be obtained:
(32)δtαeik−δxDδxeik=pik+Rik,(33)eik=0,0≤i≤M,−m2+1≤k≤0,(34)e0k=0,eMk=0,1≤k≤N.where pik=f2Uik−1−Uik−2,Uik−m1−δ1,Uik−m2−δ2,xi,tk−f2uik−1−uik−2,uik−m1−δ1,uik−m2−δ2,xi,tk.Theorem 7.
Denoteux,t to be the solution of systems (3)–(5), and (22)–(24) has numerical solution uik∣0≤i≤M,−m2+1≤k≤N. Then, we have
(35)∥ek∥∞≤Cτ2−α+h2,1≤k≤N,where C is a positive constant independent of h and τ.Proof.
Multiplying (32) with hδtαeik, one can have
(36)h∑i=1M−1δtαeikδtαeik−h∑i=1M−1δxDδxeikδtαeik=h∑i=1M−1pikδtαeik+h∑i=1M−1Rikδtαeik,1≤i≤M−1,0≤k≤N−1.Then each term of (36) will be estimated.
(37)h∑i=1M−1δtαeikδtαeik=δtαeik2.By the discrete Green formula, we have(38)−h∑i=1M−1δxDδxeikδtαeik=h∑i=0M−1Di+1/2δxei+1/2kδtαδxei+1/2k=1λh∑i=0M−1Di+1/2δxei+1/2kδxei+1/2k−∑j=1k−1dk−j−1−dk−jδxei+1/2j−dk−1δxei+1/20=1λδxek1D2−∑j=1k−1dk−j−1−dk−jδxek,δxejD−dk−1δxek,δxe0D≥1λδxek1D2−∑j=1k−1dk−j−1−dk−jδxek1D2+δxej1D22−dk−1δxek1D2+δxe01D22=1λδxek1D22−∑j=1k−1dk−j−1−dk−jδxej1D22.From the Cauchy-Schwarz inequality, we have(39)h∑i=1M−1Rikδtαeik≤h∑i=1M−1Rik22ε+ε2δtαeik2=12ε∥Rik∥2+ε2∥δtαeik∥2.From the Cauchy-Schwarz inequality and Assumption (H2), we have(40)h∑i=1M−1pikδtαeik≤ε2∥δtαeik∥2+12εh∑i=1M−1pik2≤ε2∥δtαeik∥2+12εh∑i=1M−1c32eik−1−eik−2+c4eik−m1−δ1+c5eik−m2−δ22≤ε2∥δtαeik∥2+3c322εh∑i=1M−12eik−1−eik−22+3c422εh∑i=1M−1eik−m1−δ12+3c522εh∑i=1M−1eik−m2−δ22≤ε2∥δtαeik∥2+24c322ε∥ek−1∥2+∥ek−2∥2+3c422εmax∥ek−m1∥2,∥ek−m1+1∥2+3c522εmax∥ek−m2∥2,∥ek−m2+1∥2.Inserting (37)–(40) into (36), we obtain
(41)∥δtαeik∥2+1λδxek1D22−∑j=1k−1dk−j−1−dk−jδxej1D22≤12ε∥Rik∥2+ε2∥δtαeik∥2+ε2∥δtαeik∥2+24c322ε∥ek−1∥2+∥ek−2∥2+3c422εmax∥ek−m1∥2,∥ek−m1+1∥2+3c522εmax∥ek−m2∥2,∥ek−m2+1∥2.Multiplying (41) by 2λ, and letting ε=1, we have
(42)δxek1D2≤∑j=1k−1dk−j−1−dk−jδxej1D2+λ∥Rik∥2+24λc32∥ek−1∥2+∥ek−2∥2+3λc42max∥ek−m1∥2,∥ek−m1+1∥2+3λc52max∥ek−m2∥2,∥ek−m2+1∥2.From Lemmas3 and 4, and noticing (19), we obtain
(43)δxek12≤c1c0∑j=1k−1dk−j−1−dk−jδxej12+λc0CR2τ2−α+h22+4λc32c0δxek−112+δxek−212+λc422c0maxδxek−m112,δxek−m1+112+λc522c0maxδxek−m212,δxek−m2+112,denoting
(44)Ck=1c0Γ2−αmaxCR2,4c32,c42/2,c52/2>0.Noticingλ=ταΓ2−α, one can obtain τα<1 from 0<τ<1. By (43), one can have
(45)δxek12≤c1c0∑j=1k−1dk−j−1−dk−jδxej12+Ckτ2−α+h22+δxek−112+δxek−212+maxδxek−m112,δxek−m1+112+maxδxek−m212,δxek−m2+112.From Lemmas2 and 5, we have
(46)δxek12≤Ckexp4Ck+c1c0∑j=1k−1dk−j−1−dk−jτ2−α+h22=Ckexp4Ck+c1c01−dk−1τ2−α+h22≤C1τ2−α+h22,where C1=Ckexp4Ck+c1/c0. From Lemma 4, we have
(47)∥ek∥∞≤C1τ2−α+h22=Cτ2−α+h2.We complete the proof.To discuss the stability of the difference schemes (22)–(24), we consider the following problem:
(48)∂αv∂tα−Dxvxx=fvx,t,vx,t−s1,vx,t−s2,x,t,x,t∈0,1×0,T,(49)vx,t=ϕx,t+ψx,t,x∈0,1,t∈−s,0,(50)v0,t=γt,u1,t=βt,t∈0,T,where ψx,t is the perturbation caused by ϕx,t. The following difference scheme solving for (48)–(50) can be obtained:
(51)δtαuik−δxDδxuik=f2uik−1−uik−2,uik−m1−δ1,uik−m2−δ2,xi,tk,(52)uik=ϕxi,tk+ψik,0≤i≤M,−m2+1≤k≤0,(53)u0k=γtk,uMk=βtk,1≤k≤N.Denote(54)ωik=vik−uik,0≤i≤M,−m2+1≤k≤N.where C¯ is a bounded constant independent of h and τ.Definition 8.
Assumeuik satisfy (22)–(24) and vik satisfy (51)–(53), a numerical scheme for (3)–(5) is stable if we have
(55)ωk∞≤C¯max−m2+1≤j≤0ψj1,Similar to the proof of Theorem7, the following stability result can be obtained.Theorem 9.
Under the same condition as Theorem7, numerical schemes (22)–(24) satisfy ωk∞≤C^max−m2+1≤j≤0ψj1.


## 4. Numerical Test

Now, we present the following numerical tests to testify the schemes (22)–(24). Introducing the following notations,
(56)eh,τ=max1≤k≤NUk−uk∞,Rateτ=logeh,τ1/eh,τ2logτ1/τ2,Rateh=logeh1,τ/eh2,τlogh1/h2.when considering for Rateτ, h should be fixed and small enough. While considering for Rateh, τ should be fixed and small enough.
(57)∂αu∂tα−∂u∂xx2+1∂u∂x=ux,t−0.11+u2x,t−0.05+Gx,t,x,t∈0,1×0,T,ux,t=t2+αcos2πx,x∈0,1,t∈−0.1,0,u0,t=t2+α,u1,t=t2+α,t∈0,1.Example 10.
LetDx=x2+1.Equation (57) has exact solution ux,t=t2+αcos2πx, where
(58)Gx,t=Γ3+αΓ3t2+4π2x2+1t2+αcos2πx+4πxt2+αsin2πx−t−0.12+αcos2πx1+t−0.054+2αcos22πx.Table1 presents the maximum errors in the temporal directions, where α=0.4,0.6,0.8, and h=1/2000. From the results, we can draw a conclusion that the convergence order in the temporal directions is 2−α. Table 2 presents the maximum errors in the spatial directions, where α=0.4,0.6,0.8, and τ=1/2000. From the results, we can draw a conclusion that the convergence order in the spatial directions is 2.Table 1
Simulation results of Example10, where h=1/2000 and α=0.4,0.6,0.8.
τα=0.4α=0.6α=0.8eh,τRateτeh,τRateτeh,τRateτ1/504.244e-005∗1.752e-004∗6.464e-004∗1/1001.609e-0051.3996.903e-0051.3442.851e-0041.1811/1508.950e-0061.4463.963e-0051.3691.758e-0041.1921/2006.401e-0061.1652.724e-0051.3031.253e-0041.178Table 2
Simulation results of Example10, where τ=1/2000 and α=0.4,0.6,0.8.
hα=0.4α=0.6α=0.8eh,τRateheh,τRateheh,τRateh1/106.213e-002∗6.072e-002∗5.889e-002∗1/201.534e-0022.0181.500e-0022.0171.455e-0022.0171/403.842e-0031.9973.758e-0031.9973.654e-0031.9941/809.599e-0042.0019.396e-0042.0009.185e-0041.992Figure1 gives the error planes of Example 10, where τ=h=1/100, and α=0.1,0.3,0.7,0.9. Figure 1 tells us that bigger α can bring bigger error.Example 11.
LetDx=x2+1(59)∂αu∂tα−∂u∂xx2∂u∂x=−ux,t2+ux,t−0.1+ux,t−0.2+Gx,t,x,t∈0,1×0,T,ux,t=exx21−x2t2+α,x∈0,1,t∈−0.2,0,u0,t=0,u1,t=0,t∈0,1..Figure 1
Error planes of Example10, where τ=h=1/100.Equation (59) has exact solution ux,t=exx21−x2t2+α, where
(60)Gx,t=e2xx41−x4t4+2α−exx21−x2t−0.12+α−exx21−x2t−0.22+α+Γ3+α2t2exx21−x2−ext2+α6x2−18x3+5x4+8x5+x6.Tables3 and 4 show the computational results for a different α of Example 11. When h=1/2000, the maximum errors and convergence orders are provided in Table 3. It can be shown that the 2−α order accuracy in temporal direction is verified, while when τ=1/2000, from Table 4, we can see that the 2 order accuracy in spatial direction is verified.Table 3
Simulation results of Example11, where h=1/2000 and α=0.3,0.5,0.8.
τα=0.3α=0.5α=0.8eh,τRateτeh,τRateτeh,τRateτ1/109.029e-005∗2.277e-004∗9.049e-004∗1/202.712e-0051.7357.924e-0051.5233.941e-0041.1991/408.183e-0061.7292.763e-0051.5201.713e-0041.2021/802.530e-0061.6939.707e-0061.5097.451e-0051.201Table 4
Simulation results of Example11, where τ=1/2000 and α=0.3,0.5,0.8.
hα=0.3α=0.5α=0.8eh,τRateheh,τRateheh,τRateh1/88.064e-003∗7.917e-003∗7.624e-003∗1/162.100e-0031.9412.057e-0031.9451.971e-0031.9521/325.272e-0041.9945.163e-0041.9944.957e-0041.9911/641.319e-0041.9981.293e-0041.9981.250e-0041.987


## 5. Conclusion

This paper provides a finite difference scheme for solving a type of variable coefficient time fractional subdiffusion equation with multidelay. The unconditional stability and the global convergence of the scheme in the maximum norm are proved. Numerical experiments were provided to support the theoretical results and testify the efficiency of the difference scheme.



---

*Source: 6964159-2020-08-11.xml*