A Uniform Accuracy High-Order Finite Difference and FEM for Optimal Problem Governed by Time-Fractional Diffusion Equation

Abstract

In this paper, the time fractional diffusion equations optimal control problem is solved by $3-\alpha$ order with uniform accuracy scheme in time and finite element method (FEM) in space. For the state and adjoint state equation, the piecewise linear polynomials are used to make the space variables discrete, and obtain the semidiscrete scheme of the state and adjoint state. The priori error estimates for the semidiscrete scheme for state and adjoint state equation are established. Furthermore, the $3-\alpha$ order uniform accuracy scheme is used to make the time variable discrete in the semidiscrete scheme and construct the full discrete scheme for the control problems based on the first optimal condition and ‘first optimize, then discretize’ approach. The fully discrete scheme’s stability and truncation error are analyzed. Finally, two numerical examples are denoted to show that the theoretical analysis are correct.

1. Introduction

The optimal control problems governed by differential equation usually include objective functional, control variables and state variables that need to be optimized, in which control variables and state variables are coupled in the form of differential equations, which are usually called state equations. According to the different constraints imposed on control variables or state variables, differential equation optimal control problems can be divided into unconstrained problems, control constraint problems and state constraint problems. In recent decades, the research on optimal control of integer order differential equations has made great progress, and many well-known scholars have done a lot of research work in control theory and numerical algorithms.

In the last few years, with the rapid development of the fractional calculus theory and its application, the research on fractional equation constrained optimal control has attracted extensive attention of scholars. Fractional equation constrained optimal control has been widely used in engineering fields, such as groundwater pollution control. The goal of this problem is to keep the concentration of pollutants in groundwater within an allowable range in a given region while, at the same time, minimizing the cost. Its mathematical model can be expressed as a fractional order optimal control problem with point-by-point state constraints, in which the state variable represents the concentration of pollutants and satisfies the fractional convection diffusion equation.

In [1], the authors gave an effective numerical method of the fractional optimal control problems (FOCPs) involving a singular or non-singular kernel. The distributed-order FOCPs were studied by pseudo-spectral method in [2]. The delay fractional optimal control problems was solved by fractional-order Lagrange polynomials and the collocation method with convergence analysis in [3]. In [4], the generalized shifted Chebyshev polynomials were used to construct a numerical solution for fractional optimal control problems. The distributed order FOCPs was solved by Legendre spectral collocation method with convergence analysis in [5]. Researchers can read more references in this area, such as: pseudospectral method [6], Chebyshev cardinal functions [7,8], modified hat functions [9], generalized shifted Legendre polynomials [10], fractional Birkhoff interpolation [11], weighted Jacobi polynomials [12,13], shifted Jacobi orthonormal polynomials [14], generalized fractional-order Chebyshev wavelets [15], B-spline polynomials and operational matrix [16], Laplace transform and shifted Chebyshev-Gauss collocation method [17], fractional pseudospectral method [18], generalized fractional-order Bernoulli functions [19] and so on. In [20], the second order necessary optimality condition to FOCP was given. Based on the piecewise constant functions, tensor product finite element (FE) and finite difference (FD) method, the fully discrete scheme was constructed for FOCPs in [21]. The time FOCP was solved by FEM and projected gradient algorithm in [22]. The Galerkin spectral approximation and conjugate gradient optimization algorithm was used to solve the FOCP of distributed order in [23]. The $2-\alpha$ order FD-FE scheme was applied to solve the FOCPs in [24]. The fast primal dual active set algorithm was constructed for FOCP by the finite element approximation in [25]. The authors designed and analyzed solution techniques for a linear-quadratic optimal control problems governed by fractional Laplacian using semidiscrete approach and fully discrete Approach. Others derived a priori error estimates for both solution techniques [26]. The the parallel Crank-Nicolson scheme was implemented in time and gradient projection technique to solve the FOCPs in [27]. The spectral discretization was used to solve the FOCP governed with priori error estimates in [28]. The spectral Petrov-Galerkin method was investigated for the FOCP with error estimate in [29]. The wavelets method was constructed to solve the FOCPs by using the Chebyshev polynomials of 6-th kind in [30].The authors of [31] presented an indirect low computational complexity and flexible accuracy numerical approach for FOCP by using 2nd kind Chebyshev wavelets. The authors of [32] provided an effificient numerical solution to solve two-dimensional FOCP with variable order. The fast gradient projection method for FOCP was constructed in [33]. The efficient numerical scheme to solve FOCPs based on the Hermite scaling function with $L^2$-error estimates was presented in [34]. The modified numerical scheme was devoted to solving the FOCPs of variable order in the sense of Riemann Liouville or Caputo derivatives by the non-standard FD method in [35]. The FOCP is a research hotspot; readers who are interested in FOCPs can further refer to [36,37].

According to the existing literature, the FOCP is solved by FD-FE scheme in which the difference schemes in time only use the low $2-\alpha$ order numerical schemes in time discretization. In this paper, we will construct a novel FD-FE numerical scheme for FOCP with state constraint is constructed based on the uniform accuracy $3-\alpha$ order finite difference scheme and finite elements for temporal discretization and spatial discretization, respectively. The first order optimality condition of the FOCP is analyzed. The state and adjoint state’s priori error estimate are derived. Some numerical results are used to show the theoretical result.

The remainder of this paper is organised in five sections: In Section 2, we describe the optimality condition of FOCP. We construct the semi-discrete Galerkin finite element approximate solution for the FOCP, with the convergence analysis in Section 3. In Section 4, the full discrete scheme for the FOCP is constructed and the stability and truncation error of the scheme are analyzed roughly. We describe the conjugate gradient optimization algorithm and show some numerical experiments to validate our method in Section 5. In Section 6, we give some remarks for the FOCPs’s high-order numerical scheme.

2. Optimality Condition of FOCP

Let $\mathsf{\Omega }={[0,1]}^d,I=[0,T]$, ${\mathsf{\Omega }}_T=\mathsf{\Omega }\times I,{\mathsf{\Gamma }}_T=\partial \mathsf{\Omega }\times (0,T]$, where d is the dimension of space. We consider the following FOCP

$$\begin{array}{c}\hfill \underset{(u,q)}{min}J(u,q)=\frac{1}{2}\parallel u(x,t)-u_d{(x,t)\parallel }_{L^2\left({\mathsf{\Omega }}_T\right)}^2+\frac{\gamma }{2}{\parallel q(x,t)\parallel }_{L^2\left({\mathsf{\Omega }}_T\right)}^2,\end{array}$$

(1)

here $(u,q)\in V\times \widehat{W},V=H^{\frac{\alpha }{2}}(0,T;H_0^1\left(\mathsf{\Omega }\right)\bigcap H^2\left(\mathsf{\Omega }\right)),\widehat{W}=L^2\left({\mathsf{\Omega }}_T\right)\doteq L^2(0,T;L^2\left(\mathsf{\Omega }\right)),u(x,t)$ is governed by the following time fractional diffusion equation

