Gradient Recovery Based Two-Grid Finite Element Method for Parabolic Integro-Differential Optimal Control Problems

Miao Yang

School of Mathematics, Southwest Jiaotong University, Chengdu, China.

DOI: 10.4236/jamp.2024.128170                

Abstract

In this paper, the optimal control problem of parabolic integro-differential equations is solved by gradient recovery based two-grid finite element method. Piecewise linear functions are used to approximate state and co-state variables, and piecewise constant function is used to approximate control variables. Generally, the optimal conditions for the problem are solved iteratively until the control variable reaches error tolerance. In order to calculate all the variables individually and parallelly, we introduce a gradient recovery based two-grid method. First, we solve the small scaled optimal control problem on coarse grids. Next, we use the gradient recovery technique to recover the gradients of state and co-state variables. Finally, using the recovered variables, we solve the large scaled optimal control problem for all variables independently. Moreover, we estimate priori error for the proposed scheme, and use an example to validate the theoretical results.

Keywords

Optimal Control Problem, Gradient Recovery, Two-Grid Finite Element Method

1. Introduction

The optimal control problem [1] [2] is very crucial in the field of science and engineering. For example, both population dynamics and heat conduction involve the optimal control problem [3] [4]. Solving and analyzing optimal control problems are important matters in applying the optimal control models. So far, lots of studies have been made for the optimal control problems. For instance, a new method based on optimal control theory for megawatt frequency control problem was discussed in [5]. Superconvergence analysis and error estimation of finite element method were established for convex OCP in [6]. Furthermore, a posteriori error estimation of spectral method was presented for optimal control problems governed by parabolic equations in [7]. In addition, the convergence and superconvergence of fully discrete finite elements for time fractional optimal control problems are given in [8]. We can learn more about the OCP through [9] and references therein.

When using the FEM with two grids to solve the OCP of parabolic integro-differential equation, the processes of solving all variables are interdependent. Therefore, in this article, we try to design a scheme to solve all variables independently. Based on the algorithms studied in [10] [11], we mainly study combining the high efficiency of the two-grid finite element method and the high-precision property of the gradient recovery method to design an efficient algorithm for the optimal control problem.

Gradient recovery is an important post-processing method. On the one hand, it reconstructs high-precision gradient approximations of finite element solutions [12]. On the other hand, gradient recovery can be used to construct posteriori error estimators [13]-[16]. Omar and Tristan have studied gradient recovery in adaptive finite element method for parabolic problems [17]. As early as 2003, Yan proposed a gradient recovery type a posteriori error estimator for FEM of the OCP [18]. Recently, some scholars have studied the posteriori error estimate of two grid mixed finite element methods for semilinear elliptic equations [19]. In addition, the Galerkin method for numerical solution of Volterra integro-differential equations with certain orthogonal basis function is presented in [20]. To the best of our knowledge, there is no literature about gradient recovery based two-grid finite element method for parabolic integro-differential optimal control problems. In this paper, we apply two-grid finite element method based gradient recovery method studied in [21] to study a parallel algorithm for (1)-(2). Details of the algorithm mainly include three parts. First, solving the small scaled optimal control problem for state, co-state and control variables $y_H$ , $p_H$ and $u_H$ respectively on coarse mesh with mesh size H at each time level. Next, applying the gradient recovery method to get $R_H{y}_H$ and $R_H{p}_H$ . With $R_H{y}_H$ and $R_H{p}_H$ , we solve the large scaled optimal control problem for $y_h$ , $p_h$ and $u_h$ on fine mesh with mesh size h in parallel on all time levels. Error in the approximate solution of the proposed algorithm is estimated and a numerical experiment is implemented to confirm the theoretical results. Here we provide a flowchart to illustrate the methodological approach utilized in this paper.

The framework of this article is organized as below. We introduce the discrete form of the modelling issue in Section 2. In addition, fully discrete form with intermediate variables is also given. Priori error estimate is presented in Section 3. Gradient recovery based two-grid finite element method as well as error estimation of the method are presented in Section 4. A numerical experiment is used to justify the theoretical result in Section 5. Finally, a simple summary is given in Section 6.

2. Model Problem and Its Finite Element Scheme

In this article, we study the below OCP:

$\underset{u\in K\subset U_{ad}}{\text{min}}\left\{\frac{1}{2}{\displaystyle {\int }_0^T{\displaystyle {\int }_{\Omega }{\left|y-z_d\right|}^2\text{d}x\text{d}t}}+\frac{\gamma }{2}{\displaystyle {\int }_0^T{\displaystyle {\int }_{{\Omega }_U}{\left|u\right|}^2\text{d}x\text{d}t}}\right\}$ (1)

subject to

