# Fractional Crank-Nicolson-Galerkin Finite Element Methods for Nonlinear Time Fractional Parabolic Problems with Time Delay

Authors: Lili Li; Mianfu She; Yuanling Niu
Journal: Journal of Function Spaces
 (2021)

---

## Abstract

A linearized numerical scheme is proposed to solve the nonlinear time-fractional parabolic problems with time delay. The scheme is based on the standard Galerkin finite element method in the spatial direction, the fractional Crank-Nicolson method, and extrapolation methods in the temporal direction. A novel discrete fractional Grönwall inequality is established. Thanks to the inequality, the error estimate of a fully discrete scheme is obtained. Several numerical examples are provided to verify the effectiveness of the fully discrete numerical method.


---

## Body

## 1. Introduction

In this paper, we consider the linearized fractional Crank-Nicolson-Galerkin finite element method for solving the nonlinear time-fractional parabolic problems with time delay(1)RDtαu−Δu=ft,ux,t,ux,t−τ,inΩ×0,T,ux,t=φx,t,inΩ×−τ,0,ux,t=0,on∂Ω×0,T,where Ω is a bounded convex and convex polygon in R2 (or polyhedron in R3) and τ is the delay term. RDtαu denotes the Riemann-Liouville fractional derivative, defined by
(2)RDtαu·,t=1Γ1−α∂∂t∫0tt−s−αu·,sds,0<α<1.The nonlinear fractional parabolic problems with time delay have attracted significant attention because of their wide range of applications in various fields, such as biology, physics, and engineering [1–9]. Recently, plenty of numerical methods were presented for solving the linear time-fractional diffusion equations. For instance, Chen et al. [10] used finite difference methods and the Kansa method to approximate time and space derivatives, respectively. Dehghan et al. [11] presented a fully discrete scheme based on the finite difference methods in time direction and the meshless Galerkin method in space direction and proved the scheme was unconditionally stable and convergent. Murio [12] and Zhuang and Liu [13] proposed a fully implicit finite difference numerical scheme and obtained unconditionally stability. Jin et al. [14] derived the time-fractional Crank-Nicolson scheme to approximate Riemann-Liouville fractional derivative. Li et al. [15] used a transformation to develop some new schemes for solving the time-fractional problems. The new schemes admit some advantages for both capturing the initial layer and solving the models with small parameter α. More studies can be found in [16–32].Recently, it has been one of the hot spots in the investigations of different numerical methods for the nonlinear time-fractional problems. For the analysis of the L1-type methods, we refer readers to the paper [33–40]. For the analysis of the convolution quadrature methods or the fractional Crank-Nicolson scheme, we refer to the recent papers [41–46]. The key role in the convergence analysis of the schemes is the fractional Grönwall-type inequations. However, as pointed out in [47–49], the similar fractional Grönwall-type inequations can not be directly applied to study the convergence of numerical schemes for the nonlinear time-fractional problems with delay.In this paper, we present a linearized numerical scheme for solving the nonlinear fractional parabolic problems with time delay. The time Riemann-Liouville fractional derivative is approximated by the fractional Crank-Nicolson-type time-stepping scheme, the spatial derivative is approximated by using the standard Galerkin finite element method, and the nonlinear term is approximated by the extrapolation method. To study the numerical behavior of the fully discrete scheme, we construct a novel discrete fractional type Grönwall inequality. With the inequality, we consider the convergence of the numerical methods for the nonlinear fractional parabolic problems with time delay.The rest of this article is organized as follows. In Section2, we present a linearized numerical scheme for the nonlinear time-fractional parabolic problems with delay and main convergence results. In Section 3, we present a detailed proof of the main results. In Section 4, numerical examples are given to confirm the theoretical results. Finally, the conclusions are presented in Section 5.


## 2. Fractional Crank-Nicolson-Galerkin FEMs

DenoteTh is a shape regular, quasiuniform triangulation of the Ω into d-simplexes. Let h=maxKεThdiamK. Let Xh be the finite-dimensional subspace of H01Ω consisting of continuous piecewise function on Th. Let Δt=τ/mτ be the time step size, where mτ is a positive integer. Denote N=T/Δt, tj=jΔt, j=−mτ,−mτ+1,⋯,0,1,2,⋯,N.The approximation to the Riemann-Liouville fractional derivative at pointt=tn−α/2 is given by [14]
(3)RDtn−α/2αux,t=Δt−α∑i=0nωn−iαux,ti+OΔt2≔RDtΔtαun+OΔt2,where
(4)ωiα=−1iΓα+1Γi+1Γα−i+1.For simplicity, denote∥v∥=∫Ωvx2dx1/2, ηn,α=1−α/2ηn+α/2ηn−1, η∧n,α=2−α/2ηn−1−1−α/2ηn−2, tnα=nΔtα.With the notation, the fully discrete scheme is to findUhnεXh such that
(5)RDΔtαUhn,v+∇Uhn,α,∇v=ftn−α/2,U^hn,α,Uhn−mτ,α,v,∀vεXh,n=1,2,⋯,N,and the initial condition
(6)Uhn=Rhφx,tn,n=−mτ,−mτ+1,⋯,0,where Rh:H01Ω→Xh is Ritz projection operator which satisfies the following equality [50]
(7)∇Rhu,∇v=∇u,∇v,∀uεH01Ω∩H2Ω,vεXh.We present the main convergence results here and leave their proof in the next section.Theorem 1.
Suppose the system (1) has a unique solution u satisfying
(8)∥u0∥Hr+1+∥u∥C0,T;Hr+1+∥ut∥C0,T;Hr+1+∥utt∥C0,T;H2+∥RDΔtαu∥C0,T;Hr+1≤K,and the source term ft,ux,t,ux,t−τ satisfies the Lipschitz condition
(9)ft,ux,t,ux,t−τ−ft,vx,t,vx,t−τ≤L1ux,t−vx,t+L2ux,t,τ−vx,t,τ,where K is a constant independent of n, h, and Δt and L1 and L2 are given positive constants. Then, there exists a positive constant Δt∗ such that for Δt≤Δt∗, the following estimate holds that
(10)∥un−Uhn∥≤C1∗Δt2+hr+1,n=1,2,⋯,N,where C1∗ is a positive constant independent of h and Δt.Remark 2.
The main contribution of the present study is that we obtain a discrete fractional Grönwall’s inequality. Thanks to the inequality, the convergence of the fully discrete scheme for the nonlinear time-fractional parabolic problems with delay can be obtained.Remark 3.
At present, the convergence of the proposed scheme is proved without considering the weak singularity of the solutions. In fact, if the initial layer of the problem is taken into account, there are some corrected terms at the beginning. Then, the scheme can be of order two in the temporal direction for nonsmooth initial data and some incompatible source terms. However, we still have the difficulties to get the similar discrete fractional Grönwall’s inequality. We hope to leave the challenging problems in the future.


