Numerical Analysis of Alternating Direction Implicit Orthogonal Spline Collocation Scheme for the Hyperbolic Integrodifferential Equation with a Weakly Singular Kernel

Abstract

This paper studies an alternating direction implicit orthogonal spline collocation (ADIOSC) technique for calculating the numerical solution of the hyperbolic integrodifferential problem with a weakly singular kernel in the two-dimensional domain. The integral term is approximated with the help of the second-order fractional quadrature formula introduced by Lubich. The stability and convergence analysis of the proposed strategy are proven in $L^2$-norm. Numerical results highlight the high accuracy and efficiency of the proposed strategy and clarify the theoretical prediction.

1. Introduction

In this paper, we study the hyperbolic integrodifferential equation with a weakly singular kernel [1] as:

$$\begin{array}{c}\hfill u_{tt}(x,y,t)-\mu \Delta u(x,y,t)={\int }_0^t\beta (t-s)\Delta u(x,y,s)ds+f(x,y,t),t\in (0,T],(x,y)\in \Omega .\end{array}$$

(1)

The initial conditions (ICs) and boundary condition (BC) are prescribed as:

$$\begin{array}{cc}\hfill & u(x,y,0)=v(x,y),u_t(x,y,0)=\nu (x,y),(x,y)\in \Omega ,\hfill \end{array}$$

(2)

$$\begin{array}{cc}\hfill & u(x,y,t)=0,t\in (0,T],(x,y)\in \partial \Omega ,\hfill \end{array}$$

(3)

respectively, where $\Omega =(0,L_1)\times (0,L_2)$ is a bounded convex domain with continuous boundary $\partial \Omega$, two functions $v(x,y)$ and $f(x,y,t)$ represent sufficiently smooth functions, the kernel $\beta (t)=\frac{t^{\alpha -1}}{\Gamma (\alpha )}$ denotes the weakly singular function $(0<\alpha <1),$ and $\mu$ is a positive constant. The problem (1)–(3) is used to model wave propagation subject to heat flow in materials with memory [2] and viscoelastic mechanics [3,4] and heat conduction in shape-memory materials, among other applications [5,6,7,8,9,10] and references therein.

Due to the potential applications in various fileds, considerable attention has been assigned to the progress of the theoretical and numerical solutions for the hyperbolic equation. Hlavacek [11] analyzed and proved the stability, uniqueness and existence of solutions for a class of problems, including hyperbolic problems. Heard [12] studied the problem of global solvability for the hyperbolic Volterra integrodifferential equation in a Banach space. Xu [13] studied the observability virtues of the temporal discrete approximation formula for hyperbolic integrodifferential equations. Fairweather [14] formulated and rigorously analyzed an alternating direction implicit orthogonal spline collocation (ADIOSC) formulation applied to the inhomogeneous hyperbolic equation on the unit square subject to appropriate ICs and homogeneous Dirichlet BC. Bialecki [15] proposed and analyzed an ADIOSC approach for the linear second-order hyperbolic initial-boundary value problem. Yannik [16] developed the finite element technique for certain hyperbolic and parabolic partial integrodifferential equations. Inoan et al. [17,18] studied the semi-Hyers–Ulam–Rassias stability for a class of integro-differential equations.

To the best of our knowledge, so far, there are no research papers on the hyperbolic integrodifferential equations, including a weakly singular kernel in the two-dimensional domain. The main purpose of the this paper is to introduce an efficient numerical scheme for Equation (1). The ADIOSC approach is a valuable and robust algorithm for obtaining the approximate solutions of a broad category of the ordinary and partial differential equations (ODEs and PDEs, respectively), which have been introduced in [19,20,21,22]; besides, this method also has some competitiveness and vitality compared to some classical methods [23,24,25,26,27,28,29,30].

Our main contributions are as follows: (I) We extended applications of an ADIOSC approach to the solution of the second-order hyperbolic problem containing a weakly singular kernel. The difficulty of theoretical analysis comes from the viscoelasticity term, which is essentially an integral with the time introduced in (1). For this new hyperbolic integrodifferential Equation (1), some traditional numerical quadrature rules such as the compound quadrature rules becomes difficult to analyze. Therefore, some standard approaches to treat the observability problems, such as the non-negativity skill [31], can not be used directly. We propose a new technique to handle the discrete integral term in the theoretical analysis, which is different from the method in [32,33]. (II) The ADI algorithm we considered can solve the high-dimensional problem of space into a series of one-dimensional sub-problems, which greatly reduces the computational cost [34,35,36,37,38]; in addition, the proposed OSC technology can obtain the fourth-order spatial convergence [19], which can greatly improve the numerical precision. (III) The theoretical results proved by us are validated by several numerical examples.

The layout of this paper is as follows: Section 2 introduces some basic notations and mathematical preliminaries. Section 3 proposes an ADIOSC approach for computing the solution of the problems (1)–(3). Section 4 examines the stability and convergence analysis of the full discretization. Section 5 presents two numerical examples to show high accuracy and performance for the ADIOSC approach and support the theoretical analysis. Finally, Section 6 contains some concluding remarks.

Remark 1.

Throughout this article, C represents a generic positive constant as a generic constant free of the time-space step parameters at every occurrence.

2. Notations and Auxiliary Results

In this section, we provide the basic concepts, and definitions needed to use the main results in the subsequent paper. For any bounded region $\Omega \subset {\mathbb{R}}^2,$ let us define the inner product:

$$\langle u,\upsilon \rangle ={\int }_{\Omega }u(x,y)\upsilon (x,y)dxdy,$$

with the $H^m(\Omega )$ norm on Sobolev:

$$\begin{array}{c}\hfill \begin{array}{c}\hfill {\parallel u\parallel }_{H^m}={\left(\sum _{0\le {\alpha }_1+{\alpha }_2\le m}{∥\frac{{\partial }^{{\alpha }_1+{\alpha }_2}u}{\partial x^{{\alpha }_1}\partial y^{{\alpha }_2}}∥}^2\right)}^{\frac{1}{2}},m>0.\end{array}\end{array}$$

Suppose that $I=(0,L_1)$ and $J=(0,L_2)$. Let ${\delta }_x={\{x_i\}}_{i=0}^{N_x}$ and ${\delta }_y={\{y_j\}}_{j=0}^{N_y}$ be two partitions of $\overline{I}\equiv [0,L_1]$ and $\overline{J}\equiv [0,L_2]$, so that:

$$0=x_0<x_1<\cdots <x_{N_x}=L_1,0=y_0<y_1<\cdots <y_{N_y}=L_2.$$

${\delta }_y\times {\delta }_x=\delta$ represents the partition of $\Omega$. For the sake of convenience, we denote:

$$h({\delta }_x)=\underset{1\le i\le N_x}{max}(x_i-x_{i-1}),h({\delta }_y)=\underset{1\le j\le N_y}{max}(y_j-y_{j-1}),$$

$${\gamma }_{ij}=(x_{i-1},x_i)\times (y_{j-1},y_j),h=max(h({\delta }_x),h({\delta }_y)),j=1,\dots ,N_y,i=1,\dots ,N_x.$$

Suppose that $M(r,{\delta }_y)$ and $M(r,{\delta }_x)$ are the spaces of piecewise polynomials of degree $\le r(r\ge 3)$, represented by:

$$M(r,{\delta }_x)=\{\upsilon |\upsilon \in C^1(I),\upsilon \in M_r([x_{i-1},x_i]),1\le i\le N_x,\upsilon (0)=\upsilon (L_1)=0\},$$

and,

$$M(r,{\delta }_y)=\{\upsilon |\upsilon \in C^1(I),\upsilon \in M_r([y_{j-1},y_j]),1\le j\le N_y,\upsilon (0)=\upsilon (L_2)=0\},$$

in which $M_r$ represents the set of polynomials of degree at most r.

Let:

$$M(r,{\delta }_x)\otimes M(r,{\delta }_y)=M(\delta ),$$

to be the set of all functions that represent finite linear combinations of products in the case of ${\upsilon }_1{\upsilon }_2$, such that ${\upsilon }_2\in M(r,{\delta }_y)$ and ${\upsilon }_1\in M(r,{\delta }_x)$.

Suppose that ${\{{\sigma }_l\}}_{l=1}^{r-1}$ and ${\{{\sigma }_k\}}_{k=1}^{r-1}$ are the nodal points of the two-point Gaussian quadrature rule over $\overline{I}$ and $\overline{J}$, respectively, and let:

$${\sigma }_{li}({\delta }_x)={\sigma }_l{h}_i+x_{i-1},l=1,\dots ,r-1,i=1,\dots ,N_x,$$

and,

$${\sigma }_{kj}({\delta }_y)={\sigma }_k{h}_j+y_{j-1},k=1,\dots ,r-1,j=1,\dots ,N_y,$$

be two collection of Gauss nodes in the x-direction and the y-direction, respectively.

$$\Gamma =\{({\sigma }_{li}({\delta }_x),{\sigma }_{kj}({\delta }_y))\},l,k=1,\dots ,r-1,i=1,\dots ,N_x,j=1,\dots ,N_y,$$

is the set of Gauss nodes on the domain $\Omega$.

Here, let us introduce the discrete inner product $\langle ·,·\rangle$ as:

$$\langle u,\upsilon \rangle =\sum _{i=1}^{N_x}\sum _{j=1}^{N_y}{h}_i{h}_j\sum _{l=1}^{r-1}\sum _{k=1}^{r-1}{\gamma }_l{\gamma }_k(u\upsilon )({\sigma }_{li}({\delta }_x),{\sigma }_{kj}({\delta }_y)),$$

with the norm ${|·|}_D$

$${|\upsilon |}_D^2=\langle \upsilon ,\upsilon \rangle .$$

such that u and v are defined on $\Gamma$.

