The Finite Volume Element Method for Time-Fractional Nonlinear Fourth-Order Diffusion Equation with Time Delay

Anran Li^*, Qing Yang

School of Mathematics and Statistics, Shandong Normal University, Ji’nan, China.

DOI: 10.4236/eng.2025.171004                

Abstract

In this article, a finite volume element algorithm is presented and discussed for the numerical solutions of a time-fractional nonlinear fourth-order diffusion equation with time delay. By choosing the second-order spatial derivative of the original unknown as an additional variable, the fourth-order problem is transformed into a second-order system. Then the fully discrete finite volume element scheme is formulated by using $L1$ approximation for temporal Caputo derivative and finite volume element method in spatial direction. The unique solvability and stable result of the proposed scheme are proved. A priori estimate of $L^2$ -norm with optimal order of convergence $O\left(h^2+{\tau }^{2-\alpha }\right)$ where $\tau$ and $h$ are time step length and space mesh parameter, respectively, is obtained. The efficiency of the scheme is supported by some numerical experiments.

Keywords

Time-Fractional Nonlinear Fourth-Order Diffusion Equation with Time Delay, Finite Volume Element Method, Caputo-Fractional Derivative, Optimal Priori Error Analysis

1. Introduction

Nowadays, researchers have placed more attention on the development of fractional differential equations as these equations are widely used in fractal media, mathematical biology, chemistry, statistical mechanics, engineering and so on [1]-[7]. Time delay occurs in many real-life applications such as population ecology, cell biology, control theory [8]-[13]. Therefore, development of numerical methods for fractional equations with time delay seems to be vital and essential.

In recent years, various numerical methods and theory of fractional differential equations have been studied extensively by researchers and their study comprises numerical methods such as finite difference, finite volume, finite element and so on. In [14], Danumjaya P et al. applied the mixed finite element methods to a fourth order reaction diffusion equation with different types of boundary conditions and established some priori bounds with the help of Lyapunov functional. In [15], Yang Liu et al. presented a finite difference/finite element algorithm, which is based on a finite difference approximation in time direction and finite element method in spatial direction, and discussed the numerical solutions of a time-fractional fourth-order reaction-diffusion problem with a nonlinear reaction term. Tie Zhang et al. in [16] studied the finite volume method for solving the time-fractional diffusion equations and analyzed a fully discrete numerical scheme which is based on the linear finite volume method and the L1 difference. Xinfei Liu et al. in [17] considered the nonlinear time-fractional stochastic fourth-order reaction-diffusion equation perturbed by noises based on the mixed finite element in spatial direction and the generalized BDF2-$\theta$ in temporal discretization, and obtained the semi- and fully-discrete schemes.

There have been many studies on nonlinear time delay differential equations with spatial second derivative [18]-[21]. However, limited work has been done for nonlinear fourth-order differential equations with time delay. Sarita Nandal et al. in [22] constructed a compact difference scheme for one-dimensional time-fractional fourth-order nonlinear sub-diffusion wave equation with time delay and conducted the numerical analysis of the scheme using discrete energy method. In [23] Hongxia Xie et al. constructed a compact difference scheme for two-dimensional time-fractional nonlinear fourth-order diffusion equation with time delay and proved the convergence rate in time and space.

In this article, we take into account the following time-fractional nonlinear fourth-order diffusion equation with time delay,

$\begin{array}{l}{}_0^C{D}_t^{\alpha }u\left(x\mathrm{,}y\mathrm{,}t\right)+{\Delta }^{2}u\left(x\mathrm{,}y\mathrm{,}t\right)\\ =f\left(x\mathrm{,}y\mathrm{,}t\mathrm{,}u\left(x\mathrm{,}y\mathrm{,}t\right)\mathrm{,}u\left(x\mathrm{,}y\mathrm{,}t-s\right)\right)\mathrm{,}\left(x\mathrm{,}y\mathrm{,}t\right)\in \Omega \times \left(\mathrm{0,}T\right]\mathrm{,}\end{array}$ (1.1a)

$u\left(x\mathrm{,}y\mathrm{,}t\right)=\Delta u\left(x\mathrm{,}y\mathrm{,}t\right)=\mathrm{0,}\left(x\mathrm{,}y\mathrm{,}t\right)\in \partial \Omega \times \left[\mathrm{0,}T\right]\mathrm{,}$ (1.1b)

$u\left(x\mathrm{,}y\mathrm{,}t\right)=g\left(x\mathrm{,}y\mathrm{,}t\right)\mathrm{,}\left(x\mathrm{,}y\mathrm{,}t\right)\in \overline{\Omega }\times \left[-s\mathrm{,0}\right]\mathrm{,}$ (1.1c)

where $\alpha \in \left(\mathrm{0,1}\right)$ , $\Omega =\left(\mathrm{0,}{L}_1\right)\times \left(\mathrm{0,}{L}_2\right)$ , $s>0$ is delay and $f\left(x\mathrm{,}y\mathrm{,}t\mathrm{,}u\left(x\mathrm{,}y\mathrm{,}t\right)\mathrm{,}u\left(x\mathrm{,}y\mathrm{,}t-s\right)\right)$ stands for nonlinear time delay source term, $g\left(x\mathrm{,}y\mathrm{,}t\right)$ is given and sufficiently smooth function. The fractional derivative ${}_0^c{D}_t^{\alpha }\left(x\mathrm{,}y\mathrm{,}t\right)$ is considered in Caputo sense as follows

${}_0^c{D}_t^{\alpha }u\left(x\mathrm{,}y\mathrm{,}t\right)=\frac{1}{\Gamma \left(1-\alpha \right)}{\displaystyle {\int }_0^t}\frac{\partial u\left(x\mathrm{,}y\mathrm{,}\xi \right)}{\partial \xi }\frac{1}{{\left(t-\xi \right)}^{\alpha }}\text{d}\xi \mathrm{.}$

Throughout the article, we assume that the source function $f\left(x\mathrm{,}y\mathrm{,}t\mathrm{,}\mu \mathrm{,}\nu \right)$ is sufficiently smooth likewise considered in the following sense:

The partial derivatives $f_{\mu }\left(x\mathrm{,}y\mathrm{,}t\mathrm{,}\mu \mathrm{,}\nu \right)$ and $f_{\nu }\left(x\mathrm{,}y\mathrm{,}t\mathrm{,}\mu \mathrm{,}\nu \right)$ are continuous in the ${\epsilon }_0$ neighborhood of the solution and let

$\begin{array}{l}{\overline{c}}_1=\underset{\left(x,y,t\right)\in \Omega \times \left(0,T\right],\left|{\epsilon }_1\right|\le {\epsilon }_0,\left|{\epsilon }_2\right|\le {\epsilon }_0}{\sup }\left|f_{\mu }\left(x,y,t,u\left(x,y,t\right)+{\epsilon }_1,u\left(x,y,t-s\right)+{\epsilon }_2\right)\right|,\\ {\overline{c}}_2=\underset{\left(x,y,t\right)\in \Omega \times \left(0,T\right],\left|{\epsilon }_1\right|\le {\epsilon }_0,\left|{\epsilon }_2\right|\le {\epsilon }_0}{\sup }\left|f_{\nu }\left(x,y,t,u\left(x,y,t\right)+{\epsilon }_1,u\left(x,y,t-s\right)+{\epsilon }_2\right)\right|,\\ {\overline{c}}_3=\max \left\{{\overline{c}}_1,{\overline{c}}_2\right\}.\end{array}$ (1.2)

2. Preliminary

To construct a finite volume element method for problem (1.1), we firstly divide the region $\overline{\Omega }\times \left[-s\mathrm{,}T\right]$ . Taking a positive integer $N$ , we define the temporal step size $\tau =\frac{T}{N}$ and $m=\frac{s}{\tau }$ .

Suppose $\Omega$ is a polygonal region with boundary $\partial \Omega$ . Divide $\overline{\Omega }$ into a sum of finite number of small triangles called elements that they have no overlapping internal region. All the elements constitute a triangulation of $\overline{\Omega }$ , denoted by $T_h$ , where $h$ is the maximum length of all the sides. Then we construct a dual decomposition $T_h^{\mathrm{<mo>*</mo>}}$ related to $T_h$ . Let $P_0$ be a node of a triangle, $P_i\left(i=\mathrm{1,2,}\cdots \mathrm{,6}\right)$ the adjacent nodes of $P_0$ . Choose the barycenter $Q_i$ of the triangle $\Delta P_0{P}_i{P}_{i+1}\left(P_7=P_1\right)$ and the midpoint $M_i$ of $\overline{P_0{P}_i}$ and connect successively to form a dual element $K_{P_0}^{\mathrm{<mo>*</mo>}}$ as shown in Figure 1. All the dual elements constitute the dual grid $T_h^{\mathrm{<mo>*</mo>}}$ .

Figure 1. Barycenter dual decomposition.