## 3. Proof of the Main Results

In this section, we will present a detailed proof of the main result.


### 3.1. Preliminaries and Discrete Fractional Grönwall Inequality

Firstly, we review the definition of weightsωiα and denote gnα=∑i=0nωiα. Then, we can get
(11)ω0α=g0α,ωiα=giα−gi−1α,1≤i≤n.Actually, it has been shown [51] that ωiα and gnα process following properties:
(1)
The weightsωiα can be evaluated recursively, ωiα=1−α+1/iωi−1α,i≥1,ω0α=1(2)
The sequenceωiαi=0∞ are monotone increasing −1<ωiα<ωi+1α<0,i≥0(3)
The sequencegiαi=0∞ are monotone decreasing, giα>gi+1α for i≥0 and g0α=1Noticing the definition ofgiα, RDΔtαun can be rewritten as
(12)RDΔtαun=Δt−α∑i=1ngiα−gi−1αun−i+Δt−αg0αun.In fact, rearranging this identity yields(13)RDΔtαun=Δt−α∑i=1ngn−iαδtui+Δt−αgnαu0,where δtui=ui−ui−1.Lemma 4 (see [51]).
Consider the sequenceϕn given by
(14)ϕ0=1,ϕn=∑i=1ngi−1α−giαϕn−i,n≥1.
Then,ϕn satisfies the following properties:
(i)
0<ϕn<1,∑i=jnϕn−igi−jα=1,1≤j≤n(ii)
1/Γα∑i=1nϕn−i≤nα/Γ1+α(iii)
1/ΓαΓ1+k−1α∑i=1n−1ϕn−iik−1α≤nkα/Γ1+α,k=1,2⋯Lemma 5 (see [51]).
Consider the matrix(15)W=2μΔtα0ϕ1⋯ϕn−2ϕn−100⋯ϕn−3ϕn−2⋮⋮⋱⋮⋮00⋯0ϕ100⋯00n×n.
Then,W satisfies the following properties:
(i)
Wl=0,l≥n(ii)
Wke←≤1/Γ1+kα2Γαμtnαk,2Γαμtn−1αk,⋯,2Γαμt1αk′,k=0,1,2,⋯(iii)
∑k=0lWke←=∑k=0n−1Wke←≤Eα2Γαμtnα,Eα2Γαμtn−1α,⋯,Eα2Γαμt1α′,l≥nwhere e←=1,1,⋯,1′∈ℝn, μ is a constant.Theorem 6.
Assumingun∣n=−m,−m+1,⋯,0,1,2,⋯ and fn∣n=0,1,2,⋯ are nonnegative sequence, for λi>0, i=1,2,3,4,5, if
(16)RDΔtαuj≤λ1uj+λ2uj−1+λ3uj−2+λ4uj−m+λ5uj−m−1+fj,j=1,2⋯,
then, there exists a positive constantΔt∗, for Δt<Δt∗, the following holds
(17)un≤2λ4ΓαtnαΓ1+αM+λ5ΓαtnαΓ1+αM+max1≤j≤nfjΓαtnαΓ1+α+2M+λ2MΔtα+2λ3MΔtαEα2Γαλtnα,1≤n≤N,where λ=λ1+1/g0α−g1αg0α−g1αλ2+1/g1α−g2αλ3+1/gm−1α−gmαλ4+1/gmα−gm+1αλ5, Eαz=∑k=0∞zk/Γa+kα is the Mittag-Leffler function, and M=maxu−m,u−m+1,⋯,u0.Proof.
By using the definition ofRDΔtαun in (13), we have
(18)∑k=1jgj−kαδtuk+gjαu0≤Δtαλ1uj+λ2uj−1+λ3uj−2+λ4uj−m+λ5uj−m−1+Δtαfj.
Multiplying equation (18) by ϕn−j and summing the index j from 1 to n, we get
(19)∑j=1nϕn−j∑k=1jgj−kαδtuk≤Δtα∑j=1nϕn−jλ1uj+λ2uj−1+λ3uj−2+λ4uj−m+λ5uj−m−1+Δtα∑j=1nϕn−jfj−∑j=1nϕn−jgjαu0.
We change the order of summation and make use of the definition ofϕn−j to obtain
(20)∑j=1nϕn−j∑k=1jgj−kαδtuk=∑k=1nδtuk∑j=1kϕn−jgj−kα=∑k=1nδtuk=un−u0,and using Lemma 4, we have
(21)Δtα∑j=1nϕn−jfj≤Δtαmax1≤j≤nfj∑j=1nϕn−j≤Δtαmax1≤j≤nfjΓαnαΓ1+α=max1≤j≤nfjΓαtnαΓ1+α.
Noticinggjα is monotone decreasing and using Lemma 4, we have
(22)−∑j=1nϕn−jgjαu0≤∑j=1nϕn−jgjαu0≤u0∑j=1nϕn−jgj−1α=u0.
Substituting (20), (21), and (22) into (19), we can obtain
(23)un≤Δtα∑j=1nϕn−jλ1uj+λ2uj−1+λ3uj−2+λ4uj−m+λ5uj−m−1+2u0+max1≤j≤nfjΓαtnαΓ1+α.
Applying Lemma4, we have
(24)Δtα∑j=1mϕn−juj−m≤ΓαtnαΓ1+αM,Δtα∑j=1m+1ϕn−juj−m−1≤ΓαtnαΓ1+αM.
Therefore,(25)λ4Δtα∑j=1mϕn−juj−m+λ5Δtα∑j=1m+1ϕn−juj−m−1+2u0+λ2Δtαϕn−1u0+λ3Δtαϕn−1u−1+ϕn−2u0≤λ4ΓαtnαΓ1+αM+λ5ΓαtnαΓ1+αM+2M+λ2MΔtα+2λ3MΔtα.
Denote(26)Ψn=λ4ΓαtnαΓ1+αM+λ5ΓαtnαΓ1+αM+max1≤j≤nfjΓαtnαΓ1+α+2M+λ2MΔtα+2λ3MΔtα,
Equation (23) can be rewritten as
(27)1−λ1Δtαun≤λ1Δtα∑j=1n−1ϕn−juj+λ2Δtα∑j=2nϕn−juj−1+λ3Δtα∑j=3nϕn−juj−2+λ4Δtα∑j=m+1nϕn−juj−m+λ5Δtα∑j=m+2nϕn−juj−m−1+Ψn.
LetΔt∗=1/2λ1α, when Δt≤Δt∗, we have
(28)un≤2Ψn+2Δtαλ1∑j=1n−1ϕn−juj+λ2∑j=2nϕn−juj−1+λ3∑j=3nϕn−juj−2+λ4∑j=m+1nϕn−juj−m+λ5∑j=m+2nϕn−juj−m−1.
LetV=un,un−1,⋯,u1T, then (28) can be rewritten in the following matrix form:
(29)V≤2Ψne←+λ1W1+λ2W2+λ3W3+λ4W4+λ5W5V,where
(30)W1=2Δtα0ϕ1ϕ2⋯ϕn−2ϕn−100ϕ1⋯ϕn−3ϕn−2⋮⋮⋮⋱⋮⋮000⋯ϕ1ϕ2000⋯0ϕ10000⋯0n×n,W2=2Δtα0ϕ0ϕ1⋯ϕn−3ϕn−200ϕ0⋯ϕn−4ϕn−3⋮⋮⋮⋱⋮⋮000⋯ϕ0ϕ1000⋯0ϕ00000⋯0n×n,W3=2Δtα00ϕ0⋯ϕn−4ϕn−3000⋯ϕn−5ϕn−4⋮⋮⋮⋱⋮⋮000⋯0ϕ0000⋯00000000n×n,W4=2Δtα0⋯0ϕ0ϕ1⋯ϕn−m−2ϕn−m−10⋯00ϕ0⋯ϕn−m−3ϕn−m−2⋮⋯⋮⋮⋮⋱⋮⋮0⋯000⋯ϕ0ϕ10⋯000⋯0ϕ00000⋯00⋮⋯⋮⋮⋮⋯⋮⋮0⋯000⋯00n×n,W5=2Δtα0⋯00ϕ0⋯ϕn−m−3ϕn−m−20⋯000⋯ϕn−m−4ϕn−m−3⋮⋯⋮⋮⋮⋱⋮⋮0⋯000⋯0ϕ00⋯000⋯000000⋯00⋮⋯⋮⋮⋮⋯⋮⋮0⋯000⋯00n×n.
Since the definition ofϕn, we have
(31)ϕn−j≤1gj−1α−gjαϕn.
Then,(32)W2V≤1g0α−g1αW1V,W3V≤1g1α−g2αW1V,W4V≤1gm−1α−gmαW1V,W5V≤1gmα−gm+1αW1V.
Hence, (29) can be shown as follows
(33)V≤λ1+1g0α−g1αλ2+1g1α−g2αλ3+1gm−1α−gmαλ4+1gmα−gm+1αλ5W1V+2Ψne←=WV+2Ψne←,where W=λW1.
Therefore,(34)V≤WV+2Ψne←≤WWV+2Ψne←+2Ψne←=W2V+2Ψn∑j=01Wje←≤⋯≤WnV+2Ψn∑j=0n−1Wje←.
According to Lemma5, the result can be proved.Lemma 7 (see [51]).
For any sequenceekk=0N⊂Xh, the following inequality holds:
(35)RDΔtαek,1−α2ek+α2ek−1≥12RDΔtα∥ek∥2,1≤k≤N.Lemma 8 (see [52]).
There exists a positive constantCΩ, independent of h, for any v∈HsΩ∩H01Ω, such that
(36)∥v−Rhv∥L2+h∥∇v−Rhv∥L2≤CΩhs∥v∥Hs,1≤s≤r+1.