If X represents a norm space having norm ${\parallel ·\parallel }_X$, then we introduce $C([0,T],X)$ the set of functions $\upsilon \in C({\overline{\Lambda }}_T)\equiv C^{0,0,0}({\overline{\Lambda }}_T)$ so that $\upsilon (·,t)\in X$ with $t\in [0,T]$, and

$${\parallel \upsilon \parallel }_{C([0,T],X)}=\underset{0\le t\le T}{max}{\parallel \upsilon (·,t)\parallel }_X\le \infty ,$$

where ${\Lambda }_T=\Omega \times (0,T]$.

Suppose that $C^{p,q,s}({\overline{\Lambda }}_T)$ is the set of functions $\upsilon (x,y,t)$ so that $\frac{{\partial }^{i+j+n}\upsilon }{\partial x^i\partial y^j\partial t^n}$ is continuous over ${\overline{\Lambda }}_T$ with $0\le i\le p,0\le j\le q,$ and $0\le n\le s.$ For $\upsilon \in C^{p,q,s}({\overline{\Lambda }}_T),$ we define ${\parallel \upsilon \parallel }_{C^{p,q,s}}$ by:

$${\parallel \upsilon \parallel }_{C^{p,q,s}}=\underset{0\le i\le p,0\le j\le q,0\le n\le s}{max}\underset{(x,y,t)\in {\overline{\Lambda }}_T}{max}|\frac{{\partial }^{i+j+n}\upsilon }{\partial x^i\partial y^i\partial t^n}|.$$

Then, we define $L^p(X)$ by:

$$L^p(X)=\nu :\nu (·,t)\in X,t\in {[0,T];\parallel \nu \parallel }_{L^p(X)}<\infty ,$$

in which,

$${\parallel \nu \parallel }_{L^p}(X)={\left({\int }_0^T{\parallel \nu \parallel }^p(X)dt\right)}^{1/p}{,\parallel \nu \parallel }_{L^{\infty }}(X)=\underset{0\le t\le T}{sup}{\parallel \nu \parallel }_X.$$

3. Discretization

Suppose that ${\{t_n\}}_{n=0}^K$ is a partition of temporal interval $[0,T]$ so that $t_n=nk,k=T/K$. Then, we define the convolution quadrature [39,40] to discretize the fractional integral term by:

$$\begin{array}{c}\hfill \begin{array}{c}\hfill Q_n^{(\alpha )}(\varphi )=k^{\alpha }\sum _{j=0}^n{q}_{n-j}^{(\alpha )}{\varphi }^j+k^{\alpha }{\tilde{q}}_n^{(\alpha )}{\varphi }^0\approx {\int }_0^t\beta (t_n-s)\varphi (s)ds,\end{array}\end{array}$$

(4)

where the quadrature weights $q_{n-j}^{(\alpha )}$ can be computed by their generating power series:

$$\begin{array}{c}\hfill \widehat{\beta }\left[\sigma (z)\right]={\left[\frac{1}{2}(3-4z+z^2)\right]}^{-\alpha }=\sum _{j=0}^{\infty }{q}_j^{(\alpha )}{z}^j,\end{array}$$

in which $\widehat{\beta }$ represents the Laplace transform of the convolution kernel.

Now, we employ the correction quadrature weights ${\tilde{q}}_n^{(\alpha )}$ for discretizing the integral with the second-order accuracy as follows:

$$\begin{array}{c}\hfill k^{\alpha }\sum _{j=0}^n{q}_{n-j}^{(\alpha )}+k^{\alpha }{\tilde{q}}_n^{(\alpha )}=\frac{1}{\Gamma (\alpha )}{\int }_0^{t_n}{(t_n-s)}^{\alpha -1}ds=\frac{1}{\Gamma (\alpha +1)}{t}_n^{\alpha }\le C(T),\end{array}$$

(5)

such that the quadrature formula becomes exact for constant. In the following lemma, we represent the quadrature error by:

$$E{(\varphi )}^{(\alpha )}(t_n)={\int }_0^{t_n}\beta (t_n-s)\varphi (s)ds-Q_n^{(\alpha )}(\varphi ).$$

Lemma 1

([41]). Assume that ϕ is a real, continuously differentiable function, with ${\varphi }_{tt}$ integrable and continuous over $(0,T)$. Then, the following error of the fractional quadrature rule is derived by:

$$\begin{array}{c}\hfill {|E(\varphi )}^{(\alpha )}(t_n)|\le C{k}^{1+\alpha },n\ge 1.\end{array}$$

Then, we present several lemmas required in the derivation of the error estimates. First, we define a map W, which plays an important role in the theoretical analysis. Let u be the solution of (1). Then, the differentiable map $W:[0,T]\to M(\delta )$ is defined by:

$$\Delta (u-W)=0,\mathrm{on}\Gamma \times [0,T],t>0.$$

(6)

To formulate the ADIOSC approach for the hyperbolic integrodifferential Equations (1)–(3), we restate (1) as a system of equations by letting $\varphi =\frac{\partial u}{\partial t}$. Then, the relation (1) converts to:

$$\begin{array}{c}\hfill \begin{array}{c}\hfill f(x,y,t)+{\int }_0^t\beta (t-s)\Delta u(x,y,s)=\frac{\partial \varphi }{\partial t}-\mu \Delta u(x,y,t),t\in (0,T],(x,y)\in \Omega .\end{array}\end{array}$$

(7)

Regarding the convolution quadrature defined in (4) with second-order accuracy, the OSC approach for the approximation of (1) is to find ${\{U^n\}}_{n=1}^K\in M(\delta )$ over $\Gamma$ fulfilling:

$$\begin{array}{c}\hfill \begin{array}{c}\hfill {\partial }_t{\Phi }^n-\mu \Delta U^{n+\frac{1}{2}}=k^{\alpha }\left(\sum _{p=0}^n{q}_{n-p}^{(\alpha )}\Delta U^{p+\frac{1}{2}}+{\tilde{q}}_{n+\frac{1}{2}}^{(\alpha )}\Delta U^0\right)+f^{n+\frac{1}{2}},1\le n\le K,\end{array}\end{array}$$

(8)

in which,

$$\begin{array}{c}\hfill \begin{array}{c}\hfill {\Phi }^{n+1}={\partial }_t{U}^{n+1}=\frac{U^{n+1}-U^n}{k},U^{n+\frac{1}{2}}=(U^n+U^{n+1})/2,{\partial }_t{\Phi }^n=\frac{{\Phi }^{n+1}-{\Phi }^n}{k}.\end{array}\end{array}$$

(9)

Note: ${\Phi }^n\in M(\delta )$ is an approximation of ${\varphi }^n,n=0,1,\dots ,K$ in (7). We choose the initial value $U^0=W^0$, defined by (6) with u replaced by v and ${\Phi }^0={\left(\frac{\partial W}{\partial t}\right)}_0$, obtained by differentiating (6) with respect to t and replacing $\frac{\partial u}{\partial t}$ by $\nu$. If we solve (9) for $U^{n+l}(l=\{1,\frac{1}{2}\})$ and substitute the resulting expression in (8) with $\lambda =\mu +q_0^{(\alpha )}{k}^{\alpha }$, we have:

$$\begin{array}{c}\hfill \begin{array}{c}\hfill {\partial }_t{\Phi }^n-\frac{\lambda k}{2}\Delta {\Phi }^{n+1}=k^{\alpha }\left(\sum _{p=0}^{n-1}{q}_{n-p}^{(\alpha )}\Delta U^{p+\frac{1}{2}}+{\tilde{q}}_{n+\frac{1}{2}}^{(\alpha )}\Delta U^0\right)+f^{n+\frac{1}{2}}+\lambda \Delta U^n\mathrm{on}\Gamma ,1\le n\le K.\end{array}\end{array}$$

(10)

If $E^{n+1}={\Phi }^{n+1}-{\Phi }^n$, then we can rewrite (9) and (10) as:

$$\begin{array}{c}\hfill \begin{array}{c}\hfill E^{n+1}-\frac{\lambda k^2}{4}\Delta E^{n+1}=F^{n+1}\mathrm{on}\Gamma ,1\le n\le K,\end{array}\end{array}$$

(11)

in which,

$$\begin{array}{c}\hfill \begin{array}{c}\hfill F^{n+1}=k\left[k^{\alpha }\left(\sum _{p=0}^{n-1}{q}_{n-p}^{(\alpha )}\Delta U^{p+\frac{1}{2}}+{\tilde{q}}_{n+\frac{1}{2}}^{(\alpha )}\Delta U^0\right)+f^{n+\frac{1}{2}}(\eta ({\delta }_x),\eta ({\delta }_y))+\lambda \Delta U^n+\frac{\lambda k}{2}\Delta {\Phi }^n\right].\end{array}\end{array}$$

(12)

If the small term,

$$\begin{array}{c}\hfill \frac{{\lambda }^2{k}^4}{16}\frac{{\partial }^4}{\partial x^2\partial y^2}{E}^{n+1}\end{array}$$

(13)

is added to the left-hand side of (11), then we have:

$$\begin{array}{c}\hfill \begin{array}{c}\hfill \left[1-\frac{\lambda k^2}{4}\Delta +\frac{{\lambda }^2{k}^4}{16}\frac{{\partial }^4}{\partial x^2\partial y^2}\right]E^{n+1}=F^{n+1}\mathrm{on}\Gamma ,1\le n\le K,\end{array}\end{array}$$

(14)

from which we obtain the approximation:

$$\begin{array}{c}\hfill \begin{array}{c}\hfill U^{n+1}=U^n+k{\Phi }^n+k{E}^{n+1}.\end{array}\end{array}$$

(15)

We now rewrite Equations (14) and (15) as an ADI matrix approach. For this aim, suppose that ${\{{\phi }_i\}}_{i=1}^{M_1}$ and ${\{{\psi }_j\}}_{j=1}^{M_2}$ represent bases for the subspace $M(r,{\delta }_y)$ and $M(r,{\delta }_x)$, respectively, so that $M_2=(r-1)N_y$ and $M_1=(r-1)N_x$. Make:

$$\begin{array}{c}\hfill \begin{array}{c}\hfill U^n(x,y)=\sum _{i=1}^{M_1}\sum _{j=1}^{M_2}{\alpha }_{ij}^{(n)}{\phi }_i(x){\psi }_j(y),\end{array}\end{array}$$

(16)

and,

$$\begin{array}{c}\hfill \begin{array}{c}\hfill {\Phi }^n(x,y)=\sum _{i=1}^{M_1}\sum _{j=1}^{M_2}{\beta }_{ij}^{(n)}{\phi }_i(x){\psi }_j(y),\end{array}\end{array}$$

(17)

in which,

$$\begin{array}{c}\hfill \begin{array}{c}\hfill {\alpha }^{(n)}={[{\alpha }_{11}^{(n)},{\alpha }_{12}^{(n)},\dots ,{\alpha }_{1M_2}^{(n)},{\alpha }_{21}^{(n)},\dots ,{\alpha }_{M_1{M}_2}^{(n)}]}^T,{\beta }^{(n)}={[{\beta }_{11}^{(n)},{\beta }_{12}^{(n)},\dots ,{\beta }_{1M_2}^{(n)},{\beta }_{21}^{(n)},\dots ,{\beta }_{M_1{M}_2}^{(n)}]}^T.\end{array}\end{array}$$

Define:

$$\begin{array}{c}\hfill \begin{array}{c}\hfill A_x=\{{(a_{ij}^x)}_{i,j=1}^{M_1},a_{ij}^x=-{\phi }_j^{{}^{″}}({\sigma }_i^x)\},B_x=\{{(b_{ij}^x)}_{i,j=1}^{M_1},b_{ij}^x={\phi }_j({\sigma }_i^x)\},\\ \hfill A_y=\{{(a_{ij}^y)}_{i,j=1}^{M_2},a_{ij}^y=-{\psi }_j^{{}^{″}}({\sigma }_i^y)\},B_y=\{{(b_{ij}^y)}_{i,j=1}^{M_2},b_{ij}^y={\psi }_j({\sigma }_i^y)\},\end{array}\end{array}$$

and,

$$\begin{array}{c}\hfill \begin{array}{c}\hfill \{F^{(n)}\}={[F^{(n)}({\sigma }_1^x,{\sigma }_1^y),F^{(n)}({\sigma }_1^x,{\sigma }_2^y),\dots ,F^{(n)}({\sigma }_1^x,{\sigma }_{M_2}^y),F^{(n)}({\sigma }_2^x,{\sigma }_1^y),\dots ,F^{(n)}({\sigma }_{M_1}^x,{\sigma }_{M_2}^y)]}^T.\end{array}\end{array}$$

Then, Equation (14) may be written in the following matrix-vector form with ${\gamma }_{ij}^{(n+1)}={\beta }_{ij}^{(n+1)}-{\beta }_{ij}^{(n)}$,

$$\begin{array}{c}\hfill \begin{array}{c}\hfill \left[B_1\otimes B_2+\frac{\lambda k^2}{4}(A_1\otimes B_1+B_1\otimes A_1)+\frac{{\lambda }^2{k}^4}{16}{A}_1\otimes A_2\right]{\gamma }^{(n+1)}=F^{(n+1)},\end{array}\end{array}$$

(18)

$$\begin{array}{c}\hfill \begin{array}{c}\hfill {\alpha }^{(n+1)}={\alpha }^{(n)}+k\left[{\beta }^{(n)}+{\gamma }^{(n+1)}\right],\end{array}\end{array}$$

(19)

and from Equation (12), the components of the vector $F^{(n+1)}$ in Equation (18) are given by:

$$\begin{array}{c}\hfill \begin{array}{cc}\hfill F^{(n+1)}=& k(A_1\otimes B_1+B_1\otimes A_1)\left[k^{\alpha }\left(\sum _{p=0}^{n-1}{q}_{n-p}^{(\alpha )}{\alpha }^{(p+\frac{1}{2})}+{\tilde{q}}_{n+\frac{1}{2}}^{(\alpha )}\Delta {\alpha }^{(0)}\right)+\lambda {\alpha }^{(n)}+\frac{\lambda k}{2}{\beta }^{(n)}\right]\hfill \\ & +k{f}^{n+\frac{1}{2}}(\eta ({\delta }_x),\eta ({\delta }_y)).\hfill \end{array}\end{array}$$

Then, relation (18) becomes:

$$\begin{array}{c}\hfill \begin{array}{cc}& \left[(B_x+\frac{\lambda k^2}{4}{A}_x)\otimes I_{M_2}\right]{\{{\nu }^{*}\}}^{n+1}=\{F^{(n+1)}\},\hfill \\ & \left[I_{M_1}\otimes (B_y+\frac{\lambda k^2}{4}{A}_y)\right]\{{\gamma }^{(n+1)}\}={\{{\nu }^{*}\}}^{n+1},1\le n\le K.\hfill \end{array}\end{array}$$

(20)

Therefore, we obtain ${\beta }^{(n+1)}={\gamma }^{(n+1)}+{\beta }^{(n)}$ from Equation (20) in the standard way. Then, the vector ${\alpha }^{(n+1)}$ is computed from Equation (19). Taking ${\alpha }^{(n+1)}$ in Equation (16), we achieve the desired $U^{n+1},1\le n\le K$.

4. Error Analysis of the ADIOSC Approach

In this section, we analyze the stability and convergence analysis of the ADIOSC approach in $L^2$-norm.

4.1. Stability of the ADIOSC Approach

Here, we introduce the following lemmas to examine the stability of the ADIOSC approach.

Lemma 2

([14]). For any $V,U\in M(\delta )$, we have:

$$\begin{array}{c}\hfill \langle -\Delta V,U\rangle =\langle V,-\Delta U\rangle ,\end{array}$$

(21)

$$\begin{array}{c}\hfill \langle -\Delta V,V\rangle \ge {C\parallel \nabla V\parallel }^2\ge 0,\end{array}$$

(22)

$$\begin{array}{c}\hfill \langle -\Delta V,U\rangle \le C\parallel \nabla V\parallel \parallel \nabla U\parallel .\end{array}$$

(23)

Lemma 3.

The norms ${\parallel ·\parallel }_D$ and $|·|$ are equivalent over $M(\delta ).$

Theorem 1.

Suppose $U^n\in M(\delta )$ $(1\le n\le K)$, satisfies (14) and (15) with $U^0=W^0$, then there exists a positive constant C, independent of h and k, such that:

$$|U^n{|}_D\le |U^0{|}_D+Ck\sum _{n=0}^{K-1}|f^{n+\frac{1}{2}}{|}_D+\frac{\lambda k^2}{4}\parallel \frac{{\partial }^2{\Phi }^0}{\partial x\partial y}\parallel +\parallel \nabla U^0\parallel +|{\Phi }^0{|}_D.$$

Proof. 

Let us rewrite Equation (14) in the following form:

$$\begin{array}{c}\hfill \begin{array}{cc}\hfill {\partial }_t{\Phi }^n-\mu \Delta U^{n+\frac{1}{2}}+\frac{{\lambda }^2{k}^4}{16}\frac{{\partial }^4{\partial }_t{\Phi }^n}{\partial x^2\partial y^2}=k^{\alpha }(\sum _{p=0}^n{q}_{n-p}^{(\alpha )}\Delta U^{p+\frac{1}{2}}+{\tilde{q}}_{n+\frac{1}{2}}^{(\alpha )}\Delta U^0)& +f^{n+\frac{1}{2}},1\le n\le K.\hfill \end{array}\end{array}$$

(24)

Using the inequalities:

$$\begin{array}{c}\hfill \begin{array}{c}\hfill \langle {\partial }_t{\Phi }^n,{\Phi }^{n+\frac{1}{2}}\rangle \ge \frac{1}{2k}(|{\Phi }^{n+1}{|}_D^2-|{\Phi }^n{|}_D^2),\end{array}\end{array}$$

(25)

$$\begin{array}{c}\hfill \begin{array}{c}\hfill \langle \Delta U^{n+\frac{1}{2}},{\Phi }^{n+\frac{1}{2}}\rangle =\frac{1}{2}{\partial }_t\langle \Delta U^{n+1},U^{n+1}\rangle ,\end{array}\end{array}$$

(26)

$$\begin{array}{c}\hfill \begin{array}{c}\hfill 〈\frac{{\partial }^4{\partial }_t{\Phi }^n}{\partial x^2\partial y^2},{\Phi }^{n+\frac{1}{2}}〉\ge \frac{1}{2}{\partial }_t〈\frac{{\partial }^4{\Phi }^{n+1}}{\partial x^2\partial y^2},{\Phi }^{n+1}〉,1\le n\le K,\end{array}\end{array}$$

(27)

then applying the inner product of both sides of (24) with ${\Phi }^{n+\frac{1}{2}}$, for $1\le n\le K$, we obtain:

$$\begin{array}{cc}\hfill \frac{1}{2k}(|{\Phi }^{n+1}{|}_D^2-|{\Phi }^n{|}_D^2)& -\frac{\mu }{2}{\partial }_t\langle \Delta U^{n+1},U^{n+1}\rangle +\frac{{\lambda }^2{k}^4}{32}{\partial }_t〈\frac{{\partial }^4{\Phi }^{n+1}}{\partial x^2\partial y^2},{\Phi }^{n+1}〉\hfill \\ & =〈Q_{n+\frac{1}{2}}^{(\alpha )}(\Delta U),{\Phi }^{n+\frac{1}{2}}〉+\langle f^{n+\frac{1}{2}},{\Phi }^{n+\frac{1}{2}}\rangle .\hfill \end{array}$$