In this paper, we denote by ${\overline{\Omega }}_h$ the set of the nodes of the decomposition $T_h$ , ${\dot{\Omega }}_h={\overline{\Omega }}_h\backslash \partial \Omega$ the set of the interior nodes, and ${\Omega }_h^{\mathrm{<mo>*</mo>}}$ the set of the nodes of the dual decomposition $T_h^{\mathrm{<mo>*</mo>}}$ . We assume that $T_h$ and $T_h^{\mathrm{<mo>*</mo>}}$ are quasi-uniform, i.e. let $S_{K_Q}$ and $S_{P_0}^{\mathrm{<mo>*</mo>}}$ be the areas of the triangular element $K_Q$ and dual element $S_{P_0}^{\mathrm{<mo>*</mo>}}$ respectively, there exist constants $c_1\mathrm{,}{c}_2\mathrm{,}{c}_3>0$ independent of $h$ that

$c_1{h}^2\le S_{K_Q}\le h^2\mathrm{,}Q\in {\Omega }_h^{\mathrm{<mo>*</mo>}}\mathrm{,}$

$c_2{h}^2\le S_{P_0}^{\mathrm{<mo>*</mo>}}\le c_3{h}^2\mathrm{,}{P}_0\in {\overline{\Omega }}_h\mathrm{.}$

Lemma 2.1 [20]. The $L1$ approximation formula for Caputo fractional derivative of $\alpha \left(0<\alpha <1\right)$ order is given by

$\begin{array}{l}{}_0^c{D}_t^{\alpha }f\left(t_n\right)\\ =\frac{{\tau }^{-\alpha }}{\Gamma \left(2-\alpha \right)}\left[a_0^{\left(\alpha \right)}f\left(t_n\right)-{\displaystyle \sum _{k=1}^{n-1}}\left(a_{n-k-1}^{\left(\alpha \right)}-a_{n-k}^{\left(\alpha \right)}\right)f\left(t_k\right)-a_{n-1}^{\left(\alpha \right)}f\left(t_0\right)\right]+E_0^n,\\ :=\frac{{\tau }^{-\alpha }}{\Gamma \left(2-\alpha \right)}{\displaystyle \sum _{k=1}^n}{a}_{n-k}^{\left(\alpha \right)}\left(f\left(t_k\right)-f\left(t_{k-1}\right)\right)+E_0^n,\\ ={\partial }_t^{\left(\alpha \right)}f\left(t_n\right)+E_0^n,\end{array}$ (2.1)

where $a_l^{\alpha }={\left(l+1\right)}^{1-\alpha }-l^{1-\alpha },l=0,1,\cdots$ , and

$\left|E_0^n\right|\le \frac{1}{2\Gamma \left(2-\alpha \right)}\left[\frac{1}{4}+\frac{\alpha }{\left(1-\alpha \right)\left(2-\alpha \right)}\right]\underset{t_0\le t\le t_n}{\max }\left|f_{tt}\right|{\tau }^{2-\alpha }\mathrm{.}$ (2.2)

The following statements hold for $a_l^{\left(\alpha \right)}$ :

$\begin{array}{l}1=a_0^{\alpha }>a_1^{\alpha }>a_2^{\alpha }>\cdots >a_l^{\alpha }>0;a_0^{\alpha }\to 0,\text{if}l\to \infty ,\\ \left(1-\alpha \right)l^{-\alpha }<a_{l-1}^{\left(\alpha \right)}<\left(1-\alpha \right){\left(l-1\right)}^{-\alpha },l\ge 1.\end{array}$ (2.3)

Lemma 2.2 [24]. Let $\left\{y_k\right\}$ and $\left\{z_k\right\}$ be nonnegative sequences. If

$z_k\le C+{\displaystyle \sum _{0\le k<n}}{y}_k{z}_k,k\ge 0,n\ge 1,$

then it holds that

$z_k\le C\cdot \exp \left({\displaystyle \sum _{0\le k<n}}{y}_k\right),k\ge 0,n\ge 1,$

where $C$ is a nonnegative constant.

3. Formula of the Finite Volume Element Scheme

In this section, we will present the derivation of the finite volume element scheme approximating problem (1.1). Let $v=-\Delta u$ be an auxiliary variable and the problem (1.1) can be rewritten as the following system:

$\begin{array}{l}{}_0^C{D}_t^{\alpha }u\left(x\mathrm{,}y\mathrm{,}t\right)-\Delta v\left(x\mathrm{,}y\mathrm{,}t\right)\\ =f\left(x\mathrm{,}y\mathrm{,}t\mathrm{,}u\left(x\mathrm{,}y\mathrm{,}t\right)\mathrm{,}u\left(x\mathrm{,}y\mathrm{,}t-s\right)\right)\mathrm{,}\left(x\mathrm{,}y\mathrm{,}t\right)\in \Omega \times \left(\mathrm{0,}T\right]\mathrm{,}\end{array}$ (3.1a)

$v\left(x\mathrm{,}y\mathrm{,}t\right)=-\Delta u\left(x\mathrm{,}y\mathrm{,}t\right)\mathrm{,}\left(x\mathrm{,}y\mathrm{,}t\right)\in \Omega \times \left(\mathrm{0,}T\right]\mathrm{,}$ (3.1b)

$u\left(x\mathrm{,}y\mathrm{,}t\right)=v\left(x\mathrm{,}y\mathrm{,}t\right)=\mathrm{0,}\left(x\mathrm{,}y\mathrm{,}t\right)\in \partial \Omega \times \left[\mathrm{0,}T\right]\mathrm{,}$ (3.1c)

$u\left(x,y,t\right)=g\left(x\mathrm{,}y\mathrm{,}t\right),\left(x\mathrm{,}y\mathrm{,}t\right)\in \overline{\Omega }\times \left[-s,0\right].$(3.1d)

Making use of Green’s formula, the corresponding weak formulation of (3.1) is to seek $\left(u\mathrm{,}v\right)\mathrm{<mo>:</mo>}\left[\mathrm{0,}T\right]\mapsto H_0^1\left(\Omega \right)\times H_0^1\left(\Omega \right)$ satisfying, for any $\left(\phi \mathrm{,}\psi \right)\in H_0^1\left(\Omega \right)\times H_0^1\left(\Omega \right)$

$\left({}_0^c{D}_t^{\alpha }u\mathrm{,}\phi \right)+a\left(v\mathrm{,}\phi \right)=\left(f\mathrm{,}\phi \right)\mathrm{,}\forall \phi \in H_0^1\left(\Omega \right)\mathrm{,}$ (3.2a)

$a\left(u\mathrm{,}\psi \right)=\left(v\mathrm{,}\psi \right)\mathrm{,}\forall \psi \in H_0^1\left(\Omega \right)\mathrm{,}$ (3.2b)

with $u\left(x\mathrm{,}y\mathrm{,}t\right)=g\left(x\mathrm{,}y\mathrm{,}t\right)\mathrm{,}\left(x\mathrm{,}y\mathrm{,}t\right)\in \overline{\Omega }\times \left[-s\mathrm{,0}\right]$, where $a\left(u\mathrm{,}v\right)={\displaystyle {\int }_{\Omega }}\nabla u\cdot \nabla v\text{d}x\text{d}y$ . Define the space $U_h$ as the set of piecewise-linear polynomials with respect to $T_h$ , which can be expressed by

where is the set of all linear polynomials on $K$ . It is obvious that $U_h$ is the subspace of $H_0^1\left(\Omega \right)$ . Then we choose the test function space $V_h$ as the piecewise constant function space corresponding to $T_h^{\mathrm{<mo>*</mo>}}$ :

$V_h=\left\{v_h\mathrm{<mo>:</mo>}{v}_h\text{is}\text{a}\text{constant}\text{on}\text{the}\text{interior}\text{of}\text{each}{K}^{\mathrm{<mo>*</mo>}}\in T_h^{\mathrm{<mo>*</mo>}}\right\}\mathrm{.}$

Set

$U_{0h}=\left\{u_h\in U_h\mathrm{<mo>:</mo>}{u}_h\left(P_0\right)=\mathrm{0,}\forall P_0\in {\overline{\Omega }}_h\cap \partial \Omega \right\}\mathrm{,}$

$V_{0h}=\left\{v_h\in V_h\mathrm{<mo>:</mo>}{v}_h\left(P_0\right)=\mathrm{0,}\forall P_0\in {\overline{\Omega }}_h\cap \partial \Omega \right\}\mathrm{.}$

Let ${\Pi }_h^{\mathrm{<mo>*</mo>}}$ be the interpolation project from $U_h$ to $V_h$ :

${\Pi }_h^{*}{\omega }_h={\displaystyle \sum _{P_0\in {\overline{\Omega }}_h}}{\omega }_h\left(P_0\right){\chi }_{P_0},{\omega }_h\in U_h,$(3.3)

where ${\chi }_{P_0}$ is the characteristic function of the set $K_{P_0}^{\mathrm{<mo>*</mo>}}$ .