### 3.2. Proof of Theorem1

Now, we are ready to prove our main results.Proof.
Takingt=tn−α/2 in the first equation (1), we can find that un satisfies the following equation:
(37)RDΔtαun,v+∇un,α,∇v=ftn−α/2,u∧n,α,un−mτ,α,v+Pn,v,for n=1,2,3,⋯,N and ∀v∈Xh, where
(38)Pn=RDΔtαun−RDtn−α/2αu+Δun−α/2−Δun,α+ftn−α/2,un−α/2,un−mτ−α/2−ftn−α/2,u∧n,α,un−mτ,α.
Now, we estimate the error of∥Pn∥. Actually, from the definition of un,α and u∧n,α and the regularity of the exact solution (8), we can obtain that
(39)un−α/2−un,α=1−α2un−α/2+α2un−α/2−1−α2un−α/2un−1=1−α2un−α/2−un+α2un−α/2−un−1=−1−α2+α2Δtu′ξ1+1−α2+α2Δtu′ξ2=1−α2+α2Δtu′ξ2−u′ξ1≤1−α2+α2Δt∫tn−1tnuttsds≤C1Δt2,(40)∥un−α/2−u∧n,α∥=∥un−α/2−2−α2un−1+1−α2un−2∥=∥2−α2un−α/2−2−α2un−1+1−α2un−2−1−α2un−α/2∥=∥2−α2un−α/2−un−1+1−α2un−2−un−α/2∥=∥2−α21−α2Δtu′ξ3−2−α21−α2Δtu′ξ4∥=2−α21−α2Δt∥u′ξ3−u′ξ4∥≤2−α21−α2Δt∫tn−2tn−1∥utts∥ds≤C2Δt2,where ξ1∈tn−α/2,tn,ξ2∈tn−1,tn−α/2,ξ3∈tn−α/2,tn−1,ξ4∈tn−2,tn−α/2, C1=1−α/2α/2K, C2=2−α/21−α/2K are constants.
Applying (39) and (40) and the Lipschitz condition,
(41)∥ftn−α/2,un−α/2,un−mτ−α/2−ftn−α/2,u∧n,α,un−mτ,α∥≤L1C1+L2C2Δt2,∥Δun,α−un−α/2∥≤C1Δt2,which further implies that
(42)∥Pn∥≤CKΔt2,n=1,2,3,⋯,N,here CK=L1C1+L2C2.
Denoteθhn=Rhun−Uhn,n=0,1,⋯,N.
Substituting fully scheme (5) from equation (37) and using the property in (7), we can get that
(43)RDΔtαθhn,v+∇θhn,α,v=R1n,v+Pn,v−RDΔtαun−Rhun,v,where
(44)R1n=ftn−α/2,U^hn,α,Uhn−mτ,α−ftn−α/2,u∧n,α,un−mτ,α.
Settingv=θhn,α and applying the Cauchy-Schwarz inequality, it holds that
(45)RDΔtαθhn,θhn,α+∥∇θhn,α∥2≤∥R1n∥∥θhn,α∥+∥Pn∥∥θhn,α∥+∥RDΔtαun−Rhun∥∥θhn,α∥.
Noticing the factab≤1/2a2+b2 and ∥∇θhn,α∥2≥0,
(46)RDΔtαθhn,θhn,α≤12∥R1n∥2+∥Pn∥2+∥RDΔtαun−Rhun∥2+32∥θhn,α∥2.
Together with (9) and (36), we can arrive that
(47)∥RDΔtαun−Rhun∥≤CΩhr+1∥RDΔtαun∥Hr+1≤CΩKhr+1.(48)∥u∧n,α−Rhu∧n,α∥=∥2−α2un−1−1−α2un−2−2−α2Rhun−1+1−α2Rhun−2∥≤2−α2∥un−1−Rhun−1∥+1−α2∥un−2−Rhun−2∥≤2−α2CΩhr+1∥un−1∥Hr+1+1−α2CΩhr+1∥un−2∥Hr+1≤2−α2CΩKhr+1+1−α2CΩKhr+1≤C3hr+1.Similarly, we have
(49)∥un−mτ,α−Rhun−mτ,α∥=∥1−α2un−mτ+α2un−mτ−1−1−α2Rhun−mτ−α2Rhun−mτ−1∥≤1−α2CΩKhr+1+α2CΩKhr+1≤C4hr+1,where C3=22−α/2CΩK, C4=2max1−α/2,α/2CΩK.
Therefore,(50)∥R1n∥=∥ftn−α/2,u∧n,α,un−mτ,α−ftn−α/2,U^hn,α,Uhn−mτ,α∥≤L1∥u∧n,α−U^hn,α∥+L2∥un−mτ,α−Uhn−mτ,α∥≤L1∥θ^hn,α∥+L2∥θhn−mτ,α∥+L1∥u∧n,α−Rhu∧n,α∥+L2∥un−mτ,α−Rhun−mτ,α∥≤L1∥θ^hn,α∥+L2∥θhn−mτ,α∥+L1C3+L2C4hr+1.
Substituting (42), (47), and (50) into (46) and the fact a+b+c2≤3a2+3b2+3c2, we can get
(51)RDΔtαθhn,θhn,α≤32∥θhn,α∥2+3L122∥θ^hn,α∥2+3L222∥θhn−mτ,α∥2+CK22Δt4+123L12C32+L22C42+CKK2h2r+1≤32∥θhn,α∥2+3L122∥θ^hn,α∥2+3L222∥θhn−mτ,α∥2+C42Δt2+hr+12,where C4=maxCK2,3L12C32+L22C42+CKK2.
Applying Lemma7, we have
(52)RDΔtα∥θhn∥2≤3∥θhn,α∥2+3L12∥θ^hn,α∥2+3L22∥θhn−mτ,α∥2+C4Δt2+hr+12.
In terms of the definition of∥θhn,α∥ and θ^hn,α, we obtain
(53)RDΔtα∥θhn∥2≤31−α22∥θhn∥2+3α22+3L122−α22∥θhn−1∥2+3L121−α22∥θhn−2∥2+3L221−α22∥θhn−mτ∥2+3L22α22∥θhn−mτ−1∥2+C4Δt2+hr+12.
Using Theorem6, we can find a positive constant Δt∗ such that Δt≤Δt∗, then
(54)∥θhn∥2≤C5Δt2+hr+12,where C5 is a nonnegative constant which only depends on L1,L2,C4,CK,CΩ. In terms of the definition of θhn, we have
(55)∥un−Uhn∥≤∥un−Rhun∥+∥Rhun−Uhn∥≤C1∗Δt2+hr+1.Then, we complete the proof.