$\{\begin{array}{l}\frac{\partial y}{\partial t}-A\text{Δ}{y}_t-D\text{Δ}y-{\displaystyle {\int }_0^t\mathcal{C}\left(t,s\right)\text{Δ}y\left(s\right)\text{d}s}=f+Eu\text{in}\Omega \times \left(0,T\right],\hfill \\ y=0\text{on}\partial \Omega \times \left[0,T\right],\hfill \\ {y|}_{t=0}=y^0\text{in}\Omega ,\hfill \end{array}$ (2)

where $\Omega$ , with Lipschitz boundary $\partial \Omega$ , and ${\Omega }_U$ be bounded open sets in $R^d$ , $1\le d\le 3$ . $\gamma$ is a positive regular constant. y is the state, u is the control, f, $z_d$ and $y^0$ are some given functions, $U_{ad}$ denotes a closed convex subset that includes control, where

$A={\left(a_{ij}\left(x\right)\right)}_{d\times d},D={\left(d_{ij}\left(x\right)\right)}_{d\times d},\mathcal{C}={\left(c_{ij}\left(x,t,s\right)\right)}_{d\times d},$

E is a bounded operator from $L^2\left(0,T;L^2\left({\Omega }_U\right)\right)$ to $L^2\left(0,T;L^2\left(\Omega \right)\right)$ and is independent of t.

$K=\left\{u\in U_{ad}:{\displaystyle {\int }_{\Omega }u\left(x,t\right)\text{d}x}\ge 0\right\}.$ (3)

First, we shall introduce some notations that used in this paper. Let $W^{m,p}\left(\Omega \right)$ be the Sobolev spaces on $\Omega$ and a norm ${\Vert \cdot \Vert }_{m,p}$ denoted by ${\Vert v\Vert }_{m,p}^p={\displaystyle \sum }_{\left|\alpha \right|\le m}{\Vert D^{\alpha }v\Vert }_{L^p\left(\Omega \right)}^p$ , D is differential operator, a semi-norm ${\left|\cdot \right|}_{m,p}$ denoted by ${\left|v\right|}_{m,p}^p={\displaystyle \sum }_{\left|\alpha \right|=m}{\Vert D^{\alpha }v\Vert }_{L^p\left(\Omega \right)}^p$ . We set $W_0^{m,p}\left(\Omega \right)=\left\{v\in W^{m,p}\left(\Omega \right):{v|}_{\partial \Omega }=0\right\}$ . For $p=2$ , we write as $H^m\left(\Omega \right)=W^{m,2}\left(\Omega \right)$ , $H_0^m\left(\Omega \right)=W_0^{m,2}\left(\Omega \right)$ , and ${\Vert \cdot \Vert }_m={\Vert \cdot \Vert }_{m,2}$ ,$\Vert \cdot \Vert ={\Vert \cdot \Vert }_{0,2}$ .

We denote by $L^s\left(J;W^{m,p}\left(\Omega \right)\right)$ the Banach space of all $L^s$ integrable functions from J into $W^{m,p}\left(\Omega \right)$ with norm ${\Vert v\Vert }_{L^s\left(J;W^{m,p}\left(\Omega \right)\right)}={\left({\displaystyle {\int }_0^T{\Vert v\Vert }_{W^{m,p}\left(\Omega \right)}^s\text{d}t}\right)}^{\frac{1}{s}}$ for $s\in \left[1,\infty \right)$ , $J=\left(0,T\right]$ . For simplicity of presentation, we write ${\Vert v\Vert }_{L^s\left(J;W^{m,p}\left(\Omega \right)\right)}$ by ${\Vert v\Vert }_{L^s\left(W^{m,p}\right)}$ . In the same way, we are able to define the spaces $H^1\left(J;W^{m,p}\left(\Omega \right)\right)$ and $C^k\left(J;W^{m,p}\left(\Omega \right)\right)$ . Besides, h is the spatial mesh size, C represents an arbitrary positive constant, $\text{Δ}t$ is time step, the state space $Q=L^2\left(J;V\right)$ and $V=H_0^1\left(\Omega \right)$ .

The weak form of (1)-(2) is to find $\left(y,u\right)\in Q\times K$ satisfies

$\underset{u\in K\subset U_{ad}}{\text{min}}\left\{\frac{1}{2}{\displaystyle {\int }_0^T{\displaystyle {\int }_{\Omega }{\left|y-z_d\right|}^2}}+\frac{\gamma }{2}{\displaystyle {\int }_0^T{\displaystyle {\int }_{{\Omega }_U}{\left|u\right|}^2}}\right\}$ (4)

$\begin{array}{l}\left(y_t,v\right)+\left(A\nabla y_t,\nabla v\right)+\left(D\nabla y,\nabla v\right)+{\displaystyle {\int }_0^t\left(\mathcal{C}\left(t,s\right)\nabla y\left(s\right),\nabla v\right)\text{d}s}\\ =\left(f+Eu,v\right),\forall v\in V,t\in J,\end{array}$ (5)

${y|}_{t=0}=y^0,$ (6)

From [22] and [23], we know that (4)-(6) has a unique solution $\left(y,u\right)$ only when there is a co-state $p\in Q$ such that $\left(y,p,u\right)$ meets the following optimality conditions:

$\begin{array}{l}\left(y_t,v\right)+\left(A\nabla y_t,\nabla v\right)+\left(D\nabla y,\nabla v\right)+{\displaystyle {\int }_0^t\left(\mathcal{C}\left(t,s\right)\nabla y\left(s\right),\nabla v\right)\text{d}s}\\ =\left(f+Eu,v\right),\forall v\in V,t\in J,\end{array}$ (7)

${y|}_{t=0}=y^0,$ (8)

$\begin{array}{l}-\left(p_t,q\right)-\left(A\nabla p_t,\nabla q\right)+\left(D\nabla p,\nabla q\right)+{\displaystyle {\int }_t^T\left({\mathcal{C}}^{*}\left(s,t\right)\nabla p\left(s\right),\nabla q\right)\text{d}s}\\ =\left(y-z_d,q\right),\forall q\in V,t\in J,\end{array}$ (9)

${p|}_{t=T}=0,$ (10)

$\left(\gamma u+E^{*}p,v-u\right)\ge 0,\forall v\in K,t\in J,$ (11)

where $E^{\text{*}}$ is the dual operator of E, ${\mathcal{C}}^{\text{*}}$ is the dual operator of the operator $\mathcal{C}$ . Refer to [24], inequality (11) can be expressed as

$u=-\frac{1}{\gamma }{E}^{*}p+\max \left\{0,\frac{1}{\gamma \left|{\Omega }_U\right|}{\displaystyle {\int }_{{\Omega }_U}{E}^{*}p}\right\}.$ (12)

Suppose ${\mathcal{T}}_h$ indicates a partitioning of $\Omega$ into disjoint regular triangulations $\tau$ , $h_{\tau }$ be the size of $\tau$ , $h=\underset{\tau \in {\mathcal{T}}_h}{\max }{h}_{\tau }$ . Suppose $V_h\subset V$ be defined as the following finite element space:

$V_h=\left\{v_h\in C^0\left(\overline{\Omega }\right)\cap V|{v_h|}_{\tau }\in P_1\left(\tau \right),\forall \tau \in {\mathcal{T}}_h\right\}.$(13)

The approximate space of control is written as

$U_h:=\left\{{\tilde{u}}_h\in K:\forall \tau \in {\mathcal{T}}_h,{{\tilde{u}}_h|}_{\tau }=constant\right\},$(14)

Prior to giving the fully discrete finite element equation, let’s start introducing several projection operators. The standard elliptic projection [25] $R_h:V\to V_h$ is defined as follows:

$\left(\nabla \left(\varphi -R_h\varphi \right),\nabla v_h\right)=0,\forall v_h\in V_h,\forall \varphi \in V,$ (15)

${\Vert \varphi -R_h\varphi \Vert }_s\le C{h}^{2-s}{\Vert \varphi \Vert }_2,s=0,1,\forall \varphi \in H^s\left(\Omega \right).$ (16)

The standard $L^2$ -orthogonal projection [26] $Q_h:K\to U_h$ is defined as:

$\left(u-Q_{h}u,{\tilde{u}}_h\right)=0,\forall {\tilde{u}}_h\in U_h,$(17)

${\Vert u-Q_{h}u\Vert }_{-s,r}\le C{h}^{1+s}{\Vert u\Vert }_{1,r},s=0,1,\forall u\in W^{1,r}\left(\Omega \right).$ (18)

The definition of the element average operator [27] ${\pi }_h:L^2\left(\Omega \right)\to U_h$ is given by

(19)

which has an approximation property:

${\Vert \psi -{\pi }_h\psi \Vert }_{-s,r}\le C{h}^{1+s}{\Vert \psi \Vert }_{1,r},s=0,1,\forall \psi \in W^{1,r}\left(\Omega \right).$ (20)

For the theoretical analysis in Section 3 and Section 4, we need to present a fully discrete form of the model problem and a fully discrete form with intermediate variable $\tilde{u}$ . Now, we present the fully discrete finite element approximation for (4)-(5). Suppose $\text{Δ}t>0$ , $N=T/\text{Δ}t\in Z$ , $t_n=n\text{Δ}t$ , $n\in Z$ , and we propose the following concepts:

${\psi }^n={\psi }^n\left(x\right)=\psi \left(x,t_n\right),\text{d}t{\psi }^n=\frac{{\psi }^n-{\psi }^{n-1}}{\text{Δ}t},\delta {\psi }^n={\psi }^n-{\psi }^{n-1}.$

The fully discrete finite element approximation is to find: $\left(y_h^n,u_h^n\right)\in V_h\times K_h$ ,$n=1,2,\cdots ,N$ satisfy

$\underset{u_h^n\in K_h}{\text{min}}\left\{\frac{1}{2}{\displaystyle \sum _{n=1}^N\text{Δ}t{\displaystyle {\int }_{\Omega }{\left|y_h^n-z_d^n\right|}^2}}+\frac{\gamma }{2}\text{Δ}t{\displaystyle \sum _{n=1}^N{\displaystyle {\int }_{{\Omega }_U}{\left|u_h^n\right|}^2}}\right\}$ (21)

$\begin{array}{l}\left(\text{d}t{y}_h^n,v_h\right)+\left(A\nabla \text{d}t{y}_h^n,\nabla v_h\right)+\left(D\nabla y_h^n,\nabla v_h\right)+{\displaystyle \sum }_{i=1}^n\text{Δ}t\left(\mathcal{C}\left(t_n,t_{i-1}\right)\nabla y_h^i,\nabla v_h\right)\\ =\left(f^n+E{u}_h^n,v_h\right),\forall v_h\in V_h,\end{array}$ (22)

$y_h^0=R_h{y}^0.$ (23)

where we change the integral term to the summation term and carefully set the subscript of the variables to facilitate the calculation and proof of the theoretical analysis part.

Once again, combining the optimality condition, duplet $\left(y_h^n,u_h^n\right)\in V_h\times K_h$ is the solution of (21)-(23) only when there is a co-state $p_h^{n-1}\in V_h$ such that $\left(y_h^n,p_h^{n-1},u_h^n\right)$ satisfies the following conditions:

$\begin{array}{l}\left(\text{d}t{y}_h^n,v_h\right)+\left(A\nabla \text{d}t{y}_h^n,\nabla v_h\right)+\left(D\nabla y_h^n,\nabla v_h\right)+{\displaystyle \sum }_{i=1}^n\text{Δ}t\left(\mathcal{C}\left(t_n,t_{i-1}\right)\nabla y_h^i,\nabla v_h\right)\\ =\left(f^n+E{u}_h^n,v_h\right),\forall v_h\in V_h,\end{array}$ (24)

$y_h^0=R_h{y}^0,$ (25)

$\begin{array}{l}-\left(\text{d}t{p}_h^n,q_h\right)-\left(A\nabla \text{d}t{p}_h^n,\nabla q_h\right)+\left(D\nabla p_h^{n-1},\nabla q_h\right)\\ +{\displaystyle \sum }_{i=n}^N\text{Δ}t\left({\mathcal{C}}^{*}\left(t_i,t_{n-1}\right)\nabla p_h^{i-1},\nabla q_h\right)=\left(y_h^n-z_d^n,q_h\right),\forall q_h\in V_h,\end{array}$ (26)

$p_h^N=0,$ (27)

$u_h^n=-\frac{1}{\gamma }{E}^{*}{\pi }_h{p}_h^{n-1}+\max \left\{0,\frac{1}{\gamma \left|{\Omega }_U\right|}{\displaystyle {\int }_{{\Omega }_U}{E}^{*}{p}_h^{n-1}}\right\}.$ (28)

Next, we introduce the fully discrete form with intermediate variable first. For $\forall \tilde{u}\in K$, the following is the discrete form with $\tilde{u}$ :

$\begin{array}{l}\left(\text{d}t{y}_h^n\left(\tilde{u}\right),v_h\right)+\left(A\nabla \text{d}t{y}_h^n\left(\tilde{u}\right),\nabla v_h\right)+\left(D\nabla y_h^n\left(\tilde{u}\right),\nabla v_h\right)\\ +{\displaystyle \sum }_{i=1}^n\text{Δ}t\left(\mathcal{C}\left(t_n,t_{i-1}\right)\nabla y_h^i\left(\tilde{u}\right),\nabla v_h\right)=\left(f^n+E{\tilde{u}}^n,v_h\right),\forall v_h\in V_h,\end{array}$(29)

$y_h^0\left(\tilde{u}\right)=R_h{y}^0,$ (30)

$\begin{array}{l}-\left(\text{d}t{p}_h^n\left(\tilde{u}\right),q_h\right)-\left(A\nabla \text{d}t{p}_h^n\left(\tilde{u}\right),\nabla q_h\right)+\left(D\nabla p_h^{n-1}\left(\tilde{u}\right),\nabla q_h\right)\\ +{\displaystyle \sum }_{i=n}^N\text{Δ}t\left({\mathcal{C}}^{*}\left(t_i,t_{n-1}\right)\nabla p_h^{i-1}\left(\tilde{u}\right),\nabla q_h\right)=\left(y_h^n\left(\tilde{u}\right)-z_d^n,q_h\right),\forall q_h\in V_h,\end{array}$(31)

$p_h^N\left(\tilde{u}\right)=0.$ (32)

3. A priori Error Estimates

Refer to [28], we have the following Lemmas 1-6.

Lemma 1. Suppose $\left(y_h^n\left(u\right),p_h^{n-1}\left(u\right)\right)$ satisfies (29)-(32) with $\tilde{u}=u$ and $\left(y,p\right)$ satisfies (7)-(11). For $1\le n\le N$ , we get

$\Vert y^n-y_h^n\left(u\right)\Vert +\Vert p^{n-1}-p_h^{n-1}\left(u\right)\Vert \le C\left(\text{Δ}t+h^2\right),$ (33)

$\Vert \nabla \left(y^n-y_h^n\left(u\right)\right)\Vert +\Vert \nabla \left(p^{n-1}-p_h^{n-1}\left(u\right)\right)\Vert \le C\left(\text{Δ}t+h\right).$ (34)

Lemma 2. Choosing ${\tilde{u}}^n=Q_h{u}^n$ and ${\tilde{u}}^n=u^n$ in (29)-(32) respectively. Then, for $1\le n\le N$ , we get

$\Vert \nabla \left(y_h^n\left(u\right)-y_h^n\left(Q_{h}u\right)\right)\Vert +\Vert \nabla \left(p_h^{n-1}\left(u\right)-p_h^{n-1}\left(Q_{h}u\right)\right)\Vert \le C{h}^2.$ (35)

Lemma 3. Suppose $\left(y_h^n,p_h^n\right)$ and $\left(y_h^n\left(Q_{h}u\right),p_h^n\left(Q_{h}u\right)\right)$ are the discrete solutions of (29)-(32) with ${\tilde{u}}^n=u_h^n$ , and ${\tilde{u}}^n=Q_h{u}^n$ respectively. For $1\le n\le N$ , we have

$\begin{array}{l}\Vert \nabla \left(y_h^n\left(Q_{h}u\right)-y_h^n\right)\Vert +\Vert \nabla \left(p_h^{n-1}\left(Q_{h}u\right)-p_h^{n-1}\right)\Vert \\ \le C{\left({\displaystyle \sum }_{n=1}^N{\Vert Q_h{u}^n-u_h^n\Vert }^2\text{Δ}t\right)}^{\frac{1}{2}}.\end{array}$ (36)

Lemma 4. Substituting ${\tilde{u}}^n=Q_h{u}^n$ , ${\tilde{u}}^n=u_h^n$ into (29)-(32) respectively, thus

${\displaystyle \sum }_{n=1}^N\left(Q_h{u}^n-u_h^n,p_h^{n-1}\left(Q_{h}u\right)-p_h^{n-1}\right)\text{Δ}t\ge 0.$ (37)

Lemma 5. Suppose u satisfies (7)-(11) and $u_h^n$ satisfies (24)-(28) and all the assumptions in Lemmas 1-4 are reasonable. Thus, for $1\le n\le N$ , we have

${\left({\displaystyle \sum }_{n=1}^N{\Vert Q_h{u}^n-u_h^n\Vert }^2\text{Δ}t\right)}^{\frac{1}{2}}\le C\left(h^2+\text{Δ}t\right).$ (38)

In order to prove the global superconvergence for all variables, we utilize the recovery technique on uniform meshes. We set up the recovery operator $P_h$ and $G_h$ , suppose $P_{h}v$ is a continuous quadratic function. The value of $P_{h}v$ at the nodes is defined on a patch of elements around the node by the least squares method, the details can be found in [29] [30].

For the gradients of y and p, $G_{h}v=\left(P_h{v}_x,P_h{v}_y\right)$ . For the quadratic function, which is identical to the Z-Z gradient recovery [29] [30]. We set the discrete co-state for an acceptable set

${\widehat{u}}_h^n=-\frac{1}{\gamma }{E}^{*}{p}_h^{n-1}+\max \left\{0,\frac{1}{\gamma \left|{\Omega }_U\right|}{\displaystyle {\int }_{{\Omega }_U}{E}^{*}{p}_h^{n-1}}\right\}.$ (39)

Lemma 6. Suppose $\left(y,p\right)$ and $\left(y_h^n,p_h^{n-1}\right)$ satisfies (7)-(11) and (29)-(32) respectively and all the hypotheses in Lemmas 1-5 are valid. Thus

$\Vert G_h{y}_h^n-\nabla y^n\Vert +\Vert G_h{p}_h^{n-1}-\nabla p^{n-1}\Vert \le C\left(h^2+\text{Δ}t\right).$ (40)

Next, we will deduce the result of global superconvergence for the control variable and the state variable.

Theorem 7. Suppose $\left(y,p,u\right)$ and $\left(y_h^n,p_h^{n-1},u_h^n\right)$ makes (7)-(11) and (24)-(28) valid respectively. For $1\le n\le N$ , we can derive that

$\Vert y^n-y_h^n\Vert +\Vert p^{n-1}-p_h^{n-1}\Vert \le C\left(h^2+\text{Δ}t\right),$ (41)

$\Vert \nabla \left(y^n-y_h^n\right)\Vert +\Vert \nabla \left(p^{n-1}-p_h^{n-1}\right)\Vert \le C\left(h+\text{Δ}t\right),$ (42)

$\Vert u^n-u_h^n\Vert \le C\left(h+\text{Δ}t\right).$ (43)

Proof. Apparently, using Poincaré inequality and the Lemma 1-Lemma 5, we have (41)-(42). From (12) and (28), combining (20), (41) and mean value theorem, we get

$\begin{array}{c}\Vert u^n-u_h^n\Vert \le \Vert \frac{1}{\gamma }{E}^{*}{p}^n-\frac{1}{\gamma }{E}^{*}{\pi }_h{p}_h^{n-1}\Vert +\max \left\{0,\frac{1}{\gamma \left|{\Omega }_U\right|}{\displaystyle {\int }_{{\Omega }_U}{E}^{*}{p}^n}\right\}\\ -\max \left\{0,\frac{1}{\gamma \left|{\Omega }_U\right|}{\displaystyle {\int }_{{\Omega }_U}{E}^{*}{p}_h^{n-1}}\right\}\\ \le C\Vert p^n-{\pi }_h{p}_h^{n-1}\Vert +C\Vert p^n-p_h^{n-1}\Vert \\ \le C\Vert p^n-p^{n-1}\Vert +C\Vert p^{n-1}-{\pi }_h{p}^{n-1}\Vert \\ +C\Vert {\pi }_h{p}^{n-1}-{\pi }_h{p}_h^{n-1}\Vert +C\Vert p^n-p_h^{n-1}\Vert \end{array}$

$\begin{array}{c}\le C\Vert p^n-p^{n-1}\Vert +C\Vert p^{n-1}-{\pi }_h{p}^{n-1}\Vert +C\Vert p^{n-1}-p_h^{n-1}\Vert \\ \le C\left(h+\Delta t\right).\end{array}$ (44)

□

Theorem 8. Suppose u and $u_h^n$ satisfy (7)-(11) and (29)-(32) respectively. Assuming that all the conditions in Lemmas 1-5 are valid. Thus

$\Vert u^n-{\widehat{u}}_h^n\Vert \le C\left(h^2+\text{Δ}t\right).$(45)

Proof. The proof of this theorem is similar to the proof of Theorem 7. □

4. Recovery Based Two-Grid Scheme

In this section, we propose a gradient recovery based two-grid finite element method and make a priori error estimate for the algorithm. The main idea of the scheme includes two parts corresponding to coarse and fine mesh respectively.

Step 1. Solving $\left(y_H^n,p_H^{n-1},u_H^n\right)\in V_H\times V_H\times K_H$ on the coarse grid ${\mathcal{T}}_H$ to satisfy the following optimality conditions:

$\begin{array}{l}\left(\text{d}t{y}_H^n,v_H\right)+\left(A\text{d}t\nabla y_H^n,\nabla v_H\right)+\left(D\nabla y_H^n,\nabla v_H\right)\\ +{\displaystyle \sum }_{i=1}^n\text{Δ}t\left(\mathcal{C}\left(t_n,t_{i-1}\right)\nabla y_H^i,\nabla v_H\right)=\left(f^n+E{u}_H^n,v_H\right),\forall v_H\in V_H,\end{array}$ (46)

$y_H^0=R_H{y}^0,$ (47)

$\begin{array}{l}-\left(\text{d}t{p}_H^n,q_H\right)-\left(A\text{d}t\nabla p_H^n,\nabla q_H\right)+\left(D\nabla p_H^{n-1},\nabla q_H\right)\\ +{\displaystyle \sum }_{i=n}^N\text{Δ}t\left({\mathcal{C}}^{*}\left(t_i,t_{n-1}\right)\nabla p_H^{i-1},\nabla q_H\right)=\left(y_H^n-z_d^n,q_H\right),\forall q_H\in V_H,\end{array}$ (48)

$p_H^N=0,$ (49)

${\left(\gamma u_H^n+E^{*}{p}_H^{n-1},v_H^{*}-u_H^n\right)}_U\ge 0,\forall v_H^{*}\in K_H.$ (50)

Step 2. Finding $\left({\overline{\tilde{y}}}_h^n,{\overline{\tilde{p}}}_h^{n-1},{\overline{\tilde{u}}}_h^n\right)\in V_h\times V_h\times K_h$ on the fine grid ${\mathcal{T}}_h$ to satisfy:

$\begin{array}{l}\left(\text{d}t{\overline{\tilde{y}}}_h^n,v_h\right)+\left(A\text{d}t\nabla {\overline{\tilde{y}}}_h^n,\nabla v_h\right)+\left(D\nabla {\overline{\tilde{y}}}_h^n,\nabla v_h\right)\\ +{\displaystyle \sum }_{i=1}^n\text{Δ}t\left(\mathcal{C}\left(t_n,t_{i-1}\right)G_H{y}_H^i,\nabla v_h\right)=\left(f^n+E{\widehat{u}}_H^n,v_h\right),\forall v_h\in V_h,\end{array}$(51)

${\overline{\tilde{y}}}_h^0=R_h{y}^0,$ (52)

$\begin{array}{l}-\left(\text{d}t{\overline{\tilde{p}}}_h^n,q_h\right)-\left(A\text{d}t\nabla {\overline{\tilde{p}}}_h^n,\nabla q_h\right)+\left(D\nabla {\overline{\tilde{p}}}_h^{n-1},\nabla q_h\right)\\ +{\displaystyle \sum }_{i=n}^N\text{Δ}t\left({\mathcal{C}}^{*}\left(t_i,t_{n-1}\right)G_H{p}_H^{i-1},\nabla q_h\right)=\left(y_H^n-z_d^n,q_h\right),\forall q_h\in V_h,\end{array}$(53)

${\overline{\tilde{p}}}_h^N=0,$ (54)

${\left(\gamma {\overline{\tilde{u}}}_h^n+E^{*}{p}_H^{n-1},v_h^{*}-{\overline{\tilde{u}}}_h^n\right)}_U\ge 0,\forall v_h^{*}\in K_h.$(55)

Combining the stability estimation and Theorem 7, it’s easy to get the following conclusion.

Theorem 9. Assuming that $\left(y,p,u\right)$ and $\left({\overline{\tilde{y}}}_h^n,{\overline{\tilde{p}}}_h^n,{\overline{\tilde{u}}}_h^n\right)$ make (7)-(11) and (46)-(55) correct, respectively. For $\Delta t$ sufficiently small and $1\le n\le N$ , we obtain

$\Vert \nabla \left(y^n-{\overline{\tilde{y}}}_h^n\right)\Vert +\Vert \nabla \left(p^{n-1}-{\overline{\tilde{p}}}_h^{n-1}\right)\Vert \le C\left(h+H^2+\text{Δ}t\right),$ (56)

$\Vert u^n-{\overline{\tilde{u}}}_h^n\Vert \le C\left(h+H^2+\text{Δ}t\right).$(57)

Proof. To reduce calculation, let

$y^n-{\overline{\tilde{y}}}_h^n=y^n-R_h{y}^n+R_h{y}^n-{\overline{\tilde{y}}}_h^n=:{\rho }^n+e^n,$

$p^n-{\overline{\tilde{p}}}_h^n=p^n-R_h{p}^n+R_h{p}^n-{\overline{\tilde{p}}}_h^n=:{\beta }^n+{\alpha }^n.$

Choosing $t=t_n$ in (7), subtracting (51) from (7), then using (15), thus

$\begin{array}{l}\left(\text{d}t{e}^n,v_h\right)+\left(D\nabla e^n,\nabla v_h\right)+\left(A\nabla \text{d}t{e}^n,\nabla \text{d}t{e}^n\right)\\ =-\left({\displaystyle {\int }_0^{t_n}\mathcal{C}\left(t_n,s\right)\nabla y\left(s\right)\text{d}s}-{\displaystyle \sum }_{i=1}^n\text{Δ}t\mathcal{C}\left(t_n,t_{i-1}\right)G_H{y}_H^i,\nabla v_h\right)\\ +\left(\text{d}t{y}^n-y_t^n,v_h\right)-\left(\text{d}t{\rho }^n,v_h\right)+\left(E\left(u^n-{\widehat{u}}_H^n\right),v_h\right)\\ +\left(A\nabla \text{d}t{y}^n-A\nabla y_t^n,\text{d}t\nabla e^n\right).\end{array}$(58)

Substituting $v_h=\text{d}t{e}^n$ into (58), multiplying it by $\text{Δ}t$ , then summing it over n from 1 to l ($1\le l\le N$ ) at the both sides of (58), therefore

$\begin{array}{l}\frac{1}{2}{\Vert D^{\frac{1}{2}}\nabla e^l\Vert }^2+{\displaystyle \sum }_{n=1}^l{\Vert \text{d}t{e}^n\Vert }^2\text{Δ}t+{\displaystyle \sum }_{n=1}^l{\Vert A^{\frac{1}{2}}\text{d}t\nabla e^n\Vert }^2\text{Δ}t\\ \le -{\displaystyle \sum }_{n=1}^l\left(\text{d}t{\rho }^n,\text{d}t{e}^n\right)\text{Δ}t+{\displaystyle \sum }_{n=1}^l\left(\text{d}t{y}^n-y_t^n,\text{d}t{e}^n\right)\text{Δ}t\\ -{\displaystyle \sum }_{n=1}^l\left({\displaystyle {\int }_0^{t_n}\mathcal{C}\left(t_n,s\right)\nabla y\left(s\right)\text{d}s}-{\displaystyle \sum }_{i=1}^n\text{Δ}t\mathcal{C}\left(t_n,t_i\right)\nabla y^i,\text{d}t\nabla e^n\right)\text{Δ}t\\ -{\displaystyle \sum }_{n=1}^l{\displaystyle \sum }_{i=1}^n\left(\text{Δ}t\mathcal{C}\left(t_n,t_i\right)\nabla y^i-\text{Δ}t\mathcal{C}\left(t_n,t_{i-1}\right)\nabla y^i,\text{d}t\nabla e^n\right)\text{Δ}t\\ -{\displaystyle \sum }_{n=1}^l{\displaystyle \sum }_{i=1}^n\left(\text{Δ}t\mathcal{C}\left(t_n,t_{i-1}\right)\nabla y^i-\text{Δ}t\mathcal{C}\left(t_n,t_{i-1}\right)G_H{y}_H^i,\text{d}t\nabla e^n\right)\text{Δ}t\\ +{\displaystyle \sum }_{n=1}^l\left(E\left(u^n-{\widehat{u}}_H^n\right),\text{d}t{e}^n\right)\text{Δ}t+{\displaystyle \sum }_{n=1}^l\left(A\nabla \text{d}t{y}^n-A\nabla y_t^n,\text{d}t\nabla e^n\right)\text{Δ}t\\ =:{\displaystyle \sum }_{i=1}^7{I}_i.\end{array}$(59)

For $I_1$ , using (16), we have

$I_1\le C{h}^4{\Vert y_t\Vert }_{L^2\left(H^2\right)}^2+\frac{1}{4}{\displaystyle \sum }_{n=1}^l{\Vert \text{d}t{e}^n\Vert }^2\text{Δ}t.$ (60)

For $I_2$ , according to the consequence shown in [31], we get

$\begin{array}{c}{I}_2\le C{\displaystyle \sum }_{n=1}^l{\left({\displaystyle {\int }_{t_{n-1}}^{t_n}\Vert y_{tt}\Vert \text{d}t}\right)}^2\text{Δ}t+\frac{1}{4}{\displaystyle \sum }_{n=1}^l{\Vert \text{d}t{e}^n\Vert }^2\text{Δ}t\\ \le C{\left(\text{Δ}t\right)}^2{\displaystyle {\int }_0^{t_l}{\Vert y_{tt}\Vert }^2\text{d}t}+\frac{1}{4}{\displaystyle \sum }_{n=1}^l{\Vert \text{d}t{e}^n\Vert }^2\text{Δ}t\\ \le C{\left(\text{Δ}t\right)}^2{\Vert y_{tt}\Vert }_{L^2\left(L^2\right)}^2+\frac{1}{4}{\displaystyle \sum }_{n=1}^l{\Vert \text{d}t{e}^n\Vert }^2\text{Δ}t.\end{array}$ (61)

For $I_3$ , we have

$\begin{array}{c}\left|I_3\right|\le {\displaystyle \sum }_{n=1}^l\left({\displaystyle {\int }_0^{t_n}\mathcal{C}\left(t_n,s\right)\nabla y\left(s\right)\text{d}s}-{\displaystyle \sum }_{i=1}^n\text{Δ}t\mathcal{C}\left(t_n,t_i\right)\nabla y^i,\nabla e^n\right)\\ +{\displaystyle \sum }_{n=1}^l\left({\displaystyle {\int }_0^{t_n}\mathcal{C}\left(t_n,s\right)\nabla y\left(s\right)\text{d}s}-{\displaystyle \sum }_{i=1}^n\text{Δ}t\mathcal{C}\left(t_n,t_i\right)\nabla y^i,\nabla e^{n-1}\right)\\ ={\displaystyle \sum }_{n=1}^l\left({\displaystyle {\int }_0^{t_n}\mathcal{C}\left(t_n,s\right)\nabla R_{h}y\left(s\right)\text{d}s}-{\displaystyle \sum }_{i=1}^n\text{Δ}t\mathcal{C}\left(t_n,t_i\right)\nabla R_h{y}^i,\nabla e^n\right)\\ +{\displaystyle \sum }_{n=1}^l\left({\displaystyle {\int }_0^{t_n}\mathcal{C}\left(t_n,s\right)\nabla R_{h}y\left(s\right)\text{d}s}-{\displaystyle \sum }_{i=1}^n\text{Δ}t\mathcal{C}\left(t_n,t_i\right)\nabla R_h{y}^i,\nabla e^{n-1}\right)\\ \le C{\left(\text{Δ}t\right)}^2\left({\Vert \nabla R_h{y}_t\Vert }_{L^2\left(L^2\right)}^2+{\Vert \nabla R_{h}y\Vert }_{L^2\left(L^2\right)}^2\right)+C{\displaystyle \sum }_{n=1}^l{\Vert \nabla e^n\Vert }^2\text{Δ}t\\ +C{\displaystyle \sum }_{n=1}^l{\Vert \nabla e^{n-1}\Vert }^2\text{Δ}t.\end{array}$ (62)

where we also used (15), and

$\begin{array}{l}{\displaystyle {\int }_0^{t_n}\mathcal{C}\left(t_n,s\right)\nabla R_{h}y\left(s\right)\text{d}s}-{\displaystyle \sum }_{i=1}^n\mathcal{C}\left(t_n,t_i\right)\nabla R_h{y}^i\text{Δ}t\\ \le C\text{Δ}t\left({\Vert \nabla R_h{y}_t\Vert }_{L^2\left(L^2\right)}+{\Vert \nabla R_{h}y\Vert }_{L^2\left(L^2\right)}\right).\end{array}$

For $I_4$ , we have

$\left|I_4\right|\le C{\displaystyle \sum }_{n=1}^l{\displaystyle \sum }_{i=1}^n{\Vert \nabla y^i\Vert }^2\text{Δ}t+C{\displaystyle \sum }_{n=1}^l{\Vert \nabla e^n\Vert }^2\text{Δ}t+C{\displaystyle \sum }_{n=1}^l{\Vert \nabla e^{n-1}\Vert }^2\text{Δ}t.$ (63)

For $I_5$ , from Lemma 6, we have

$\left|I_5\right|\le C\left(H^4+{\left(\text{Δ}t\right)}^2\right)+\frac{1}{8}{\Vert \nabla e^l\Vert }^2+C{\displaystyle \sum }_{n=1}^l{\Vert \nabla e^n\Vert }^2\text{Δ}t.$ (64)

For $I_6$ , using Theorem 8, we have

$I_6\le C\left(H^4+{\left(\text{Δ}t\right)}^2\right)+\frac{1}{4}{\displaystyle \sum }_{n=1}^l{\Vert \text{d}t{e}^n\Vert }^2\text{Δ}t.$ (65)

Let’s estimate $I_7$ , it’s similar to $I_2$

$\left|I_7\right|\le C{\left(\text{Δ}t\right)}^2{\Vert \nabla y_{tt}\Vert }_{L^2\left(L^2\right)}^2+C{\displaystyle \sum }_{n=1}^l{\Vert \text{d}t\nabla e^n\Vert }^2\text{Δ}t.$ (66)

Adding up $I_1,I_2,I_3,I_4,I_5,I_6,I_7$ , combining (16), the triangle inequality and discrete Gronwall’s inequality, thus

$\Vert \nabla \left(y^n-{\overline{\tilde{y}}}_h^n\right)\Vert \le C\left(h+H^2+\text{Δ}t\right).$(67)

Putting $t=t_{n-1}$ into (9), subtracting (53) from (9), then using (15), we get

$\begin{array}{l}-\left(\text{d}t{\alpha }^n,q_h\right)+\left(D\nabla {\alpha }^{n-1},\nabla q_h\right)\\ =-\left(\text{d}t{p}^n-p_t^{n-1},q_h\right)+\left(\text{d}t{\beta }^n,q_h\right)+{\displaystyle \sum }_{i=n}^N\left(\text{Δ}t{\mathcal{C}}^{*}\left(t_i,t_{n-1}\right)G_H{p}_H^{i-1},\nabla q_h\right)\\ +\left(A\nabla p_t^{n-1},\nabla q_h\right)-\left(A\nabla \text{d}t{\overline{\tilde{p}}}_h^n,\nabla q_h\right)-{\displaystyle {\int }_{t_{n-1}}^T\left({\mathcal{C}}^{*}\left(s,t_{n-1}\right)\nabla p\left(s\right),\nabla q_h\right)\text{d}s}\\ +\left(\delta z_d^n-\delta y^n,q_h\right)+\left(y^n-y_H^n,q_h\right).\end{array}$(68)

Taking $q_h=-\text{d}t{\alpha }^n$ in (68), multiplying it by $\text{Δ}t$ and summing it over n from $l+1$ to N ($0\le l\le N-1$ ) at both sides of (68), making use of (16), (42) and the triangle inequality, thus

$\begin{array}{l}{\displaystyle \sum }_{n=l+1}^N{\Vert \text{d}t{\alpha }^n\Vert }^2\text{Δ}t+{\displaystyle \sum }_{n=l+1}^N{\Vert A^{\frac{1}{2}}\text{d}t\nabla {\alpha }^n\Vert }^2\text{Δ}t+\frac{1}{2}{\Vert D^{\frac{1}{2}}\nabla {\alpha }^l\Vert }^2\\ \le {\displaystyle \sum }_{n=l+1}^N\left(\text{d}t{p}^n-p_t^{n-1},\text{d}t{\alpha }^n\right)\text{Δ}t-{\displaystyle \sum }_{n=l+1}^N\left(\text{d}t{\beta }^n,\text{d}t{\alpha }^n\right)\text{Δ}t\\ +{\displaystyle \sum }_{n=l+1}^N\left({\displaystyle {\int }_{t_{n-1}}^T{\mathcal{C}}^{*}\left(s,t_{n-1}\right)\nabla p\left(s\right)\text{d}s}-{\displaystyle \sum }_{i=n}^N\text{Δ}t{\mathcal{C}}^{*}\left(t_i,t_{n-1}\right)G_H{p}_H^{i-1},\text{d}t\nabla {\alpha }^n\right)\text{Δ}t\\ +{\displaystyle \sum }_{n=l+1}^N\left(A\nabla \text{d}t{p}^n-A\nabla p_t^n,\text{d}t\nabla {\alpha }^n\right)\text{Δ}t+{\displaystyle \sum }_{n=l+1}^N\left(A\nabla \left(p_t^n-p_t^{n-1}\right),\text{d}t\nabla {\alpha }^n\right)\text{Δ}t\\ -{\displaystyle \sum }_{n=l+1}^N\left(\delta z_d^n-\delta y^n,\text{d}t{\alpha }^n\right)\text{Δ}t-{\displaystyle \sum }_{n=l+1}^N\left(y^n-y_H^n,\text{d}t{\alpha }^n\right)\text{Δ}t\\ =:{\displaystyle \sum }_{j=1}^7{J}_j.\end{array}$ (69)

Notice that

$-\left(D\nabla {\alpha }^{n-1},\text{d}t\nabla {\alpha }^n\right)\ge \frac{1}{2\text{Δ}t}\left({\Vert D^{\frac{1}{2}}\nabla {\alpha }^{n-1}\Vert }^2-{\Vert D^{\frac{1}{2}}\nabla {\alpha }^n\Vert }^2\right).$ (70)

Next, we estimate the right sides of (69). Similar to (61), we have

$J_1\le C{\left(\text{Δ}t\right)}^2{\Vert p_{tt}\Vert }_{L^2\left(L^2\right)}^2+\frac{1}{4}{\displaystyle \sum }_{n=l+1}^N{\Vert \text{d}t{\alpha }^n\Vert }^2\text{Δ}t.$ (71)

For $J_2$ , using Cauchy inequality and (16), we get

$J_2\le C{h}^4{\Vert p_t\Vert }_{L^2\left(H^2\right)}^2+\frac{1}{4}{\displaystyle \sum }_{n=l+1}^N{\Vert \text{d}t{\alpha }^n\Vert }^2\text{Δ}t.$ (72)

For $J_3$ , we have

(73)

For $L_1$ , similar to $I_3$ , we get

$\begin{array}{c}{L}_1\le C{\left(\text{Δ}t\right)}^2\left({\Vert \nabla R_h{p}_t\Vert }_{L^2\left(L^2\right)}^2+{\Vert \nabla R_{h}p\Vert }_{L^2\left(L^2\right)}^2\right)\\ +C{\displaystyle \sum }_{n=l+1}^N{\Vert \nabla {\alpha }^n\Vert }^2\text{Δ}t+C{\displaystyle \sum }_{n=l+1}^N{\Vert \nabla {\alpha }^{n-1}\Vert }^2\text{Δ}t.\end{array}$ (74)

For $L_2$ , it’s easy to get

$L_2\le C{\displaystyle \sum }_{n=l+1}^N{\displaystyle \sum }_{i=n}^N{\left(\text{Δ}t\right)}^2{\Vert \nabla p_t\Vert }^2+C{\displaystyle \sum }_{n=l+1}^N{\Vert \nabla {\alpha }^n\Vert }^2\text{Δ}t+C{\displaystyle \sum }_{n=l+1}^N{\Vert \nabla {\alpha }^{n-1}\Vert }^2\text{Δ}t.$ (75)

For $L_3$ , it’s similar to $I_5$

$L_3\le C\left(H^4+{\left(\text{Δ}t\right)}^2\right)+\frac{1}{8}{\Vert \nabla {\alpha }^N\Vert }^2+C{\displaystyle \sum }_{n=l+1}^N{\Vert \nabla {\alpha }^n\Vert }^2\text{Δ}t.$ (76)

For $J_4$ , it’s similar to (61), we get

$J_4\le C{\left(\text{Δ}t\right)}^2{\Vert \nabla p_{tt}\Vert }_{L^2\left(L^2\right)}^2+\frac{1}{4}{\displaystyle \sum }_{n=l+1}^N{\Vert \text{d}t\nabla {\alpha }^n\Vert }^2\text{Δ}t.$ (77)

For $J_5$ , it’s easy to get

$J_5\le C{\left(\text{Δ}t\right)}^2{\Vert \nabla p_{tt}\Vert }_{L^2\left(L^2\right)}^2+\frac{1}{4}{\displaystyle \sum }_{n=l+1}^N{\Vert \text{d}t\nabla {\alpha }^n\Vert }^2\text{Δ}t.$ (78)

For $J_6$ , using Cauchy inequality and the smoothness of y and $z_d$ , we get

$J_6\le C{\left(\text{Δ}t\right)}^2\left({\Vert y_t\Vert }_{L^2\left(L^2\right)}^2+{\Vert {\left(z_d\right)}_t\Vert }_{L^2\left(L^2\right)}^2\right)+\frac{1}{8}{\displaystyle \sum }_{n=l+1}^N\Vert \text{d}t{\alpha }^n{}^2\Vert \text{Δ}t.$ (79)

For $J_7$ , it’s easy to get

$\begin{array}{c}{J}_7\le C{\Vert y^n-y_H^n\Vert }_{L^2\left(L^2\right)}^2+\frac{1}{8}{\displaystyle \sum }_{n=l+1}^N{\Vert \text{d}t{\alpha }^n\Vert }^2\text{Δ}t\\ \le C(H^4+\left(\text{Δ}t{)}^2\right)+\frac{1}{8}{\displaystyle \sum }_{n=l+1}^N{\Vert \text{d}t{\alpha }^n\Vert }^2\text{Δ}t.\end{array}$ (80)

Suming up $J_1,J_2,J_3,J_4,J_5,J_6,J_7$ , combining (16), the triangle inequality and discrete Gronwall’s inequality, we can get

$\Vert \nabla \left(p^{n-1}-{\overline{\tilde{p}}}_h^{n-1}\right)\Vert \le C\left(h+H^2+\text{Δ}t\right).$(81)

Note that

$u^n=-\frac{1}{\gamma }{E}^{*}{p}^n+\max \left\{0,\frac{1}{\gamma \left|{\Omega }_U\right|}{\displaystyle {\int }_{{\Omega }_U}{E}^{*}{p}^n}\right\},$

${\overline{\tilde{u}}}_h^n=-\frac{1}{\gamma }{E}^{*}{\pi }_h{\overline{\tilde{p}}}_H^{n-1}+\max \left\{0,\frac{1}{\gamma \left|{\Omega }_U\right|}{\displaystyle {\int }_{{\Omega }_U}{E}^{*}{\overline{\tilde{p}}}_H^{n-1}}\right\}.$

It’s similar to (44)

$\begin{array}{c}\Vert u^n-{\overline{\tilde{u}}}_h^n\Vert \le \frac{1}{\gamma }\Vert E^{*}{p}^n-E^{*}{\pi }_h{\overline{\tilde{p}}}_H^{n-1}\Vert +\max \left\{0,\frac{1}{\gamma \left|{\Omega }_U\right|}{\displaystyle {\int }_{{\Omega }_U}{E}^{*}{p}^n}\right\}\\ -\max \left\{0,\frac{1}{\gamma \left|{\Omega }_U\right|}{\displaystyle {\int }_{{\Omega }_U}{E}^{*}{\overline{\tilde{p}}}_H^{n-1}}\right\}\\ \le C\Vert p^n-{\pi }_h{\overline{\tilde{p}}}_H^{n-1}\Vert +\frac{1}{\gamma \left|{\text{Ω}}_U\right|}\Vert {\displaystyle {\int }_{{\Omega }_U}\left(E^{*}{p}^n-E^{*}{\overline{\tilde{p}}}_H^{n-1}\right)}\Vert \\ \le C\Vert p^n-p^{n-1}\Vert +C\Vert p^{n-1}-{\pi }_h{p}^{n-1}\Vert \\ +C\Vert {\pi }_h{p}^{n-1}-{\pi }_h{\overline{\tilde{p}}}_H^{n-1}\Vert +C\Vert p^n-{\overline{\tilde{p}}}_H^{n-1}\Vert \\ \le C\Vert p^n-p^{n-1}\Vert +C\Vert p^{n-1}-{\pi }_h{p}^{n-1}\Vert +C\Vert p^{n-1}-{\overline{\tilde{p}}}_H^{n-1}\Vert \\ \le C\left(h+H^2+\text{Δ}t\right).\end{array}$(82)

□

5. Numerical Experiment

In this section, we implement an experiment to check the theoretical results studied in Section 3 and Section 4. The numerical experiment was done by utilizing MATLAB finite element package iFEM [32].

In the numerical experiment, we take $\Omega ={\Omega }_U={\left[0,1\right]}^2$ . The stopping criterion is that the control variable satisfies $\Vert u^{n+1}-u^n\Vert <{10}^{-5}$ . We mainly show the error of the FEM, and the gradient recovery based two-grid finite element method.

We solved the following control problem:

$\underset{K}{\text{min}}\frac{1}{2}{\displaystyle {\int }_0^T\left({\displaystyle {\int }_{\Omega }{\left(y-z_d\right)}^2}+{\displaystyle {\int }_{\Omega }{u}^2}\right)}$ (83)

subject to

$\{\begin{array}{l}\frac{\partial y}{\partial t}-\text{Δ}{y}_t-\text{Δ}y-{\displaystyle {\int }_0^t\left(t-s\right)\text{Δ}y}=f+u \text{in}\Omega ,0<t\le 0.01;\hfill \\ {y|}_{\partial \Omega }=0,\hfill \end{array}$ (84)

where $K=\left\{u\in U_{ad}:{\displaystyle {\int }_{\Omega }u\left(x,t\right)\text{d}x}\ge 0\right\}$ , with data and solutions are given as:

$\{\begin{array}{l}p=-\left(T-t\right)\sin \pi x_1\sin \pi x_2,\hfill \\ u=\max \left(-p,0\right),\hfill \\ y=e^{2t}{x}_1\left(1-x_1\right)x_2\left(1-x_2\right),\hfill \\ z_d=y+\frac{\partial p}{\partial t}-\Delta p_t+\Delta p+{\displaystyle {\int }_t^T\left(s-t\right)\Delta p},\hfill \\ f=\frac{\partial y}{\partial t}-\Delta y_t-\Delta y-{\displaystyle {\int }_0^t\left(t-s\right)\Delta y}-u.\hfill \end{array}$ (85)

In the numerical experiment, we use $\text{Δ}t=\frac{h}{1000}$ so that the error in time direction does not influence the error in spatial direction, and choose $T=0.01$ to compute time levels 400, 1600, 6400 for $h=\frac{1}{4},\frac{1}{16},\frac{1}{64}$ . In Table 1, we show the error of finite element method for y, p in both $L^2$ norm and $H^1$ norm and control variable u in $L^2$ norm. In Table 2, we show the error of gradient recovery based two-grid finite element solutions. Table 3 gives the calculation time of numerical example using FEM and gradient recovery based two-grid finite element method respectively.

Table 1. Error of finite element method with $h=\frac{1}{4},\frac{1}{16},\frac{1}{64}$ , $\text{Δ}t=\frac{h}{1000}$ at $t=0.005$ .

Freedom number	$L^2\left(H^1\left(\Omega \right)\right)$		$L^2\left(L^2\left(\Omega \right)\right)$
h	$y-y_h$	$p-p_h$	$u-u_h$	$y-y_h$	$p-p_h$
$\frac{1}{4}$	0.059117	0.003992	0.000885	0.003970	0.000349
$\frac{1}{16}$	0.015191	0.001074	0.000186	0.000265	2.45275e−05
$\frac{1}{64}$	0.003805	0.000272	4.44124e−05	2.77247e−05	1.49681e−06

Table 2. Error of gradient recovery based two-grid method with $h=H^2=\frac{1}{4},\frac{1}{16},\frac{1}{64}$ , $\text{Δ}t=\frac{h}{1000}$ at $t=0.005$ .

Freedom number	$L^2\left(H^1\left(\Omega \right)\right)$		$L^2\left(L^2\left(\Omega \right)\right)$
h	$y-y_h$	$p-p_h$	$u-u_h$	$y-y_h$	$p-p_h$
$\frac{1}{4}$	0.059117	0.003992	0.000886	0.003971	0.000349
$\frac{1}{16}$	0.015192	0.001075	0.000187	0.000257	2.46157e−05
$\frac{1}{64}$	0.003805	0.000272	4.44501e−05	1.61259e−05	1.58670e−06

Table 3. Calculation time.

Freedom number	Finite element method	Gradient recovery based two-grid method
h	Cputime(s)	Cputime(s)
$\frac{1}{4}$	0.7790	0.9920
$\frac{1}{16}$	4.7140	5.4150
$\frac{1}{64}$	342.0670	197.9580

The convergence order graphs of three variables calculated by finite element method and gradient recovery based on two-grid finite element method are also given, see Figure 1 and Figure 2 respectively. From the figures, we can know that the numerical results are consistent with the theory discussed in Section 3-4.

Figure 1. The convergence order of finite element method at $t=0.005$ .

Figure 2. The convergence order of doing gradient recovery based two-grid finite element method at $t=0.005$ .

6. Conclusion

In this paper, we present a gradient recovery based two-grid finite element method for parabolic integro-differential optimal control problem, which is of innovative significance. Combining the high efficiency of two-grid finite element method and the high precision of gradient recovery, we estimate the priori error of state variable, co-state variable and control variable. Finally, a numerical example is used to illustrate the correctness of the theoretical results. In our future work, we will study a posteriori error estimation as well as an adaptive method for the constraint optimal control problem.

Conflicts of Interest

The author declares no conflicts of interest regarding the publication of this paper.

References

[1]	Barcikowski, A., Goeppner, G., Norem, J., Rotela, E., Rusthoven, B., Sharma, S., et al. (2001). An E/sup +/e/sup-/ Collider in the VLHC Tunnel. PACS2001. Proceedings of the 2001 Particle Accelerator Conference, Chicago, 18-22 June 2001, 3493-3495. https://doi.org/10.1109/pac.2001.988155
[2]	Locatelli, A. and Sieniutycz, S. (2002) Optimal Control: An Introduction. Applied Mechanics Reviews, 55, B48-B49. https://doi.org/10.1115/1.1470675
[3]	Zhao, C., Wang, M. and He, Z. (2004) Well Posedness and Optimal Control Problems for a Class of Periodic Population Dynamics. Mathematica Applicata, 17, 551-556.
[4]	İskender Eroğlu, B.B., Avcı, D. and Özdemir, N. (2017) Optimal Control Problem for a Conformable Fractional Heat Conduction Equation. Acta Physica Polonica A, 132, 658-662. https://doi.org/10.12693/aphyspola.132.658
[5]	Fosha, C. and Elgerd, O. (1970) The Megawatt-Frequency Control Problem: A New Approach via Optimal Control Theory. IEEE Transactions on Power Apparatus and Systems, 89, 563-577. https://doi.org/10.1109/tpas.1970.292603
[6]	Kambe, T. (1984) Axisymmetric Vortex Solution of Navier-Stokes Equation. Journal of the Physical Society of Japan, 53, 13-15. https://doi.org/10.1143/jpsj.53.13
[7]	Chen, Y., Huang, Y. and Liu, W. (2010) Error Estimates and Superconvergence of Mixed Finite Element Methods for Convex Optimal Control Problems. Journal of Scientific Computing, 42, 382-403.
[8]	Tang, Y. (2021) Convergence and Superconvergence of Fully Discrete Finite Element for Time Fractional Optimal Control Problems. American Journal of Computational Mathematics, 11, 53-63. https://doi.org/10.4236/ajcm.2021.111005
[9]	Lewis, F.L., Vrabie, D.L. and Syrmos, V.L. (2012) Optimal Control. Wiley. https://doi.org/10.1002/9781118122631
[10]	Xu, J. (1994) A Novel Two-Grid Method for Semilinear Elliptic Equations. SIAM Journal on Scientific Computing, 15, 231-237. https://doi.org/10.1137/0915016
[11]	Xu, J. (2000) Some Two-Grid Finite Element Methods. Contemporary Mathematics, 157, 79.
[12]	Huang, Y., Liang, Q. and Yi, N. (2010) High Order Compact Schemes for Gradient Approximation. Science China Mathematics, 53, 1903-1918. https://doi.org/10.1007/s11425-010-3081-0
[13]	Xu, J. and Zhang, Z. (2003) Analysis of Recovery Type a Posteriori Error Estimators for Mildly Structured Grids. Mathematics of Computation, 73, 1139-1153. https://doi.org/10.1090/s0025-5718-03-01600-4
[14]	Yan, N. and Zhou, A. (2001) Gradient Recovery Type a Posteriori Error Estimates for Finite Element Approximations on Irregular Meshes. Computer Methods in Applied Mechanics and Engineering, 190, 4289-4299. https://doi.org/10.1016/s0045-7825(00)00319-4
[15]	Liu, W. and Yan, N. (2001) Some a Posteriori Error Estimators for p-Laplacian Based on Residual Estimation or Gradient Recovery. Journal of Scientific Computing, 16, 435-477. https://doi.org/10.1023/a:1013246424707
[16]	Cao, H., Huang, Y. and Yi, N. (2023) Gradient Recovery Based a Posteriori Error Estimator for the Adaptive Direct Discontinuous Galerkin Method. Calcolo, 60, Article No. 18. https://doi.org/10.1007/s10092-023-00513-9
[17]	Lakkis, O. and Pryer, T. (2011) Gradient Recovery in Adaptive Finite-Element Methods for Parabolic Problems. IMA Journal of Numerical Analysis, 32, 246-278. https://doi.org/10.1093/imanum/drq019
[18]	Yan, N. (2003) A Posteriori Error Estimators of Gradient Recovery Type for Fem of a Model Optimal Control Problem. Advances in Computational Mathematics, 19, 323-336. https://doi.org/10.1023/a:1022800401298
[19]	Wen, Y., Chen, L. and Dai, J. (2023) A Posteriori Error Estimate of Two Grid Mixed Finite Element Methods for Semilinear Elliptic Equations. Journal of Applied Mathematics and Physics, 11, 361-376. https://doi.org/10.4236/jamp.2023.112020
[20]	Taiwo, O.A., Alhassan, L.K., Odetunde, O.S. and Alabi, O.O. (2023) Galerkin Method for Numerical Solution of Volterra Integro-Differential Equations with Certain Orthogonal Basis Function. International Journal of Modern Nonlinear Theory and Application, 12, 68-80. https://doi.org/10.4236/ijmnta.2023.122005
[21]	Yi, N. (2011) A Posteriori Error Estimates Based on Gradient Recovery and Adaptive Finite Element Methods. Xiangtan University Press, 30-44.
[22]	Russell, D.L., Lions, J.L. and Mitter, S.K. (1972) Optimal Control of Systems Governed by Partial Differential Equations. The American Mathematical Monthly, 79, 1049-1050. https://doi.org/10.2307/2318095
[23]	Shen, W., Ge, L. and Yang, D. (2013) Finite Element Methods for Optimal Control Problems Governed by Linear Quasi-Parabolic Integro-Differential Equations. International Journal of Numerical Analysis and Modeling, 10, 536-550.
[24]	Meyer, C. and Rösch, A. (2004) Superconvergence Properties of Optimal Control Problems. SIAM Journal on Control and Optimization, 43, 970-985. https://doi.org/10.1137/s0363012903431608
[25]	Ciarlet, P.G. and Oden, J.T. (1978) The Finite Element Method for Elliptic Problems. Journal of Applied Mechanics, 45, 968-969. https://doi.org/10.1115/1.3424474
[26]	Douglas, J. and Roberts, J.E. (1985) Global Estimates for Mixed Methods for Second Order Elliptic Equations. Mathematics of Computation, 44, 39-52. https://doi.org/10.1090/s0025-5718-1985-0771029-9
[27]	Liu, H. and Yan, N. (2007) Recovery Type Superconvergence and a Posteriori Error Estimates for Control Problems Governed by Stokes Equations. Journal of Computational and Applied Mathematics, 209, 187-207. https://doi.org/10.1016/j.cam.2006.10.083
[28]	Xu, C. and Li, H. (2023) Two-Grid Methods of Finite Element Approximation for Parabolic Integro-Differential Optimal Control Problems. Electronic Research Archive, 31, 4818-4842. https://doi.org/10.3934/era.2023247
[29]	Zienkiewicz, O.C. and Zhu, J.Z. (1992) The Superconvergent Patch Recovery and a Posteriori Error Estimates. Part 2: Error Estimates and Adaptivity. International Journal for Numerical Methods in Engineering, 33, 1365-1382. https://doi.org/10.1002/nme.1620330703
[30]	Zienkiewicz, O.C. and Zhu, J.Z. (1992) The Superconvergent Patch Recovery (SPR) and Adaptive Finite Element Refinement. Computer Methods in Applied Mechanics and Engineering, 101, 207-224. https://doi.org/10.1016/0045-7825(92)90023-d
[31]	Russell, T.F. (1985) Time Stepping along Characteristics with Incomplete Iteration for a Galerkin Approximation of Miscible Displacement in Porous Media. SIAM Journal on Numerical Analysis, 22, 970-1013. https://doi.org/10.1137/0722059
[32]	Chen, L. (2009) iFEM: An Integrated Finite Element Methods Package in MATLAB. Technical Report.