We define $u^k=u\left(\cdot \mathrm{,}{t}_k\right)$ . Using Taylor’s series, the following equations can be easily obtained

$u^k=2u^{k-1}-u^{k-2}+O\left({\tau }^2\right)\mathrm{,}k\in N_{+}\mathrm{.}$

For $u^{k-m}$ with $m=\frac{s}{\tau }$ when $m\notin N_{+}$ :

Case 1: $m\in \left(\mathrm{0,1}\right)$

$u^{k-m}=\left(2-m\right)u^{k-1}-\left(1-m\right)u^{k-2}+O\left({\tau }^2\right)\mathrm{<mo>;</mo>}$

Case 2: $m>1$

$u^{k-m}=\left(\lceil m\rceil -m\right)u^{k-\lfloor m\rfloor }+\left(m-\lfloor m\rfloor \right)u^{k-\lceil m\rceil }+O\left({\tau }^2\right)\mathrm{.}$

Denote

${\widehat{u}}^k=2u^{k-1}-u^{k-2}\mathrm{,}$

and

Linearization of the non-linear source term $f\left(x\mathrm{,}y\mathrm{,}t\mathrm{,}u\left(x\mathrm{,}y\mathrm{,}t\right)\mathrm{,}u\left(x\mathrm{,}y\mathrm{,}t-s\right)\right)$ by Taylor’s series yields

$f\left(x_i\mathrm{,}{y}_j\mathrm{,}{t}_k\mathrm{,}{u}_{ij}^k\mathrm{,}{u}_{ij}^{k-m}\right)=f\left(x_i\mathrm{,}{y}_j\mathrm{,}{t}_k\mathrm{,}{\widehat{u}}_{ij}^k\mathrm{,}{\tilde{u}}_{ij}^{k-m}\right)+O\left({\tau }^2\right)\mathrm{.}$

We denote ${\overline{f}}_h=f\left(\cdot \mathrm{,}{\widehat{u}}_h\mathrm{,}{\tilde{u}}_h\right)$ and $E_1=O\left({\tau }^2\right)$ .

The semidiscrete finite volume scheme is to find a pair $\left(u_h\mathrm{,}{v}_h\right)\mathrm{<mo>:</mo>}\left[\mathrm{0,}T\right]\mapsto U_{0h}\times U_{0h}$ such that

$\left({}_0^c{D}_t^{\alpha }{u}_h\mathrm{,}{\Pi }_h^{\mathrm{<mo>*</mo>}}{\phi }_h\right)+a\left(v_h\mathrm{,}{\Pi }_h^{\mathrm{<mo>*</mo>}}{\phi }_h\right)=\left({\overline{f}}_h\mathrm{,}{\Pi }_h^{\mathrm{<mo>*</mo>}}{\phi }_h\right)\mathrm{,}\forall {\phi }_h\in U_{0h}\mathrm{,}$(3.4a)

$a\left(u_h\mathrm{,}{\Pi }_h^{\mathrm{<mo>*</mo>}}{\psi }_h\right)=\left(v_h\mathrm{,}{\Pi }_h^{\mathrm{<mo>*</mo>}}{\psi }_h\right)\mathrm{,}\forall {\psi }_h\in U_{0h}\mathrm{,}$ (3.4b)

$u_h\left(\cdot \mathrm{,}t\right)=g\left(\cdot \mathrm{,}t\right)\mathrm{,}t\in \left[-s\mathrm{,0}\right]\mathrm{,}$ (3.4c)

where $a\left(u_h,{\Pi }_h^{*}{\psi }_h\right)=-{\displaystyle {\sum }_{P_0\in {\overline{\Omega }}_h}}{\psi }_h\left(P_0\right){\displaystyle {\int }_{\partial K_{P_0}^{*}}}\frac{\partial u_h}{\partial n}\text{d}s$.

Denote $\mu ={\tau }^{\alpha }\Gamma \left(2-\alpha \right)$ . Applying Lemme 2.1, consider the completely discretization finite volume scheme as follows: find $\left(u_h^n\mathrm{,}{v}_h^n\right)\in U_{0h}\times U_{0h}$ , such that:

$\begin{array}{l}\left(u_h^n,{\Pi }_h^{*}{\phi }_h\right)+\mu a\left(v_h^n,{\Pi }_h^{*}{\phi }_h\right)\\ =a_{n-1}^{\left(\alpha \right)}\left(u_h^0,{\Pi }_h^{*}{\phi }_h\right)+{\displaystyle \sum _{k=1}^{n-1}}\left(a_{n-k-1}^{\left(\alpha \right)}-a_{n-k}^{\left(\alpha \right)}\right)\left(u_h^k,{\Pi }_h^{*}{\phi }_h\right)\\ +\mu \left({\overline{f}}_h^n,{\Pi }_h^{*}{\phi }_h\right),n=1,2,\cdots ,\forall {\phi }_h\in U_{0h},\end{array}$(3.5a)

$a\left(u_h^n\mathrm{,}{\Pi }_h^{\mathrm{<mo>*</mo>}}{\psi }_h\right)=\left(v_h^n\mathrm{,}{\Pi }_h^{\mathrm{<mo>*</mo>}}{\psi }_h\right)\mathrm{,}n=\mathrm{1,2,}\cdots \mathrm{,}\forall {\psi }_h\in U_{0h}\mathrm{,}$ (3.5b)

$u_h\left(x\mathrm{,}y\mathrm{,}t\right)=P_{h}g\left(x\mathrm{,}y\mathrm{,}t\right)\mathrm{,}t\in \left[-s\mathrm{,0}\right]\mathrm{,}$ (3.5c)

where $u_h^n$ and $v_h^n$ are approximation of $u\left(\cdot \mathrm{,}{t}_n\right)$ and $v\left(\cdot \mathrm{,}{t}_n\right)$ , respectively. The projection $P_h$ will be defined lately.

Lemma 3.1 [25]. The bilinear form $a\left(\cdot \mathrm{,}{\Pi }_h^{\mathrm{<mo>*</mo>}}\cdot \right)$ is symmetric and positive definite:

$a\left(u_h\mathrm{,}{\Pi }_h^{\mathrm{<mo>*</mo>}}{\omega }_h\right)=a\left({\omega }_h\mathrm{,}{\Pi }_h^{\mathrm{<mo>*</mo>}}{u}_h\right)\mathrm{,}\forall u_h\mathrm{,}{\omega }_h\in U_{0h}\mathrm{,}$

$a\left(u_h\mathrm{,}{\Pi }_h^{\mathrm{<mo>*</mo>}}{u}_h\right)={\left|u_h\right|}_1^2\mathrm{,}\forall u_h\in U_{0h}\mathrm{.}$

Lemma 3.2 [25]. There hold the following statements:

(i) $\left(u_h\mathrm{,}{\Pi }_h^{\mathrm{<mo>*</mo>}}{\overline{u}}_h\right)=\left({\overline{u}}_h\mathrm{,}{\Pi }_h^{\mathrm{<mo>*</mo>}}{u}_h\right)\mathrm{,}\forall u_h\mathrm{,}{\overline{u}}_h\in U_{0h}$.

(ii) Set $\left|\Vert u_h\Vert \right|={\left(u_h\mathrm{,}{\Pi }_h^{\mathrm{<mo>*</mo>}}{u}_h\right)}^{\frac{1}{2}}$ . Then $\left|\Vert \cdot \Vert \right|$ is equivalent to $\Vert \cdot \Vert$ on $U_{0h}$ , that is, there exist positive constants $c_1$ and $c_2$ such that

$c_1\Vert u_h\Vert \le \left|\Vert u_h\Vert \right|\le c_2\Vert u_h\Vert \mathrm{,}\forall u_h\in U_{0h}\mathrm{.}$

Applying Hölder inequality, it’s obvious that

$\left|\left(u_h\mathrm{,}{\Pi }_h^{\mathrm{<mo>*</mo>}}{\overline{u}}_h\right)\right|\le \left|\Vert u_h\Vert \right|\left|\Vert {\overline{u}}_h\Vert \right|\mathrm{,}\forall u_h\mathrm{,}{\overline{u}}_h\in U_{0h}\mathrm{.}$

Theorem 3.1. The finite volume scheme (3) is uniquely solvable.

proof. Since the finite volume scheme is linear, we can obtain the unique solvability by proving that the relevant homogeneous problem:

$\left(u_h^n\mathrm{,}{\Pi }_h^{\mathrm{<mo>*</mo>}}{\phi }_h\right)+\mu a\left(v_h^n\mathrm{,}{\Pi }_h^{\mathrm{<mo>*</mo>}}{\phi }_h\right)=\mathrm{0,}$ (3.6a)

$a\left(u_h^n\mathrm{,}{\Pi }_h^{\mathrm{<mo>*</mo>}}{\psi }_h\right)=\left(v_h^n\mathrm{,}{\Pi }_h^{\mathrm{<mo>*</mo>}}{\psi }_h\right)\mathrm{,}$ (3.6b)