## 3.1. Preliminaries and Discrete Fractional Grönwall Inequality

Firstly, we review the definition of weightsωiα and denote gnα=∑i=0nωiα. Then, we can get
(11)ω0α=g0α,ωiα=giα−gi−1α,1≤i≤n.Actually, it has been shown [51] that ωiα and gnα process following properties:
(1)
The weightsωiα can be evaluated recursively, ωiα=1−α+1/iωi−1α,i≥1,ω0α=1(2)
The sequenceωiαi=0∞ are monotone increasing −1<ωiα<ωi+1α<0,i≥0(3)
The sequencegiαi=0∞ are monotone decreasing, giα>gi+1α for i≥0 and g0α=1Noticing the definition ofgiα, RDΔtαun can be rewritten as
(12)RDΔtαun=Δt−α∑i=1ngiα−gi−1αun−i+Δt−αg0αun.In fact, rearranging this identity yields(13)RDΔtαun=Δt−α∑i=1ngn−iαδtui+Δt−αgnαu0,where δtui=ui−ui−1.Lemma 4 (see [51]).
Consider the sequenceϕn given by
(14)ϕ0=1,ϕn=∑i=1ngi−1α−giαϕn−i,n≥1.
Then,ϕn satisfies the following properties:
(i)
0<ϕn<1,∑i=jnϕn−igi−jα=1,1≤j≤n(ii)
1/Γα∑i=1nϕn−i≤nα/Γ1+α(iii)
1/ΓαΓ1+k−1α∑i=1n−1ϕn−iik−1α≤nkα/Γ1+α,k=1,2⋯Lemma 5 (see [51]).
Consider the matrix(15)W=2μΔtα0ϕ1⋯ϕn−2ϕn−100⋯ϕn−3ϕn−2⋮⋮⋱⋮⋮00⋯0ϕ100⋯00n×n.
Then,W satisfies the following properties:
(i)
Wl=0,l≥n(ii)
Wke←≤1/Γ1+kα2Γαμtnαk,2Γαμtn−1αk,⋯,2Γαμt1αk′,k=0,1,2,⋯(iii)
∑k=0lWke←=∑k=0n−1Wke←≤Eα2Γαμtnα,Eα2Γαμtn−1α,⋯,Eα2Γαμt1α′,l≥nwhere e←=1,1,⋯,1′∈ℝn, μ is a constant.Theorem 6.
Assumingun∣n=−m,−m+1,⋯,0,1,2,⋯ and fn∣n=0,1,2,⋯ are nonnegative sequence, for λi>0, i=1,2,3,4,5, if
(16)RDΔtαuj≤λ1uj+λ2uj−1+λ3uj−2+λ4uj−m+λ5uj−m−1+fj,j=1,2⋯,
then, there exists a positive constantΔt∗, for Δt<Δt∗, the following holds
(17)un≤2λ4ΓαtnαΓ1+αM+λ5ΓαtnαΓ1+αM+max1≤j≤nfjΓαtnαΓ1+α+2M+λ2MΔtα+2λ3MΔtαEα2Γαλtnα,1≤n≤N,where λ=λ1+1/g0α−g1αg0α−g1αλ2+1/g1α−g2αλ3+1/gm−1α−gmαλ4+1/gmα−gm+1αλ5, Eαz=∑k=0∞zk/Γa+kα is the Mittag-Leffler function, and M=maxu−m,u−m+1,⋯,u0.Proof.
By using the definition ofRDΔtαun in (13), we have
(18)∑k=1jgj−kαδtuk+gjαu0≤Δtαλ1uj+λ2uj−1+λ3uj−2+λ4uj−m+λ5uj−m−1+Δtαfj.
Multiplying equation (18) by ϕn−j and summing the index j from 1 to n, we get
(19)∑j=1nϕn−j∑k=1jgj−kαδtuk≤Δtα∑j=1nϕn−jλ1uj+λ2uj−1+λ3uj−2+λ4uj−m+λ5uj−m−1+Δtα∑j=1nϕn−jfj−∑j=1nϕn−jgjαu0.
We change the order of summation and make use of the definition ofϕn−j to obtain
(20)∑j=1nϕn−j∑k=1jgj−kαδtuk=∑k=1nδtuk∑j=1kϕn−jgj−kα=∑k=1nδtuk=un−u0,and using Lemma 4, we have
(21)Δtα∑j=1nϕn−jfj≤Δtαmax1≤j≤nfj∑j=1nϕn−j≤Δtαmax1≤j≤nfjΓαnαΓ1+α=max1≤j≤nfjΓαtnαΓ1+α.
Noticinggjα is monotone decreasing and using Lemma 4, we have
(22)−∑j=1nϕn−jgjαu0≤∑j=1nϕn−jgjαu0≤u0∑j=1nϕn−jgj−1α=u0.
Substituting (20), (21), and (22) into (19), we can obtain
(23)un≤Δtα∑j=1nϕn−jλ1uj+λ2uj−1+λ3uj−2+λ4uj−m+λ5uj−m−1+2u0+max1≤j≤nfjΓαtnαΓ1+α.
Applying Lemma4, we have
(24)Δtα∑j=1mϕn−juj−m≤ΓαtnαΓ1+αM,Δtα∑j=1m+1ϕn−juj−m−1≤ΓαtnαΓ1+αM.
Therefore,(25)λ4Δtα∑j=1mϕn−juj−m+λ5Δtα∑j=1m+1ϕn−juj−m−1+2u0+λ2Δtαϕn−1u0+λ3Δtαϕn−1u−1+ϕn−2u0≤λ4ΓαtnαΓ1+αM+λ5ΓαtnαΓ1+αM+2M+λ2MΔtα+2λ3MΔtα.
Denote(26)Ψn=λ4ΓαtnαΓ1+αM+λ5ΓαtnαΓ1+αM+max1≤j≤nfjΓαtnαΓ1+α+2M+λ2MΔtα+2λ3MΔtα,
Equation (23) can be rewritten as
(27)1−λ1Δtαun≤λ1Δtα∑j=1n−1ϕn−juj+λ2Δtα∑j=2nϕn−juj−1+λ3Δtα∑j=3nϕn−juj−2+λ4Δtα∑j=m+1nϕn−juj−m+λ5Δtα∑j=m+2nϕn−juj−m−1+Ψn.
LetΔt∗=1/2λ1α, when Δt≤Δt∗, we have
(28)un≤2Ψn+2Δtαλ1∑j=1n−1ϕn−juj+λ2∑j=2nϕn−juj−1+λ3∑j=3nϕn−juj−2+λ4∑j=m+1nϕn−juj−m+λ5∑j=m+2nϕn−juj−m−1.
LetV=un,un−1,⋯,u1T, then (28) can be rewritten in the following matrix form:
(29)V≤2Ψne←+λ1W1+λ2W2+λ3W3+λ4W4+λ5W5V,where
(30)W1=2Δtα0ϕ1ϕ2⋯ϕn−2ϕn−100ϕ1⋯ϕn−3ϕn−2⋮⋮⋮⋱⋮⋮000⋯ϕ1ϕ2000⋯0ϕ10000⋯0n×n,W2=2Δtα0ϕ0ϕ1⋯ϕn−3ϕn−200ϕ0⋯ϕn−4ϕn−3⋮⋮⋮⋱⋮⋮000⋯ϕ0ϕ1000⋯0ϕ00000⋯0n×n,W3=2Δtα00ϕ0⋯ϕn−4ϕn−3000⋯ϕn−5ϕn−4⋮⋮⋮⋱⋮⋮000⋯0ϕ0000⋯00000000n×n,W4=2Δtα0⋯0ϕ0ϕ1⋯ϕn−m−2ϕn−m−10⋯00ϕ0⋯ϕn−m−3ϕn−m−2⋮⋯⋮⋮⋮⋱⋮⋮0⋯000⋯ϕ0ϕ10⋯000⋯0ϕ00000⋯00⋮⋯⋮⋮⋮⋯⋮⋮0⋯000⋯00n×n,W5=2Δtα0⋯00ϕ0⋯ϕn−m−3ϕn−m−20⋯000⋯ϕn−m−4ϕn−m−3⋮⋯⋮⋮⋮⋱⋮⋮0⋯000⋯0ϕ00⋯000⋯000000⋯00⋮⋯⋮⋮⋮⋯⋮⋮0⋯000⋯00n×n.
Since the definition ofϕn, we have
(31)ϕn−j≤1gj−1α−gjαϕn.
Then,(32)W2V≤1g0α−g1αW1V,W3V≤1g1α−g2αW1V,W4V≤1gm−1α−gmαW1V,W5V≤1gmα−gm+1αW1V.
Hence, (29) can be shown as follows
(33)V≤λ1+1g0α−g1αλ2+1g1α−g2αλ3+1gm−1α−gmαλ4+1gmα−gm+1αλ5W1V+2Ψne←=WV+2Ψne←,where W=λW1.
Therefore,(34)V≤WV+2Ψne←≤WWV+2Ψne←+2Ψne←=W2V+2Ψn∑j=01Wje←≤⋯≤WnV+2Ψn∑j=0n−1Wje←.
According to Lemma5, the result can be proved.Lemma 7 (see [51]).
For any sequenceekk=0N⊂Xh, the following inequality holds:
(35)RDΔtαek,1−α2ek+α2ek−1≥12RDΔtα∥ek∥2,1≤k≤N.Lemma 8 (see [52]).
There exists a positive constantCΩ, independent of h, for any v∈HsΩ∩H01Ω, such that
(36)∥v−Rhv∥L2+h∥∇v−Rhv∥L2≤CΩhs∥v∥Hs,1≤s≤r+1.