(28)

Now pay attention to the first term on the right-hand side of (28), and with the definition in (4), we have:

$$\begin{array}{cc}\hfill 〈Q_{n+\frac{1}{2}}^{(\alpha )}(\Delta U),{\partial }_t{U}^{n+1}〉& ={\partial }_t〈Q_{n+1}^{(\alpha )}(\Delta U),U^{n+1}〉-〈{\partial }_t{Q}_{n+\frac{1}{2}}^{(\alpha )}(\Delta U),U^{n+\frac{1}{2}}〉\hfill \\ & -\langle k^{\alpha -1}{q}_0\Delta U^{n+1},U^{n+1}\rangle ,\hfill \end{array}$$

(29)

where,

$$\begin{array}{c}\hfill \begin{array}{cc}\hfill 〈{\partial }_t{Q}_{n+\frac{1}{2}}^{(\alpha )}(\Delta U),U^{n+\frac{1}{2}}〉& =〈k^{\alpha }(\sum _{p=0}^n{\partial }_t{q}_{n-p}^{(\alpha )}\Delta U^{p+\frac{1}{2}}+{\partial }_t{\tilde{q}}_{n+\frac{1}{2}}^{(\alpha )}\Delta U^0),U^{n+\frac{1}{2}}〉,\hfill \end{array}\end{array}$$

(30)

and,

$$\begin{array}{c}\hfill \begin{array}{c}\hfill {\partial }_t〈Q_{n+1}^{(\alpha )}(\Delta U),U^{n+1}〉=\frac{1}{k}\left[〈Q_{n+1}^{(\alpha )}(\Delta U),U^{n+1}〉-〈Q_n^{(\alpha )}(\Delta U),U^n〉\right].\end{array}\end{array}$$

(31)

From Lemma 2, we know:

$$\begin{array}{c}\hfill \begin{array}{c}\hfill -\langle k^{\alpha -1}{q}_0\Delta U^{n+1},U^{n+1}\rangle \ge 0.\end{array}\end{array}$$

(32)

Using Equations (30) and (31), dropping the positive term (32) in (29), we have:

$$\begin{array}{c}\hfill \begin{array}{cc}& 〈Q_{n+\frac{1}{2}}^{(\alpha )}(\Delta U),{\partial }_t{U}^{n+1}〉\le {\partial }_t〈Q_{n+1}^{(\alpha )}(\Delta U),U^{n+1}〉-〈{\partial }_t{Q}_{n+\frac{1}{2}}^{(\alpha )}(\Delta U),U^{n+\frac{1}{2}}〉.\hfill \end{array}\end{array}$$

(33)

Multiplying (28) by $2k$ and summing over $n=0,\dots ,N-1,$ we obtain:

$$\begin{array}{c}\hfill \begin{array}{cc}& |{\Phi }^N{|}_D^2-\mu \langle \Delta U^N,U^N\rangle +\frac{{\lambda }^2{k}^4}{16}\parallel \frac{{\partial }^2{\Phi }^N}{\partial x\partial y}{\parallel }^2-2k\sum _{n=0}^{N-1}{\partial }_t\langle Q_{n+1}^{(\alpha )}(\Delta U),U^{n+1}\rangle \hfill \\ & +2k\sum _{n=0}^{N-1}\langle {\partial }_t{Q}_{n+\frac{1}{2}}^{(\alpha )}(\Delta U),U^{n+\frac{1}{2}}\rangle \le Ck\sum _{n=0}^{N-1}(|{\Phi }^{n+\frac{1}{2}}{|}_D^2+|f^{n+\frac{1}{2}}{|}_D^2)\hfill \\ & +\frac{{\lambda }^2{k}^4}{16}\parallel \frac{{\partial }^2{\Phi }^0}{\partial x\partial y}{\parallel }^2+\parallel \nabla U^0{\parallel }^2+{|{\Phi }^0|}_D^2.\hfill \end{array}\end{array}$$

(34)

Again, noting that $\frac{{\lambda }^2{k}^4}{16}\parallel \frac{{\partial }^2{\Phi }^N}{\partial x\partial y}{\parallel }^2\ge 0$, dropping this positive term in (34), using Lemma 2 that $-\mu \langle \Delta U^N,U^N\rangle \ge C{\parallel \nabla U^N\parallel }^2$ and the Grönwall inequality, we have:

$$\begin{array}{c}\hfill \begin{array}{cc}\hfill |{\Phi }^N{|}_D^2+C{\parallel \nabla U^N\parallel }^2\le & |2k\sum _{n=0}^{N-1}{\partial }_t\langle Q_{n+1}^{(\alpha )}(\Delta U),U^{n+1}\rangle |+2k|\sum _{n=0}^{N-1}\langle {\partial }_t{Q}_{n+\frac{1}{2}}^{(\alpha )}(\Delta U),U^{n+\frac{1}{2}}\rangle |\hfill \\ & +Ck\sum _{n=0}^{N-1}(|{\Phi }^{n+\frac{1}{2}}{|}_D^2+|f^{n+\frac{1}{2}}{|}_D^2)+\frac{{\lambda }^2{k}^4}{16}\parallel \frac{{\partial }^2{\Phi }^0}{\partial x\partial y}{\parallel }^2+\parallel \nabla U^0{\parallel }^2+{|{\Phi }^0|}_D^2.\hfill \end{array}\end{array}$$

(35)

For the first term on the right-hand side of (35), we have:

$$\begin{array}{c}\hfill \begin{array}{cc}\hfill \left|k\sum _{n=0}^{N-1}{\partial }_t\langle Q_{n+1}^{(\alpha )}(\Delta U),U^n\rangle \right|& =\left|\langle Q_N^{(\alpha )}(\Delta U),U^N\rangle \right|=\left|k^{\alpha }\sum _{p=0}^{N-1}{q}_{N-p}^{(\alpha )}\langle \Delta U^p,U^N\rangle \right|\hfill \\ & \le C{k}^{\alpha }\sum _{p=0}^{N-1}{q}_{N-p}^{(\alpha )}(\frac{1}{2}|\nabla U^p{|}_D^2+\frac{1}{2}|\nabla U^N{|}_D^2).\hfill \end{array}\end{array}$$

(36)

For $q\in C^1,$ with (5), we know:

$$\begin{array}{c}\hfill \begin{array}{c}\hfill k^{\alpha }\sum _{j=0}^n{q}_{n-j}^{(\alpha )}+k^{\alpha }{\tilde{q}}_n^{(\alpha )}=\frac{1}{\Gamma (\alpha )}{\int }_0^{t_n}{(t_n-s)}^{\alpha -1}ds=\frac{1}{\Gamma (\alpha +1)}{t}_n^{\alpha }\le C(T),\end{array}\end{array}$$

(37)

then we get:

$$\begin{array}{c}\hfill \begin{array}{c}\hfill {\partial }_t\left(k^{\alpha }\sum _{j=0}^n{q}_{n-j}^{(\alpha )}+k^{\alpha }{\tilde{q}}_n^{(\alpha )}\right)\le C.\end{array}\end{array}$$

(38)

By estimating the second term on the left-hand side of (35) and using Equations (4) and (38), we can write:

$$\begin{array}{c}\hfill \begin{array}{ccc}\hfill \left|k\sum _{n=0}^{N-1}〈{\partial }_t{Q}_{n-\frac{1}{2}}^{(\alpha )}(\Delta U),U^{n-\frac{1}{2}}〉\right|& =\left|k^{\alpha +1}〈\sum _{n=0}^{N-1}\sum _{p=0}^n{\partial }_t{q}_{n-p+\frac{1}{2}}^{(\alpha )}\Delta U^p,U^{n-\frac{1}{2}}〉\right|\hfill & \\ \hfill \le k\sum _{p=0}^{N-1}{|\nabla U^{p-\frac{1}{2}}|}_D^2.\end{array}\end{array}$$

(39)

Letting (36) and (39) in (35), we have:

$$\begin{array}{c}\hfill \begin{array}{cc}\hfill |{\Phi }^N{|}_D^2& +C\parallel \nabla U^N{\parallel }^2\le C{k}^{\alpha }\sum _{p=0}^{N-1}{q}_{N-p}^{(\alpha )}(|\nabla U^p{|}_D^2+|\nabla U^N{|}_D^2)+2k\sum _{p=0}^{N-1}{|\nabla U^{p-\frac{1}{2}}|}_D^2\hfill \\ & +Ck\sum _{n=0}^{N-1}(|{\Phi }^{n+\frac{1}{2}}{|}_D^2+|f^{n+\frac{1}{2}}{|}_D^2)+\frac{{\lambda }^2{k}^4}{16}\parallel \frac{{\partial }^2{\Phi }^0}{\partial x\partial y}{\parallel }^2+\parallel \nabla U^0{\parallel }^2+{|{\Phi }^0|}_D^2.\hfill \end{array}\end{array}$$

(40)

We now consider the term $|{\Phi }^N{|}_D^2$ in (40) that:

$$\begin{array}{c}\hfill |{\Phi }^N{|}_D^2=\frac{1}{k^2}|U^{N+1}-U^N{|}_D^2\ge \frac{1}{k^2}(|U^{N+1}{|}_D-|U^N{{|}_D)}^2,1\le N\le K.\end{array}$$

(41)

Using the Grönwall inequality for (40), with (41), we have:

$$\begin{array}{c}\hfill \begin{array}{cc}\hfill \frac{1}{k^2}(|U^{N+1}{|}_D-|U^N{|}_D{)}^2+C{\parallel \nabla U^N\parallel }^2\le & Ck\sum _{n=0}^{N-1}{|f^{n+\frac{1}{2}}|}_D^2+\frac{{\lambda }^2{k}^4}{16}{\parallel \frac{{\partial }^2{\Phi }^0}{\partial x\partial y}\parallel }^2\hfill \\ & +\parallel \nabla U^0{\parallel }^2+{|{\Phi }^0|}_D^2,\hfill \end{array}\end{array}$$