admits solely trivial solution. Setting ${\phi }_h=u_h^n$ and ${\psi }_h=v_h^n$ in (3.6a) and (3.6b), respectively. Using Lemma 2.3, we have

$a\left(u_h^n\mathrm{,}{\Pi }_h^{\mathrm{<mo>*</mo>}}{v}_h^n\right)=a\left(v_h^n\mathrm{,}{\Pi }_h^{\mathrm{<mo>*</mo>}}{u}_h^n\right)\mathrm{.}$

Multiplying (3.6b) by $\mu$ , then we substract the resulting equation from (3.6a) to obtain

$\left(u_h^n\mathrm{,}{\Pi }_h^{\mathrm{<mo>*</mo>}}{u}_h^n\right)=-\mu \left(v_h^n\mathrm{,}{\Pi }_h^{\mathrm{<mo>*</mo>}}{v}_h^n\right)\mathrm{.}$

Applying Lemma 2.4, we arrive at $-\left(v_h^n\mathrm{,}{\Pi }_h^{\mathrm{<mo>*</mo>}}{v}_h^n\right)\le 0$ , then it can be easily obtained that

$\left(u_h^n\mathrm{,}{\Pi }_h^{\mathrm{<mo>*</mo>}}{u}_h^n\right)={\left|\Vert u_h^n\Vert \right|}^2=0.$

From the above equality, it is obvious that $u_h^n=v_h^n=0$ . Therefore we show that the solution of (3.6a)-(3.6b) is zero which implies that the scheme (3.5) is uniquely solvable. This proves the theorem.

4. Convergence and Stability Analysis

In this subsection, to analyze and discuss fully discrete a priori error results, we need to introduce an auxiliary projection $P_h\mathrm{<mo>:</mo>}{H}_0^1\left(\Omega \right)\cap H^2\left(\Omega \right)\mapsto U_{0h}$ defined by

$a\left(P_{h}u\mathrm{,}{\chi }_h\right)=a\left(u\mathrm{,}{\chi }_h\right)\mathrm{,}\forall {\chi }_h\in V_{0h}\mathrm{,}$ (4.1)

Lemma 4.1 [25]. Let $P_{h}u$ be the auxiliary projection of $u$ defined by (4.1) and $u\in W^{\mathrm{3,}p}\left(\Omega \right)\cap H_0^1\left(\Omega \right)$ then

$\Vert u-P_{h}u\Vert \le C{h}^2{\Vert u\Vert }_{3,p},\left(p>1\right).$

Theorem 4.1. Let $u\in C^2\left(\left[\mathrm{0,}T\right]\mathrm{,}{H}_0^1\left(\Omega \right)\cap H^2\left(\Omega \right)\right)$ be the solution of (1) and $u_h^n$ be the solution of the finite volume scheme (3) with $u_h^0=P_h{u}^0$ respectively, then the optimal error result in $L^2$ -norm hold

$\Vert u\left(t_n\right)-u_h^n\Vert \le C\left(h^2+{\tau }^{2-\alpha }\right)\mathrm{,}$ (4.2)

where $C$ is independent of $h$ and $\tau$ .

proof. To simplify the process of writing in the proof, we now split the errors as

$u^n-u_h^n=\left(u^n-P_h{u}^n\right)+\left(P_h{u}^n-u_h^n\right)={\rho }^n+{\eta }^n,$

$v^n-v_h^n=\left(v^n-P_h{v}^n\right)+\left(P_h{v}^n-v_h^n\right)={\xi }^n+{\theta }^n,$

$f^n-{\overline{f}}_h^n=\left(f^n-{\overline{f}}^n\right)+\left({\overline{f}}^n-{\overline{f}}_h^n\right)=E_1+\left({\overline{f}}^n-{\overline{f}}_h^n\right),$

${\sigma }^n={\partial }_t^{\left(\alpha \right)}{P}_h{u}^n-{}_0^C{D}_t^{\alpha }{u}^n\mathrm{.}$

Using (4.1), we note that

$a\left(\rho ,{\Pi }_h^{*}{\chi }_h\right)=0,a\left(\xi ,{\Pi }_h^{*}{\chi }_h\right)=0,\forall {\chi }_h\in U_h.$

Substracting the Equation (3.5) from (3.1), we obtain the error equations as follows

$\begin{array}{l}\left({}_0^c{D}_t^{\alpha }{u}^n-{\partial }_t^{\left(\alpha \right)}{u}_h^n\mathrm{,}{\Pi }_h^{\mathrm{<mo>*</mo>}}{\phi }_h\right)+a\left({\theta }^n\mathrm{,}{\Pi }_h^{\mathrm{<mo>*</mo>}}{\phi }_h\right)\\ =\left({\overline{f}}^n-{\overline{f}}_h^n\mathrm{,}{\Pi }_h^{\mathrm{<mo>*</mo>}}{\phi }_h\right)+\left(E_1\mathrm{,}{\Pi }_h^{\mathrm{<mo>*</mo>}}{\phi }_h\right)\mathrm{,}\forall {\phi }_h\in U_{0h}\mathrm{,}\end{array}$(4.3a)

$a\left({\eta }^n\mathrm{,}{\Pi }_h^{\mathrm{<mo>*</mo>}}{\psi }_h\right)=\left({\xi }^n+{\theta }^n\mathrm{,}{\Pi }_h^{\mathrm{<mo>*</mo>}}{\psi }_h\right)\mathrm{,}\forall {\psi }_h\in U_{0h}\mathrm{.}$ (4.3b)

Setting ${\phi }_h={\eta }^n$ and ${\psi }_h={\theta }^n$ , substracting the Equations (4.3b) from (4.3a), we arrive at

$\begin{array}{l}\left({\partial }_t^{\left(\alpha \right)}{\eta }^n\mathrm{,}{\Pi }_h^{\mathrm{<mo>*</mo>}}{\eta }^n\right)+\left({\theta }^n\mathrm{,}{\Pi }_h^{\mathrm{<mo>*</mo>}}{\theta }^n\right)\\ =\left({\sigma }^n\mathrm{,}{\Pi }_h^{\mathrm{<mo>*</mo>}}{\eta }^n\right)-\left({\xi }^n\mathrm{,}{\Pi }_h^{\mathrm{<mo>*</mo>}}{\theta }^n\right)+\left({\overline{f}}^n-{\overline{f}}_h^n\mathrm{,}{\Pi }_h^{\mathrm{<mo>*</mo>}}{\eta }^n\right)+\left(E_1\mathrm{,}{\Pi }_h^{\mathrm{<mo>*</mo>}}{\eta }^n\right)\mathrm{.}\end{array}$ (4.4)

based on the result $a\left({\theta }^n\mathrm{,}{\Pi }_h^{\mathrm{<mo>*</mo>}}{\eta }^n\right)=a\left({\eta }^n\mathrm{,}{\Pi }_h^{\mathrm{<mo>*</mo>}}{\theta }^n\right)$ obtained by Lemma 3.1.

Substituting the definition of ${\partial }_t^{\left(\alpha \right)}{\eta }^n$ and multiplying the equations by $\mu$ , we obtain that

$\begin{array}{l}\left({\eta }^n,{\Pi }_h^{*}{\eta }^n\right)+\mu \left({\theta }^n,{\Pi }_h^{*}{\theta }^n\right)\\ ={\displaystyle \sum _{k=1}^{n-1}}\left(a_{n-k-1}^{\left(\alpha \right)}-a_{n-k}^{\left(\alpha \right)}\right)\left({\eta }^k,{\Pi }_h^{*}{\eta }^n\right)+a_{n-1}^{\left(\alpha \right)}\left({\eta }^0,{\Pi }_h^{*}{\eta }^n\right)\\ -\mu \left({\xi }^n,{\Pi }_h^{*}{\theta }^n\right)+\mu \left({\sigma }^n,{\Pi }_h^{*}{\eta }^n\right)+\mu \left({\overline{f}}^n-{\overline{f}}_h^n,{\Pi }_h^{*}{\eta }^n\right)+\mu \left(E_1,{\Pi }_h^{*}{\eta }^n\right).\end{array}$(4.5)

Making use of Cauchy-Schwarz inequality and Lemma 3.2 and multiplying the equations by 4, we get