## 3.2. Proof of Theorem1

Now, we are ready to prove our main results.Proof.
Takingt=tn−α/2 in the first equation (1), we can find that un satisfies the following equation:
(37)RDΔtαun,v+∇un,α,∇v=ftn−α/2,u∧n,α,un−mτ,α,v+Pn,v,for n=1,2,3,⋯,N and ∀v∈Xh, where
(38)Pn=RDΔtαun−RDtn−α/2αu+Δun−α/2−Δun,α+ftn−α/2,un−α/2,un−mτ−α/2−ftn−α/2,u∧n,α,un−mτ,α.
Now, we estimate the error of∥Pn∥. Actually, from the definition of un,α and u∧n,α and the regularity of the exact solution (8), we can obtain that
(39)un−α/2−un,α=1−α2un−α/2+α2un−α/2−1−α2un−α/2un−1=1−α2un−α/2−un+α2un−α/2−un−1=−1−α2+α2Δtu′ξ1+1−α2+α2Δtu′ξ2=1−α2+α2Δtu′ξ2−u′ξ1≤1−α2+α2Δt∫tn−1tnuttsds≤C1Δt2,(40)∥un−α/2−u∧n,α∥=∥un−α/2−2−α2un−1+1−α2un−2∥=∥2−α2un−α/2−2−α2un−1+1−α2un−2−1−α2un−α/2∥=∥2−α2un−α/2−un−1+1−α2un−2−un−α/2∥=∥2−α21−α2Δtu′ξ3−2−α21−α2Δtu′ξ4∥=2−α21−α2Δt∥u′ξ3−u′ξ4∥≤2−α21−α2Δt∫tn−2tn−1∥utts∥ds≤C2Δt2,where ξ1∈tn−α/2,tn,ξ2∈tn−1,tn−α/2,ξ3∈tn−α/2,tn−1,ξ4∈tn−2,tn−α/2, C1=1−α/2α/2K, C2=2−α/21−α/2K are constants.
Applying (39) and (40) and the Lipschitz condition,
(41)∥ftn−α/2,un−α/2,un−mτ−α/2−ftn−α/2,u∧n,α,un−mτ,α∥≤L1C1+L2C2Δt2,∥Δun,α−un−α/2∥≤C1Δt2,which further implies that
(42)∥Pn∥≤CKΔt2,n=1,2,3,⋯,N,here CK=L1C1+L2C2.
Denoteθhn=Rhun−Uhn,n=0,1,⋯,N.
Substituting fully scheme (5) from equation (37) and using the property in (7), we can get that
(43)RDΔtαθhn,v+∇θhn,α,v=R1n,v+Pn,v−RDΔtαun−Rhun,v,where
(44)R1n=ftn−α/2,U^hn,α,Uhn−mτ,α−ftn−α/2,u∧n,α,un−mτ,α.
Settingv=θhn,α and applying the Cauchy-Schwarz inequality, it holds that
(45)RDΔtαθhn,θhn,α+∥∇θhn,α∥2≤∥R1n∥∥θhn,α∥+∥Pn∥∥θhn,α∥+∥RDΔtαun−Rhun∥∥θhn,α∥.
Noticing the factab≤1/2a2+b2 and ∥∇θhn,α∥2≥0,
(46)RDΔtαθhn,θhn,α≤12∥R1n∥2+∥Pn∥2+∥RDΔtαun−Rhun∥2+32∥θhn,α∥2.
Together with (9) and (36), we can arrive that
(47)∥RDΔtαun−Rhun∥≤CΩhr+1∥RDΔtαun∥Hr+1≤CΩKhr+1.(48)∥u∧n,α−Rhu∧n,α∥=∥2−α2un−1−1−α2un−2−2−α2Rhun−1+1−α2Rhun−2∥≤2−α2∥un−1−Rhun−1∥+1−α2∥un−2−Rhun−2∥≤2−α2CΩhr+1∥un−1∥Hr+1+1−α2CΩhr+1∥un−2∥Hr+1≤2−α2CΩKhr+1+1−α2CΩKhr+1≤C3hr+1.Similarly, we have
(49)∥un−mτ,α−Rhun−mτ,α∥=∥1−α2un−mτ+α2un−mτ−1−1−α2Rhun−mτ−α2Rhun−mτ−1∥≤1−α2CΩKhr+1+α2CΩKhr+1≤C4hr+1,where C3=22−α/2CΩK, C4=2max1−α/2,α/2CΩK.
Therefore,(50)∥R1n∥=∥ftn−α/2,u∧n,α,un−mτ,α−ftn−α/2,U^hn,α,Uhn−mτ,α∥≤L1∥u∧n,α−U^hn,α∥+L2∥un−mτ,α−Uhn−mτ,α∥≤L1∥θ^hn,α∥+L2∥θhn−mτ,α∥+L1∥u∧n,α−Rhu∧n,α∥+L2∥un−mτ,α−Rhun−mτ,α∥≤L1∥θ^hn,α∥+L2∥θhn−mτ,α∥+L1C3+L2C4hr+1.
Substituting (42), (47), and (50) into (46) and the fact a+b+c2≤3a2+3b2+3c2, we can get
(51)RDΔtαθhn,θhn,α≤32∥θhn,α∥2+3L122∥θ^hn,α∥2+3L222∥θhn−mτ,α∥2+CK22Δt4+123L12C32+L22C42+CKK2h2r+1≤32∥θhn,α∥2+3L122∥θ^hn,α∥2+3L222∥θhn−mτ,α∥2+C42Δt2+hr+12,where C4=maxCK2,3L12C32+L22C42+CKK2.
Applying Lemma7, we have
(52)RDΔtα∥θhn∥2≤3∥θhn,α∥2+3L12∥θ^hn,α∥2+3L22∥θhn−mτ,α∥2+C4Δt2+hr+12.
In terms of the definition of∥θhn,α∥ and θ^hn,α, we obtain
(53)RDΔtα∥θhn∥2≤31−α22∥θhn∥2+3α22+3L122−α22∥θhn−1∥2+3L121−α22∥θhn−2∥2+3L221−α22∥θhn−mτ∥2+3L22α22∥θhn−mτ−1∥2+C4Δt2+hr+12.
Using Theorem6, we can find a positive constant Δt∗ such that Δt≤Δt∗, then
(54)∥θhn∥2≤C5Δt2+hr+12,where C5 is a nonnegative constant which only depends on L1,L2,C4,CK,CΩ. In terms of the definition of θhn, we have
(55)∥un−Uhn∥≤∥un−Rhun∥+∥Rhun−Uhn∥≤C1∗Δt2+hr+1.Then, we complete the proof.