(42)

where $C\parallel \nabla U^N{\parallel }^2\ge 0,$ and we obtain:

$$\begin{array}{c}\hfill \begin{array}{cc}\hfill \frac{1}{k}(|U^{N+1}{|}_D-|U^N{|}_D)\le & Ck\sum _{n=0}^{N-1}|f^{n+\frac{1}{2}}{|}_D+\frac{\lambda k^2}{4}\parallel \frac{{\partial }^2{\Phi }^0}{\partial x\partial y}\parallel +\parallel \nabla U^0\parallel +|{\Phi }^0{|}_D.\hfill \end{array}\end{array}$$

(43)

Changing index N to n, summing for n from 0 to $K-1$, then we obtain:

$$\begin{array}{c}\hfill \begin{array}{cc}\hfill |U^K{|}_D\le & |U^0{|}_D+Ck\sum _{n=0}^{K-1}|f^{n+\frac{1}{2}}{|}_D+\frac{\lambda k^2}{4}\parallel \frac{{\partial }^2{\Phi }^0}{\partial x\partial y}\parallel +\parallel \nabla U^0\parallel +|{\Phi }^0{|}_D,\hfill \end{array}\end{array}$$

(44)

which finishes the proof of the theorem. □

4.2. Convergence of the ADIOSC Approach

Let us introduce the intermediate projections as differentiable maps $W:[0,T]\to M(\delta )$ in order to achieve an optimal convergence order satisfying [42] Equation (2.19).

$$\begin{array}{c}\hfill \begin{array}{c}\hfill \langle \Delta (u-W),w\rangle =0,w\in M(\delta ),\end{array}\end{array}$$

(45)

in which u represents the solution of (1) over $\Gamma \times [0,T]$. In the following, we will obtain the estimates on $u-W$ and its temporal derivatives.

For obtaining an optimal order $L^2$ error estimate for the ADIOSC approximation, let us remember the following two lemmas.

Lemma 4

([20]). For ${\partial }^{l}u/\partial t^l\in H^{r-j+3}$, with $0\le t\le T,0\le l,j\le 1,$ and W introduced in the relation (21), we have:

$$\begin{array}{c}\hfill \begin{array}{c}\hfill {∥\frac{{\partial }^l(u-W)}{\partial t^l}∥}_{H^j}\le C{h}^{r-j+1}{∥\frac{{\partial }^{l}u}{\partial t^l}∥}_{H^{r-j+3}},\end{array}\end{array}$$

where the constant C is independent of h.

Lemma 5

([19]). For any ${\partial }^{i}u/\partial t^i\in H^{r+3}$, with $0\le t\le T,$ for $i=0,1$, we have:

$$\begin{array}{c}\hfill \begin{array}{c}\hfill {\left|\frac{{\partial }^{l+i}(u-W)}{\partial x^{l_1}\partial y^{l_2}\partial t^i}\right|}_D\le C{h}^{r+1-l}{∥\frac{{\partial }^{i}u}{\partial t^i}∥}_{H^{r+3}},\end{array}\end{array}$$

in which $0\le l=l_1+l_2\le 4.$

The discrete problem (14) generates a system of algebraic equations at the nodes $t=t_n$. Therefore, we can indicate that the obtained system includes a unique solution at every temporal step $t=t_n$ based on the a priori bound for $U^n$. In the following, we obtain an error estimate.

Theorem 2.

Suppose u is the solution of (1)–(3), if $u\in C^{2,0,3}\cap C^{0,2,3}\cap C^{0,0,4},$ and $\frac{\partial u}{\partial t}\in C([0,T],H^{r+3}),$ $U^n$ $(1\le n\le K)$, satisfies (14) and (15) with $U^0=W^0$. Then, we have:

$$|u(t_n)-U^n{|}_D\le C(h^{r+1}+k^2).$$

Proof. 

For $1\le n\le K$, it is easy to show that (14) is equivalent to:

$$\begin{array}{c}\hfill \begin{array}{c}\hfill 〈\left[1-\frac{\lambda k^2}{4}\Delta +\frac{{\lambda }^2{k}^4}{16}\frac{{\partial }^4}{\partial x^2\partial y^2}\right]E^{n+1},v〉=\langle F^{n+1},v\rangle ,v\in M(\delta ),1\le n\le K,\end{array}\end{array}$$

(46)

where:

$$\begin{array}{c}\hfill \begin{array}{c}\hfill F^{n+1}=k\left[k^{\alpha }\left(\sum _{p=0}^{n-1}{q}_{n-p}^{(\alpha )}\Delta U^{p+\frac{1}{2}}+{\tilde{q}}_{n+\frac{1}{2}}^{(\alpha )}\Delta U^0\right)+f^{n+\frac{1}{2}}({\sigma }_{li}({\delta }_x),{\sigma }_{kj}({\delta }_y))+\lambda \Delta U^n+\frac{\lambda k}{2}\Delta {\Phi }^n\right].\end{array}\end{array}$$

(47)

On substituting for $E^{n+1}$ and dividing by k, we take the form:

$$\begin{array}{c}\hfill \begin{array}{cc}\hfill 〈{\partial }_t{\Phi }^n,v〉-〈\mu \Delta U^{n+\frac{1}{2}},v〉+〈\frac{{\lambda }^2{k}^4}{16}\frac{{\partial }^4{\partial }_t{\Phi }^n}{\partial x^2\partial y^2},v〉=& 〈k^{\alpha }(\sum _{p=0}^n{q}_{n-p}^{(\alpha )}\Delta U^{p+\frac{1}{2}}+{\tilde{q}}_{n+\frac{1}{2}}^{(\alpha )}\Delta U^0),v〉\hfill \\ & +〈f^{n+\frac{1}{2}},v〉,v\in M(\delta ),1\le n\le K,\hfill \end{array}\end{array}$$

(48)

$$\begin{array}{c}\hfill \begin{array}{c}\hfill {\partial }_t{U}^{n+1}={\Phi }^{n+1}=\frac{U^{n+1}-U^n}{k},1\le n\le K.\end{array}\end{array}$$

(49)

From (1), it follows that:

$$\begin{array}{c}\hfill \begin{array}{cc}& 〈{\partial }_t{\varphi }^n,v〉-〈\mu \Delta u^{n+\frac{1}{2}},v〉+〈\frac{{\lambda }^2{k}^4}{16}\frac{{\partial }^4{\partial }_t{\varphi }^n}{\partial x^2\partial y^2},v〉=〈k^{\alpha }(\sum _{p=0}^n{q}_{n-p}^{(\alpha )}\Delta u^{p+\frac{1}{2}}+{\tilde{q}}_{n+\frac{1}{2}}^{(\alpha )}\Delta u^0),v〉\hfill \\ & +〈f^{n+\frac{1}{2}},v〉+〈{\partial }_t{\varphi }^n-{\left(\frac{\partial u}{\partial t}\right)}_{n+\frac{1}{2}},v〉+\langle E{(\Delta u)}^{(\alpha )}(t_{n+1}),v\rangle ,v\in M(\delta ),1\le n\le K,\hfill \end{array}\end{array}$$

(50)

$$\begin{array}{c}\hfill \begin{array}{c}\hfill {\partial }_t{u}^{n+1}={\varphi }^{n+\frac{1}{2}}+{\partial }_t{u}^{n+1}-{\left(\frac{\partial u}{\partial t}\right)}_{n+1/2},1\le n\le K.\end{array}\end{array}$$

(51)

Based on W introduced in Equation (45), we can write:

$$\begin{array}{cc}\hfill & u(t_n)-U^n=(u(t_n)-W^n)-(U^n-W^n)={\eta }^n-{\xi }^n,1\le n\le K,\end{array}$$

(52)

$$\begin{array}{cc}& {\varphi }^n-{\Phi }^n=\left({\varphi }^n-{\left(\frac{\partial W}{\partial t}\right)}_n\right)-\left({\Phi }^n-{\left(\frac{\partial W}{\partial t}\right)}_n\right)={\widehat{\eta }}^n-{\widehat{\xi }}^n,1\le n\le K.\end{array}$$

(53)

Subtracting Equation (50) from Equation (48) and Equation (51) from Equation (49), and using Equations (45) and (52) yields:

$$\begin{array}{c}\hfill \begin{array}{cc}& 〈{\partial }_t{\widehat{\xi }}^n,v〉-〈\mu \Delta {\xi }^{n+\frac{1}{2}},v〉+〈\frac{{\lambda }^2{k}^4}{16}\frac{{\partial }^4{\partial }_t{\widehat{\xi }}^n}{\partial x^2\partial y^2},v〉-〈k^{\alpha }(\sum _{p=0}^n{q}_{n-p}^{(\alpha )}\Delta {\xi }^{p+\frac{1}{2}}+{\tilde{q}}_{n+\frac{1}{2}}^{(\alpha )}\Delta {\xi }^0),v〉\hfill \\ & =\langle {\zeta }^n,v\rangle ,v\in M(\delta ),1\le n\le K,\hfill \end{array}\end{array}$$

(54)

$$\begin{array}{c}\hfill \begin{array}{c}\hfill {\partial }_t{\xi }^n={\widehat{\xi }}^{n+\frac{1}{2}}-{\rho }^n,1\le n\le K,\end{array}\end{array}$$

(55)

where,

$$\begin{array}{c}\hfill \begin{array}{c}\hfill {\zeta }^n={\partial }_t{\widehat{\eta }}^n-{\partial }_t{\Phi }^n+{\left(\frac{\partial \varphi }{\partial t}\right)}_{n+1/2}+\frac{{\lambda }^2{k}^4}{16}\frac{{\partial }^4}{\partial x^2\partial y^2}({\widehat{\eta }}^n-{\varphi }^n)+E{(\Delta \eta )}^{(\alpha )}(t_{n+1}),\end{array}\end{array}$$