$\begin{array}{l}4{\left|\Vert {\eta }^n\Vert \right|}^2+4\mu {\left|\Vert {\theta }^n\Vert \right|}^2\\ \le 2{\displaystyle \sum _{k=1}^{n-1}}\left(a_{n-k-1}^{\left(\alpha \right)}-a_{n-k}^{\left(\alpha \right)}\right)\left({\left|\Vert {\eta }^k\Vert \right|}^2+{\left|\Vert {\eta }^n\Vert \right|}^2\right)+2a_{n-1}^{\left(\alpha \right)}\left({\left|\Vert {\eta }^0\Vert \right|}^2+{\left|\Vert {\eta }^n\Vert \right|}^2\right)\\ +2\mu \left({\left|\Vert {\xi }^n\Vert \right|}^2+{\left|\Vert {\theta }^n\Vert \right|}^2\right)+4\mu \left({\sigma }^n,{\Pi }_h^{*}{\eta }^n\right)+4\mu \left({\overline{f}}^n-{\overline{f}}_h^n,{\Pi }_h^{*}{\eta }^n\right)\\ +4\mu \left(E_1,{\Pi }_h^{*}{\eta }^n\right).\end{array}$(4.6)

From Lemma 2.1, we see that

${\displaystyle \sum _{k=1}^{n-1}}\left(a_{n-k-1}^{\left(\alpha \right)}-a_{n-k}^{\left(\alpha \right)}\right)=1-a_{n-1}^{\left(\alpha \right)}.$

Hence we can rewrite (4.6) as follows

$\begin{array}{c}2{\left|\Vert {\eta }^n\Vert \right|}^2\le 2{\displaystyle \sum _{k=1}^{n-1}}\left(a_{n-k-1}^{\left(\alpha \right)}-a_{n-k}^{\left(\alpha \right)}\right){\left|\Vert {\eta }^k\Vert \right|}^2+2a_{n-1}^{\left(\alpha \right)}{\left|\Vert {\eta }^0\Vert \right|}^2+2\mu {\left|\Vert {\xi }^n\Vert \right|}^2\\ +4\mu \left({\sigma }^n,{\Pi }_h^{*}{\eta }^n\right)+4\mu \left({\overline{f}}^n-{\overline{f}}_h^n,{\Pi }_h^{*}{\eta }^n\right)+4\mu \left(E_1,{\Pi }_h^{*}{\eta }^n\right).\end{array}$(4.7)

Choosing $\epsilon =\frac{3}{4}$ , using $\epsilon$ -Cauchy inequality and (4.7), we obtain

$\begin{array}{l}{\left|\Vert {\eta }^n\Vert \right|}^2\le 2{\displaystyle \sum _{k=1}^{n-1}}\left(a_{n-k-1}^{\left(\alpha \right)}-a_{n-k}^{\left(\alpha \right)}\right){\left|\Vert {\eta }^k\Vert \right|}^2+2a_{n-1}^{\left(\alpha \right)}{\left|\Vert {\eta }^0\Vert \right|}^2+2\mu {\left|\Vert {\xi }^n\Vert \right|}^2\\ +3\mu \left({\left|\Vert {\sigma }^n\Vert \right|}^2+{\left|\Vert {\overline{f}}^n-{\overline{f}}_h^n\Vert \right|}^2+{\left|\Vert E_1\Vert \right|}^2\right)\mathrm{.}\end{array}$(4.8)

Applying Lemma 3.2 yields

$\begin{array}{l}{\Vert {\eta }^n\Vert }^2\le 2{\displaystyle \sum _{k=1}^{n-1}}\left(a_{n-k-1}^{\left(\alpha \right)}-a_{n-k}^{\left(\alpha \right)}\right){\Vert {\eta }^k\Vert }^2+2a_{n-1}^{\left(\alpha \right)}{\Vert {\eta }^0\Vert }^2+2\mu {\Vert {\xi }^n\Vert }^2\\ +3\mu \left({\Vert {\sigma }^n\Vert }^2+{\Vert {\overline{f}}^n-{\overline{f}}_h^n\Vert }^2+{\Vert E_1\Vert }^2\right),\\ :=2{\displaystyle \sum _{k=1}^{n-1}}\left(a_{n-k-1}^{\left(\alpha \right)}-a_{n-k}^{\left(\alpha \right)}\right){\Vert {\eta }^k\Vert }^2+e_1+e_2+e_3+e_4+e_5.\end{array}$(4.9)

For the next process, now we need to consider $e_1$ - $e_5$ .

Considering $e_1$ , based on the initial condition, it’s clear that

${\eta }^0=P_h{u}^0-u_h^0=0.$ (4.10)

For $e_2$ , we can easily get the following inequality by Lemma 4.1,

$\begin{array}{l}\Vert {\xi }^n\Vert \le C{h}^2\underset{0\le t\le t_n}{\max }{\Vert v\Vert }_{\mathrm{3,}p}\mathrm{.}\hfill \end{array}$ (4.11)

For $e_3$ , we rewrite ${\sigma }^n$ as the following form

$\begin{array}{l}{\sigma }^n=\left({\partial }_t^{\left(\alpha \right)}{P}_h{u}^n-{\partial }_t^{\left(\alpha \right)}{u}^n\right)+\left({\partial }_t^{\left(\alpha \right)}{u}^n-{}_0^c{D}_t^{\alpha }{u}^n\right)\\ \mathrm{<mo>:</mo>}={\sigma }_1^n+{\sigma }_2^n\mathrm{.}\end{array}$

Using Lemma 2.1, we obtain

$\begin{array}{c}\left|{\sigma }_1^n\right|=\left|\frac{{\tau }^{-\alpha }}{\Gamma \left(2-\alpha \right)}{\displaystyle \sum _{k=1}^n}{a}_{n-k}^{\left(\alpha \right)}{\displaystyle {\int }_{t_{k-1}}^{t_k}}\left(P_h-I\right)u_t\text{d}t\right|\\ \le \frac{{\tau }^{-\alpha }}{\Gamma \left(2-\alpha \right)}{\displaystyle \sum _{k=1}^n}{a}_{n-k}^{\left(\alpha \right)}\tau \underset{0\le t\le t_n}{\max }\left|\left(P_h-I\right)u_t\right|\\ \le C{\tau }^{1-\alpha }{\displaystyle \sum _{k=1}^n}{a}_{n-k}^{\left(\alpha \right)}\underset{0\le t\le t_n}{\max }\left|\left(P_h-I\right)u_t\right|\\ =C{\left(n\tau \right)}^{1-\alpha }\underset{0\le t\le t_m}{\max }\left|\left(P_h-I\right)u_t\right|,\end{array}$

and

$\left|{\sigma }_2^n\right|\le \frac{1}{2\Gamma \left(1-\alpha \right)}\left[\frac{1}{4}+\frac{\alpha }{\left(1-\alpha \right)\left(2-\alpha \right)}\right]\underset{0\le t\le t_n}{\max }\left|u_{tt}\right|{\tau }^{2-\alpha }\mathrm{.}$

Applying Lemma 4.1, it yields

$\Vert {\sigma }_1^n\Vert ={\left({\displaystyle {\int }_{\Omega }}{\left|{\sigma }_1^n\right|}^2\text{d}x\text{d}y\right)}^{\frac{1}{2}}\le C{\left(n\tau \right)}^{1-\alpha }\underset{0\le t\le t_n}{\max }\Vert \left(P_h-I\right)u_t\Vert \le C{h}^2\underset{0\le t\le t_n}{\max }{\Vert u_t\Vert }_{\mathrm{3,}p}\mathrm{,}$

and

$\Vert {\sigma }_2^n\Vert ={\left({\displaystyle {\int }_{\Omega }}{\left|{\sigma }_2^n\right|}^2\text{d}x\text{d}y\right)}^{\frac{1}{2}}\le C{\tau }^{2-\alpha }\underset{0\le t\le t_n}{\max }\Vert u_{tt}\Vert \mathrm{.}$

Then we use the triangle inequality to derive

$\begin{array}{l}{\Vert {\sigma }^n\Vert }^2\le 2\left({\Vert {\sigma }_1^n\Vert }^2+{\Vert {\sigma }_2^n\Vert }^2\right)\le C{h}^4\underset{0\le t\le t_n}{\max }{\Vert u_t\Vert }_{\mathrm{3,}p}^2+C{\tau }^{4-2\alpha }\underset{0\le t\le t_n}{\max }{\Vert u_{tt}\Vert }^2\mathrm{.}\hfill \end{array}$ (4.12)

For $e_4$ , it follows from the definition of ${\widehat{u}}^n$ :

$\begin{array}{c}\left|{\widehat{u}}^n-{\widehat{u}}_h^n\right|=\left|\left(2u^{n-1}-u^{n-2}\right)-\left(2u_h^{n-1}-u_h^{n-2}\right)\right|\\ \le 2\left|u^{n-1}-u_h^{n-1}\right|+\left|u^{n-1}-u_h^{n-2}\right|\mathrm{.}\end{array}$

Similarly,