$$\begin{array}{c}\hfill \left\{\begin{array}{c}{}_0{D}_t^{\alpha }u(x,t)-\Delta u(x,t)=q(x,t)+f(x,t),(x,t)\in {\mathsf{\Omega }}_T,\hfill \\ u(x,t)=0,(x,t)\in {\mathsf{\Gamma }}_T,\hfill \\ u(x,0)=0,x\in \mathsf{\Omega },\hfill \end{array}\right.\end{array}$$

(2)

where ${}_0{D}_t^{\alpha }u(x,t)$ is $\alpha (0<\alpha <1)$ order fractional derivative of left Caputo to the state $u(x,t)$ with regard to time variable t defined as following [38]

$$\begin{array}{c}\hfill {}_0{D}_t^{\alpha }u(x,t)=\frac{1}{\mathsf{\Gamma }(1-\alpha )}{\int }_0^t\frac{{\partial }_{\tau }u(x,t)}{{(t-\tau )}^{\alpha }}d\tau .\end{array}$$

(3)

According to the existence and uniqueness of solution in [39] to Equation (2), it can be see that there exists a mapping $q\to u=u\left(q\right)$ defined by (2). It is easy to see that the cost function become $J\left(q\right)\doteq J(q,u\left(q\right)),q\in \widehat{W}$.

Then the FOCP (1) is equivalent to find $q^{*}\in \widehat{W}$, satisfy that

$$\begin{array}{c}\hfill J\left(q^{*}\right):=\underset{q}{min}J\left(q\right),\end{array}$$

(4)

where $q\in \widehat{W}$. The problem (4)’s first order necessary optimality condition is defined as follows

$$\begin{array}{c}\hfill J^{{}^{\prime }}\left(q^{*}\right)\left(\delta q\right)=0,\forall \delta q\in L^2\left({\mathsf{\Omega }}_T\right),\end{array}$$

(5)

where $J^{{}^{\prime }}\left(q^{*}\right)$ is the G$\widehat{a}$teaux differential of $J\left(q\right)$ at $q^{*}$ in the direction $\delta q$.

In the following Lemma 1, we will give the calculation of the G$\widehat{a}$teaux differential of $J\left(q\right)$. Based on the idea of Lemma 2.1 in [23], the following lemma is easy to prove. For the convenience of the readers, we give the details of proof as follows.

Lemma 1.

The gradient of $J\left(q\right)$ is determined by the following equation

$$\begin{array}{c}\hfill J^{{}^{\prime }}\left(q\right)\left(\delta q\right)={(Z\left(q\right)+\gamma q,\delta q)}_{{\mathsf{\Omega }}_T},\forall \delta q\in L^2\left({\mathsf{\Omega }}_T\right),\end{array}$$

(6)

where $Z\left(q\right)$ is defined by the adjoint state equation as follows

$$\begin{array}{c}\hfill \left\{\begin{array}{c}{}_t{D}_T^{\alpha }Z(x,t)-\Delta Z(x,t)=u(x,t)-u_d(x,t),(x,t)\in {\mathsf{\Omega }}_T,\hfill \\ Z(x,t)=0,(x,t)\in {\mathsf{\Gamma }}_T,\hfill \\ Z(x,T)=0,x\in \mathsf{\Omega },\hfill \end{array}\right.\end{array}$$

(7)

here ${}_t{D}_T^{\alpha }Z(x,t)$ is denoted by the α order right Caputo fractional derivative as follows

$$\begin{array}{c}\hfill {}_t{D}_T^{\alpha }Z(x,t)=-\frac{1}{\mathsf{\Gamma }(1-\alpha )}{\int }_t^T\frac{{\partial }_{\tau }Z(x,\tau )}{{(\tau -t)}^{\alpha }}d\tau .\end{array}$$

(8)

Proof. 

Using the chain rule and direct calculation, we can obtain that

$$\begin{array}{ccc}\hfill J^{{}^{\prime }}\left(q\right)\left(\delta q\right)& =& \gamma q\left(\delta q\right)+(u-u_d)u^{{}^{\prime }}\left(q\right)\left(\delta q\right)\hfill \\ \hfill & =& {\int }_{{\mathsf{\Omega }}_T}\gamma q·\gamma qdxdt+{\int }_{{\mathsf{\Omega }}_T}(u-u_d)u^{{}^{\prime }}\left(q\right)\left(\delta q\right)dxdt.\hfill \end{array}$$

(9)

For the $u^{{}^{\prime }}\left(q\right)\left(\delta q\right)$, we will give detailed calculation process in the following parts. Firstly, we denote as the $\delta u(x,t)$ as the derivative of $u\doteq u\left(q\right)$ in the direction $\delta q$ as following

$$\begin{array}{c}\hfill \delta u(x,t):=u^{{}^{\prime }}\left(q\right)\left(\delta q\right)=\underset{\overline{\eta }\to 0}{lim}\frac{u(q+\overline{\eta }\delta q)-u\left(q\right)}{\overline{\eta }}.\end{array}$$

It is easy to check that $\delta u$ satisfy the following problem

$$\begin{array}{c}\hfill \left\{\begin{array}{c}{}_0{D}_t^{\alpha }\delta u(x,t)-\Delta \delta u(x,t)=\delta q(x,t),\forall (x,t)\in {\mathsf{\Omega }}_T,\hfill \\ \delta u(x,t)=0,(x,t)\in {\mathsf{\Gamma }}_T,\hfill \\ \delta u(x,0)=0,x\in \mathsf{\Omega }.\hfill \end{array}\right.\end{array}$$

(10)

In order to obtain the first order necessary optimality condition for (5), we multiply both sides of the first Equation (7) by $\delta u$, and then integrate the first equation of (7) on domain ${\mathsf{\Omega }}_T$ to get the following equation

$$\begin{array}{c}\hfill {\int }_{{\mathsf{\Omega }}_T}(u-u_d)\delta udxdt={\int }_{{\mathsf{\Omega }}_T}{(}_t{D}_T^{\alpha }Z(x,t)-\Delta Z(x,t))\delta udxdt.\end{array}$$

(11)

Using the boundary conditions in (7) and (10), we can easily prove that

$$\begin{array}{c}\hfill {\int }_{{\mathsf{\Omega }}_T}\Delta Z(x,t)\delta udxdt={\int }_{{\mathsf{\Omega }}_T}Z(x,t)\Delta \delta udxdt.\end{array}$$

(12)

Based on the fractional integration by parts of [40], it is easy to check that

$$\begin{array}{ccc}\hfill & & {\int }_{{\mathsf{\Omega }}_T}{}_t{D}_T^{\alpha }Z(x,t)\delta udxdt={\int }_{{\mathsf{\Omega }}_T}{(}_t{D}_T^{\alpha }Z(x,t)-\frac{Z(x,T)}{\mathsf{\Gamma }(1-\alpha ){(T-t)}^{\alpha }})\delta udxdt\hfill \\ & =& {\int }_{{\mathsf{\Omega }}_T}{}_t^R{D}_T^{\alpha }Z(x,t)\delta udxdt={\int }_{{\mathsf{\Omega }}_T}Z(x,t){}_0^R{D}_t^{\alpha }\delta udxdt\hfill \\ & =& {\int }_{{\mathsf{\Omega }}_T}Z(x,t){}_0{D}_t^{\alpha }\delta udxdt+{\int }_{{\mathsf{\Omega }}_T}\frac{Z\delta u(x,t_0)}{\mathsf{\Gamma }(1-\alpha )t^{\alpha }}dxdt={\int }_{{\mathsf{\Omega }}_T}Z(x,t){}_0{D}_t^{\alpha }\delta udxdt,\hfill \end{array}$$

(13)

where ${}_0^R{D}_t^{\alpha }Z$ and ${}_t^R{D}_T^{\alpha }Z$ are left and right Riemann Liouville fractional derivative, respectively.

Based on the results of (12) and (13), and simple calculation, it is easy to obtain the following

$$\begin{array}{c}\hfill {\int }_{{\mathsf{\Omega }}_T}(u-u_d)\delta udxdt={\int }_{{\mathsf{\Omega }}_T}({}_0{D}_t^{\alpha }\delta u-\Delta \delta u)Zdxdt={\int }_{\mathsf{\Omega }}\delta q·Zdxdt.\end{array}$$

(14)

This together with (9), lead to (6). □

3. Semidiscrete Scheme for FOCP

In this part, we will construct the semidiscrete Galerkin FEM solution of FOCP (1) and (2) and analysis the error analysis of the semidiscrete scheme.

The state equation’s weak formulation is defined as following

$$\begin{array}{c}\hfill {(}_0{D}_t^{\alpha }u\left(q\right),v)+(\nabla u\left(q\right),\nabla v)=(q+f,v),\forall v\in H_0^1\left(\mathsf{\Omega }\right).\end{array}$$

(15)

Divide the domain $\mathsf{\Omega }$ into quasi-uniform FEM partitions ${\mathcal{T}}_h$, where h denote the maximum diameter. The FE space ${\mathcal{V}}_h$ on ${\mathcal{T}}_h$ is defined by.

$$\begin{array}{c}\hfill {\mathcal{V}}_h=\{v_h\in H_0^1\left(\mathsf{\Omega }\right)\cap C\left(\mathsf{\Omega }\right)|v_h\mathrm{is}\mathrm{a}\mathrm{linear}\mathrm{function}\mathrm{over}K,\forall K\in {\mathcal{T}}_h\}.\end{array}$$

The semidiscrete scheme to (15) reads: find $u_h\left(q\right)\in {\mathcal{V}}_h$ satisfying the following equation

$$\begin{array}{c}\hfill {(}_0{D}_t^{\alpha }{u}_h\left(q\right),v_h)+(\nabla u_h\left(q\right),\nabla v_h)=(q+f,v_h),\forall v_h\in {\mathcal{V}}_h.\end{array}$$

(16)

In the following Lemma 2, we will give the error of $u\left(q\right)-u_h\left(q\right)$. Because the following lemma can be easily proved by the idea of [41], we only give the conclusion of the following lemma and omit the tedious proof process.

Lemma 2.

Let $u\left(q\right)$ be the solution of (15) and $u_h\left(q\right)$ be the solution of (16), then

$$\begin{array}{c}\hfill \parallel u\left(q\right)-u_h{\left(q\right)\parallel }_{L^2\left({\mathsf{\Omega }}_T\right)}+h\parallel \nabla (u\left(q\right)-u_h\left(q\right)){\parallel }_{L^2\left({\mathsf{\Omega }}_T\right)}\le C{h}^2{\parallel q+f\parallel }_{L^2\left({\mathsf{\Omega }}_T\right)}.\end{array}$$

The semidiscrete scheme to (1) and (2) can be charactierized as

$$\begin{array}{c}\hfill \underset{(u_h,q_h)}{min}J(u_h,q_h)\end{array}$$

(17)

subject to

$$\begin{array}{c}\hfill \left\{\begin{array}{c}{(}_0{D}_t^{\alpha }{u}_h,v_h)+(\nabla u_h,\nabla v_h)=(q_h+f,v_h),\forall v_h\in {\mathcal{V}}_h,\hfill \\ u_h\left(0\right)=0.\hfill \end{array}\right.\end{array}$$

(18)

where $u_h,q_h$ are the finite element solution to u and q, respectively. Similar to problem (1), it is easy to check that

$$\begin{array}{c}\hfill J^{{}^{\prime }}\left(\delta q_h\right)={(\gamma q_h+Z_h,\delta q_h)}_{L^2\left(\mathsf{\Omega }\right)}=0,\forall \delta q_h\in {\mathcal{V}}_h.\end{array}$$

(19)

where $Z_h$ is the following discrete adjoint equation

$$\begin{array}{c}\hfill \left\{\begin{array}{c}{(}_t{D}_T^{\alpha }{Z}_h,v_h)+(\nabla Z_h,\nabla v_h)=(u_h-u_d,v_h),\forall v_h\in {\mathcal{V}}_h,\hfill \\ Z_h\left(T\right)=0.\hfill \end{array}\right.\end{array}$$

(20)

In order to analysis the convergence (18) and (19), we introduce $L^2$ projection ${\mathcal{P}}_h$ and Ritz projection ${\mathcal{R}}_h$ defined by

$$\begin{array}{c}\hfill \left\{\begin{array}{c}(\varphi -{\mathcal{P}}_h\varphi ,w_h)=0,\forall w_h\in {\mathcal{V}}_h,\hfill \\ (\nabla (\varphi -{\mathcal{R}}_h\varphi ),\nabla w_h)=0,\forall w_h\in {\mathcal{V}}_h.\hfill \end{array}\right.\end{array}$$

(21)

Next, we will give the error of ${\mathcal{P}}_h\varphi -\varphi$ and ${\mathcal{R}}_h\varphi -\varphi$ in Lemma 3. Based on the method in [41], we can easily prove the results in the following Lemma 3. Therefore, we omit the proof details and directly give the results as follows.

Lemma 3.

The projections ${\mathcal{P}}_h$, ${\mathcal{R}}_h$ are $L^2$-projection, Ritz-projection, respectively, and satisfy the following inequalities

$$\begin{array}{ccc}& & \parallel {\mathcal{P}}_h{\varphi -\varphi \parallel }_{L^2\left(\mathsf{\Omega }\right)}\le C{h}^2{\parallel \varphi \parallel }_{H^2\left(\mathsf{\Omega }\right)},\hfill \\ & & \parallel {\mathcal{R}}_h{\varphi -\varphi \parallel }_{L^2\left(\mathsf{\Omega }\right)}+h\parallel \nabla ({\mathcal{R}}_h\varphi -\varphi ){\parallel }_{L^2\left(\mathsf{\Omega }\right)}\le C{h}^2{\parallel \varphi \parallel }_{H^2\left(\mathsf{\Omega }\right)}.\hfill \end{array}$$

In order to obtain the semidiscrete scheme’s priori error estimates for (18) and (20), we define the following two auxiliary problems

$$\begin{array}{c}\hfill \left\{\begin{array}{c}{(}_0{D}_t^{\alpha }u\left(q_h\right),v)+(\nabla u\left(q_h\right),\nabla v)=(q_h+f,v),\forall v\in H_0^1\left(\mathsf{\Omega }\right),\hfill \\ u\left(q_h\right)\left(0\right)=0,\hfill \end{array}\right.\end{array}$$

(22)

and

$$\begin{array}{c}\hfill \left\{\begin{array}{c}{(}_t{D}_T^{\alpha }Z\left(q_h\right),v)+(\nabla Z\left(q_h\right),\nabla v)=(u\left(q_h\right)-u_d,v),\forall v\in H_0^1\left(\mathsf{\Omega }\right),\hfill \\ Z\left(q_h\right)\left(T\right)=0.\hfill \end{array}\right.\end{array}$$

(23)

It is obvious to find that $u_h$ is the semidiscrete scheme of $u\left(q_h\right)$. Based on the Lemma 2, one can immediately obtain the following

$$\begin{array}{c}\hfill \parallel u\left(q_h\right)-u_h{\parallel }_{L^2\left({\mathsf{\Omega }}_T\right)}+h\parallel \nabla (u\left(q_h\right)-u_h){\parallel }_{L^2\left({\mathsf{\Omega }}_T\right)}\le C{h}^2{\parallel f+q_h\parallel }_{L^2\left({\mathsf{\Omega }}_T\right)}.\end{array}$$

(24)

Compare (15) and (22), it is available to get

$$\begin{array}{c}\hfill {(}_0{D}_t^{\alpha }(u\left(q_h\right)-u\left(q\right)),v_h)+(\nabla (u\left(q_h\right)-u\left(q\right)),\nabla v_h)=(q_h-q,v_h),\forall v\in H_0^1\left(\mathsf{\Omega }\right).\end{array}$$

Using the stability estimates of the state equation, we have

$$\begin{array}{c}\hfill \parallel u\left(q_h\right)-u\left(q\right){\parallel }_{L^2\left({\mathsf{\Omega }}_T\right)}+\parallel \nabla (u\left(q_h\right)-u\left(q\right)){\parallel }_{L^2\left({\mathsf{\Omega }}_T\right)}\le Ch{\parallel q_h-q\parallel }_{L^2\left({\mathsf{\Omega }}_T\right)}.\end{array}$$

(25)

From (24) and (25), it is easy to find that the estimates of the state variable is dependent on control variable.

In the next Lemma 4, we will estimate the error of $q-q_h$.

Lemma 4.

Let $(u,Z,q)$ and $(u_h,Z_h,q_h)$ be the solutions of (2), (5), (7) and (18), (19), (20), respectively. Then the estimate following as

$$\begin{array}{c}\hfill \parallel q-q_h{\parallel }_{L^2\left({\mathsf{\Omega }}_T\right)}\le C{\parallel Z\left(q_h\right)-Z_h\parallel }_{L^2\left({\mathsf{\Omega }}_T\right)}\end{array}$$

holds.

Proof. 

Based on (6) and (19), it is obvious that

$$\begin{array}{ccc}\hfill \gamma \parallel q-q_h{\parallel }_{L^2\left({\mathsf{\Omega }}_T\right)}^2& =& {\int }_{{\mathsf{\Omega }}_T}\gamma q(q-q_h)dxdt-{\int }_{{\mathsf{\Omega }}_T}\gamma q_h(q-q_h)dxdt\hfill \\ & =& {\int }_{{\mathsf{\Omega }}_T}Z(q_h-q)dxdt+{\int }_{{\mathsf{\Omega }}_T}{Z}_h(q-q_h)dxdt.\hfill \end{array}$$

Based on (2) and (22) and using Green formulation of [40], we can easily obtain that

$$\begin{array}{ccc}& & {\int }_{{\mathsf{\Omega }}_T}(Z-Z\left(q_h\right))(q_h-q)dxdt\hfill \\ & =& {\int }_{{\mathsf{\Omega }}_T}{}_0{D}_t^{\alpha }(u\left(q_h\right)-u)(Z-Z\left(q_h\right))dxdt+{\int }_{{\mathsf{\Omega }}_T}\nabla (u\left(q_h\right)-u)·\nabla (Z-Z\left(q_h\right)dxdt\hfill \\ & =& {\int }_{{\mathsf{\Omega }}_T}{}_t{D}_T^{\alpha }(Z-Z\left(q_h\right))(u\left(q_h\right)-u)dxdt+{\int }_{{\mathsf{\Omega }}_T}\nabla (Z-Z\left(q_h\right)·\nabla (u\left(q_h\right)-u)dxdt\hfill \\ & =& -{\int }_{{\mathsf{\Omega }}_T}{(u\left(q_h\right)-u)}^{2}dxdt\le 0.\hfill \end{array}$$

Thus we carries at

$$\begin{array}{ccc}& & \gamma \parallel q-q_h{\parallel }_{L^2\left({\mathsf{\Omega }}_T\right)}^2\le {\int }_{{\mathsf{\Omega }}_T}(Z\left(q_h\right)-Z_h)(q_h-q)dxdt\hfill \\ & \le & C\left(\delta \right)\parallel Z\left(q_h\right)-Z_h{\parallel }_{L^2\left({\mathsf{\Omega }}_T\right)}^2+C\delta {\parallel q-q_h\parallel }_{L^2\left({\mathsf{\Omega }}_T\right)}^2.\hfill \end{array}$$

Choosing $\delta =\frac{\gamma }{2C}$ yields the final result. □

Next, we will obtain the estimate of $\parallel q-q_h{\parallel }_{L^2\left({\mathsf{\Omega }}_T\right)}$. Firstly, we will estimate of $\parallel Z\left(q_h\right)-Z_h{\parallel }_{L^2\left({\mathsf{\Omega }}_T\right)}$. In order to estimate $\parallel Z\left(q_h\right)-Z_h{\parallel }_{L^2\left({\mathsf{\Omega }}_T\right)}$, we introduce another auxiliary problem

$$\begin{array}{c}\hfill \left\{\begin{array}{c}{(}_t{D}_T^{\alpha }Z\left(u_h\right),v)+(Z\left(u_h\right),\nabla v)=(u_h-u_d,v),\forall v\in H_0^1\left(\mathsf{\Omega }\right),\hfill \\ Z\left(u_h\right)\left(T\right)=0.\hfill \end{array}\right.\end{array}$$

(26)

It is easy to check that $Z_h$ is the semidiscrete FEM solution of $Z\left(u_h\right)$. Based on the method of [41], we have the following error estimate of $Z\left(u_h\right)-Z_h$.

Lemma 5.

If $Z_h$ and $Z\left(u_h\right)$ are the solution of (20) and (23), respectively. Then the estimate is as follows

$$\begin{array}{c}\hfill \parallel Z_h-Z\left(u_h\right){\parallel }_{L^2\left({\mathsf{\Omega }}_T\right)}+h\parallel \nabla (Z_h-Z\left(u_h\right)){\parallel }_{L^2\left({\mathsf{\Omega }}_T\right)}\le C{h}^2{\parallel u_h-u_d\parallel }_{L^2\left({\mathsf{\Omega }}_T\right)}\end{array}$$

holds.

Based on the idea of [41], we split the error $Z_h-Z\left(u_h\right)$ into

$$\begin{array}{c}\hfill Z_h-Z\left(u_h\right)=Z_h-{\mathcal{P}}_{h}Z\left(u_h\right)+{\mathcal{P}}_{h}Z\left(u_h\right)-Z\left(u_h\right)\doteq \eta +\theta .\end{array}$$

Based on the Lemma 3, one can immediately obtain the error estimate of $\theta$. In the next, we only need to estimate $\eta$. First of all, one can infer from ${\Delta }_h{\mathcal{R}}_h={\mathcal{P}}_h\Delta$, where ${\Delta }_h$ is the discrete Laplace operator ${\Delta }_h:{\mathcal{V}}_h\to {\mathcal{V}}_h$:

$$\begin{array}{c}\hfill -({\Delta }_h\varphi ,v_h)=(\nabla \varphi ,\nabla v_h),\forall \varphi ,v_h\in {\mathcal{V}}_h.\end{array}$$

Therefore,

$$\begin{array}{ccc}& & {}_t{D}_T^{\alpha }\eta -{\Delta }_h\eta ={\mathcal{P}}_h(u_h-u_d){-}_t{D}_T^{\alpha }{\mathcal{P}}_{h}Z\left(u_h\right)-{\Delta }_h{\mathcal{P}}_{h}Z\left(u_h\right)\hfill \\ & =& {\Delta }_h({\mathcal{P}}_{h}Z\left(u_h\right)-{\mathcal{R}}_{h}Z\left(u_h\right))+{\Delta }_h{\mathcal{R}}_{h}Z\left(u_h\right)+{\mathcal{P}}_h(u_h-u_d)-{\mathcal{P}}_h{(}_t{D}_T^{\alpha }Z\left(u_h\right))\hfill \\ & =& {\Delta }_h({\mathcal{P}}_{h}Z\left(u_h\right)-{\mathcal{R}}_{h}Z\left(u_h\right))+{\mathcal{P}}_h\Delta Z\left(u_h\right)+{\mathcal{P}}_h(u_h-u_d)-{\mathcal{P}}_h{(}_t{D}_T^{\alpha }Z\left(u_h\right))\hfill \\ & =& {\Delta }_h({\mathcal{P}}_{h}Z\left(u_h\right)-{\mathcal{R}}_{h}Z\left(u_h\right)).\hfill \end{array}$$

Therefore, we obtain that

$$\begin{array}{c}\hfill \eta ={\int }_t^T{E}_h(\tau -t){\Delta }_h({\mathcal{P}}_{h}Z\left(u_h\right)-{\mathcal{R}}_{h}Z\left(u_h\right))d\tau ,\end{array}$$

(27)

where

$$\begin{array}{c}\hfill E_h\left(t\right)v={\displaystyle \sum _{i=1}^M}{t}^{\alpha -1}{\Xi }_{\alpha ,\alpha }(-{\lambda }_i^h{t}^{\alpha })(v,{\varphi }_i^h){\varphi }_i^h,\end{array}$$

here ${\left\{({\lambda }_i^h,{\varphi }_i^h)\right\}}_{i=1}^M$ are the eigenvalues and eigenfunctions of ${\Delta }_h$.

Using Lemma 3, for $l=0,1$, we have

$$\begin{array}{ccc}\hfill & & {\int }_0^T{\parallel \eta \parallel }_{H^1\left(\mathsf{\Omega }\right)}^{2}dt\le {\int }_0^T\parallel {\Delta }_h{({\mathcal{P}}_{h}Z\left(u_h\right)-{\mathcal{R}}_{h}Z\left(u_h\right)\parallel }_{H^{l-2}}^{2}ds\hfill \\ & \le & {\int }_0^T\parallel {\Delta }_h({\mathcal{P}}_{h}Z\left(u_h\right)-{\mathcal{R}}_{h}Z\left(u_h\right){\parallel }_{H^l}^{2}ds\le C{h}^{4-2l}{\parallel Z\left(u_h\right)\parallel }_{L^2(0,T;H^{l+2}\left(\mathsf{\Omega }\right))}^2\hfill \\ & \le & C{h}^{4-2l}{\parallel u_h-u_d\parallel }_{L^2(0,T;H^{l+1}\left(\mathsf{\Omega }\right))}^2.\hfill \end{array}$$

(28)

Therefore, we can obtain the theorem result by using triangle inequality.

Moreover, (23) and (26), we deduce

$$\begin{array}{c}\hfill {(}_t{D}_T^{\alpha }(Z\left(u_h\right)-Z\left(q_h\right)),v)+(\nabla (Z\left(u_h\right)-Z\left(q_h\right)),\nabla v)=(u_h-u\left(q_h\right),v),\forall v\in H_0^1\left(\mathsf{\Omega }\right).\end{array}$$

According to the existence and uniqueness of solution, the adjoint state equation implies

$$\begin{array}{c}\hfill \parallel Z\left(u_h\right)-Z\left(q_h\right){\parallel }_{L^2\left({\mathsf{\Omega }}_T\right)}+\parallel \nabla (Z\left(u_h\right)-Z\left(q_h\right)){\parallel }_{L^2\left({\mathsf{\Omega }}_T\right)}\le C{\parallel u_h-u\left(q_h\right)\parallel }_{L^2\left({\mathsf{\Omega }}_T\right)}.\end{array}$$

(29)

Therefore, by (25), (29) and Lemma 5, we can obtain

$$\begin{array}{ccc}\hfill \parallel Z_h-Z\left(q_h\right){\parallel }_{L^2\left({\mathsf{\Omega }}_T\right)}& \le & \parallel Z\left(q_h\right)-Z\left(u_h\right){\parallel }_{L^2\left({\mathsf{\Omega }}_T\right)}+{\parallel Z\left(u_h\right)-Z_h\parallel }_{L^2\left({\mathsf{\Omega }}_T\right)}\hfill \\ & \le & C\parallel u\left(q_h\right)-u_h{\parallel }_{L^2\left({\mathsf{\Omega }}_T\right)}+C{h}^2{\parallel u_h-u_d\parallel }_{L^2\left({\mathsf{\Omega }}_T\right)}\hfill \\ & \le & C{h}^2\parallel f+q_h{\parallel }_{L^2\left({\mathsf{\Omega }}_T\right)}+C{h}^2{\parallel u_h-u_d\parallel }_{L^2\left({\mathsf{\Omega }}_T\right)}\le C{h}^2,\hfill \end{array}$$

(30)

and

$$\begin{array}{ccc}\hfill \parallel \nabla (Z_h-Z\left(q_h\right)){\parallel }_{L^2\left({\mathsf{\Omega }}_T\right)}& \le & \parallel \nabla (Z\left(q_h\right)-Z\left(u_h\right)){\parallel }_{L^2\left({\mathsf{\Omega }}_T\right)}+{\parallel \nabla (Z\left(u_h\right)-Z_h)\parallel }_{L^2\left({\mathsf{\Omega }}_T\right)}\hfill \\ & \le & C\parallel u\left(q_h\right)-u_h{\parallel }_{L^2\left({\mathsf{\Omega }}_T\right)}+Ch{\parallel u_h-u_d\parallel }_{L^2\left({\mathsf{\Omega }}_T\right)}\hfill \\ & \le & C{h}^2\parallel f+q_h{\parallel }_{L^2\left({\mathsf{\Omega }}_T\right)}+Ch{\parallel u_h-u_d\parallel }_{L^2\left({\mathsf{\Omega }}_T\right)}\le Ch.\hfill \end{array}$$

(31)

Based on the above Lemmas 2–5, we will give the error estimation for $q-q_h,u-u_h,Z-Z_h$ in the following Theorem 1.

Theorem 1.

Let $(u,Z,q)$ and $(u_h,Z_h,q_h)$ be the solution of (2), (5), (7) and (18)–(20), respectively. Then the estimate follows as

$$\begin{array}{c}\hfill \parallel q-q_h{\parallel }_{L^2\left({\mathsf{\Omega }}_T\right)}+\parallel u-u_h{\parallel }_{L^2\left({\mathsf{\Omega }}_T\right)}+{\parallel Z-Z_h\parallel }_{L^2\left({\mathsf{\Omega }}_T\right)}\\ \hfill +h\parallel \nabla (u-u_h){\parallel }_{L^2\left({\mathsf{\Omega }}_T\right)}+h{\parallel \nabla (Z-Z_h)\parallel }_{L^2\left({\mathsf{\Omega }}_T\right)}\le C{h}^2\end{array}$$

holds.

Proof. 

Based on Lemma 4 and inequality (30), one can obtain

$$\begin{array}{c}\hfill \parallel q-q_h{\parallel }_{L^2\left({\mathsf{\Omega }}_T\right)}\le C{h}^2.\end{array}$$

(32)

By (24), (25) and (32), we can obtain

$$\begin{array}{c}\hfill \parallel u\left(q\right)-u_h{\parallel }_{L^2\left({\mathsf{\Omega }}_T\right)}+h{\parallel \nabla (u\left(q\right)-u_h)\parallel }_{L^2\left({\mathsf{\Omega }}_T\right)}\le C{h}^2.\end{array}$$

(33)

It follows from the adjoint (7) and auxiliary problem (23) that

$$\begin{array}{c}\hfill {(}_t{D}_T^{\alpha }(Z-Z\left(q_h\right)),v)+(\nabla (Z-Z\left(q_h\right)),\nabla v)=(u-u\left(q_h\right),v),\forall v\in H_0^1\left(\mathsf{\Omega }\right).\end{array}$$

Combining with the above equation, the adjoint state’s stability estimate and (25), we can obtain

$$\begin{array}{c}\hfill \parallel Z\left(q_h\right){-Z\parallel }_{L^2\left({\mathsf{\Omega }}_T\right)}+\parallel \nabla (Z\left(q_h\right)-Z){\parallel }_{L^2\left({\mathsf{\Omega }}_T\right)}\le C{\parallel q_h-q\parallel }_{L^2\left({\mathsf{\Omega }}_T\right)}.\end{array}$$

(34)

Further, using (32)–(34) gives

$$\begin{array}{c}\hfill \parallel Z_h{-Z\parallel }_{L^2\left({\mathsf{\Omega }}_T\right)}+h{\parallel \nabla (Z_h-Z)\parallel }_{L^2\left({\mathsf{\Omega }}_T\right)}\le C{h}^2.\end{array}$$

(35)

Based on (32), (33) and (35), we complete the prove of the theorem. □

4. Fully Discrete Scheme for the FOCP

4.1. The FD-FE Scheme for the State Equation

In order to construct the FD-FE scheme, we divide the time domain I into subdomins with $\tau =T/N,t_k=k\tau ,0\le k\le N$. Following [42,43], we approximate the fractional derivative as follows

$$\begin{array}{c}\hfill {}_0{D}_t^{\alpha }u(x,t_1)=\frac{1}{\mathsf{\Gamma }(1-\alpha )}{\int }_0^{t_1}\frac{{\partial }_{s}u(x,s)}{{(t_1-s)}^{\alpha }}ds=\frac{1}{\mathsf{\Gamma }(1-\alpha )}{\int }_0^{t_1}\frac{{\partial }_s{J}_{[t_0,t_2]}u(x,s)}{{(t_1-s)}^{\alpha }}ds+R_{\tau }^1,\end{array}$$

(36)

where $J_{[t_0,t_2]}u(x,t)$ is the quadratic interpolation on $[t_0,t_2]$ with respect to time variable, i.e.,

$$\begin{array}{c}\hfill J_{[t_0,t_2]}u(x,t)={\varpi }_{0,0}\left(t\right)u^0\left(x\right)+{\varpi }_{1,0}\left(t\right)u^1\left(x\right)+{\varpi }_{2,0}\left(t\right)u^2\left(x\right),\end{array}$$

(37)

with $u^i\left(x\right)$ is the numerical solution at $t_i$ of $u(x,t)$, and ${\varpi }_{i,0}\left(t\right),i=0,1,2,$ are the quadratic interpolation function at points $t_0,t_1$ and $t_2$ defined by

$$\begin{array}{c}\hfill {\varpi }_{0,0}\left(t\right)=\frac{(t-t_1)(t-t_2)}{2{\tau }^2},{\varpi }_{1,0}\left(t\right)=\frac{(t-t_2)(t-t_0)}{-{\tau }^2},{\varpi }_{2,0}\left(t\right)=\frac{(t-t_1)(t-t_0)}{2{\tau }^2}.\end{array}$$

(38)

Bringing (37) and (38) into (36), we have

$$\begin{array}{c}\hfill {}_0{D}_t^{\alpha }u(x,t_1)=B_1^{0,0}{u}^0\left(x\right)+B_1^{1,0}{u}^1\left(x\right)+B_1^{2,0}{u}^2\left(x\right)+R_{\tau }^1,\end{array}$$

(39)

with

$$\begin{array}{ccc}\hfill B_1^{i,0}& =& \frac{1}{\mathsf{\Gamma }(1-\alpha )}{\int }_0^{t_1}{(t_1-s)}^{-\alpha }{\varpi }_{i,0}^{{}^{\prime }}\left(s\right)ds,i=0,1,2,\hfill \\ \hfill R_{\tau }^1& =& \frac{1}{\mathsf{\Gamma }(1-\alpha )}{\int }_0^{t_1}{(t_1-s)}^{-\alpha }{r}^1(x,s)ds,\hfill \end{array}$$

here $r^1(x,t)\doteq u(x,t)-J_{[t_0,t_2]}u(x,t)$.

Similar to (39), we can obtain the approximation solution for $u_2\left(x\right)$ as follows

$$\begin{array}{c}\hfill {}_0{D}_t^{\alpha }u(x,t_2)=B_2^{0,0}{u}^0\left(x\right)+B_2^{1,0}{u}^1\left(x\right)+B_2^{2,0}{u}^2\left(x\right)+R_{\tau }^2,\end{array}$$

(40)

where

$$\begin{array}{ccc}\hfill B_2^{i,0}& =& \frac{1}{\mathsf{\Gamma }(1-\alpha )}{\int }_0^{t_2}{(t_2-s)}^{-\alpha }{\varpi }_{i,0}^{{}^{\prime }}\left(s\right)ds,i=0,1,2,\hfill \\ \hfill R_{\tau }^2& =& \frac{1}{\mathsf{\Gamma }(1-\alpha )}{\int }_0^{t_2}{(t_2-s)}^{-\alpha }{r}^2(x,s)ds,\hfill \end{array}$$

with $r^2(x,t)\doteq u(x,t)-J_{[t_0,t_2]}u(x,t)$.

For $k\ge 3$, we have

$$\begin{array}{ccc}\hfill {}_0{D}_t^{\alpha }u(x,t_k)& =& \frac{1}{\mathsf{\Gamma }(1-\alpha )}{\int }_0^{t_k}\frac{{\partial }_{s}u(x,s)}{{(t_k-s)}^{\alpha }}ds\hfill \\ & =& \frac{1}{\mathsf{\Gamma }(1-\alpha )}\left[{\int }_0^{t_1}\frac{{\partial }_{s}u(x,s)}{{(t_k-s)}^{\alpha }}ds+{\displaystyle \sum _{j=1}^{k-1}}{\int }_{t_j}^{t_{j+1}}\frac{{\partial }_{s}u(x,s)}{{(t_k-s)}^{\alpha }}ds\right].\hfill \end{array}$$

(41)

On $[t_j,t_{j+1}]$, the approximated solution to $u(x,t)$ can be defined by

$$\begin{array}{ccc}\hfill u(x,t)& =& {\omega }_{0,j}\left(t\right)u^{j-1}\left(x\right)+{\omega }_{1,j}\left(t\right)u^j\left(x\right)+{\omega }_{2,j}\left(t\right)u^{j+1}\left(x\right)+r^j(x,t)\hfill \\ & \doteq & J_{[t_j,t_{j+1}]}u(x,t)+r^j(x,t),\hfill \end{array}$$

(42)

where ${\omega }_{i,j}\left(t\right),i=0,1,2;j=1,\cdots ,k-1$, are the quadratic interpolation basis function at points $t_{j-1},t_j,t_{j+1}$ defined by

$$\begin{array}{c}\hfill {\omega }_{0,j}\left(t\right)=\frac{(t-t_{j+1})(t-t_j)}{2{\tau }^2},{\omega }_{1,j}\left(t\right)=\frac{(t-t_{j+1})(t-t_{j-1})}{-{\tau }^2},{\omega }_{2,j}\left(t\right)=\frac{(t-t_j)(t-t_{j-1})}{2{\tau }^2}.\end{array}$$

Through careful calculation, one can obtain

$$\begin{array}{ccc}\hfill {}_0{D}_t^{\alpha }u(x,t_k)& =& \frac{1}{\mathsf{\Gamma }(1-\alpha )}\{{\int }_0^{t_1}\frac{{\partial }_s\left[J_{[t_0,t_2]}u(x,s)\right]}{{(t_k-s)}^{\alpha }}ds+{\displaystyle \sum _{j=1}^{k-1}}{\int }_{t_j}^{t_{j+1}}\frac{{\partial }_s\left[J_{[t_j,t_{j+1}]}u(x,s)\right]}{{(t_k-s)}^{\alpha }}ds\hfill \\ \hfill & & +{\int }_0^{t_1}{(t_k-s)}^{-\alpha }{r}^2(x,s)ds+{\displaystyle \sum _{j=1}^{k-1}}{\int }_{t_j}^{t_{j+1}}{(t_k-s)}^{-\alpha }{r}^{j+1}(x,s)ds\}\hfill \\ & =& B_k^{0,0}{u}^0\left(x\right)+B_k^{1,0}{u}^1\left(x\right)+B_k^{2,0}{u}^2\left(x\right)\hfill \\ & & +{\displaystyle \sum _{j=1}^{k-1}}[B_k^{0,j}{u}^{j-1}\left(x\right)+B_k^{1,j}{u}^j\left(x\right)+B_k^{2,j}{u}^{j+1}\left(x\right)]+R_{\tau }^k,\hfill \end{array}$$

(43)

where

$$\begin{array}{ccc}\hfill B_k^{i,0}& =& \frac{1}{\mathsf{\Gamma }(1-\alpha )}{\int }_0^{t_1}{(t_k-s)}^{-\alpha }{\varpi }_{i,0}^{{}^{\prime }}\left(s\right)ds,i=0,1,2,\hfill \\ \hfill B_k^{i,j}& =& \frac{1}{\mathsf{\Gamma }(1-\alpha )}{\int }_{t_j}^{t_{j+1}}{(t_k-s)}^{-\alpha }{\varpi }_{i,j}^{{}^{\prime }}\left(s\right)ds,i=0,1,2;j=1,\cdots ,k-1,\hfill \end{array}$$

and

$$\begin{array}{c}\hfill R_{\tau }^k={\int }_0^{t_1}{(t_k-s)}^{-\alpha }{r}^2(x,s)ds+{\displaystyle \sum _{j=1}^{k-1}}{\int }_{t_j}^{t_{j+1}}{r}^{j+1}(x,s)ds,\end{array}$$

with $r^{j+1}(x,t)$ is defined by $r^{j+1}(x,t)\doteq u(x,t)-J_{[t_j,t_{j+1}]}u(x,t)$ for $j\ge 2$ and $r^2(x,t)\doteq u(x,t)-J_{[t_0,t_2]}u(x,t)$.

In all cases, the left Caputo fractional derivative ${}_0{D}_t^{\alpha }u(x,t_k)$ can be determined by a linear combination of $u^j\left(x\right)$. Furthermore, by simplifying the calculation, it is easy to find that all $B_k^{i,j}$ are proportional to ${\tau }^{-\alpha }$. Therefore, we summarize (39), (40) and (43) to write down them uniformly as

$$\begin{array}{c}\hfill {}_0{D}_t^{\alpha }u(x,t_k)\approx {}_0{D}_{\tau }^{\alpha }u(x,t_k),\end{array}$$

(44)

where the newly introduced operator ${}_0{D}_{\tau }^{\alpha }$ is the discrete Caputo derivative defined by

$$\begin{array}{c}\hfill {}_0{D}_{\tau }^{\alpha }u(x,t_k)=\left\{\begin{array}{c}{\displaystyle \frac{1}{\mathsf{\Gamma }(3-\alpha ){\tau }^{\alpha }}[{\widehat{D}}_0{u}^0\left(x\right)+{\widehat{D}}_1{u}^1\left(x\right)+{\widehat{D}}_2{u}^2\left(x\right)],k=1,}\hfill \\ {\displaystyle \frac{1}{\mathsf{\Gamma }(3-\alpha ){\tau }^{\alpha }}[{\tilde{D}}_0{u}^0\left(x\right)+{\tilde{D}}_1{u}^1\left(x\right)+{\tilde{D}}_2{u}^2\left(x\right)],k=2,}\hfill \\ {\displaystyle \frac{1}{\mathsf{\Gamma }(3-\alpha ){\tau }^{\alpha }}\{\overline{A_k}{u}^0\left(x\right)+{\overline{B}}_k{u}^1\left(x\right)+{\overline{C}}_k{u}^2\left(x\right)}\hfill \\ {\displaystyle +{\displaystyle \sum _{j=1}^{k-1}}(A_j{u}^{k-j-1}\left(x\right)+B_j{u}^{k-j}\left(x\right)+C_j{u}^{k-j+1}\left(x\right))\},k\ge 3,}\hfill \end{array}\right.\end{array}$$

(45)

where all the coefficients in (45) are constants and can be computed analytically as follows

$$\begin{array}{ccc}& & {\widehat{D}}_0=\frac{1}{2}(3\alpha -4),{\widehat{D}}_1=2(1-\alpha ),{\widehat{D}}_2=\frac{\alpha }{2},\hfill \\ & & {\tilde{D}}_0=\frac{1}{2^{\alpha }}(3\alpha -2),{\tilde{D}}_1=-\frac{4\alpha }{2^{\alpha }},{\tilde{D}}_2=\frac{1}{2^{\alpha }}(\alpha +2),\hfill \\ & & {\displaystyle {\overline{A}}_k=-\frac{3}{2}(2-\alpha )k^{1-\alpha }+\frac{1}{2}(2-\alpha ){(k-1)}^{1-\alpha }+k^{2-\alpha }-{(k-1)}^{2-\alpha },}\hfill \\ & & {\displaystyle {\overline{B}}_k=2{(k-1)}^{2-\alpha }+2(2-\alpha )k^{1-\alpha }-2k^{2-\alpha },}\hfill \\ & & {\displaystyle {\overline{C}}_k=-\frac{2-\alpha }{2}\left[k^{1-\alpha }+{(k-1)}^{1-\alpha }\right]+k^{2-\alpha }-{(k-1)}^{2-\alpha },}\hfill \\ & & {\displaystyle A_j=-\frac{1}{2}(2-\alpha )\left[j^{1-\alpha }+{(j-1)}^{1-\alpha }\right]+j^{2-\alpha }-{(j-1)}^{2-\alpha },}\hfill \\ & & {\displaystyle B_j=2(2-\alpha ){(j-1)}^{1-\alpha }-2j^{2-\alpha }+2{(j-1)}^{2-\alpha },}\hfill \\ & & {\displaystyle C_j=\frac{1}{2}(2-\alpha )j^{1-\alpha }-\frac{3}{2}(2-\alpha ){(j-1)}^{1-\alpha }+j^{2-\alpha }-{(j-1)}^{2-\alpha }.}\hfill \end{array}$$

If the solution is sufficiently smooth on time variable, based on the idea of [42,43,44], we have

$$\begin{array}{c}\hfill \parallel R_{\tau }^{k+1}{\parallel }_0\le C{\tau }^{3-\alpha }{\parallel u\parallel }_{W^{3,\infty }(0,T;L^2\left(\mathsf{\Omega }\right))},k=1,2,\cdots ,N,\end{array}$$

(46)

where $W^{3,\infty }(I;\mathsf{\Omega })={\left\{\varphi (x,t)\right|\parallel \varphi (x,t)\parallel }_{L^2\left(\mathsf{\Omega }\right)}\in W^{3,\infty }\left(I\right)\}$ and $W^{3,\infty }\left(I\right)$ is the norm Sobolev space.

Then, the full discrete FD-FE scheme for (15) reads: find $u_h^k\in {\mathcal{V}}_h,k=1,2,\cdots ,N$ satisfying

$$\begin{array}{c}\hfill \left\{\begin{array}{c}{(}_0{D}_{\tau }^{\alpha }{u}_h^k,v_h)+(\nabla u_h^k,\nabla v_h)=(f^k,v_h),\forall v_h\in {\mathcal{V}}_h,\hfill \\ u_h^k\left(0\right)=0.\hfill \end{array}\right.\end{array}$$

(47)

The purpose is that analyzing the stability of the scheme (47), we firstly introduce some notations to reformulate (45) for $k\ge 3$. We denote

$$\begin{array}{c}\hfill {\alpha }_0=\mathsf{\Gamma }(3-\alpha ){\tau }^{\alpha },{\beta }_0=C_1=\frac{4-\alpha }{2}.\end{array}$$

(48)

When $k=3$, we take

$$\begin{array}{c}\hfill \begin{array}{c}{\displaystyle D_0^3=-({\overline{A}}_3+A_2){\beta }_0^{-1},D_1^3=-({\overline{B}}_3+A_1+B_2){\beta }_0^{-1},D_2^3=-({\overline{C}}_3+B_1+C_2){\beta }_0^{-1},}\hfill \end{array}\end{array}$$

(49)

For $k\ge 4$,

$$\begin{array}{c}\hfill \begin{array}{c}{D}_{k-1}^k=-(B_1+C_2){\beta }_0^{-1},d_{k-2}^k=-(A_1+B_2+C_3){\beta }_0^{-1},\hfill \\ D_{k-i}^k=-(A_{i-1}+B_i+C_{i+1}){\beta }_0^{-1},i=3,4,\cdots ,k-3,\hfill \\ D_2^k=-({\overline{C}}_k+A_{k-3}+B_{k-2}+C_{k-1}){\beta }_0^{-1},D_1^k=-({\overline{B}}_k+A_{k-2}+B_{k-1}){\beta }_0^{-1},\hfill \\ D_0^k=-({\overline{A}}_k+A_{k-1}){\beta }_0^{-1}.\hfill \end{array}\end{array}$$

(50)

Therefore, (47) can be rewritten into the equivalent form as follows

$$\begin{array}{c}\hfill \left\{\begin{array}{c}{\displaystyle (u_h^k,v_h)+{\alpha }_0{\beta }_0^{-1}(\nabla u_h^k,\nabla v_h)=({\displaystyle \sum _{i=1}^k}{D}_{k-i}^k{u}_h^{k-i},v_h)+{\alpha }_0{\beta }_0^{-1}(f^k,v_h),\forall v_h\in {\mathcal{V}}_h,}\hfill \\ u_h^k\left(0\right)=0.\hfill \end{array}\right.\end{array}$$

(51)

In order to analyze the stability analysis of the scheme (51) we firstly show the properties of the coefficients $D_{k-i}^k$. See the following Lemma 6.

Lemma 6.

For given $0<\alpha <1$, $k\ge 4$, the scheme coefficients of (47) satisfy

(1) ${\beta }_0=C_1=\frac{4-\alpha }{2}\in (\frac{3}{2},2)$. (2) ${\displaystyle {\displaystyle \sum _{i=0}^k}{D}_{k-i}^k=1}$.

(3) $D_{k-i}^k>0$, $3\le i\le k$. (4) $D_{k-1}^k>0$.

(5) There exists ${\alpha }_0\in (0,1)$ such that $D_{k-2}^k>0$ if $\alpha \in (0,{\alpha }_0)$, and $D_{k-2}^k<0$ if $\alpha \in ({\alpha }_0,1)$.

(6) ${\displaystyle \frac{1}{4}}{\left(D_{k-1}^k\right)}^2+D_{k-2}^k>0$.

Proof. 

For a detailed proof, see Appendix A. □

Next, as $k\ge 4$, denote

$$\begin{array}{c}\hfill \rho =\frac{1}{2}{D}_{k-1}^k.\end{array}$$

(52)

Recombining of the terms in (51), we obtain

$$\begin{array}{c}\hfill \begin{array}{c}{\displaystyle (u_h^k-\rho u_h^{k-1},v_h)+{\alpha }_0{\beta }_0^{-1}(\nabla u_h^k,\nabla v_h)}\hfill \\ {\displaystyle =\rho (u_h^{k-1}-\rho u_h^{k-2},v_h)+({\rho }^2+D_{k-2}^k)(u_h^{k-2}-\rho u_h^{k-3},v_h)}\hfill \\ {\displaystyle +({\rho }^3+\rho D_{k-2}^k+D_{k-3}^k)(u_h^{k-3},v_h)+D_{k-4}^k(u_h^{k-4},v_h)+\cdots +D_0^k(u_h^0,v_h)+{\alpha }_0{\beta }_0^{-1}(f^k,v_h)}\hfill \\ {\displaystyle =\rho (u_h^{k-1}-\rho u_h^{k-2},v_h)+({\rho }^2+D_{k-2}^k)(u_h^{k-2}-\rho u_h^{k-3},v_h)}\hfill \\ {\displaystyle +({\rho }^3+\rho D_{k-2}^k+D_{k-3}^k)(u_h^{k-3},v_h)}\hfill \\ {\displaystyle +\cdots +({\rho }^{k-2}+{\rho }^{k-4}{D}_{k-2}^k+\cdots +\rho D_3^k+D_2^k)(u_h^2-\rho u_h^1,v_h)}\hfill \\ {\displaystyle +({\rho }^{k-1}+{\rho }^{k-3}{D}_{k-2}^k+\cdots +\rho D_2^k+D_1^k)(u_h^1-\rho u_h^0,v_h)}\hfill \\ {\displaystyle +({\rho }^k+{\rho }^{k-2}{D}_{k-2}^k+\cdots +\rho D_1^k+D_0^k)(u_h^0,v_h)+{\alpha }_0{\beta }_0^{-1}(f^k,v_h).}\hfill \end{array}\end{array}$$

Now we denote

$$\begin{array}{c}\hfill {\overline{D}}_{k-i}^k={\rho }^i+{\displaystyle \sum _{j=2}^i}{\rho }^{i-j}{D}_{k-j}^k,i=2,3,4,\cdots ,k,\end{array}$$

(53)

$$\begin{array}{c}\hfill {\overline{u}}_h^0=u_h^0,{\overline{u}}_h^i=u_h^i-\rho u_h^{i-1},i=1,2,\cdots ,k.\end{array}$$

(54)

As $k=3$, we have

$$\begin{array}{c}\hfill \begin{array}{c}{\displaystyle (u_h^3-\rho u_h^2,v_h)+{\alpha }_0{\beta }_0^{-1}(\nabla u_h^3,\nabla v_h)}\hfill \\ {\displaystyle ={\overline{D}}_2^3({\overline{u}}_h^2,v_h)+{\overline{D}}_1^3({\overline{u}}_h^1,v_h)+{\overline{D}}_0^3({\overline{u}}_h^0,v_h)+{\alpha }_0{\beta }_0^{-1}(f^3,v_h),}\hfill \end{array}\end{array}$$

where

$$\begin{array}{c}\hfill {\overline{D}}_2^3=D_2^3-\rho ,{\overline{D}}_1^3={\overline{D}}_2^3\rho +D_1^3,{\overline{D}}_0^3={\overline{D}}_1^3\rho +D_0^3.\end{array}$$

(55)

Then the equivalent form to (51) can be determined by

$$\begin{array}{c}\hfill \left\{\begin{array}{c}{\displaystyle ({\overline{u}}_h^3,v_h)+{\alpha }_0{\beta }_0^{-1}(\nabla u_h^3,\nabla v_h)=({\displaystyle \sum _{i=1}^3}{\overline{D}}_{3-i}^3{\overline{u}}_h^{3-i},v_h)+{\alpha }_0{\beta }_0^{-1}(f^3,v_h),}\hfill \\ {\displaystyle ({\overline{u}}_h^k,v_h)+{\alpha }_0{\beta }_0^{-1}(\nabla u_h^k,\nabla v_h)=(\rho {\overline{u}}_h^{k-1}+{\displaystyle \sum _{i=2}^{k-1}}{\overline{D}}_{k-i}^k{\overline{u}}_h^{k-i}}\hfill \\ +{\overline{D}}_0^k{\overline{u}}_h^0,v_h)+{\alpha }_0{\beta }_0^{-1}(f^k,v_h),k\ge 4.\hfill \end{array}\right.\end{array}$$

(56)

The new coefficients of (56) have some good properties for $k\ge 4$ which are given by the following Lemma 7.

Lemma 7.

For given $0<\alpha <1$, $k\ge 4$, the coefficients of (56) satisfy

(1) $0<\rho <\frac{2}{3}$.          (2) ${\overline{D}}_{k-i}^k>0,2\le i\le k$.

(3) ${\displaystyle \rho +{\displaystyle \sum _{i=2}^{k-1}}{\overline{D}}_{k-i}^k+{\overline{D}}_0^k\le 1}$.    (4) ${\displaystyle \frac{1}{{\overline{D}}_0^k}<\frac{1}{D_0^k}<\frac{k^{\alpha }}{(1-\alpha )(2-\alpha )}}$.

Proof. 

For a detailed proof, see Appendix B. □

For $k=3$, the new coefficients of (56) have some good properties by the following Lemma 8.

Lemma 8.

For given $0<\alpha <1$, for $k=3$, the coefficients of the first equation in the scheme (56) satisfy

(1) ${\overline{D}}_{3-i}^3>0,i=1,2,3$.   (2) ${\overline{D}}_2^3+{\overline{D}}_1^3+{\overline{D}}_0^3\le 1$.   (3) ${\displaystyle {\overline{D}}_2^3-\rho \le 0}$.

Proof. 

For a detailed proof, see Appendix C. □

Next, we will analysis the stability of (47) in the following Theorem 2.

Theorem 2.

Let ${\left\{u_h^k\right\}}_{k=1}^N$ be the numerical solution of (47). Then the estimate is following as

$$\begin{array}{c}\hfill \parallel u_h^k{\parallel }_0\le 3\sqrt{20T^{\alpha }\mathsf{\Gamma }(1-\alpha )d}\underset{i}{max}{\parallel f^i\parallel }_0,1\le k\le N.\end{array}$$

(57)

Proof. 

Firstly, we analysis the case for $k=1,2$. Based on the fact that $u_h^0=0$, therefore, we have

$$\begin{array}{c}\hfill {\widehat{D}}_1(u_h^1,v_h)+{\widehat{D}}_2(u_h^2,v_h)+{\alpha }_0(\nabla u_h^1,\nabla v_h)={\alpha }_0(f^1,v_h),\end{array}$$

(58)

$$\begin{array}{c}\hfill {\tilde{D}}_1(u_h^1,v_h)+{\tilde{D}}_2(u_h^2,v_h)+{\alpha }_0(\nabla u_h^2,\nabla v_h)={\alpha }_0(f^2,v_h).\end{array}$$

(59)

We choose $v_h=-{\tilde{D}}_1{u}_h^1$ in (58), $v_h={\widehat{D}}_2{u}_h^2$ in (59) and add them together, we get

$$\begin{array}{ccc}\hfill & & -{\widehat{D}}_1{\tilde{D}}_1\parallel u_h^1{\parallel }_0^2+{\widehat{D}}_2{\tilde{D}}_2\parallel u_h^2{\parallel }_0^2-{\alpha }_0{\tilde{D}}_1\parallel \nabla u_h^1{\parallel }_0^2+{\alpha }_0{\widehat{D}}_2{\parallel \nabla u_h^2\parallel }_0^2\hfill \\ \hfill & =& -{\alpha }_0{\tilde{D}}_1(f^1,u_h^1)+{\alpha }_0{\widehat{D}}_2(f^2,u_h^2)\le \frac{1}{2}\frac{{(-{\alpha }_0{\tilde{D}}_1)}^2}{-{\widehat{D}}_1{\tilde{D}}_1}\parallel f^1{\parallel }_0^2+\frac{1}{2}(-{\widehat{D}}_1{\tilde{D}}_1){\parallel u_h^1\parallel }_0^2\hfill \\ & & +\frac{1}{2}\frac{{\left({\alpha }_0{\widehat{D}}_2\right)}^2}{{\widehat{D}}_2{\tilde{D}}_2}\parallel f^2{\parallel }_0^2+\frac{1}{2}{\widehat{D}}_2{\tilde{D}}_2{\parallel u_h^2\parallel }_0^2.\hfill \end{array}$$

(60)

Simply the (60), we obtain that

$$\begin{array}{ccc}\hfill & & -{\widehat{D}}_1{\tilde{D}}_1\parallel u_h^1{\parallel }_0^2+{\widehat{D}}_2{\tilde{D}}_2\parallel u_h^2{\parallel }_0^2-2{\alpha }_0{\tilde{D}}_1\parallel \nabla u_h^1{\parallel }_0^2+2{\alpha }_0{\widehat{D}}_2{\parallel \nabla u_h^2\parallel }_0^2\hfill \\ & \le & \frac{{(-{\alpha }_0{\tilde{D}}_1)}^2}{-{\widehat{D}}_1{\tilde{D}}_1}\parallel f^1{\parallel }_0^2+\frac{{\left({\alpha }_0{\widehat{D}}_2\right)}^2}{{\widehat{D}}_2{\tilde{D}}_2}{\parallel f^2\parallel }_0^2.\hfill \end{array}$$

(61)

According to (61), we have

$$\begin{array}{ccc}& & \parallel u_h^1{\parallel }_0^2+{\alpha }_0{\beta }_0^{-1}\parallel \nabla u_h^1{\parallel }_0^2\le \parallel u_h^1{\parallel }_0^2+\frac{2{\alpha }_0}{{\widehat{D}}_1}{\parallel \nabla u_h^1\parallel }_0^2\hfill \\ & \le & \left[\frac{{(-{\alpha }_0{\tilde{D}}_1)}^2}{-{\widehat{D}}_1{\tilde{D}}_1}+\frac{{\left({\alpha }_0{\widehat{D}}_2\right)}^2}{{\widehat{D}}_2{\tilde{D}}_2}\right]/(-{\widehat{D}}_1{\tilde{D}}_1)\underset{i}{max}{\parallel f^i\parallel }_0^2\hfill \\ & =& {(2-\alpha )}^2{(1-\alpha )}^2{\left[\mathsf{\Gamma }(1-\alpha )\right]}^2{\tau }^{2\alpha }\left[\frac{1}{4{(1-\alpha )}^2}+\frac{4^{\alpha }}{16(1-\alpha )(2+\alpha )}\right]\underset{i}{max}{\parallel f^i\parallel }_0^2\hfill \\ & \le & 4{\left[\mathsf{\Gamma }(1-\alpha )\right]}^2{T}^{\alpha }\underset{i}{max}\parallel f^i{\parallel }_0^2\le 12d{\left[\mathsf{\Gamma }(1-\alpha )\right]}^2{T}^{\alpha }\underset{i}{max}{\parallel f^i\parallel }_0^2,\hfill \end{array}$$

where d is the dimension of space.

Therefore, ${\overline{u}}_h^1=u_h^1-\rho u_h^0$ and $0<\rho <\frac{2}{3}$, as $k=1$, we have

$$\begin{array}{c}\hfill \parallel {\overline{u}}_h^1{\parallel }_0^2+{\alpha }_0{\beta }_0^{-1}\parallel \nabla u_h^1{\parallel }_0^2\le 20d{\left[\mathsf{\Gamma }(1-\alpha )\right]}^2{T}^{\alpha }\underset{i}{max}{\parallel f^i\parallel }_0^2.\end{array}$$

(62)

Similar, according to (61), we have

$$\begin{array}{ccc}& & \parallel u_h^2{\parallel }_0^2+{\alpha }_0{\beta }_0^{-1}\parallel \nabla u_h^2{\parallel }_0^2\le \parallel u_h^2{\parallel }_0^2+\frac{2{\alpha }_0}{{\tilde{D}}_2}{\parallel \nabla u_h^2\parallel }_0^2\hfill \\ & \le & \frac{1}{{\widehat{D}}_2{\tilde{D}}_2}\left[\frac{{(-{\alpha }_0{\tilde{D}}_1)}^2}{-{\widehat{D}}_1{\tilde{D}}_1}+\frac{{\left({\alpha }_0{\widehat{D}}_2\right)}^2}{{\widehat{D}}_2{\tilde{D}}_2}\right]\underset{i}{max}{\parallel f^i\parallel }_0^2\hfill \\ & =& {(2-\alpha )}^2{(1-\alpha )}^2{\left[\mathsf{\Gamma }(1-\alpha )\right]}^2{\tau }^{2\alpha }\left[\frac{4}{(1-\alpha )(2+\alpha )}+\frac{4^{\alpha }}{{(2+\alpha )}^2}\right]\underset{i}{max}{\parallel f^i\parallel }_0^2\hfill \\ & \le & 12{\left[\mathsf{\Gamma }(1-\alpha )\right]}^2{T}^{\alpha }\underset{i}{max}\parallel f^i{\parallel }_0^2\le 12d{\left[\mathsf{\Gamma }(1-\alpha )\right]}^2{T}^{\alpha }\underset{i}{max}{\parallel f^i\parallel }_0^2.\hfill \end{array}$$

Then $k=2$, ${\overline{u}}_h^2=u_h^2-\rho u_h^1$, and

$$\begin{array}{c}\hfill \parallel {\overline{u}}_h^2{\parallel }_0^2+{\alpha }_0{\beta }_0^{-1}\parallel \nabla u_h^2{\parallel }_0^2\le 20d{\left[\mathsf{\Gamma }(1-\alpha )\right]}^2{T}^{\alpha }\underset{i}{max}{\parallel f^i\parallel }_0^2.\end{array}$$

(63)

When $k\ge 3$, letting $v_h=2{\overline{u}}_h^k$ in (56), we get

$$\begin{array}{ccc}\hfill & & 2\parallel {\overline{u}}_h^k{\parallel }_0^2+2{\alpha }_0{\beta }_0^{-1}(\nabla u_h^k,\nabla {\overline{u}}_h^k)\hfill \\ & =& 2\rho ({\overline{u}}_h^{k-1},{\overline{u}}_h^k)+2{\displaystyle \sum _{i=2}^{k-1}}{\overline{D}}_{k-i}^k({\overline{u}}_h^{k-i},{\overline{u}}_h^k)+2{\overline{D}}_0^k(u_h^0,{\overline{u}}^k)+2{\alpha }_0{\beta }_0^{-1}(f^k,{\overline{u}}_h^k).\hfill \end{array}$$

(64)

Using $(\nabla u_h^k,\nabla {\overline{u}}_h^k)=\parallel \nabla {\overline{u}}_h^k{\parallel }_0^2+\parallel \nabla u_h^k{\parallel }_0^2-{\rho }^2{\parallel \nabla u_h^{k-1}\parallel }_0^2$ and the integration by parts, we have

$$\begin{array}{ccc}& & 2\parallel {\overline{u}}_h^k{\parallel }_0^2+{\alpha }_0{\beta }_0^{-1}\parallel \nabla {\overline{u}}_h^k{\parallel }_0^2+{\alpha }_0{\beta }_0^{-1}\parallel \nabla u_h^k{\parallel }_0^2-{\alpha }_0{\beta }_0^{-1}{\rho }^2{\parallel \nabla u_h^{k-1}\parallel }_0^2\hfill \\ & =& 2\rho ({\overline{u}}_h^{k-1},{\overline{u}}_h^k)+2{\displaystyle \sum _{i=2}^{k-1}}{\overline{D}}_{k-i}^k({\overline{u}}_h^{k-i},{\overline{u}}_h^k)+2{\overline{D}}_0^k(u_h^0,{\overline{u}}_h^k)-2{\alpha }_0{\beta }_0^{-1}(I{f}^k,\nabla {\overline{u}}_h^k),\hfill \end{array}$$

where

$$\begin{array}{ccc}\hfill I{f}^k& =& ({\int }_0^{x_1}f(\tau ,x_2,\cdots ,x_d,t_k)d\tau ,{\int }_0^{x_2}f(x_1,\tau ,x_3,\cdots ,x_d,t_k)d\tau ,\hfill \\ & & \cdots ,{\int }_0^{x_d}f(x_1,\cdots ,x_{d-1},\tau ,t_k)d\tau )\hfill \end{array}$$

is a integral vector and $x_1,x_2,\cdots ,x_d\in (0,1)$.

According to Lemma 7, we know the $k-1$ coefficients are positive of the right hand side in (58), we have

$$\begin{array}{c}\hfill \begin{array}{ccc}& & {\displaystyle 2\parallel {\overline{u}}_h^k{\parallel }_0^2+{\alpha }_0{\beta }_0^{-1}\parallel \nabla {\overline{u}}_h^k{\parallel }_0^2+{\alpha }_0{\beta }_0^{-1}\parallel \nabla u_h^k{\parallel }_0^2-{\alpha }_0{\beta }_0^{-1}{\rho }^2{\parallel \nabla u_h^{k-1}\parallel }_0^2}\hfill \\ & \le \hfill & {\displaystyle \rho (\parallel {\overline{u}}_h^{k-1}{\parallel }_0^2+\parallel {\overline{u}}_h^k{\parallel }_0^2)+{\displaystyle \sum _{i=2}^{k-1}}{\overline{D}}_{k-i}^k(\parallel {\overline{u}}_h^{k-i}{\parallel }_0^2+\parallel {\overline{u}}_h^k{\parallel }_0^2)+{\overline{D}}_0^k(\parallel u_h^0{\parallel }_0^2+\parallel {\overline{u}}_h^k{\parallel }_0^2)}\hfill \\ & & {\displaystyle +{\alpha }_0{\beta }_0^{-1}\parallel I{f}^k{\parallel }_0^2+{\alpha }_0{\beta }_0^{-1}{\parallel \nabla {\overline{u}}_h^k\parallel }_0^2.}\hfill \end{array}\end{array}$$

(65)

According to (1)–(3) in the Lemma 8, we can find (65) is still satisfy for $k=3$. By directly computation, it can deduce that $\rho +{\overline{D}}_1^3+{\overline{D}}_2^3\le 1$. According to (3) in the Lemma 7, we have ${\displaystyle \rho +{\overline{D}}_0^k+{\displaystyle \sum _{i=2}^{k-1}}{\overline{D}}_{k-i}^k\le 1}$, then

$$\begin{array}{ccc}& & \parallel {\overline{u}}_h^k{\parallel }_0^2+{\alpha }_0{\beta }_0^{-1}\parallel \nabla u_h^k{\parallel }_0^2-{\alpha }_0{\beta }_0^{-1}{\rho }^2{\parallel \nabla u_h^{k-1}\parallel }_0^2\hfill \\ & \le & \rho \parallel {\overline{u}}_h^{k-1}{\parallel }_0^2+{\displaystyle \sum _{i=2}^{k-1}}{\overline{D}}_{k-i}^k\parallel {\overline{u}}_h^{k-i}{\parallel }_0^2+{\overline{D}}_0^k\parallel u_h^0{\parallel }_0^2+{\alpha }_0{\beta }_0^{-1}{\parallel I{f}^k\parallel }_0^2.\hfill \end{array}$$

We have

$$\begin{array}{ccc}& & \parallel {\overline{u}}_h^k{\parallel }_0^2+{\alpha }_0{\beta }_0^{-1}{\parallel \nabla u_h^k\parallel }_0^2\hfill \\ & \le & \rho (\parallel {\overline{u}}_h^{k-1}{\parallel }_0^2+{\alpha }_0{\beta }_0^{-1}\rho \parallel \nabla u_h^{k-1}{\parallel }_0^2)+{\displaystyle \sum _{i=2}^{k-1}}{\overline{D}}_{k-i}^k\parallel {\overline{u}}_h^{k-i}{\parallel }_0^2+{\overline{D}}_0^k\parallel u_h^0{\parallel }_0^2+{\alpha }_0{\beta }_0^{-1}{\parallel I{f}^k\parallel }_0^2.\hfill \end{array}$$

According to (1) and (4) in the Lemma 7, we have

$$\begin{array}{ccc}\hfill & & \parallel {\overline{u}}_h^k{\parallel }_0^2+{\alpha }_0{\beta }_0^{-1}{\parallel \nabla u_h^k\parallel }_0^2\hfill \\ \hfill & \le & \rho (\parallel {\overline{u}}_h^{k-1}{\parallel }_0^2+{\alpha }_0{\beta }_0^{-1}\rho \parallel \nabla u_h^{k-1}{\parallel }_0^2)+{\displaystyle \sum _{i=2}^{k-1}}{\overline{D}}_{k-i}^k(\parallel {\overline{u}}_h^{k-i}{\parallel }_0^2+{\alpha }_0{\beta }_0^{-1}\rho \parallel \nabla u_h^{k-i}{\parallel }_0^2)\hfill \\ & & +{\overline{D}}_0^k(\parallel u_h^0{\parallel }_0^2+\frac{T^{\alpha }\mathsf{\Gamma }(1-\alpha )}{{\beta }_0}\parallel I{f}^k{\parallel }_0^2).\hfill \end{array}$$

(66)

For $\parallel I{f}^k{\parallel }_0^2$ in the (66), we use Cauchy–Schwarz inequality, then

$$\begin{array}{c}\hfill \parallel I{f}^k{\parallel }_0^2\le \parallel [{\displaystyle \sum _{i=1}^d}{\int }_0^1|f(x_1,x_2,\cdots ,x_d,t_k){|}^{2}d{x}_i{]}^{\frac{1}{2}}{\parallel }_0^2=d{\parallel f(x_1,x_2,\cdots ,x_d,t_k)\parallel }_0^2.\end{array}$$

(67)

According to (1) in the Lemma 6 and (66), (67) becomes

$$\begin{array}{ccc}\hfill & & \parallel {\overline{u}}_h^k{\parallel }_0^2+{\alpha }_0{\beta }_0^{-1}{\parallel \nabla u_h^k\parallel }_0^2\hfill \\ \hfill & \le & \rho (\parallel {\overline{u}}_h^{k-1}{\parallel }_0^2+{\alpha }_0{\beta }_0^{-1}\rho \parallel \nabla u_h^{k-1}{\parallel }_0^2)+{\displaystyle \sum _{i=2}^{k-1}}{\overline{D}}_{k-i}^k(\parallel {\overline{u}}_h^{k-i}{\parallel }_0^2+{\alpha }_0{\beta }_0^{-1}\rho \parallel \nabla u_h^{k-i}{\parallel }_0^2)\hfill \\ \hfill & & +{\overline{D}}_0^k(\parallel u_h^0{\parallel }_0^2+\frac{T^{\alpha }\mathsf{\Gamma }(1-\alpha )d}{{\beta }_0}\parallel f^k{\parallel }_0^2)\hfill \\ \hfill & \le & \rho (\parallel {\overline{u}}_h^{k-1}{\parallel }_0^2+{\alpha }_0{\beta }_0^{-1}\parallel \nabla u_h^{k-1}{\parallel }_0^2)+{\displaystyle \sum _{i=2}^{k-1}}{\overline{D}}_{k-i}^k(\parallel {\overline{u}}_h^{k-i}{\parallel }_0^2+{\alpha }_0{\beta }_0^{-1}\parallel \nabla u_h^{k-i}{\parallel }_0^2)\hfill \\ & & +{\overline{D}}_0^k(\parallel u_h^0{\parallel }_0^2+T^{\alpha }\mathsf{\Gamma }(1-\alpha )d\underset{i}{max}\parallel f^i{\parallel }_0^2).\hfill \end{array}$$

(68)

Next, we will prove the following estimate using mathematics induction:

$$\begin{array}{c}\hfill \parallel {\overline{u}}_h^k{\parallel }_0^2+{\alpha }_0{\beta }_0^{-1}\parallel \nabla u_h^k{\parallel }_0^2\le 20T^{\alpha }\mathsf{\Gamma }(1-\alpha )d\underset{1\le i\le k}{max}{\parallel f^i\parallel }_0^2.\end{array}$$

(69)

According to (62) and (63), we can easily check (69) for $k=1,2$. Therefore assuming (69) is through for $j=1,2,\cdots ,k-1$:

$$\begin{array}{c}\hfill \parallel {\overline{u}}_h^j{\parallel }_0^2+{\alpha }_0{\beta }_0^{-1}\parallel \nabla u_h^j{\parallel }_0^2\le 20T^{\alpha }\mathsf{\Gamma }(1-\alpha )d\underset{1\le i\le k}{max}{\parallel f^i\parallel }_0^2,1\le j\le k-1.\end{array}$$

We deduce from (68),

$$\begin{array}{c}\hfill \parallel {\overline{u}}_h^k{\parallel }_0^2+{\alpha }_0{\beta }_0^{-1}\parallel \nabla u_h^k{\parallel }_0^2\le (\rho +{\displaystyle \sum _{i=2}^{k-1}}{\overline{D}}_{k-i}^k+{\overline{D}}_0^k)(20T^{\alpha }\mathsf{\Gamma }(1-\alpha )d\underset{1\le i\le k}{max}\parallel f^i{\parallel }_0^2).\end{array}$$

Then (69) is proven, i.e.,

$$\begin{array}{ccc}\hfill \parallel {\overline{u}}_h^k\parallel & =& \parallel u_h^k-\rho u_h^{k-1}\parallel \le (20T^{\alpha }\mathsf{\Gamma }(1-\alpha )d\underset{1\le i\le k}{max}\parallel f^i{{\parallel }_0^2)}^{\frac{1}{2}}\hfill \\ & =& \sqrt{20T^{\alpha }\mathsf{\Gamma }(1-\alpha )d}\underset{1\le i\le k}{max}{\parallel f^i\parallel }_0.\hfill \end{array}$$

(70)

Finally, we turn to estimate $\parallel u_h^k\parallel$. Applying the triangle inequality and (70) yields

$$\begin{array}{ccc}& & \parallel u_h^k\parallel =\parallel {\overline{u}}_h^k+\rho u_h^{k-1}\parallel \le \rho \parallel u_h^{k-1}\parallel +\sqrt{20T^{\alpha }\mathsf{\Gamma }(1-\alpha )d}\underset{1\le i\le k}{max}{\parallel f^i\parallel }_0\hfill \\ & \le & \rho (\parallel u_h^{k-2}\parallel +\sqrt{20d{T}^{\alpha }\mathsf{\Gamma }(1-\alpha )}\underset{1\le i\le k}{max}\parallel f^i{\parallel }_0)+\sqrt{20d{T}^{\alpha }\mathsf{\Gamma }(1-\alpha )}\underset{1\le i\le k}{max}{\parallel f^i\parallel }_0\hfill \\ & \le & \cdots \le (1+\rho +\cdots +{\rho }^{k-2}+{\rho }^{k-1})\sqrt{20d{T}^{\alpha }\mathsf{\Gamma }(1-\alpha )}\underset{1\le i\le k}{max}{\parallel f^i\parallel }_0\hfill \\ & \le & \frac{1}{1-\rho }\sqrt{20d{T}^{\alpha }\mathsf{\Gamma }(1-\alpha )}\underset{1\le i\le k}{max}\parallel f^i{\parallel }_0\le \sqrt{20d{T}^{\alpha }\mathsf{\Gamma }(1-\alpha )}\underset{1\le i\le k}{max}{\parallel f^i\parallel }_0.\hfill \end{array}$$

The proof is completed. □

4.2. The Adjoint Equation’s FD-FE Scheme

In this part, we will analysis the FOCP’s full discretization. For the cost functional discretization form is defined by

$$\begin{array}{ccc}\hfill & & J_{h,\tau }({\mathbf{u}}_h,{\mathbf{q}}_h)=\frac{2\tau }{6}(\parallel u_h^0-u_d^0{\parallel }_{L^2\left(\mathsf{\Omega }\right)}^2+{\displaystyle \sum _{l=0}^{N-1}}4{\parallel u_h^{2l+1}-u_d^{2l+1}\parallel }_{L^2\left(\mathsf{\Omega }\right)}^2\hfill \\ & & +{\displaystyle \sum _{l=0}^{N-1}}2\parallel u_h^{2l}-u_d^{2l}{\parallel }_{L^2\left(\mathsf{\Omega }\right)}^2+\parallel u_h^{2N}-u_d^{2N}{\parallel }_{L^2\left(\mathsf{\Omega }\right)}^2+\gamma \parallel q_h^0{\parallel }_{L^2\left(\mathsf{\Omega }\right)}^2+\gamma {\displaystyle \sum _{l=0}^{N-1}}4{\parallel q_h^{2l+1}\parallel }_{L^2\left(\mathsf{\Omega }\right)}^2\hfill \\ & & +\gamma {\displaystyle \sum _{l=0}^{N-1}}2\parallel q_h^{2l}{\parallel }_{L^2\left(\mathsf{\Omega }\right)}^2+\gamma {\parallel q_h^{2N}\parallel }_{L^2\left(\mathsf{\Omega }\right)}^2)\doteq {({\mathbf{u}}_h,{\mathbf{q}}_h)}_{{\mathsf{\Omega }}_{T,\tau }},\hfill \end{array}$$

(71)

here the Simpson rule was used to make the time integral of the cost function discrete, and ${\mathbf{u}}_h=(u_h^1,\cdots ,u_h^N),{\mathbf{q}}_h=(q_h^1,\cdots ,q_h^N)$.

Using (71), one obtain the full discretization of the FOCP (1) and (2), finding $({\mathbf{u}}_h,{\mathbf{q}}_h)\in {\mathcal{V}}_h^N\times {\mathcal{V}}_h^N$, such that

$$\begin{array}{c}\hfill \underset{{\mathbf{q}}_h\in {\mathcal{V}}_h^N}{min}{J}_{h,\tau }({\mathbf{u}}_h,{\mathbf{q}}_h)\end{array}$$

(72)

subject to

$$\begin{array}{c}\hfill \left\{\begin{array}{c}{(}_0{D}_{\tau }^{\alpha }{u}_h^k,v_h)+(\nabla u_h^k,\nabla v_h)=(f^k+q_h^k,v_h),\forall v_h\in {\mathcal{V}}_h,\hfill \\ u^0=0,k=1,2,\cdots ,N,\hfill \end{array}\right.\end{array}$$

(73)

where the control variable was implicity discretized by variational discretization concept.

Similar to (73), we construct the numerical scheme for (7) as follows

$$\begin{array}{c}\hfill \left\{\begin{array}{c}{(}_{\tau }{D}_T^{\alpha }{Z}_h^k,v_h)+(\nabla Z_h^k,\nabla v_h)=(u_h^{k+1}-u_d^{k+1},v_h),\forall v_h\in {\mathcal{V}}_h,\hfill \\ Z_h^N=0,k=0,\cdots ,N-1,\hfill \end{array}\right.\end{array}$$

(74)

where ${}_{\tau }{D}_T^{\alpha }{Z}_h^k$ is the discrete right Caputo derivative as follows

$$\begin{array}{c}\hfill \begin{array}{c}{\displaystyle {}_{\tau }{D}_T^{\alpha }{Z}^k=\left\{\begin{array}{c}{\displaystyle \frac{{\tau }^{-\alpha }}{\mathsf{\Gamma }(3-\alpha )}({\overline{E}}_0{Z}_h^{N-2}+{\overline{E}}_1{Z}_h^{N-1}+{\overline{E}}_2{Z}_h^N),j=N-1,}\hfill \\ {\displaystyle \frac{{\tau }^{-\alpha }}{\mathsf{\Gamma }(3-\alpha )}(E_0{Z}_h^{N-2}+E_1{Z}_h^{N-1}+E_2{Z}_h^N),j=N-2,}\hfill \\ {\displaystyle \frac{{\tau }^{-\alpha }}{\mathsf{\Gamma }(3-\alpha )}[{\overline{F}}_j{Z}_h^{N-2}+{\overline{G}}_j{Z}_h^{N-1}+{\overline{H}}_j{Z}_h^N}\hfill \\ +{\displaystyle {\displaystyle \sum _{k=j+1}^{N-1}}}(F_k{Z}_h^{k-1-j}+G_k{Z}_h^{k-j}+H_k{Z}_h^{k+1-j})],\hfill \\ j\le N-3;0\le k\le N-1,\hfill \end{array}\right.}\hfill \end{array}\end{array}$$

(75)

where

$$\begin{array}{c}\hfill \begin{array}{c}{\displaystyle {\overline{E}}_0=\frac{\alpha }{2},{\overline{E}}_1=2-2\alpha ,{\overline{E}}_2=\frac{3\alpha -4}{2},E_0={\frac{\alpha +2}{2}}^{\alpha },E_1=-\frac{4\alpha }{2^{\alpha }},E_2=\frac{3\alpha -2}{2^{\alpha }},}\hfill \\ {\displaystyle {\overline{F}}_j=-\frac{2-\alpha }{2}[{(N-j-1)}^{1-\alpha }+{(N-j)}^{1-\alpha }]+{(N-j)}^{2-\alpha }-{(N-j-1)}^{2-\alpha },}\hfill \\ {\displaystyle {\overline{G}}_j=-2{(N-j)}^{2-\alpha }+2(2-\alpha ){(N-j)}^{1-\alpha }+2{(N-j-1)}^{2-\alpha },}\hfill \\ {\displaystyle {\overline{H}}_j=-\frac{3(2-\alpha )}{2}{(N-j)}^{1-\alpha }+\frac{2-\alpha }{2}{(N-j-1)}^{1-\alpha }+{(N-j)}^{2-\alpha }-{(N-j-1)}^{2-\alpha },}\hfill \\ {\displaystyle F_k=-\frac{3}{2}(2-\alpha ){(k-1)}^{1-\alpha }+\frac{1}{2}(2-\alpha )k^{1-\alpha }+k^{2-\alpha }-{(k-1)}^{2-\alpha },}\hfill \\ {\displaystyle G_k=2{(k-1)}^{2-\alpha }+2(2-\alpha ){(k-1)}^{1-\alpha }-2k^{2-\alpha },}\hfill \\ {\displaystyle H_k=-\frac{1}{2}(2-\alpha )[k^{1-\alpha }+{(k-1)}^{1-\alpha }]+k^{2-\alpha }-{(k-1)}^{2-\alpha }.}\hfill \end{array}\end{array}$$

(76)

For the Equation (19), we have the following discretization scheme

$$\begin{array}{c}\hfill \gamma q_h^{k+1}+Z_h^k=0,k=0,1,\cdots ,N-1.\end{array}$$

(77)

Therefore, the discrete optimality conditions of FOCP (1) and (2) are given by

$$\begin{array}{c}\hfill \left\{\begin{array}{c}{(}_0{D}_{\tau }^{\alpha }{u}_h^{k+1},v_h)+(\nabla u_h^{k+1},\nabla v_h)=(f^{k+1}+q_h^{k+1},v_h),\forall v_h\in {\mathcal{V}}_h,\hfill \\ {(}_{\tau }{D}_T^{\alpha }{Z}_h^k,v_h)+(\nabla Z_h^k,\nabla v_h)=(u_h^{k+1}-u_d^{k+1},v_h),\forall v_h\in {\mathcal{V}}_h,\hfill \\ (\gamma q_h^{k+1}+Z_h^k,v_h)=0,\forall v_h\in {\mathcal{V}}_h,\hfill \\ u_h^0=0,Z_h^N=0,0\le k\le N-1.\hfill \end{array}\right.\end{array}$$

(78)

Next, we are going to give the error of the FD-FE scheme for the FOCP based on the idea of [23] and the above results.

Theorem 3.

Let $(u,Z,q)$ and $(u_h^k,Z_h^k,q_h^k)$ be the solution of (2), (5), (7) and (78), respectively. Then the estimate as follows

$$\begin{array}{c}\hfill \parallel q(x,t_k)-q_h^k{\parallel }_{L^2\left(\mathsf{\Omega }\right)}+\parallel u(x,t_k)-u_h^k{\parallel }_{L^2\left(\mathsf{\Omega }\right)}+{\parallel Z(x,t_k)-Z_h^k\parallel }_{L^2\left(\mathsf{\Omega }\right)}\\ \hfill +h\parallel \nabla (u(x,t_k)-u_h^k){\parallel }_{L^2\left(\mathsf{\Omega }\right)}+h{\parallel \nabla (Z(x,t_k)-Z_h^k)\parallel }_{L^2\left(\mathsf{\Omega }\right)}\\ \hfill \le C(h^2+{\tau }^{3-\alpha }),\forall k\in \{1,2,\cdots ,N\}\end{array}$$

hold for $u,Z\in W^{3,\infty }(0,T;H_0^1\left(\mathsf{\Omega }\right)),q\in L^{\infty }(0,T;H_0^1\left(\mathsf{\Omega }\right)\bigcap H^2\left(\mathsf{\Omega }\right))$, and the constant C is independent of $h,\tau$.

5. Numerical Examples

5.1. The Conjugate Gradient (CG) of Optimization Algorithm

In this part, we will carry out two numerical examples to show the prior error estimates of numerical scheme to the FOCP (1) and (2) in previous section. The details of CG algorithm for the (1) and (2) optimization problem is described as follows.

We hereafter denote $l=1,2,\cdots ,N$ without explanation. Let ${\mathbf{q}}_h^{\left(0\right)}=(q_h^{\left(0\right),1},\cdots ,q_h^{\left(0\right),N})$ be the initial value of control variable and ${\mathbf{u}}_h\left({\mathbf{q}}_h^{\left(0\right)}\right)=(u_h^1\left(q_h^{\left(0\right),1}\right),\cdots ,u_h^N\left(q_h^{\left(0\right),N}\right))$ be the corresponding state variable defined by (18) which is semidiscretization of the FOCP. Let $\parallel J^{{}^{\prime }}\left({\mathbf{q}}_h^{\left(0\right)}\right){\parallel }_{{\mathsf{\Omega }}_{T,\tau }}\le ϵ$ be the stopping criterion with $ϵ$ being a tolerance. The adjoint state ${\mathbf{Z}}_h\left({\mathbf{q}}_h^{\left(0\right)}\right)=(Z_h^1\left(q_h^{\left(0\right),1}\right),\cdots ,Z_h^N\left(q_h^{\left(0\right),N}\right))$ can be obtained from the adjoint state Equation (20) when ${\mathbf{u}}_h\left({\mathbf{q}}_h^{\left(0\right)}\right)$ and ${\mathbf{q}}_h^{\left(0\right)}$ are given. The objective function’s gradient at ${\mathbf{q}}_h^{\left(0\right)}$ is defined by

$$\begin{array}{c}\hfill {\mathbf{d}}_h^{\left(0\right)}=J^{{}^{\prime }}\left({\mathbf{q}}_h^{\left(0\right)}\right)={\mathbf{Z}}_h\left(q_h^{\left(0\right)}\right)+\gamma {\mathbf{q}}_h^{\left(0\right)}.\end{array}$$

We choose the initial conjugate direction such that it is the same as the gradient direction, namely

$$\begin{array}{c}\hfill {\mathbf{s}}_h^{\left(0\right)}={\mathbf{d}}_h^{\left(0\right)}.\end{array}$$

Supposing the $k-th$ iteration ${\mathbf{q}}_h^{\left(k\right)}$, ${\mathbf{d}}_h^{\left(k\right)}$ and ${\mathbf{s}}_h^{\left(k\right)}$ are known, we update ${\mathbf{q}}_h^{(k+1)}$ via

$$\begin{array}{c}\hfill {\mathbf{q}}_h^{(k+1)}={\mathbf{q}}_h^{\left(k\right)}-{\rho }_k{\mathbf{s}}_h^{\left(k\right)},\end{array}$$

where ${\mathbf{d}}_h^{\left(k\right)}=(d_h^{\left(k\right),1},\cdots ,d_h^{\left(k\right),N})$, ${\mathbf{s}}_h^{\left(k\right)}=(s_h^{\left(k\right),1},\cdots ,s_h^{\left(k\right),N})$ and the iteration step size ${\rho }_k$ is determined by

$$\begin{array}{c}\hfill J_{h,\tau }({\mathbf{q}}_h^{\left(k\right)}-{\rho }_k{\mathbf{d}}_h^{\left(k\right)})=\underset{\rho >0}{min}{J}_{h,\tau }({\mathbf{q}}_h^{\left(k\right)}-\rho {\mathbf{d}}_h^{\left(k\right)}).\end{array}$$

Due to ${({\mathbf{d}}_h^{(k+1)},{\mathbf{s}}_h^{\left(k\right)})}_{{\mathsf{\Omega }}_{T,\tau }}=0$ and

$$\begin{array}{c}\hfill {\mathbf{d}}_h^{(k+1)}={\mathbf{Z}}_h\left(q_h^{(k+1)}\right)+\gamma {\mathbf{q}}_h^{(k+1)}={\mathbf{Z}}_h\left(q_h^{(k+1)}\right)+\gamma ({\mathbf{q}}_h^{\left(k\right)}-{\rho }_k{\mathbf{s}}_h^{\left(k\right)}),\end{array}$$

where ${\rho }_k$ is characterized as

$$\begin{array}{c}\hfill {({\mathbf{Z}}_h\left({\mathbf{q}}_h^{(k+1)}\right)+\gamma ({\mathbf{q}}_h^{\left(k\right)}-{\rho }_k{\mathbf{s}}_h^{\left(k\right)}),{\mathbf{s}}_h^{\left(k\right)})}_{{\mathsf{\Omega }}_{T,\tau }}=0,\end{array}$$

(79)

here ${\mathbf{Z}}_h\left({\mathbf{q}}_h^{(k+1)}\right)=(Z_h^1\left(q_h^{(k+1),1}\right),\cdots ,Z_h^N\left(q_h^{(k+1),N}\right))$ is the finite element solution of (74) and ${\mathbf{u}}_h^{(k+1)}\in {\mathcal{V}}_h^N$ given by (73).

The optimal $k-th$ iteration ${\rho }_k$ can be solved efficiently by (79). Indeed, the adjoint state ${\mathbf{Z}}_h\left({\mathbf{q}}_h^{(k+1)}\right)$ dependent on ${\rho }_k$, let ${\tilde{\mathbf{u}}}_h^{\left(k\right)}=({\tilde{u}}_h^{\left(k\right),1},\cdots ,{\tilde{u}}_h^{\left(k\right),N})$, ${\tilde{\mathbf{Z}}}_h^{\left(k\right)}=({\tilde{Z}}_h^{\left(k\right),1},\cdots ,{\tilde{Z}}_h^{\left(k\right),N})$ denote respectively, the solution of

$$\begin{array}{c}\hfill {(}_0{D}_{\tau }^{\alpha }{\tilde{u}}_h^{\left(k\right),l},v_h)+(\nabla {\tilde{u}}_h^{\left(k\right),l},\nabla v_h)=(s_h^{\left(k\right),l},v_h),\forall v_h\in {\mathcal{V}}_h,\end{array}$$

(80)

$$\begin{array}{c}\hfill {(}_{\tau }{D}_T^{\alpha }{\tilde{Z}}_h^{\left(k\right),l},v_h)+(\nabla {\tilde{Z}}_h^{\left(k\right),l},\nabla v_h)=({\tilde{u}}_h^{\left(k\right),l},v_h).\end{array}$$

(81)

$u_h^l\left(q_h^{\left(k\right),l}\right)$ and $Z_h^l\left(q_h^{\left(k\right),l}\right)$ are, respectively, the solutions of the following equations

$$\begin{array}{c}\hfill {(}_0{D}_{\tau }^{\alpha }{u}_h^l\left(q_h^{\left(k\right),l}\right),v_h)+(\nabla u_h^l\left(q_h^{\left(k\right),l}\right),\nabla v_h)=(f^l+q_h^{\left(k\right),l},v_h),\forall v_h\in {\mathcal{V}}_h,\end{array}$$

(82)

$$\begin{array}{c}\hfill {(}_{\tau }{D}_T^{\alpha }{Z}_h^l\left(q_h^{\left(k\right),l}\right),v_h)+(\nabla Z_h^l\left(q_h^{\left(k\right),l}\right),\nabla v_h)=(u_h^l\left(q_h^{\left(k\right),l}\right)-u_d^l,v_h),\forall v_h\in {\mathcal{V}}_h.\end{array}$$

(83)

Then it can be checked that ${\mathbf{Z}}_h\left({\mathbf{q}}_h^{\left(k\right)}\right)-{\rho }_k{\tilde{\mathbf{Z}}}_h^{\left(k\right)}$ solves (80)–(83), that is

$$\begin{array}{c}\hfill {\mathbf{Z}}_h\left({\mathbf{q}}_h^{(k+1)}\right)={\mathbf{Z}}_h\left({\mathbf{q}}_h^{\left(k\right)}\right)-{\rho }_k{\tilde{\mathbf{Z}}}_h^{\left(k\right)}.\end{array}$$

Putting this expression into (79) gives

$$\begin{array}{c}\hfill {({\mathbf{Z}}_h\left({\mathbf{q}}_h^{\left(k\right)}\right)-{\rho }_k{\tilde{\mathbf{Z}}}_h^{\left(k\right)}+\gamma ({\mathbf{q}}_h^{\left(k\right)}-{\rho }_k{\mathbf{s}}_h^{\left(k\right)}),{\mathbf{s}}_h^{\left(k\right)})}_{{\mathsf{\Omega }}_{T,\tau }}=0.\end{array}$$

Obviously, ${\mathbf{d}}_h^{\left(k\right)}={\mathbf{Z}}_h\left({\mathbf{q}}_h^{\left(k\right)}\right)+\gamma {\mathbf{q}}_h^{\left(k\right)}$ holds. Let ${\tilde{\mathbf{d}}}_h^{\left(k\right)}={\tilde{\mathbf{Z}}}_h^{\left(k\right)}+{\mathbf{s}}_h^{\left(k\right)}$, then we obtain

$$\begin{array}{c}\hfill {\rho }_k=\frac{{({\mathbf{d}}_h^{\left(k\right)},{\mathbf{s}}_h^{\left(k\right)})}_{{\mathsf{\Omega }}_{T,\tau }}}{{({\tilde{\mathbf{d}}}_h^{\left(k\right)},{\mathbf{s}}_h^{\left(k\right)})}_{{\mathsf{\Omega }}_{T,\tau }}}.\end{array}$$

(84)

Furthermore, we can check that ${\mathbf{d}}_h^{(k+1)}={\mathbf{d}}_h^{\left(k\right)}-{\rho }_k·{\tilde{\mathbf{d}}}_h^{\left(k\right)}$ holds, denote

$$\begin{array}{c}\hfill {\beta }_k=\frac{\parallel {\mathbf{d}}_h^{(k+1)}{\parallel }_{0,{\mathsf{\Omega }}_{T,\tau }}^2}{\parallel {\mathbf{d}}_h^{\left(k\right)}{\parallel }_{0,{\mathsf{\Omega }}_{T,\tau }}^2}\end{array}$$

and we choose the ${\mathbf{s}}_h^{(k+1)}$ defined by the following equation

$$\begin{array}{c}\hfill {\mathbf{s}}_h^{(k+1)}={\mathbf{d}}_h^{(k+1)}+{\beta }_k{\mathbf{s}}_h^{\left(k\right)}.\end{array}$$

Based on the above fact, we have ${({\mathbf{d}}_h^{\left(k\right)},{\mathbf{s}}_h^{(k-1)})}_{{\mathsf{\Omega }}_{T,\tau }}=0.$ It can improve the optimimum k-th iterative step size ${\rho }_k$, which is defined as follows

$$\begin{array}{c}\hfill {\rho }_k=\frac{{({\mathbf{d}}_h^{\left(k\right)},{\mathbf{s}}_h^{\left(k\right)})}_{{\mathsf{\Omega }}_{T,\tau }}}{{({\tilde{\mathbf{d}}}_h^{\left(k\right)},{\mathbf{s}}_h^{\left(k\right)})}_{{\mathsf{\Omega }}_{T,\tau }}}=\frac{{({\mathbf{d}}_h^{\left(k\right)},{\mathbf{d}}_h^{\left(k\right)}+{\beta }_{k-1}{\mathbf{s}}_h^{(k-1)})}_{{\mathsf{\Omega }}_{T,\tau }}}{{({\tilde{\mathbf{d}}}_h^{\left(k\right)},{\mathbf{s}}_h^{\left(k\right)})}_{{\mathsf{\Omega }}_{T,\tau }}}=\frac{{({\mathbf{d}}_h^{\left(k\right)},{\mathbf{d}}_h^{\left(k\right)})}_{{\mathsf{\Omega }}_{T,\tau }}}{{({\tilde{\mathbf{d}}}_h^{\left(k\right)},{\mathbf{s}}_h^{\left(k\right)})}_{{\mathsf{\Omega }}_{T,\tau }}}.\end{array}$$

(85)

The overall process of CG algorithm as summarized below.

The CG optimization algorithm. Choosing ${\mathbf{q}}_h^{\left(0\right)}$ for the initial value of control variable.

(I)

Solving problems (82), (83), let ${\mathbf{d}}_h^{\left(0\right)}={\mathbf{Z}}_h\left({\mathbf{q}}_h^{\left(0\right)}\right)+\gamma {\mathbf{q}}_h^{\left(0\right)},{\mathbf{s}}_h^{\left(0\right)}={\mathbf{d}}_h^{\left(0\right)}$. Set $k=0$.

(II)

Solving problems (80), (81), and set ${\tilde{\mathbf{d}}}_h^{\left(k\right)}={\tilde{\mathbf{Z}}}_h^{\left(k\right)}+{\mathbf{s}}_h^{\left(k\right)}$, ${\rho }_k=\frac{{({\mathbf{d}}_h^{\left(k\right)},{\mathbf{d}}_h^{\left(k\right)})}_{{\mathsf{\Omega }}_{T,\tau }}}{{({\tilde{\mathbf{d}}}_h^{\left(k\right)},{\mathbf{s}}_h^{\left(k\right)})}_{{\mathsf{\Omega }}_{T,\tau }}}.$

(III)

Update ${\mathbf{q}}_h^{(k+1)}={\mathbf{q}}_h^{\left(k\right)}-{\rho }_k{\mathbf{s}}_h^{\left(k\right)}$, ${\mathbf{d}}_h^{(k+1)}={\mathbf{d}}_h^{\left(k\right)}-{\rho }_k{\tilde{\mathbf{d}}}_h^{\left(k\right)}$.

(IV)

If $\parallel {\mathbf{d}}_h^{(k+1)}{\parallel }_{{\mathsf{\Omega }}_{T,\tau }}\le$ tolerence, then ${\mathbf{q}}_h^{*}={\mathbf{q}}_h^{(k+1)}$, and solve problems (73) and (74) to get ${\mathbf{u}}_h\left({\mathbf{q}}_h^{*}\right)$ and ${\mathbf{Z}}_h\left({\mathbf{q}}_h^{*}\right)$. Else, let ${\beta }_k=\frac{\parallel {\mathbf{d}}_h^{(k+1)}{\parallel }_{{\mathsf{\Omega }}_{T,\tau }}^2}{\parallel {\mathbf{d}}_h^{\left(k\right)}{\parallel }_{{\mathsf{\Omega }}_{T,\tau }}^2}$, ${\mathbf{s}}_h^{(k+1)}={\mathbf{d}}_h^{(k+1)}+{\beta }_k{\mathbf{s}}_h^{\left(k\right)}$. Set $k\leftarrow k+1$, repeat (II)–(IV).

5.2. Numerical Results

In this part, we give two numerical examples to show that the theorem results are correct, which are 1D and 2D FOCP, respectively. In all the following examples, we take $T=1$ and $\gamma =1$. The following two examples were implemented on a LAPTOP-H91AOQNL computer with a Intel(R) Core(TM) i7-10510U CPU @ 1.80GHz and 12.00 GB of RAM by using MATLAB.

Example 1.

Considering the problem(1) and (2) with the desired state $u_d(x,t)$ and the right function $f(x,t)$ be defined as follows

$$\begin{array}{ccc}& & u_d(x,t)=-4{\pi }^2{(1-t)}^{4}sin\left(2\pi x\right)+t^{4}sin\left(2\pi x\right)-\frac{24{(1-t)}^{4-\alpha }}{\mathsf{\Gamma }(5-\alpha )}sin\left(2\pi x\right),\hfill \\ & & f(x,t)={(1-t)}^{4}sin\left(2\pi x\right)+t^{4}4{\pi }^{2}sin\left(2\pi x\right)+\frac{24t^{4-\alpha }}{\mathsf{\Gamma }(5-\alpha )}sin\left(2\pi x\right).\hfill \end{array}$$

After direct calculation, we obtain exact analytical solution state variable and control variable:

$$\begin{array}{c}\hfill u(x,t)=t^{4}sin\left(2\pi x\right),q\left(t\right)=-{(1-t)}^{4}sin\left(2\pi x\right).\end{array}$$

Let the error of u be defined as following: ${\displaystyle e_u=\underset{i}{max}\left|\right|u_h^i\left(x\right)-u(x,t_i){\left|\right|}_{L^2\left(\mathsf{\Omega }\right)}}$, and the error of q be ${\displaystyle e_q=\underset{i}{max}\left|\right|q_h^i\left(x\right)-q(x,t_i){\left|\right|}_{L^2\left(\mathsf{\Omega }\right)}.}$

From

Table 1. Error $e_q$ and $e_u$ for the spatial convergence rates with $\alpha =0.3$.

h	${\mathit{e}}_{\mathit{q}}$	Rate	${\mathit{e}}_{\mathit{u}}$	Rate
$\frac{1}{8}$	5.09132788$\times {10}^{-2}$	-	5.09260622$\times {10}^{-2}$	-
$\frac{1}{16}$	1.24529070$\times {10}^{-2}$	2.03155941	1.24543702$\times {10}^{-2}$	2.03175209
$\frac{1}{32}$	3.09634055$\times {10}^{-3}$	2.00784650	3.09650807$\times {10}^{-3}$	2.00793795
$\frac{1}{64}$	7.73008270$\times {10}^{-4}$	2.00200840	7.73051611$\times {10}^{-4}$	2.00200556
$\frac{1}{128}$	1.93158032$\times {10}^{-4}$	2.00070217	1.93195589$\times {10}^{-4}$	2.00050257
$\frac{1}{256}$	4.82568312$\times {10}^{-5}$	2.00097658	4.82961118$\times {10}^{-4}$	2.00008320

,

Table 2. Error $e_q$ and $e_u$ for the spatial convergence rates with $\alpha =0.5$.

h	${\mathit{e}}_{\mathit{q}}$	Rate	${\mathit{e}}_{\mathit{u}}$	Rate
$\frac{1}{8}$	5.02474475$\times {10}^{-2}$	-	5.02718614$\times {10}^{-2}$	-
$\frac{1}{16}$	1.22959862$\times {10}^{-2}$	2.03086285	1.22994227$\times {10}^{-2}$	2.03116050
$\frac{1}{32}$	3.05765198$\times {10}^{-3}$	2.00769134	3.05818033$\times {10}^{-3}$	2.00784521
$\frac{1}{64}$	7.63342242$\times {10}^{-4}$	2.00202227	7.63490535$\times {10}^{-4}$	2.00199130
$\frac{1}{128}$	1.90714809$\times {10}^{-4}$	2.00091315	1.90812917$\times {10}^{-4}$	2.00045143
$\frac{1}{256}$	4.76172991$\times {10}^{-5}$	2.00185918	4.77085028$\times {10}^{-5}$	1.99984051

and

Table 3. Error $e_q$ and $e_u$ for the spatial convergence rates with $\alpha =0.7$.

h	${\mathit{e}}_{\mathit{q}}$	Rate	${\mathit{e}}_{\mathit{u}}$	Rate
$\frac{1}{16}$	4.94152479$\times {10}^{-2}$	-	4.94428033$\times {10}^{-2}$	-
$\frac{1}{32}$	1.20997787$\times {10}^{-2}$	2.02997560	1.21038816$\times {10}^{-2}$	2.03029075
$\frac{1}{64}$	3.00932932$\times {10}^{-3}$	2.00746677	3.00999385$\times {10}^{-3}$	2.00763734
$\frac{1}{128}$	7.51326636$\times {10}^{-4}$	2.00192983	7.51521618$\times {10}^{-4}$	2.00187403
$\frac{1}{256}$	1.87735003$\times {10}^{-4}$	2.00074258	1.87864365$\times {10}^{-4}$	2.00012316
$\frac{1}{512}$	4.68939021$\times {10}^{-5}$	2.00122542	4.70134232$\times {10}^{-5}$	1.99854680

, we take $\tau =2^{-9};h=2^{-l},l=4,5,6,7,8,9$, $\alpha =0.3,0.5,0.7$, respectively. From Table 1, Table 2 and Table 3, we find the spatial accuracy is 2, this result is accord with the theoretical analysis obtained in Theorem 3.

Next, we check the temporal convergence rate. In

Table 4. Error $e_q$ and $e_u$ for the time convergence rates with $\alpha =0.3$.

$\mathit{\tau }$	${\mathit{e}}_{\mathit{q}}$	Rate	${\mathit{e}}_{\mathit{u}}$	Rate
$\frac{1}{32}$	1.81302869$\times {10}^{-3}$	-	1.81311070$\times {10}^{-3}$	-
$\frac{1}{64}$	2.75010932$\times {10}^{-4}$	2.72084088	2.75021268$\times {10}^{-4}$	2.72085192
$\frac{1}{128}$	4.27236424$\times {10}^{-5}$	2.68638241	4.27420180$\times {10}^{-5}$	2.68581626
$\frac{1}{256}$	6.54546760$\times {10}^{-6}$	2.70646648	6.57125488$\times {10}^{-6}$	2.70141421
$\frac{1}{512}$	9.77513735$\times {10}^{-7}$	2.74330738	1.01038698$\times {10}^{-6}$	2.70126093

,

Table 5. Error $e_q$ and $e_u$ for the time convergence rates with $\alpha =0.5$.

$\mathit{\tau }$	${\mathit{e}}_{\mathit{q}}$	Rate	${\mathit{e}}_{\mathit{u}}$	Rate
$\frac{1}{32}$	3.12780360$\times {10}^{-3}$	-	3.12813403$\times {10}^{-3}$	-
$\frac{1}{64}$	5.54169271$\times {10}^{-4}$	2.49675130	5.54220128$\times {10}^{-4}$	2.49677131
$\frac{1}{128}$	9.77009694$\times {10}^{-5}$	2.50388193	9.77492654$\times {10}^{-5}$	2.50330134
$\frac{1}{256}$	1.7184216$\times {10}^{-5}$	2.50728883	1.72467119$\times {10}^{-5}$	2.50276451
$\frac{1}{512}$	2.97915718$\times {10}^{-6}$	2.52810789	3.05597106$\times {10}^{-6}$	2.49661855

and

Table 6. Error $e_q$ and $e_u$ for the time convergence rates with $\alpha =0.7$.

$\mathit{\tau }$	${\mathit{e}}_{\mathit{q}}$	Rate	${\mathit{e}}_{\mathit{u}}$	Rate
$\frac{1}{32}$	5.54807629$\times {10}^{-3}$	-	5.54989898$\times {10}^{-3}$	-
$\frac{1}{64}$	1.09332414$\times {10}^{-3}$	2.34326644	1.09360389$\times {10}^{-3}$	2.34337123
$\frac{1}{128}$	2.25329722$\times {10}^{-4}$	2.27861164	2.25445955$\times {10}^{-4}$	2.27823674
$\frac{1}{256}$	4.53794774$\times {10}^{-5}$	2.31192572	4.54816054$\times {10}^{-5}$	2.30942654
$\frac{1}{512}$	9.12804809$\times {10}^{-6}$	2.31366169	9.24005981$\times {10}^{-6}$	2.29930908

, we take $\alpha =0.3,0.5,0.7$, respectively. We let $\tau =2^{-l},l=5,6,7,8,9$ and list the value of $e_{\tau ,h}$ and the corresponding order when $\alpha$, $\tau$ takes a series of different values, where $h=O\left({\tau }^{\frac{3-\alpha }{2}}\right)$ is taken. When $\alpha$ takes $0.3,0.5$, and $0.7$, the convergence order tends to $2.7,2.5$, and $2.3$, respectively, this can show that the time’s convergence rate is about $3-\alpha$.

Next, we plot the numerical solution for u and q with the conditions of $N=32,h=\frac{1}{76}$ in Figure 1, where the numerical solution of u is on the left and numerical solution of q is on the right with $\alpha =0.5$.

In next example, we will use the scheme for the FOCP (1) and (2) for the two dimensional in space. Denote $x=(x_1,x_2)\in R^2$ in the following example.

Example 2.

Considering the problem (1) and (2) with 2D problem and the desired state $u_d(x,t)$ and the right function $f(x,t)$ be defined as follows

$$\begin{array}{ccc}\hfill u_d(x,t)& =& t^{4}sin\left(2\pi x_1\right)sin\left(2\pi x_2\right)-\frac{24{(1-t)}^{4-\alpha }}{\mathsf{\Gamma }(5-\alpha )}sin\left(2\pi x_1\right)sin\left(2\pi x_2\right)\hfill \\ & & -8{\pi }^2{(1-t)}^{4}sin\left(2\pi x_1\right)sin\left(2\pi x_2\right),\hfill \\ \hfill f(x,t)& =& \frac{24t^{4-\alpha }}{\mathsf{\Gamma }(5-\alpha )}sin\left(2\pi x_1\right)sin\left(2\pi x_2\right)+{(1-t)}^{4}sin\left(2\pi x_1\right)sin\left(2\pi x_2\right)\hfill \\ & & +t^{4}8{\pi }^{2}sin\left(2\pi x_1\right)sin\left(2\pi x_2\right).\hfill \end{array}$$

By direct calculation, we obtain exact analytical solution state variable and control variable:

$$\begin{array}{c}\hfill u(x,t)=t^{4}sin\left(2\pi x_1\right)sin\left(2\pi x_2\right),q\left(t\right)=-{(1-t)}^{4}sin\left(2\pi x_1\right)sin\left(2\pi x_2\right).\end{array}$$

Let the error of u be defined as follows: ${\displaystyle e_u=\underset{i}{max}\left|\right|u_h^i\left(x\right)-u(x,t_i){\left|\right|}_{L^2\left(\mathsf{\Omega }\right)}}$, and the error of q be ${\displaystyle e_q=\underset{i}{max}\left|\right|q_h^i\left(x\right)-q(x,t_i){\left|\right|}_{L^2\left(\mathsf{\Omega }\right)}.}$

In

Table 7. Error $e_q$ and $e_u$ for the spatial convergence rates with $\alpha =0.3$.

h	${\mathit{e}}_{\mathit{q}}$	Rate	${\mathit{e}}_{\mathit{u}}$	Rate
$\frac{1}{8}$	1.11534867$\times {10}^{-2}$	-	1.11532795$\times {10}^{-2}$	-
$\frac{1}{16}$	2.96296656$\times {10}^{-3}$	1.91238053	2.96291062$\times {10}^{-3}$	1.91238097
$\frac{1}{32}$	7.63819538$\times {10}^{-4}$	1.95573861	7.63805481$\times {10}^{-4}$	1.95573793
$\frac{1}{64}$	1.93916779$\times {10}^{-4}$	1.97779417	1.93913963$\times {10}^{-4}$	1.97778857
$\frac{1}{128}$	4.88530203$\times {10}^{-5}$	1.98891798	4.88531596$\times {10}^{-5}$	1.98889291

,

Table 8. Error $e_q$ and $e_u$ for the spatial convergence rates with $\alpha =0.5$.

h	${\mathit{e}}_{\mathit{q}}$	Rate	${\mathit{e}}_{\mathit{u}}$	Rate
$\frac{1}{8}$	1.11068234$\times {10}^{-2}$	-	1.11065513$\times {10}^{-2}$	-
$\frac{1}{16}$	2.94892688$\times {10}^{-3}$	1.91318430	2.94885444$\times {10}^{-3}$	1.91318440
$\frac{1}{32}$	7.59885002$\times {10}^{-4}$	1.95633704	7.59867696$\times {10}^{-4}$	1.95633445
$\frac{1}{64}$	1.92856611$\times {10}^{-4}$	1.97825249	1.92854388$\times {10}^{-4}$	1.97823627
$\frac{1}{128}$	4.85725519$\times {10}^{-5}$	1.98931541	4.85748893$\times {10}^{-5}$	1.98922935

and

Table 9. Error $e_q$ and $e_u$ for the spatial convergence rates with $\alpha =0.7$.

h	${\mathit{e}}_{\mathit{q}}$	Rate	${\mathit{e}}_{\mathit{u}}$	Rate
$\frac{1}{8}$	7.17426167$\times {10}^{-1}$	-	3.98340005$\times {10}^{-2}$	-
$\frac{1}{16}$	4.90137294$\times {10}^{-1}$	5.49644442	1.12084191$\times {10}^{-2}$	1.82941756
$\frac{1}{32}$	2.97665744$\times {10}^{-3}$	7.36334890	2.97659460$\times {10}^{-3}$	1.91284816
$\frac{1}{64}$	1.94692375$\times {10}^{-4}$	1.95640433	7.66980380$\times {10}^{-4}$	1.95640117
$\frac{1}{128}$	1.94692375$\times {10}^{-4}$	1.97802058	1.94691821$\times {10}^{-4}$	1.97799739

, we take $\tau =h/2;h=2^{-l},l=4,5,6,7,8,9$, $\alpha =0.3,0.5,0.7$, respectively. From Table 7, Table 8 and Table 9, we find the spatial accuracy is 2, this result is accord with the theoretical analysis obtained in Theorem 3.

Next, we will study the convergence order of time. In the following, we choose $h=2^{-l},$$l=3,4,5,6,7.$ and $\tau =O\left(h^{\frac{2}{3-\alpha }}\right)$. In

Table 10. Error $e_q$ and $e_u$ for the time convergence rates with $\alpha =0.3$.

h	${\mathit{e}}_{\mathit{q}}$	Rate	${\mathit{e}}_{\mathit{u}}$	Rate
$\frac{1}{8}$	4.65123837$\times {10}^{-1}$	-	3.96008205$\times {10}^{-2}$	-
$\frac{1}{16}$	1.11534867$\times {10}^{-2}$	5.38204819	1.11532795$\times {10}^{-2}$	1.82806233
$\frac{1}{32}$	2.96429450$\times {10}^{-3}$	1.91173409	2.96423899$\times {10}^{-3}$	1.91173431
$\frac{1}{64}$	7.64321219$\times {10}^{-4}$	1.95543779	7.64307151$\times {10}^{-4}$	1.95543733
$\frac{1}{128}$	1.94087839$\times {10}^{-4}$	1.97746935	1.94084832$\times {10}^{-4}$	1.97746515

,

Table 11. Error $e_q$ and $e_u$ for the time convergence rates with $\alpha =0.5$.

h	${\mathit{e}}_{\mathit{q}}$	Rate	${\mathit{e}}_{\mathit{u}}$	Rate
$\frac{1}{8}$	4.65214791$\times {10}^{-1}$	-	3.94176400$\times {10}^{-2}$	-
$\frac{1}{16}$	1.10979951$\times {10}^{-2}$	5.38952598	1.10977217$\times {10}^{-2}$	1.82857786
$\frac{1}{32}$	2.94892688$\times {10}^{-3}$	1.91203712	2.94885444$\times {10}^{-3}$	1.91203702
$\frac{1}{64}$	7.60347877$\times {10}^{-4}$	1.95545851	7.60330462$\times {10}^{-4}$	1.95545611
$\frac{1}{128}$	1.93057130$\times {10}^{-4}$	1.97763179	1.93054617$\times {10}^{-4}$	1.97761752

and

Table 12. Error $e_q$ and $e_u$ for the time convergence rates with $\alpha =0.7$.

h	${\mathit{e}}_{\mathit{q}}$	Rate	${\mathit{e}}_{\mathit{u}}$	Rate
$\frac{1}{8}$	4.65337865$\times {10}^{-1}$	-	3.91653829$\times {10}^{-2}$	-
$\frac{1}{16}$	1.10285894$\times {10}^{-2}$	5.39895839	1.10283434$\times {10}^{-2}$	1.82836296
$\frac{1}{32}$	2.93078394$\times {10}^{-3}$	1.91188976	2.93072027$\times {10}^{-3}$	1.91188891
$\frac{1}{64}$	7.55565458$\times {10}^{-4}$	1.95565796	7.55552247$\times {10}^{-4}$	1.95565184
$\frac{1}{128}$	1.91829133$\times {10}^{-4}$	1.97773490	1.91829807$\times {10}^{-4}$	1.97770461

, the $\alpha$ choose the value of $0.3,0.5$, and $0.7$, respectively. From Table 10, Table 11 and Table 12, it is easy to check that the time convergence order is $3-\alpha$, based on the fact that the rate of convergence is close to 2 under the condition $\tau =O\left(h^{\frac{2}{3-\alpha }}\right)$.

In the end of Example 2, we show the CPU time variation about $M=h^{-1}$ and N in Figure 2 under the conduction of $N=2^8$ and $M=h^{-1}=2^8$, respectively. From Figure 2, we obtain that the CPU time increases with the increase of $M=h^{-1}$ or N and almost not affected by the change of $\alpha$.

6. Conclusions

In this paper, we constructed a novel FD-FE scheme for the time FOCP based on the uniform accuracy $(3-\alpha )$ order FD scheme in time and FE scheme in space. We firstly give the stability analysis for the full discrete scheme for the time fractional partial differential equation based on the uniform accuracy $(3-\alpha )$ order FD scheme in time and FE scheme in space. The priori error estimates of the semidiscrete scheme underwent rigorous theoretical analysis. Some numerical examples are devoted to verify the correctness of the theoretical analysis. Due to the nonlocality of the time fractional derivative, the discrete high-order numerical scheme has a large amount of computation and storage, which is difficult to calculate. Especially for three-dimensional practical engineering problems, the computational efficiency of the algorithm needs to be further improved.

In our future work, we will investigate Structural optimization of viscoelastic materials and structural optimization design of viscoelastic composite plates based the idea of [45]. In the future, it is expected that the construction of the higher-order efficient scheme for time FOCP can be applied to the structural optimization design of practical engineering materials based on the ideas of [46,47,48]. In particular, we are going to use the above efficient higher-order scheme for the optimization of a composite structure of viscoelastic materials or structure with memory.