(56)

$$\begin{array}{c}\hfill \begin{array}{c}\hfill {\rho }^n={\widehat{\eta }}^{n+1/2}-{\partial }_t{\eta }^{n+1}+{\partial }_t{u}^{n+1}-{\left(\frac{\partial u}{\partial t}\right)}_{n+1/2}.\end{array}\end{array}$$

(57)

Note that the bound of $E{(\Delta \eta )}^{(\alpha )}(t_{n+1})$ has been given in Lemma 1.

In Equation (54), we choose $v={\widehat{\xi }}^{n+\frac{1}{2}}={\partial }_t{\xi }^n+{\rho }^n,$ (from (55)), and use the definition in (4), then:

$$\begin{array}{c}\hfill \begin{array}{cc}& 〈{\partial }_t{\widehat{\xi }}^n,{\widehat{\xi }}^{n+\frac{1}{2}}〉-\mu 〈\Delta {\xi }^{n+\frac{1}{2}},{\partial }_t{\xi }^n〉+\frac{{\lambda }^2{k}^4}{16}〈\frac{{\partial }^4{\partial }_t{\widehat{\xi }}^n}{\partial x^2\partial y^2},{\widehat{\xi }}^{n+\frac{1}{2}}〉-〈Q_{n+\frac{1}{2}}^{(\alpha )}(\Delta \xi ),{\widehat{\xi }}^{n+\frac{1}{2}}〉\hfill \\ & =\langle {\zeta }^n,{\widehat{\xi }}^{n+\frac{1}{2}}\rangle +\mu 〈\Delta {\xi }^{n+\frac{1}{2}},{\rho }^n〉,1\le n\le K.\hfill \end{array}\end{array}$$

(58)

Now, for $1\le n\le K$, with ${\xi }^0=0$, using (4), the fourth term on the left-hand side of (58) can be stated as:

$$\begin{array}{c}\hfill \begin{array}{c}\hfill 〈Q_{n+\frac{1}{2}}^{(\alpha )}(\Delta \xi ),{\partial }_t{\xi }^n〉={\partial }_t〈Q_{n+1}^{(\alpha )}(\Delta \xi ),{\xi }^n〉-〈{\partial }_t{Q}_{n-\frac{1}{2}}^{(\alpha )}(\Delta \xi ),{\xi }^{n-\frac{1}{2}}〉-\langle k^{\alpha -1}{q}_0\Delta {\xi }^n,{\xi }^n\rangle ,\end{array}\end{array}$$

(59)

where,

$$\begin{array}{c}\hfill \begin{array}{cc}\hfill 〈{\partial }_t{Q}_{n-\frac{1}{2}}^{(\alpha )}(\Delta \xi ),{\xi }^{n-\frac{1}{2}}〉& =〈k^{\alpha }\left(\sum _{p=0}^n{\partial }_t{q}_{n-p}^{(\alpha )}\Delta {\xi }^{p-\frac{1}{2}}+{\partial }_t{\tilde{q}}_{n-\frac{1}{2}}^{(\alpha )}\Delta {\xi }^0\right),{\xi }^{n-\frac{1}{2}}〉\hfill \\ & =〈k^{\alpha }\sum _{p=0}^n{\partial }_t{q}_{n-p}^{(\alpha )}\Delta {\xi }^{p-\frac{1}{2}},{\xi }^{n-\frac{1}{2}}〉,\hfill \end{array}\end{array}$$

(60)

and,

$$\begin{array}{c}\hfill \begin{array}{c}\hfill {\partial }_t〈Q_{n+1}^{(\alpha )}(\Delta \xi ),{\xi }^n〉=\frac{1}{\tau }\left[〈Q_{n+1}^{(\alpha )}(\Delta \xi ),{\xi }^{n+1}〉-〈Q_n^{(\alpha )}(\Delta \xi ),{\xi }^n〉\right].\end{array}\end{array}$$

(61)

From a simple calculation, it is easy to know that:

$$\begin{array}{c}\hfill \begin{array}{c}\hfill 〈\frac{{\partial }^4{\partial }_t{\widehat{\xi }}^n}{\partial x^2\partial y^2},{\widehat{\xi }}^{n+\frac{1}{2}}〉\ge \frac{1}{2}{\partial }_t〈\frac{{\partial }^4{\widehat{\xi }}^n}{\partial x^2\partial y^2},{\widehat{\xi }}^n〉,〈{\partial }_t{\widehat{\xi }}^n,{\widehat{\xi }}^{n+\frac{1}{2}}〉\ge \frac{1}{2}{\partial }_t\langle {\widehat{\xi }}^n,{\widehat{\xi }}^n\rangle ,\end{array}\end{array}$$

(62)

and,

$$\begin{array}{c}\hfill \begin{array}{c}\hfill -〈\Delta {\xi }^{n+\frac{1}{2}},{\partial }_t{\xi }^n〉=-\frac{1}{2}{\partial }_t\langle \Delta {\xi }^n,{\xi }^n\rangle ,1\le n\le K.\end{array}\end{array}$$

(63)

Using (59)–(63) for the left-hand side of (58), Lemma 2 for the second term on the right-hand side, and Young’s inequality that $ab\le \frac{a^2+b^2}{2}(a,b\in R),$ it is easy to see that:

$$\begin{array}{c}\hfill \begin{array}{cc}\hfill & {\partial }_t\left[〈{\widehat{\xi }}^n,{\widehat{\xi }}^n〉-\mu 〈\Delta {\xi }^n,{\xi }^n〉+\frac{{\lambda }^2{k}^4}{16}〈\frac{{\partial }^4{\widehat{\xi }}^n}{\partial x^2\partial y^2},{\widehat{\xi }}^n〉-〈Q_{n+1}^{(\alpha )}(\Delta \xi ),{\xi }^n〉\right]-\langle k^{\alpha -1}{q}_0\Delta {\xi }^{n+\frac{1}{2}},{\xi }^{n+\frac{1}{2}}\rangle \hfill \\ & +〈{\partial }_t{Q}_{n-\frac{1}{2}}^{(\alpha )}(\Delta \xi ),{\xi }^{n-\frac{1}{2}}〉\le C(|{\zeta }^n{|}_D^2+|{\widehat{\xi }}^{n+\frac{1}{2}}{|}_D^2+\mu |\nabla {\xi }^{n+\frac{1}{2}}{|}_D^2+|\nabla {\rho }^n{|}_D^2),1\le n\le K.\hfill \end{array}\end{array}$$

(64)

In view of (22) in Lemma 2, we conclude that the fifth term in the left-hand side form (64) is the non-negative term:

$$\begin{array}{c}\hfill \begin{array}{c}\hfill k^{\alpha -1}{q}_0\langle \Delta {\xi }^{n+\frac{1}{2}},{\xi }^{n+\frac{1}{2}}\rangle \ge k^{\alpha -1}{q}_0{\parallel \nabla {\xi }^{n+\frac{1}{2}}\parallel }^2\ge 0.\end{array}\end{array}$$

(65)

Dropping the positive term (65) in (64), multiplying (64) by k, and summing from $n=0$ to $n=K-1,$ and using the fact that ${\xi }^0={\widehat{\xi }}^0=0$, we obtain:

$$\begin{array}{c}\hfill \begin{array}{cc}& |{\widehat{\xi }}^K{|}_D^2-\mu \langle \Delta {\xi }^K,{\xi }^K\rangle +\frac{{\lambda }^2{k}^4}{16}\parallel \frac{{\partial }^2{\widehat{\xi }}^K}{\partial x\partial y}{\parallel }^2\le \left|k\sum _{n=0}^{K-1}{\partial }_t\langle Q_{n+1}^{(\alpha )}(\Delta \xi ),{\xi }^n\rangle \right|\hfill \\ & +k\sum _{n=0}^{K-1}\left|\langle {\partial }_t{Q}_{n-\frac{1}{2}}^{(\alpha )}(\Delta \xi ),{\xi }^{n-\frac{1}{2}}\rangle \right|+Ck\sum _{n=0}^{K-1}(|{\widehat{\xi }}^{n+\frac{1}{2}}{|}_D^2+\mu |\nabla {\xi }^{n+\frac{1}{2}}{|}_D^2)+k\sum _{n=0}^{K-1}(|{\zeta }^n{|}_D^2+|\nabla {\rho }^n{|}_D^2).\hfill \end{array}\end{array}$$

(66)

In view of (4) and (61), the first term in the right-hand side form (66) can be restated as:

$$\begin{array}{c}\hfill \begin{array}{cc}& k\sum _{n=0}^{K-1}{\partial }_t\langle Q_{n+1}^{(\alpha )}(\Delta \xi ),{\xi }^n\rangle =\langle Q_K^{(\alpha )}(\Delta \xi ),{\xi }^K\rangle =k^{\alpha }\sum _{p=0}^{K-1}{q}_{K-p}^{(\alpha )}\langle \Delta {\xi }^p,{\xi }^K\rangle ,\hfill \end{array}\end{array}$$

(67)

By virtue of Young’s inequality used above, we obtain:

$$\begin{array}{c}\hfill \begin{array}{c}\hfill k^{\alpha }\sum _{p=0}^{K-1}{q}_{K-p}^{(\alpha )}\langle \Delta {\xi }^p,{\xi }^K\rangle \le C{k}^{\alpha }\sum _{p=0}^{K-1}{q}_{K-p}^{(\alpha )}(\frac{1}{2}|\nabla {\xi }^p{|}_D^2+\frac{1}{2}|\nabla {\xi }^K{|}_D^2).\end{array}\end{array}$$

(68)

For $q\in C^1,$ with (5), we know that:

$$\begin{array}{c}\hfill \begin{array}{c}\hfill k^{\alpha }\sum _{j=0}^n{q}_{n-j}^{(\alpha )}+k^{\alpha }{\tilde{q}}_n^{(\alpha )}=\frac{1}{\Gamma (\alpha )}{\int }_0^{t_n}{(t_n-s)}^{\alpha -1}ds=\frac{1}{\Gamma (\alpha +1)}{t}_n^{\alpha }\le C(T),\end{array}\end{array}$$

(69)

then we have,

$$\begin{array}{c}\hfill \begin{array}{c}\hfill {\partial }_t(k^{\alpha }\sum _{j=0}^n{q}_{n-j}^{(\alpha )}+k^{\alpha }{\tilde{q}}_n^{(\alpha )})\le C.\end{array}\end{array}$$

(70)

Meanwhile, estimating the second term on the right-hand side of (66) and using (4), (70), and (60), we obtain:

$$\begin{array}{c}\hfill \begin{array}{cc}\hfill k\sum _{n=0}^{K-1}〈{\partial }_t{Q}_{n-\frac{1}{2}}^{(\alpha )}(\Delta \xi ),{\xi }^{n-\frac{1}{2}}〉& =k^{\alpha +1}〈\sum _{n=0}^{K-1}\sum _{p=0}^n{\partial }_t{q}_{n-p+\frac{1}{2}}^{(\alpha )}\Delta {\xi }^p,{\xi }^{n-\frac{1}{2}}〉\hfill \\ & \le k\sum _{p=0}^{K-1}{|\nabla {\xi }^{p-\frac{1}{2}}|}_D^2.\hfill \end{array}\end{array}$$

(71)

From (22) in Lemma 2, we note that:

$$\begin{array}{c}\hfill \begin{array}{c}\hfill -\langle \Delta {\xi }^K,{\xi }^K\rangle \ge {\parallel \nabla {\xi }^K\parallel }^2\ge 0,\end{array}\end{array}$$

(72)

and,

$$\begin{array}{c}\hfill \begin{array}{c}\hfill \frac{{\lambda }^2{k}^4}{16}k\sum _{n=0}^K\parallel \frac{{\partial }^2{\widehat{\xi }}^{n+\frac{1}{2}}}{\partial x\partial y}{\parallel }^2\ge 0.\end{array}\end{array}$$

(73)

Dropping two positive terms (72) and (73) in (66), we have:

$$\begin{array}{c}\hfill \begin{array}{cc}& |{\widehat{\xi }}^K{|}_D^2+|\nabla {\xi }^K{|}_D^2\le Ck\sum _{n=0}^{K-1}(|{\widehat{\xi }}^{n+\frac{1}{2}}{|}_D^2+\frac{\lambda }{2}|\nabla {\xi }^{n+\frac{1}{2}}{|}_D)+k\sum _{n=0}^{K-1}(|{\zeta }^n{|}_D^2+|\nabla {\rho }^n{|}_D)\hfill \\ & +Ck\sum _{p=0}^{K-1}{|\nabla {\xi }^{p-\frac{1}{2}}|}_D^2.\hfill \end{array}\end{array}$$

(74)

Then, (74) reduces to:

$$\begin{array}{c}\hfill \begin{array}{cc}\hfill |{\widehat{\xi }}^K{|}_D^2+(1-\frac{C\lambda }{4}){|\nabla {\xi }^K|}_D^2& \le Ck\sum _{n=0}^{K-1}(|{\widehat{\xi }}^{n+\frac{1}{2}}{|}_D^2+\frac{\lambda }{2}|\nabla {\xi }^{n+\frac{1}{2}}{|}_D)\hfill \\ & +k\sum _{n=0}^{K-1}(|{\zeta }^n{|}_D^2+|\nabla {\rho }^n{|}_D)+C{k}^{\alpha }\sum _{p=0}^{K-1}{q}_{K-p+\frac{1}{2}}^{(\alpha )}{|\nabla {\xi }^{p-\frac{1}{2}}|}_D^2.\hfill \end{array}\end{array}$$

(75)

Using the discrete Grönwall inequality, we obtain:

$$\begin{array}{c}\hfill \begin{array}{c}\hfill |{\widehat{\xi }}^n{|}_D\le Ck\sum _{n=0}^{K-1}(|{\zeta }^n{|}_D^2+|\nabla {\rho }^n{|}_D^2),\end{array}\end{array}$$

(76)

for k sufficiently small. We use the following bounds:

$$\begin{array}{c}\hfill \begin{array}{c}\hfill \parallel {\left(\frac{\partial \varphi }{\partial t}\right)}_{n+\frac{1}{2}}-{\partial }_t{\varphi }_n\parallel \le C{k}^2,\end{array}\end{array}$$

(77)

and,

$$\begin{array}{c}\hfill \begin{array}{c}\hfill \parallel \frac{{\partial }^4{\delta }_t{\varphi }^n}{\partial x^2\partial y^2}\parallel \le C.\end{array}\end{array}$$

(78)

Then, we obtain:

$$\begin{array}{c}\hfill \begin{array}{c}\hfill k\sum _{n=0}^{K-1}{|{\zeta }^n|}_D^2\le C{k}^2+k\sum _{n=0}^{K-1}{\partial }_t{\widehat{\eta }}^n+\frac{{\lambda }^2{k}^4}{16}k\sum _{n=0}^{K-1}\frac{{\partial }^4{\delta }_t{\widehat{\eta }}^n}{\partial x^2\partial y^2}+E{(\eta )}^{(\alpha )}(t_{n+1}).\end{array}\end{array}$$

(79)

Since bound of ${\eta }^n$ is known from Lemmas 4 and 5, with ${\widehat{\eta }}^n$ replaced by ${\eta }^n$, we obtain:

$$\begin{array}{c}\hfill \begin{array}{c}\hfill k\sum _{n=0}^{K-1}{|{\zeta }^n|}_D^2\le C(k^4+h^{2r+2}),\end{array}\end{array}$$

(80)

then it holds that,

$$\begin{array}{c}\hfill \begin{array}{cc}\hfill k\sum _{n=0}^{K-1}{|\nabla {\rho }^n|}_D^2& \le C\left[k^4+k\sum _{n=0}^{K-1}|{\widehat{\eta }}^{n+1/2}{|}_D^2+k\sum _{n=0}^{K-1}{|{\partial }_t{\eta }^n|}_D^2\right]\hfill \\ & \le C\left[k^4+\parallel \frac{\partial \eta }{\partial t}{\parallel }_{L^2(H_0^1)}^2+\parallel \frac{\partial \eta }{\partial t}{\parallel }_{L^{\infty }(H_0^1)}^2\right].\hfill \end{array}\end{array}$$

(81)

Then we have:

$$\begin{array}{c}\hfill \begin{array}{c}\hfill k\sum _{n=0}^K{|\nabla {\rho }^n|}_D\le C(k^2+h^{r+1}).\end{array}\end{array}$$

(82)

With (80) and (82) in (76), we arrive at the desired result. □

Note that Theorem 2 is first proved in this paper, which shows the competitiveness of OSC methods.

5. Numerical Results

This section considers two numerical examples to verify that the proposed strategy is accurate, applicable, and effective. In the numerical simulations, we use the space of the piecewise Hermite bicubics, for which $r=3$, with the standard value and scaled slope basis functions on identical uniform partitions in both the x and y variable with $N_x=N_y$. The IC is approximated by choosing $U^0={\widehat{U}}^0,$ the OSC elliptic projection of $u^0$. To check the proposed method, we define temporal and spatial convergence orders:

$$\begin{array}{cc}\hfill & C_k={log}_2\left(\frac{\parallel L^{\infty }(2N)\parallel }{\parallel L^{\infty }(N)\parallel }\right),C_h={log}_2\left(\frac{\parallel L^{\infty }(2N_x)\parallel }{\parallel L^{\infty }(N_x)\parallel }\right),\hfill \end{array}$$

respectively, where $L^{\infty }=\underset{1\le j\le N-1}{max}|U(x_j,T)-u(x_j,T)|$ stands for the maximum norm (i.e., absolute error). Furthermore, all numerical results are calculated by MATLAB R2016b on a personal computer with 4GB RAM.

Example 1.

Let us consider the hyperbolic integrodifferential Equations (1)–(3) for $(x,y)\in (0,1)\times (0,1)$ at $T=1$ with the exact solution:

$$\begin{array}{c}\hfill u(x,y,t)=t^{\alpha +3}sin(\pi x)sin(\pi y),\end{array}$$

thus, the associated forcing term is selected as:

$$\begin{array}{c}\hfill f(x,y,t)=sin(\pi y)sin(\pi x)\left(2{\pi }^2{t}^{2\alpha +3}\frac{\Gamma (\alpha +4)}{\Gamma (2\alpha +4)}+2{\pi }^2{t}^{\alpha +3}+(\alpha +3)(\alpha +2)t^{\alpha +1}\right).\end{array}$$

The ICs $v(x,y)=\nu (x,y)=0$ are obtained from the exact solution.

Table 1. The maximum norm errors $L^{\infty }$ and spatial convergence rates $C_h$ with temporal step $k=\frac{1}{1000}$.

${\mathit{N}}_{\mathit{x}}={\mathit{N}}_{\mathit{y}}$	$\mathit{\alpha }=0.25$			$\mathit{\alpha }=0.5$			$\mathit{\alpha }=0.75$
	${\mathit{L}}^{\mathbf{\infty }}$	${\mathit{C}}_{\mathit{h}}$	CPU(s)	${\mathit{L}}^{\mathbf{\infty }}$	${\mathit{C}}_{\mathit{h}}$	CPU(s)	${\mathit{L}}^{\mathbf{\infty }}$	${\mathit{C}}_{\mathit{h}}$	CPU(s)