## 4. Numerical Examples

In this section, we give two examples to verify our theoretical results. The errors are all calculated in L2-norm.Example 1.
Consider the nonlinear time-fractional Mackey-Glass-type equation(56)RDtαux,y,t=Δux,y,t−2ux,y,t+ux,y,t−0.11+u2x,y,t−0.1+fx,y,t,x,y∈0,12,t∈0,1,ux,y,t=t2sinπxsinπy,x,y∈0,12,t∈−0.1,0,where
(57)fx,y,t=2t2−αΓ3−αsinπxsinπy+2t2π2sinπxsinπy−2t2sinπxsinπy−t−0.12sinπxsinπy1+t−0.12sinπxsinπy2.
The exact solution is given as(58)ux,t=t2sinπxsinπy.In order to test the convergence order in temporal direction, we fixedM=40 for α=0.4,α=0.6 and different N. Similarly, to obtain the convergence order in spatial direction, we fixed N=100 for α=0.4,α=0.6, and different M. Table 1 gives the errors and convergence orders in temporal direction by using the Q-FEM. Table 1 shows that the convergence order in temporal direction is 2. Similarly, Tables 2 and 3 give the errors and convergence orders in spatial direction by using the L-FEM and Q-FEM, respectively. These numerical results correspond to our theoretical convergence order.Example 2.
Consider the following nonlinear time-fractional Nicholson’s blowflies equation(59)RDtαux,y,z,t=Δux,y,z,t−2ux,y,z,t+ux,y,z,t−0.1exp−ux,y,z,t−0.1+fx,y,z,t,x,y,z∈0,13,t∈0,1,ux,y,z,t=t2sinπxsinπysinπz,x,y,z∈0,13,t∈−0.1,0,where
(60)fx,y,z,t=2t2−α/Γ3−αsinπxsinπysinπz+2t2π2−1sinπxsinπysinπz−−0.12sinπxsinπysinπzexp−t−0.12sinπxsinπysinπz,the exact solution is given as
(61)ux,t=t2sinπxsinπysinπz.Table 1
The errors and convergence orders in temporal direction by using Q-FEM.
Mα=0.4α=0.6ErrorsOrdersErrorsOrders51.6856e−03∗5.3999e−03∗102.9420e−042.51841.2503e−032.1106205.9619e−052.30303.0266e−042.0465401.3851e−052.10587.4700e−052.0185Table 2
The errors and convergence orders in spatial direction by using L-FEM.
Mα=0.4α=0.6ErrorsOrdersErrorsOrders57.2603e−02∗7.2065e−02∗101.9449e−021.90031.9297e−021.9009208.7594e−031.96738.6948e−031.9662404.9508e−031.98344.9180e−031.9807Table 3
The errors and convergence orders in spatial direction by using Q-FEM.
Mα=0.4α=0.6ErrorsOrdersErrorsOrders52.0750e−03∗2.0746e−03∗102.4888e−043.05962.5148e−043.0443207.3251e−053.01657.5802e−052.9577403.0946e−052.99523.4200e−052.7666In this example, in order to test the convergence order in temporal and spatial direction, we solve this problem by using the L-FEM withM=N and the Q-FEM with N=M3/2, respectively. Tables 4 and 5 show that the convergence orders in temporal and spatial direction are 2 and 3, respectively. The numerical results confirm our theoretical convergence order.Table 4
The errors and orders in temporal and spatial direction by using L-FEM.
Mα=0.4α=0.6ErrorsOrdersErrorsOrders58.3275e−02∗8.3375e−02∗102.2615e−021.88062.2732e−021.8749205.8356e−031.95435.8662e−031.9542401.4707e−031.98841.4784e−031.9884Table 5
The errors and orders in temporal direction and spatial direction by using Q-FEM.
Mα=0.4α=0.6ErrorsOrdersErrorsOrdersN=M3/286.7379e−04∗6.9141e−04∗103.1416e−043.04593.4945e−043.0579121.9415e−043.09681.9787e−043.1196141.1891e−043.18061.1992e−043.2485


## 5. Conclusions

We proposed a linearized fractional Crank-Nicolson-Galerkin FEM for the nonlinear fractional parabolic equations with time delay. A novel fractional Grönwall-type inequality is developed. With the help of the inequality, we prove convergence of the numerical scheme. Numerical examples confirm our theoretical results.



---

*Source: 9981211-2021-03-31.xml*