$\begin{array}{l}\left|{\tilde{u}}^{n-m}-{\tilde{u}}_h^{n-m}\right|\\ \le \{\begin{array}{l}\left|u^{n-m}-u_h^{n-m}\right|\mathrm{,}\text{if}m\in N_{+}\\ \left(2-m\right)\left|u^{n-1}-u_h^{n-1}\right|+\left(1-m\right)\left|u^{n-2}-u_h^{n-2}\right|\mathrm{,}\text{if}m\in \left(\mathrm{0,1}\right)\\ \left(\lceil m\rceil -m\right)\left|u^{k-\lfloor m\rfloor }-u_h^{k-\lfloor m\rfloor }\right|+\left(m-\lfloor m\rfloor \right)\left|u^{k-\lceil m\rceil }-u_h^{k-\lceil m\rceil }\right|\mathrm{,}\text{if}m>1\text{and}m\notin N_{+}\mathrm{.}\end{array}\end{array}$

Noting that $m-\lfloor m\rfloor <2$ and $\lceil m\rceil -m<2$ , by the triangle inequality and Taylor’s series, we can check that the following inequalities hold :

Case 1: if $m>1$ and $m\notin N_{+}$

$\begin{array}{c}\Vert {\overline{f}}^n-{\overline{f}}_h^n\Vert =\Vert {\overline{c}}_1\left({\widehat{u}}^n-{\widehat{u}}_h^n\right)+{\overline{c}}_2\left({\tilde{u}}^{n-m}-{\tilde{u}}_h^{n-m}\right)\Vert \\ \le 2{\overline{c}}_3\left(\Vert u^{n-1}-u_h^{n-1}\Vert +\Vert u^{n-2}-u_h^{n-2}\Vert +\Vert u^{n-\lceil m\rceil }-u_h^{n-\lceil m\rceil }\Vert +\Vert u^{n-\lfloor m\rfloor }-u_h^{n-\lfloor m\rfloor }\Vert \right)\\ \le 2{\overline{c}}_3(\Vert {\rho }^{n-1}\Vert +\Vert {\rho }^{n-2}\Vert +\Vert {\rho }^{n-\lceil m\rceil }\Vert +\Vert {\rho }^{n-\lfloor m\rfloor }\Vert \\ +\Vert {\eta }^{n-1}\Vert +\Vert {\eta }^{n-2}\Vert +\Vert {\eta }^{n-\lceil m\rceil }\Vert +\Vert {\eta }^{n-\lfloor m\rfloor }\Vert ).\end{array}$

Thus

$\begin{array}{c}{\Vert {\overline{f}}^n-{\overline{f}}_h^n\Vert }^2\le 16{\overline{c}}_3{({\Vert {\rho }^{n-1}\Vert }^2+{\Vert {\rho }^{n-2}\Vert }^2+{\Vert {\rho }^{n-\lceil m\rceil }\Vert }^2+\Vert {\rho }^{n-\lfloor m\rfloor }\Vert }^2\\ +{\Vert {\eta }^{n-1}\Vert }^2+{\Vert {\eta }^{n-2}\Vert }^2+{\Vert {\eta }^{n-\lceil m\rceil }\Vert }^2+{\Vert {\eta }^{n-\lfloor m\rfloor }\Vert }^2)\mathrm{<mo>;</mo>}\end{array}$

Case 2: if $m\in \left(\mathrm{0,1}\right)$

${\Vert {\overline{f}}^n-{\overline{f}}_h^n\Vert }^2\le 16{\overline{c}}_3\left({\Vert {\rho }^{n-1}\Vert }^2+{\Vert {\rho }^{n-2}\Vert }^2+{\Vert {\eta }^{n-1}\Vert }^2+{\Vert {\eta }^{n-2}\Vert }^2\right)\mathrm{<mo>;</mo>}$

Case 3: if $m\in N_{+}$

${\Vert {\overline{f}}^n-{\overline{f}}_h^n\Vert }^2\le 16{\overline{c}}_3\left({\Vert {\rho }^{n-1}\Vert }^2+{\Vert {\rho }^{n-2}\Vert }^2+{\Vert {\rho }^{n-m}\Vert }^2+{\Vert {\eta }^{n-1}\Vert }^2+{\Vert {\eta }^{n-2}\Vert }^2+{\Vert {\eta }^{n-m}\Vert }^2\right)\mathrm{.}$

Applying Lemma 4.1, we obtain that

$\Vert {\rho }^k\Vert \le C{h}^2\underset{0\le t\le t_k}{\max }{\Vert u\Vert }_{\mathrm{3,}p}\mathrm{,}$

for $k=n-\mathrm{1,}n-\mathrm{2,}n-\lceil m\rceil \mathrm{,}n-\lfloor m\rfloor$ .

Based on the above derivation, we reach that

$\begin{array}{l}{\Vert {\overline{f}}^n-{\overline{f}}_h^n\Vert }^2\\ \le C{h}^4\underset{0\le t\le t_n}{\max }{\Vert u\Vert }_{\mathrm{3,}p}^2+16{\overline{c}}_3\left({\Vert {\eta }^{n-1}\Vert }^2+{\Vert {\eta }^{n-2}\Vert }^2+{\Vert {\eta }^{n-m}\Vert }^2+{\Vert {\eta }^{n-\lceil m\rceil }\Vert }^2+{\Vert {\eta }^{n-\lfloor m\rfloor }\Vert }^2\right)\mathrm{.}\end{array}$ (4.13)

For $e_5$ , we have

${\Vert E_1\Vert }^2\le C{\tau }^4\mathrm{.}$ (4.14)

Substitute the results for $e_1$ - $e_5$ into (4.9) to obtain

$\begin{array}{c}{\Vert {\eta }^n\Vert }^2\le 2{\displaystyle \sum _{k=1}^{n-1}}\left(a_{n-k-1}^{\left(\alpha \right)}-a_{n-k}^{\left(\alpha \right)}\right){\Vert {\eta }^k\Vert }^2\\ +48{\overline{c}}_3{\tau }^{\alpha }\left({\Vert {\eta }^{n-1}\Vert }^2+{\Vert {\eta }^{n-2}\Vert }^2+{\Vert {\eta }^{n-m}\Vert }^2+{\Vert {\eta }^{n-\lceil m\rceil }\Vert }^2+{\Vert {\eta }^{n-\lfloor m\rfloor }\Vert }^2\right)\\ +C{\tau }^{\alpha }(h^4\underset{0\le t\le t_n}{\max }{\Vert v\Vert }_{3,p}^2+h^4\underset{0\le t\le t_n}{\max }{\Vert u\Vert }_{3,p}^2+h^4\underset{0\le t\le t_n}{\max }{\Vert u_t\Vert }_{3,p}^2\\ +{\tau }^{4-2\alpha }\underset{0\le t\le t_n}{\max }{\Vert u_{tt}\Vert }^2+{\tau }^4).\end{array}$(4.15)

Using Lemma 2.2, we have

$\begin{array}{c}{\Vert {\eta }^n\Vert }^2\le C\left(h^4+{\tau }^{4-2\alpha }\right)\exp \left(2{\displaystyle \sum _{k=1}^{n-1}}\left(a_{n-k-1}^{\left(\alpha \right)}-a_{n-k}^{\left(\alpha \right)}\right)+240{\overline{c}}_3\right)\\ =C\left(h^4+{\tau }^{4-2\alpha }\right)\exp \left(2\left(1-a_{n-1}^{\left(\alpha \right)}\right)+240{\overline{c}}_3\right)\\ \le C\left(h^4+{\tau }^{4-2\alpha }\right)\exp \left(2+240{\overline{c}}_3\right).\end{array}$ (4.16)

Denoting $\overline{C}={\left[C\exp \left(2+240{\overline{c}}_3\right)\right]}^{\frac{1}{2}}$ , we obtain that

$\begin{array}{l}\Vert {\eta }^n\Vert \le \overline{C}\left(h^2+{\tau }^{2-\alpha }\right)\mathrm{.}\hfill \end{array}$ (4.17)

Applying Lemma 4.1, we conclude that

$\Vert u\left(t_n\right)-u_h^n\Vert \le \Vert {\rho }^n\Vert +\Vert {\eta }^n\Vert \le C\left(h^2+{\tau }^{2-\alpha }\right)\mathrm{.}$ (4.18)

Thus the proof of the theorem is completed.

Next, we analyze the numerical stability of the finite volume scheme (3.5). The numerical stability means that a small perturbation of the initial value implies a small perturbation of the numerical solution.

Theorem 4.2. we suppose that $\left(U^n\mathrm{,}{V}^n\right)$ is the solution of perturbation equation and $\iota$ is a small perturbation of $g$ , denoting ${\lambda }_h^k=u_h^k-U^k$ , ${\gamma }_h^k=v_h^k-V^k$ , the following stability result hold

${\Vert {\lambda }_h^n\Vert }^2+\mu {\Vert {\gamma }_h^n\Vert }^2\le C\underset{-m\le k\le 0}{\max }{\Vert {\iota }^k\Vert }^2\mathrm{.}$ (4.19)

This theorem can be proved by using the same way as the proof of Theorem 4.1.

5. Numerical Experiments

We now present some numerical experiments to verify our theoretical statements. For the purpose of manifesting the stability and convergence rate of the proposed scheme, we first consider an example in which the exact solutions are known. Then in the second example, we consider a more realistic problem for which the exact solution is not given beforehand. In addition, we choose $\Omega =\left(\mathrm{0,1}\right)\times \left(\mathrm{0,1}\right)$ and $T=1$ for all examples.

Example 1. Choose the exact solution for the problem (1.1) to be

$u=t^3\sin \left(\pi x\right)\sin \left(\pi y\right)\mathrm{,}\left(x\mathrm{,}y\mathrm{,}t\right)\in \overline{\Omega }\times \left[-s\mathrm{,}T\right]\mathrm{.}$

associated with the following source function

$\begin{array}{c}f=u\left(x\mathrm{,}y\mathrm{,}t-s\right)+u^2\left(x\mathrm{,}y\mathrm{,}t\right)\\ +\left[\frac{\Gamma \left(4\right)}{\Gamma \left(4-\alpha \right)}{t}^{3-\alpha }-{\left(t-s\right)}^3+4{\pi }^4{t}^3\right]\sin \left(\pi x\right)\sin \left(\pi y\right)\\ -t^6{\sin }^2\left(\pi x\right){\sin }^2\left(\pi y\right)\end{array}$

In Table 1 and Table 2, the error in $L^2$ -norms are listed with delay parameter $s=0.2$ . To testify the convergence order in spatial direction, we keep varying $h$ and $\tau =h^2$ . Similarly, to testify the convergence order in temporal direction, we keep varying $\tau$ and $h=\tau$ . From the numerical results we can see that the scheme (3) is stable and has the convergence rate of $O\left(h^2+{\tau }^{2-\alpha }\right)$ . In order to present a comparison with our scheme, this example is also numerically solved by a central difference scheme and the results of it are listed in Table 3 and Table 4. By comparison, we can find that our scheme is more advantageous in accuracy, the comparison between the two methods is presented in Figure 2. Furthermore, The numerical solutions and exact solutions are plotted in Figure 3.

Table 1. Errors and spatial convergence orders of finite volume scheme with $h\mathrm{<mo>=</mo>}\sqrt{\tau }$ at $T=1$ .

s	$\alpha$	h	$\Vert u-u_h\Vert$	order
0.2	0.3	$\frac{1}{10}$	3.6198E−03
		$\frac{1}{20}$	9.0193E−04	2.00
		$\frac{1}{30}$	4.0061E−04	2.00
		$\frac{1}{40}$	2.2529E−04	2.00
	0.5	$\frac{1}{10}$	3.6134E−03
		$\frac{1}{20}$	9.0018E−04	2.01
		$\frac{1}{30}$	3.9981E−04	2.00
		$\frac{1}{40}$	2.2485E−04	2.00
	0.7	$\frac{1}{10}$	3.6089E−03
		$\frac{1}{20}$	8.9874E−04	2.01
		$\frac{1}{30}$	3.9910E−04	2.00
		$\frac{1}{40}$	2.2449E−04	2.00

Table 2. Errors and time convergence orders of finite volume scheme with $h=\tau$ at $T=1$ .

s	$\alpha$	$\tau$	$\Vert u-u_h\Vert$	order
0.2	0.3	$\frac{1}{10}$	3.5423E−03
		$\frac{1}{20}$	8.8172E−04	2.01
		$\frac{1}{30}$	3.9178E−04	2.00
		$\frac{1}{40}$	2.2045E−04	2.00
	0.5	$\frac{1}{10}$	3.5672E−03
		$\frac{1}{20}$	8.9237E−04	2.00
		$\frac{1}{30}$	3.9809E−04	1.99
		$\frac{1}{40}$	2.2476E−04	1.99
	0.9	$\frac{1}{10}$	3.7612E−03
		$\frac{1}{20}$	9.9265E−04	1.92
		$\frac{1}{30}$	4.6512E−04	1.87
		$\frac{1}{40}$	2.7484E−04	1.83

Table 3. Errors and spatial convergence orders of central difference scheme with $h=\sqrt{\tau }$ at $T=1$ .

s	$\alpha$	h	$\Vert u-u_h\Vert$	order
0.2	0.3	$\frac{1}{10}$	8.3102E−03
		$\frac{1}{20}$	2.0635E−03	2.01
		$\frac{1}{30}$	9.1595E−04	2.00
		$\frac{1}{40}$	5.1498E−04	2.00
	0.5	$\frac{1}{10}$	8.3036E−03
		$\frac{1}{20}$	2.0618E−03	2.01
		$\frac{1}{30}$	9.1516E−04	2.00
		$\frac{1}{40}$	5.1456E−04	2.00
	0.7	$\frac{1}{10}$	8.2989E−03
		$\frac{1}{20}$	2.0603E−03	2.01
		$\frac{1}{30}$	9.1444E−04	2.00
		$\frac{1}{40}$	5.1419E−04	2.00

Table 4. Errors and time convergence orders of central difference scheme with $h=\tau$ at $T=1$ .

s	$\alpha$	$\tau$	$\Vert u-u_h\Vert$	order
0.2	0.3	$\frac{1}{10}$	8.2269E−03
		$\frac{1}{20}$	2.0429E−03	2.01
		$\frac{1}{30}$	9.0707E−04	2.00
		$\frac{1}{40}$	5.1014E−04	2.00
	0.5	$\frac{1}{10}$	8.2519E−03
		$\frac{1}{20}$	2.0536E−03	2.01
		$\frac{1}{30}$	9.1338E−04	2.00
		$\frac{1}{40}$	5.1445E−04	2.00
	0.9	$\frac{1}{10}$	8.4475E−03
		$\frac{1}{20}$	2.1541E−03	1.97
		$\frac{1}{30}$	9.8047E−04	1.94
		$\frac{1}{40}$	5.6455E−04	1.92

Figure 2. Error comparison between the two methods when $\alpha =0.3$ . (a) Spatial errors for $u_h$ ; (b) Temporal errors for $u_h$ .

Figure 3. Numerical solution and exact solution in Example 1 at $T=1$ with $\alpha =0.3$ , $\tau =h^2=\frac{1}{400}$ . (a) Numerical solution $u_h$ ; (b) Numerical solution $v_h$ ; (c) Exact solution $u$ ; (d) Exact solution $v$ .

Example 2. We consider the problem (1.1) with the following initial and boundary conditions

$u\left(x\mathrm{,}y\mathrm{,}t\right)=\mathrm{0,}\left(x\mathrm{,}y\mathrm{,}t\right)\in \partial \Omega \times \left[\mathrm{0,}T\right]\mathrm{,}$

$\Delta u\left(x\mathrm{,}y\mathrm{,}t\right)=\mathrm{0,}\left(x\mathrm{,}y\mathrm{,}t\right)\in \partial \Omega \times \left[\mathrm{0,}T\right]\mathrm{,}$

$u\left(x\mathrm{,}y\mathrm{,}t\right)=2\left(1-t^4\right)\sin \left(5\pi x\right)\sin \left(5\pi y\right)\mathrm{,}\left(x\mathrm{,}y\mathrm{,}t\right)\in \overline{\Omega }\times \left[-s\mathrm{,0}\right]\mathrm{,}$

$\Delta u\left(x\mathrm{,}y\mathrm{,}t\right)=100{\pi }^2\left(1-t^4\right)\sin \left(5\pi x\right)\sin \left(5\pi y\right)\mathrm{,}\left(x\mathrm{,}y\mathrm{,}t\right)\in \overline{\Omega }\times \left[-s\mathrm{,0}\right]\mathrm{,}$

and source function

$f\left(x\mathrm{,}y\mathrm{,}t\right)=u^2\left(x\mathrm{,}y\mathrm{,}t-s\right)+u\left(x\mathrm{,}y\mathrm{,}t\right)\mathrm{.}$

We choose the numerical solution with $h=\frac{1}{120}$ and $\tau =h^2$ as the approximating exact solution.

In Table 5 and Table 6, the error in $L^2$ -norms are listed with delay parameter

Table 5. Errors and spatial convergence orders of finite volume scheme with $h=\sqrt{\tau }$ at $T=1$ .

s	$\alpha$	h	$\Vert u-u_h\Vert$	order
0.3	0.2	$\frac{1}{10}$	1.5841E−06
		$\frac{1}{20}$	6.8489E−07	1.21
		$\frac{1}{30}$	3.2876E−07	1.81
		$\frac{1}{40}$	1.8257E−07	2.04
0.5	0.2	$\frac{1}{10}$	1.1614E−06
		$\frac{1}{20}$	5.0213E−07	1.21
		$\frac{1}{30}$	2.4103E−07	1.81
		$\frac{1}{40}$	1.3385E−07	2.04
0.7	0.2	$\frac{1}{10}$	6.8891E−07
		$\frac{1}{20}$	2.9784E−07	1.21
		$\frac{1}{30}$	1.4297E−07	1.81
		$\frac{1}{40}$	7.9394E−08	2.04

Table 6. Errors and time convergence orders of finite volume scheme with $h=\tau$ at $T=1$ .

s	$\alpha$	$\tau$	$\Vert u-u_h\Vert$	order
0.3	0.2	$\frac{1}{10}$	1.5708E−06
		$\frac{1}{20}$	6.7389E−07	1.22
		$\frac{1}{30}$	3.2117E−07	1.83
		$\frac{1}{40}$	1.7739E−07	2.06
0.5	0.2	$\frac{1}{10}$	1.1452E−06
		$\frac{1}{20}$	4.8876E−07	1.23
		$\frac{1}{30}$	2.3183E−07	1.84
		$\frac{1}{40}$	1.2759E−07	2.08
0.7	0.2	$\frac{1}{10}$	6.7562E−07
		$\frac{1}{20}$	2.8683E−07	1.24
		$\frac{1}{30}$	1.3542E−07	1.85
		$\frac{1}{40}$	7.4271E−08	2.09

Figure 4. Numerical solution in Example 2 with $\alpha =0.3$ , $s=0.2$ . (a) t = 0.2; (b) t = 0.4; (c) t = 0.6; (d) t = 0.8.

Figure 5. Numerical solution in Example 2 with $\alpha =0.3$ , $s=0.5$ . (a) t = 0.2; (b) t = 0.4; (c) t = 0.6; (d) t = 0.8.

$s=0.2$ . To testify the convergence order in spatial direction, we keep varying $h$ and $\tau =h^2$ . Similarly, to testify the convergence order in temporal direction, we keep varying $\tau$ and $h=\tau$ . From tables we can see that numerical results are consistent with the theoretical results. In Figure 4 and Figure 5, the numerical solution with the time evolution are plotted when the delay s = 0.2, 0.5, respectively. It indicate that the delay effect on the behavior of the numerical solution.

6. Conclusion

In this article, a finite volume element scheme, which can achieve the convergence rate of $O\left(h^2+{\tau }^{2-\alpha }\right)$ , has been derived for the two-dimensional time-fractional nonlinear fourth-order diffusion equation with time delay. The stability and convergence analyses of our scheme are proved by Gronwall lemma. Then, the numerical experiments are given to verify the effectiveness of the proposed scheme.

Conflicts of Interest

The authors declare no conflicts of interest regarding the publication of this paper.

References

[1]	Arenas, A.J., González-Parra, G. and Chen-Charpentier, B.M. (2016) Construction of Nonstandard Finite Difference Schemes for the SI and SIR Epidemic Models of Fractional Order. Mathematics and Computers in Simulation, 121, 48-63. https://doi.org/10.1016/j.matcom.2015.09.001
[2]	Machado, J.A.T. (2015) Fractional Order Description of DNA. Applied Mathematical Modelling, 39, 4095-4102. https://doi.org/10.1016/j.apm.2014.12.037
[3]	Debnath, L. (2003) Recent Applications of Fractional Calculus to Science and Engineering. International Journal of Mathematics and Mathematical Sciences, 2003, 3413-3442. https://doi.org/10.1155/s0161171203301486
[4]	Kumar, S. (2014) A New Analytical Modelling for Fractional Telegraph Equation via Laplace Transform. Applied Mathematical Modelling, 38, 3154-3163. https://doi.org/10.1016/j.apm.2013.11.035
[5]	Lotfy, K. (2017) A Novel Solution of Fractional Order Heat Equation for Photothermal Waves in a Semiconductor Medium with a Spherical Cavity. Chaos, Solitons & Fractals, 99, 233-242. https://doi.org/10.1016/j.chaos.2017.04.017
[6]	Nikan, O., Jafari, H. and Golbabai, A. (2020) Numerical Analysis of the Fractional Evolution Model for Heat Flow in Materials with Memory. Alexandria Engineering Journal, 59, 2627-2637. https://doi.org/10.1016/j.aej.2020.04.026
[7]	Tuan, N.H., Aghdam, Y.E., Jafari, H. and Mesgarani, H. (2020) A Novel Numerical Manner for Two‐dimensional Space Fractional Diffusion Equation Arising in Transport Phenomena. Numerical Methods for Partial Differential Equations, 37, 1397-1406. https://doi.org/10.1002/num.22586
[8]	Beta, C., Bertram, M., Mikhailov, A.S., Rotermund, H.H. and Ertl, G. (2003) Controlling Turbulence in a Surface Chemical Reaction by Time-Delay Autosynchronization. Physical Review E, 67, Article ID: 046224. https://doi.org/10.1103/physreve.67.046224
[9]	Benchohra, M. and Berhoun, F. (2010) Impulsive Fractional Differential Equations with State-Dependent Delay. Communications in Applied Analysis, 14, 213.
[10]	Li, D., Zhang, C. and Wang, W. (2012) Long Time Behavior of Non-Fickian Delay Reaction-Diffusion Equations. Nonlinear Analysis: Real World Applications, 13, 1401-1415. https://doi.org/10.1016/j.nonrwa.2011.11.005
[11]	Tenreiro Machado, J.A. and Lopes, A.M. (2017) A Fractional Perspective on the Trajectory Control of Redundant and Hyper-Redundant Robot Manipulators. Applied Mathematical Modelling, 46, 716-726. https://doi.org/10.1016/j.apm.2016.11.005
[12]	Boukal, Y., Zasadzinski, M., Darouach, M. and Radhy, N.E. (2016) Robust Functional Observer Design for Uncertain Fractional-Order Time-Varying Delay Systems. 2016 American Control Conference (ACC), Boston, 6-8 July 2016, 2741-2746. https://doi.org/10.1109/acc.2016.7525333
[13]	Rezounenko, A.V. and Wu, J. (2006) A Non-Local PDE Model for Population Dynamics with State-Selective Delay: Local Theory and Global Attractors. Journal of Computational and Applied Mathematics, 190, 99-113. https://doi.org/10.1016/j.cam.2005.01.047
[14]	Danumjaya, P. and Pani, A.K. (2011) Mixed Finite Element Methods for a Fourth Order Reaction Diffusion Equation. Numerical Methods for Partial Differential Equations, 28, 1227-1251. https://doi.org/10.1002/num.20679
[15]	Liu, Y., Du, Y., Li, H., He, S. and Gao, W. (2015) Finite Difference/Finite Element Method for a Nonlinear Time-Fractional Fourth-Order Reaction-Diffusion Problem. Computers & Mathematics with Applications, 70, 573-591. https://doi.org/10.1016/j.camwa.2015.05.015
[16]	Zhang, T. and Guo, Q. (2018) The Finite Difference/Finite Volume Method for Solving the Fractional Diffusion Equation. Journal of Computational Physics, 375, 120-134. https://doi.org/10.1016/j.jcp.2018.08.033
[17]	Liu, X. and Yang, X. (2021) Mixed Finite Element Method for the Nonlinear Time-Fractional Stochastic Fourth-Order Reaction-Diffusion Equation. Computers & Mathematics with Applications, 84, 39-55. https://doi.org/10.1016/j.camwa.2020.12.004
[18]	Zhang, Q., Ran, M. and Xu, D. (2016) Analysis of the Compact Difference Scheme for the Semilinear Fractional Partial Differential Equation with Time Delay. Applicable Analysis, 96, 1867-1884. https://doi.org/10.1080/00036811.2016.1197914
[19]	Pimenov, V.G., Hendy, A.S. and De Staelen, R.H. (2017) On a Class of Non-Linear Delay Distributed Order Fractional Diffusion Equations. Journal of Computational and Applied Mathematics, 318, 433-443. https://doi.org/10.1016/j.cam.2016.02.039
[20]	Meis, T. and Marcowitz, U. (2012) Numerical Solution of Partial Differential Equations, volume 32. Springer.
[21]	Ferreira, J.A. and DaSilva, P.M. (2008) Energy Estimates for Delay Diffusion-Reaction Equations. Journal of Computational Mathematics, 26, 536-553.
[22]	Nandal, S. and Narain Pandey, D. (2020) Numerical Solution of Non-Linear Fourth Order Fractional Sub-Diffusion Wave Equation with Time Delay. Applied Mathematics and Computation, 369, Article ID: 124900. https://doi.org/10.1016/j.amc.2019.124900
[23]	Xie, H. and Yang, Q. (2022) Compact Difference Scheme for Time-Fractional Nonlinear Fourth-Order Diffusion Equation with Time Delay. Results in Applied Mathematics, 16, Article ID: 100339. https://doi.org/10.1016/j.rinam.2022.100339
[24]	Holte, J.M. (2009) Discrete Gronwall Lemma and Applications. MAA-NCS Meeting 2009, New York, 1-7.
[25]	Li, R.H., Chen, Z.Y. and Wu, W. (2000) Generalized Difference Methods for Differential Equations: Numerical Analysis of Finite Volume Methods. CRC Press.