2	$2.9832\times {10}^{-3}$	−	$1.0776$	$2.8761\times {10}^{-3}$	−	$0.9874$	$2.7865\times {10}^{-3}$	−	$0.8468$
4	$1.8401\times {10}^{-4}$	$4.0190$	$0.9796$	$1.7336\times {10}^{-4}$	$4.0523$	$2.2420$	$1.6336\times {10}^{-4}$	$4.0923$	$0.8571$
8	$8.7871\times {10}^{-6}$	$4.3883$	$1.3376$	$7.9503\times {10}^{-6}$	$4.4466$	1.1716	$6.8184\times {10}^{-6}$	$4.5825$	$1.2096$
16	$5.1878\times {10}^{-7}$	$4.0822$	$3.7028$	$5.7066\times {10}^{-7}$	$3.8003$	3.5713	$8.1472\times {10}^{-7}$	$3.0651$	$3.5900$

displays the maximum norm errors $L^{\infty }$, spatial convergence orders $C_h,$ and CPU run times (in seconds) for temporal steps as $k=\frac{1}{1000}$, which accords to the fourth-order convergence in the space variable.

Table 2. The maximum norm errors $L^{\infty }$ and temporal convergence rates $C_k$ with $h=k$.

N	$\mathit{\alpha }=0.25$			$\mathit{\alpha }=0.5$			$\mathit{\alpha }=0.75$
	${\mathit{L}}^{\mathbf{\infty }}$	${\mathit{C}}_{\mathit{k}}$	CPU(s)	${\mathit{L}}^{\mathbf{\infty }}$	${\mathit{C}}_{\mathit{k}}$	CPU(s)	${\mathit{L}}^{\mathbf{\infty }}$	${\mathit{C}}_{\mathit{k}}$	CPU(s)
16	$4.0700\times {10}^{-3}$	−	$3.4033\times {10}^{-2}$	$2.9282\times {10}^{-3}$	−	$2.8789\times {10}^{-2}$	$2.9780\times {10}^{-3}$	−	$2.1668\times {10}^{-2}$
32	$1.0134\times {10}^{-3}$	$2.0059$	$1.7277\times {10}^{-2}$	$7.9555\times {10}^{-4}$	$1.8800$	$1.7437\times {10}^{-2}$	$9.2241\times {10}^{-4}$	$1.6909$	$1.8405\times {10}^{-2}$
64	$3.6099\times {10}^{-4}$	$1.4891$	$4.1779\times {10}^{-2}$	$1.9709\times {10}^{-4}$	$2.0131$	$3.8693\times {10}^{-2}$	$2.4596\times {10}^{-4}$	$1.9070$	$3.1814\times {10}^{-2}$
128	$5.6850\times {10}^{-5}$	$2.6667$	$1.0558\times {10}^{-1}$	$4.7291\times {10}^{-5}$	$2.0592$	$1.1978\times {10}^{-1}$	$6.2457\times {10}^{-5}$	$1.9775$	$9.5469\times {10}^{-2}$
256	$1.3096\times {10}^{-5}$	$2.1181$	$2.7117\times {10}^{-1}$	$1.1132\times {10}^{-5}$	$2.0868$	$3.1955\times {10}^{-1}$	$1.5427\times {10}^{-5}$	$2.0174$	$2.8704\times {10}^{-1}$

reports the maximum norm errors $L^{\infty }$, temporal convergence orders $C_k,$ and CPU run times (in seconds) for $h=k,$ which supports the second-order accuracy in the time variable. We also give the CPU time for the ADIOSC method to simulate the hyperbolic integrodifferential problem, which shows using a shorter time to calculate desired results. Figure 1 and Figure 2 visually show the second-order temporal and fourth-order spatial accuracy of the fully discrete scheme by fixing related parameters, which is consistent with the theoretical analysis.

Example 2.

We consider the hyperbolic integrodifferential Equations (1)–(3) for $(x,y)\in (0,1)\times (0,1)$ at $T=1$ with the exact solution:

$$\begin{array}{c}\hfill u(x,y,t)=\left(t+t^{\alpha +3}\right)sin(\pi x)sin(\pi y),\end{array}$$

so that the associated source term is obtained by:

$$\begin{array}{cc}\hfill f(x,y,t)=sin(\pi y)sin(\pi x)& (\frac{2{\pi }^2{t}^{2\alpha +3}\Gamma (\alpha +4)}{\Gamma (2\alpha +4)}+2{\pi }^2{t}^{\alpha +3}\hfill \\ & +2{\pi }^{2}t+2{\pi }^2{t}^{\alpha +1}\frac{\Gamma (2)}{\Gamma (\alpha +2)}+(\alpha +3)(\alpha +2)t^{\alpha +1}).\hfill \end{array}$$

The ICs $v(x,y)=0,\nu (x,y)=\sin (\pi x)\sin (\pi y)$ are computed from the exact solution.

Table 3. The maximum norm errors $L^{\infty }$ and spatial convergence rates $C_h$ with temporal step $k=\frac{1}{2000}$.

${\mathit{N}}_{\mathit{x}}={\mathit{N}}_{\mathit{y}}$	$\mathit{\alpha }=0.25$			$\mathit{\alpha }=0.5$			$\mathit{\alpha }=0.75$
	${\mathit{L}}^{\mathbf{\infty }}$	${\mathit{C}}_{\mathit{h}}$	CPU(s)	${\mathit{L}}^{\mathbf{\infty }}$	${\mathit{C}}_{\mathit{h}}$	CPU(s)	${\mathit{L}}^{\mathbf{\infty }}$	${\mathit{C}}_{\mathit{h}}$	CPU(s)
2	$6.1678\times {10}^{-3}$	−	$4.0305$	$6.7022\times {10}^{-3}$	−	$4.1241$	$6.7399\times {10}^{-3}$	−	$3.9646$
4	$3.9801\times {10}^{-4}$	$3.9539$	$4.3535$	$4.1981\times {10}^{-4}$	$3.9968$	$4.2544$	$4.2331\times {10}^{-4}$	$3.9929$	$4.2929$
8	$2.3045\times {10}^{-5}$	$4.1103$	$5.4502$	$2.2954\times {10}^{-5}$	$4.1929$	$5.4879$	$2.2744\times {10}^{-5}$	$4.2182$	$5.3585$
16	$1.1177\times {10}^{-6}$	$4.3658$	$14.343$	$1.1482\times {10}^{-6}$	$4.3213$	$14.256$	$1.0555\times {10}^{-6}$	$4.4295$	$14.077$

and

Table 4. The maximum norm errors $L^{\infty }$ and temporal convergence rates $C_k$ with $h=k$.

N	$\mathit{\alpha }=0.25$			$\mathit{\alpha }=0.5$			$\mathit{\alpha }=0.75$
	${\mathit{L}}^{\mathbf{\infty }}$	${\mathit{C}}_{\mathit{k}}$	CPU(s)	${\mathit{L}}^{\mathbf{\infty }}$	${\mathit{C}}_{\mathit{k}}$	CPU(s)	${\mathit{L}}^{\mathbf{\infty }}$	${\mathit{C}}_{\mathit{k}}$	CPU(s)
6	$1.0631\times {10}^{-2}$	−	$1.6929\times {10}^{-2}$	$1.3011\times {10}^{-2}$	−	$1.5590\times {10}^{-2}$	$1.6680\times {10}^{-2}$	−	$1.6576\times {10}^{-3}$
12	$7.4165\times {10}^{-3}$	$0.5195$	$2.6292\times {10}^{-3}$	$5.6240\times {10}^{-3}$	$1.2100$	$2.8282\times {10}^{-3}$	$4.4871\times {10}^{-3}$	$1.8943$	$2.8819\times {10}^{-3}$
24	$1.9494\times {10}^{-3}$	$1.9277$	$4.1437\times {10}^{-3}$	$1.6935\times {10}^{-3}$	$1.7316$	$4.3295\times {10}^{-3}$	$1.6529\times {10}^{-3}$	$1.4408$	$4.3514\times {10}^{-3}$
48	$3.9667\times {10}^{-4}$	$2.2970$	$7.9410\times {10}^{-3}$	$3.9026\times {10}^{-4}$	$2.1175$	$9.5093\times {10}^{-3}$	$4.3759\times {10}^{-4}$	$1.9173$	$8.7784\times {10}^{-3}$
96	$5.7792\times {10}^{-5}$	$2.7790$	$1.8284\times {10}^{-2}$	$5.6877\times {10}^{-5}$	$2.7785$	$1.9089\times {10}^{-2}$	$7.2056\times {10}^{-5}$	$2.6024$	$1.8466\times {10}^{-2}$

show the maximum norm errors $L^{\infty }$, spatial and temporal convergence, and CPU times (in seconds) of the ADIOSC technique for $\alpha =0.25,0.5,0.75$. Table 3 lists the maximum norm errors $L^{\infty }$, spatial convergence orders $C_h$, and CPU run times (in seconds) for the temporal step as $k=\frac{1}{2000}$, which accords to the fourth-order convergence in the space variable. Table 4 discusses the maximum norm errors $L^{\infty }$, temporal convergence orders $C_k$, and CPU run times (in seconds) for $h=k,$ which supports the second-order accuracy in the time variable. We also present the CPU running time of the ADIOSC scheme to simulate the hyperbolic integrodifferential problem, which indicates using a shorter time to compute a desired result.

6. Conclusions

This paper formulated an ADIOSC numerical method for the hyperbolic integrodifferential equation with a weakly singular kernel in two-dimension domain, from which we have used the second-order convolution quadrature. It was shown that the proposed method is stable and convergent with order two for time, and order four for space. Numerical results illustrated efficiency and performance of the proposed ADIOSC method. In our future work, we will consider compatible wavelet techniques, such as the Sinc methods [43,44,45,46] for solving Problems (1)–(3), which may obtain the exponential convergence accuracy in the space direction.