Fractional Crank-Nicolson-Galerkin finite element methods for nonlinear time fractional parabolic problems with time delay

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\"{o}nwall inequality is established. Thanks to the inequality, the error estimate of fully discrete scheme is obtained. Several numerical examples are provided to verify the effectiveness of the fully discrete numerical method.


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 u(x, t), u(x, t − τ )), in Ω × (0, T ], u(x, t) = ϕ(x, t), in Ω × (−τ, 0], u(x, t) = 0, on ∂Ω × (0, T ], (1.1) where Ω is a bounded convex and convex polygon in R 2 (or polyhedron in R 3 ), τ is the delay term. R D α t u denotes the Riemann-Liouville fractional derivative, defined by The nonlinear fractional parabolic problems with time delay have attracted significant attention because of their widely range of applications in various fields, such as biology, physics and engineering [1,2,3,4,5,6,7,8,9], etc. 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 full 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 [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,17,18,19,20,21,22,23,24,25,26,27,28,29,30,31,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,34,35,36,37,38,39,40]. For the analysis of the convolution quadrature methods or the fractional Crank-Nicolson scheme, we refer to the recent papers [41,42,43,44,45,46]. The key role in the convergence analysis of the schemes is the fractional Grönwall type inequations. However, as pointed out in [47,48,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 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 Section 2, 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) and the initial condition where R h : H 1 0 (Ω) → X h is Ritz projection operator which satisfies following equality [50] ∇R h u, ∇v = ∇u, ∇v , We present the main convergence results here and leave its proof in the next section.
Theorem 1 Suppose the system (1.1) has a unique solution u satisfying 5) and the source term f (t, u(x, t), u(x, t − τ )) satisfies the Lipschitz condition where K is a constant independent of n, h and ∆t, L 1 and L 2 are given positive constants.
Then there exists a positive constant ∆t * such that for ∆t ≤ ∆t * , the following estimate holds that u n − U n h ≤ C * 1 (∆t 2 + h r+1 ), n = 1, 2, · · · , N, where C * 1 is a positive constant independent of h and ∆t.

Remark 1
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 2 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, 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 term. However, we still have the difficulties to get the similar discrete fractional Grönwall's inequality. We hope to leave the challenging problems in future.

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

Preliminaries and discrete fractional Grönwall inequality
Firstly, we review the definition of weights ω Actually, it has been shown [51] that ω (α) i and g (α) n process following properties i , R D α ∆t u n can be rewritten as In fact, rearranging this identity yields Lemma 1 ([51]) Consider the sequence {φ n } given by Then {φ n } satisfies the following properties Then, W satisfies the following properties then there exists a positive constant ∆t * , for ∆t < ∆t * , the following holds is the Mittag-Leffler function, and M = max{u −m , u −m+1 , . . . , u 0 } Proof. By using the definition of R D α ∆t u n in (3.2), we have Multiplying the equation (3.3) by φ n−j and summing the index j from 1 to n, we get We change the order of summation and make use of the definition of φ n−j to obtain and using the Lemma 1, we have (3.6) Noticing g (α) j is monotone decreasing and using Lemma 1, we have  Therefore (3.8) can be rewritten as Let ∆t * = α 1 2λ 1 , when ∆t ≤ ∆t * , we have Let V = (u n , u n−1 , · · · , u 1 ) T , then (3.9) can be rewritten in the following matrix form where Since the definition of φ n , we have Then, Hence, (3.11) can be shown as follows where W = λW 1 . Therefore, According to Lemma 2, the result can be proved.

Proof of Theorem 1
Now, we are ready to prove our main results.
Proof. Taking t = t n− α 2 in the first equation (1.1) we can find that u n satisfies the following equation for n = 1, 2, 3, . . . , N and ∀v ∈ X h , where (3.14) Now, we estimate the error of P n . Actually, from the definition of u n,α andû n,α and the regularity of the exact solution (2.5), we can obtain that and (3.15) and (3.16) and the Lipschitz condition and which further implies that P n ≤ C K (∆t) 2 , n = 1, 2, 3, . . . , N, Denote θ n h = R h u n − U n h , n = 0, 1, . . . , N. Substituting fully scheme (2.2) from equation (3.13) and using the property in (2.4), we can get that where ,û n,α , u n−mτ ,α ). Setting v = θ n,α h and applying Cauchy-Schwarz inequality, it holds that Noticing the fact ab ≤ 1 2 (a 2 + b 2 ) and ∇θ n,α h 2 ≥ 0, (3.20) Together with (2.6) and (3.12), we can arrive that and û n,α − R hû n,α similarly, we have In terms of the definition of θ n,α h andθ n,α h , we obtain Using Theorem 2, we can find a positive constant ∆t * such that ∆t ≤ ∆t * , then where C 5 is a nonnegative constant which only dependents on L 1 , L 2 , C 4 , C K , C Ω . In terms of the definition of θ n h , we have Then, we complete the proof.

Numerical examples
In this section, we give two examples to verify our theoretical results. The errors are all calculated in L2-norm.
In this example, in order to test the convergence order in temporal and spatial direction, we solve this problem by using the L-FEM with M = N and the Q-FEM with N = M (3/2) , respectively. Table4 and Table5 show that the convergence orders in temporal and spatial direction are 2 and 3, respectively. The numerical results confirm our theoretical convergence order.