A Compact Difference Method for Viscoelastic Plate Vibration Equation

Cailian Wu, Congcong Wei, Ailing Zhu^*

School of Mathematics and Statistics, Shandong Normal University, Jinan, China.

DOI: 10.4236/eng.2021.1311045                

Abstract

In this paper, a fourth-order viscoelastic plate vibration equation is transformed into a set of two second-order differential equations by introducing an intermediate variable. A three-layer compact difference scheme for the initial-boundary value problem of the viscoelastic plate vibration equation is established. Then the stability and convergence of the difference scheme are analyzed by the energy method, and the convergence order is . Finally, some numerical examples are given of which results verify the accuracy and validity of the scheme.

Keywords

Viscoelastic Plate Vibration Equation, Compact Difference Method, Stability, Convergence

1. Introduction

Consider the following initial-boundary value problem of the viscoelastic plate vibration equation

$\{\begin{array}{ll}(\text{a})\mathrm{<mtext>\; \hspace{0.17em}\; </mtext>}\mu \frac{\partial u}{\partial t}+\frac{{\partial }^{2}u}{\partial t^2}+a{\Delta }^{2}u+ku=f\left(x,y,t\right),\hfill & \left(x,y,t\right)\in \Omega \times \left(0,T\right],\hfill \\ (\text{b})\mathrm{<mtext>\; \hspace{0.17em}\; </mtext>}{u|}_{t=0}={\phi }_1\left(x,y\right),{\frac{\partial u}{\partial t}|}_{t=0}={\phi }_2\left(x,y\right),\hfill & \left(x,y\right)\in \Omega ,\hfill \\ (\text{c})\mathrm{<mtext>\; \hspace{0.17em}\; </mtext>}{u|}_{x=0}={\psi }_1\left(y,t\right),{u|}_{x=b}={\psi }_2\left(y,t\right),\hfill & \hfill \\ \mathrm{<mtext>\; \hspace{0.17em}\; </mtext>\; <mtext>\; \hspace{0.17em}\; </mtext>\; <mtext>\; \hspace{0.17em}\; </mtext>\; <mtext>\; \hspace{0.17em}\; </mtext>\; <mtext>\; \hspace{0.17em}\; </mtext>}{u|}_{y=0}={\psi }_3\left(x,t\right),{u|}_{y=c}={\psi }_4\left(x,t\right),\hfill & t\in \left(0,T\right],\hfill \\ (\text{d})\mathrm{<mtext>\; \hspace{0.17em}\; </mtext>}{\Delta u|}_{x=0}=g_1\left(y,t\right),{\Delta u|}_{x=b}=g_2\left(y,t\right),\hfill & \hfill \\ \mathrm{<mtext>\; \hspace{0.17em}\; </mtext>\; <mtext>\; \hspace{0.17em}\; </mtext>\; <mtext>\; \hspace{0.17em}\; </mtext>\; <mtext>\; \hspace{0.17em}\; </mtext>\; <mtext>\; \hspace{0.17em}\; </mtext>}{\Delta u|}_{y=0}=g_3\left(x,t\right),{\Delta u|}_{y=c}=g_4\left(x,t\right),\hfill & t\in \left(0,T\right],\hfill \end{array}$ (1.1)

where $\Omega =\left(\mathrm{0,}b\right)\times \left(\mathrm{0,}c\right)\subset R^2$, ${\Delta }^{2}u=\frac{{\partial }^{4}u}{\partial x^4}+2\frac{{\partial }^{4}u}{\partial x^2\partial y^2}+\frac{{\partial }^{4}u}{\partial y^4}$, $\mu \left(\mu >0\right)$ is the viscosity coefficient, $k\left(k>0\right)$ is the elasticity coefficient, $a=\frac{D}{\rho l}$, D, $\rho$ and l are bending rigidity, board density and sheet thickness, respectively, and are constants. The source term $f\left(x\mathrm{,}y\mathrm{,}t\right)$ is known smooth functions, ${\phi }_1\left(x\mathrm{,}y\right)$ and ${\phi }_2\left(x\mathrm{,}y\right)$ is initial value functions, ${\psi }_i\left(y\mathrm{,}t\right)$ and $g_i\left(y\mathrm{,}t\right)$ $\left(i=1,2,3,4\right)$ are boundary value functions.

Vibration is a common physical phenomenon, but vibration can cause serious harm in some cases, such as damage to the mechanical structure and buildings, noise that causes serious environmental pollution and so on. Studying the vibration of structural systems such as beams and plates is of important significance in lowering the damage rate of the mechanical structure, increasing safety and reducing the maintenance cost. And the use of viscoelastic damping can reduce vibration hazards. Viscoelastic plate vibration equations can describe and solve bridge and aerospace vibration and noise control issues, etc., theoretical analysis and numerical calculations have attracted more and more attention of scientific research and engineering personnel. In [1] [2], Cheng et al. investigate Rayleigh wave propagation and the solvability of fractional order damped vibration in viscoelastic media. In [3], the plate vibration equation is derived in detail. In [4], Wang proves existence and uniqueness of the local solution refer to the fourth-order damped wave equation. The compact finite difference method which usually approximates the second-order derivatives has the advantages of using fewer and compactly arranged nodes to achieve high precision, avoiding introducing virtual nodes to handle boundary conditions easier and keeping desirable tridiagonal nature of the finite difference equation. In [5] [6] [7] [8], the compact finite difference scheme for wave equations is developed. Deng et al. proposed a fourth-order compact alternating direction implicit format for solving two-dimensional linear hyperbolic equations in [9]. For fourth-order parabolic equations, a compact finite difference scheme is established by introducing intermediate variables in [10]. In [11] [12] [13] [14] [15], for a class of two-dimensional fourth-order differential equations such as plate vibration equation, a compact finite difference scheme with fourth-order spatial accuracy is developed by introducing a second-order derivative as an intermediate variable. However, we have not seen any articles that solve the initial-boundary value problem of viscoelastic plate vibration equation by compact difference method until now. The goal of this paper is to construct a compact implicit difference scheme for the problem (1.1).

The outline of this paper is as follows. In Section 2, the fourth-order viscoelastic plate vibration equation is transformed into a second-order system of equations. The space derivative terms of the equations are discretized by the fourth-order compact difference scheme, while the time derivative terms are discretized by the second-order central difference scheme. Finally, the compact implicit difference scheme of (1.1) is constructed. In Section 3, we prove that the presented compact difference scheme is convergent and unconditionally stable by the energy method. In Section 4, some numerical examples are given to verify the accuracy of the scheme, which show that the scheme has high practicability.

2. Compact Difference Scheme

In this section, we develop a compact difference scheme of the problem (1.1), and introduce intermediate function $v=-a\Delta u$, (1.1) is equivalent to

$\{\begin{array}{ll}(\text{a})\mathrm{<mtext>\; \hspace{0.17em}\; </mtext>}\mu \frac{\partial u}{\partial t}+\frac{{\partial }^{2}u}{\partial t^2}-\Delta v+ku=f\left(x,y,t\right),\hfill & \left(x,y,t\right)\in \Omega \times \left(0,T\right],\hfill \\ (\text{b})\mathrm{<mtext>\; \hspace{0.17em}\; </mtext>}a\Delta u+v=0,\mathrm{<mtext>\; \hspace{0.17em}\; </mtext>\; <mtext>\; \hspace{0.17em}\; </mtext>\; <mtext>\; \hspace{0.17em}\; </mtext>}\hfill & \left(x,y,t\right)\in \Omega \times \left(0,T\right],\hfill \\ (\text{c})\mathrm{<mtext>\; \hspace{0.17em}\; </mtext>}{u|}_{t=0}={\phi }_1\left(x,y\right),{\frac{\partial u}{\partial t}|}_{t=0}={\phi }_2\left(x,y\right),\mathrm{<mtext>\; \hspace{0.17em}\; </mtext>\; <mtext>\; \hspace{0.17em}\; </mtext>}\hfill & \left(x,y\right)\in \Omega ,\hfill \\ (\text{d})\mathrm{<mtext>\; \hspace{0.17em}\; </mtext>}{u|}_{x=0}={\psi }_1\left(y,t\right),{u|}_{x=b}={\psi }_2\left(y,t\right),\hfill & \hfill \\ \mathrm{<mtext>\; \hspace{0.17em}\; </mtext>\; <mtext>\; \hspace{0.17em}\; </mtext>\; <mtext>\; \hspace{0.17em}\; </mtext>\; <mtext>\; \hspace{0.17em}\; </mtext>\; <mtext>\; \hspace{0.17em}\; </mtext>}{u|}_{y=0}={\psi }_3\left(x,t\right),{u|}_{y=c}={\psi }_4\left(x,t\right),\hfill & t\in \left(0,T\right],\hfill \\ (\text{e})\mathrm{<mtext>\; \hspace{0.17em}\; </mtext>}{v|}_{x=0}=-a{g}_1\left(y,t\right),{v|}_{x=b}=-a{g}_2\left(y,t\right),\hfill & \hfill \\ \mathrm{<mtext>\; \hspace{0.17em}\; </mtext>\; <mtext>\; \hspace{0.17em}\; </mtext>\; <mtext>\; \hspace{0.17em}\; </mtext>\; <mtext>\; \hspace{0.17em}\; </mtext>\; <mtext>\; \hspace{0.17em}\; </mtext>}{v|}_{y=0}=-a{g}_3\left(x,t\right),{v|}_{y=c}=-a{g}_4\left(x,t\right),\mathrm{<mtext>\; \hspace{0.17em}\; </mtext>\; <mtext>\; \hspace{0.17em}\; </mtext>}\hfill & t\in \left(0,T\right].\hfill \end{array}$ (2.1)

Take positive integers $M_x$, $M_y$ and N, let $h_x=\frac{b}{M_x}$, $h_y=\frac{c}{M_y}$, $h=\max \left\{h_x\mathrm{,}{h}_y\right\}$, and $\tau =\frac{T}{N}$. Denote $x_i=i{h}_x\left(i=0,1,\cdots ,M_x\right)$, $y_j=j{h}_y\left(j=0,1,\cdots ,M_y\right)$, $t_n=n\tau \left(n=0,1,\cdots ,N\right)$. Defining the grid function $u_{ij}^n$ on node $\left(x_i\mathrm{,}{y}_j\mathrm{,}{t}_n\right)$, we introduce the following notations

$\begin{array}{l}{u}_{ij}^{n+\frac{1}{2}}=\frac{1}{2}\left(u_{ij}^n+u_{ij}^{n+1}\right),{\delta }_{\hat{t}}{u}_{ij}^n=\frac{1}{2\tau }\left(u_{ij}^{n+1}-u_{ij}^{n-1}\right),\\ {\delta }_t{u}_{ij}^{n+\frac{1}{2}}=\frac{1}{\tau }\left(u_{ij}^{n+1}-u_{ij}^n\right),{\delta }_t^2{u}_{ij}^n=\frac{1}{{\tau }^2}\left(u_{ij}^{n+1}-2u_{ij}^n+u_{ij}^{n-1}\right),\\ {\delta }_x^2{u}_{ij}^n=\frac{1}{h_x^2}\left(u_{i+1,j}^n-2u_{ij}^n+u_{i-1,j}^n\right),{\delta }_y^2{u}_{ij}^n=\frac{1}{h_y^2}\left(u_{i,j+1}^n-2u_{ij}^n+u_{i,j-1}^n\right),\end{array}$

and the compact operators

$A_x{u}_{ij}=\{\begin{array}{ll}\frac{1}{12}\left(u_{i-1,j}+10u_{ij}+u_{i+1,j}\right),\hfill & 1\le i\le M_x-1,0\le j\le M_y,\hfill \\ u_{ij},\hfill & i=0,M_x,0\le j\le M_y,\hfill \end{array}$

$A_y{u}_{ij}=\{\begin{array}{ll}\frac{1}{12}\left(u_{i,j-1}+10u_{ij}+u_{i,j+1}\right),\hfill & 1\le j\le M_y-1,0\le i\le M_x,\hfill \\ u_{ij},\hfill & j=0,M_y,0\le i\le M_x.\hfill \end{array}$

Let $A_h=A_x{A}_y$, $B_h=A_y{\delta }_x^2+A_x{\delta }_y^2$.

Lemma 2.1 [16] If $g\left(x\right)\in C^6\left[c-h\mathrm{,}c+h\right]$ , it holds

$\begin{array}{l}\frac{1}{12}\left[g^{″}\left(c-h\right)+10g^{″}\left(c\right)+g^{″}\left(c+h\right)\right]\\ =\frac{1}{h^2}\left[g\left(c+h\right)-2g\left(c\right)+g\left(c-h\right)\right]+\frac{h^4}{240}{g}^{\left(6\right)}\left(\xi \right),\end{array}$

where $\xi \in \left(c-h\mathrm{,}c+h\right)$.

From now on, the derivation of the compact difference scheme is presented carefully.

Considering the Equation (2.1a), for $1\le i\le M_x-1$, $1\le j\le M_y-1$, $1\le n\le N-1$, on the node $\left(x_i\mathrm{,}{y}_j\mathrm{,}{t}_n\right)$, we have

$\begin{array}{l}\mu \frac{\partial u}{\partial t}\left(x_i,y_j,t_n\right)+\frac{{\partial }^{2}u}{\partial t^2}\left(x_i,y_j,t_n\right)-\frac{{\partial }^{2}v}{\partial x^2}\left(x_i,y_j,t_n\right)-\frac{{\partial }^{2}v}{\partial y^2}\left(x_i,y_j,t_n\right)+ku\left(x_i,y_j,t_n\right)\\ =f\left(x_i,y_j,t_n\right).\end{array}$ (2.2)

From the Taylor expansion, we obtain

$\frac{\partial u}{\partial t}\left(x_i,y_j,t_n\right)={\delta }_{\hat{t}}{u}_{ij}^n-\frac{{\tau }^2}{6}\frac{{\partial }^{3}u}{\partial t^3}\left(x_i,y_j,t_n\right)+O\left({\tau }^4\right),$

$\frac{{\partial }^{2}u}{\partial t^2}\left(x_i,y_j,t_n\right)={\delta }_t^2{u}_{ij}^n-\frac{{\tau }^2}{12}\frac{{\partial }^{4}u}{\partial t^4}\left(x_i,y_j,t_n\right)+O\left({\tau }^4\right),$

$\begin{array}{c}\frac{{\partial }^{2}v}{\partial x^2}\left(x_i,y_j,t_n\right)=\frac{1}{4}\left[\frac{{\partial }^{2}v}{\partial x^2}\left(x_i,y_j,t_{n+1}\right)+2\frac{{\partial }^{2}v}{\partial x^2}\left(x_i,y_j,t_n\right)+\frac{{\partial }^{2}v}{\partial x^2}\left(x_i,y_j,t_{n-1}\right)\right]\\ -\frac{{\tau }^2}{4}\frac{{\partial }^{4}v}{\partial x^2\partial t^2}\left(x_i,y_j,t_n\right)+O\left({\tau }^4\right),\end{array}$

$\begin{array}{c}\frac{{\partial }^{2}v}{\partial y^2}\left(x_i\mathrm{,}{y}_j\mathrm{,}{t}_n\right)=\frac{1}{4}\left[\frac{{\partial }^{2}v}{\partial y^2}\left(x_i\mathrm{,}{y}_j\mathrm{,}{t}_{n+1}\right)+2\frac{{\partial }^{2}v}{\partial y^2}\left(x_i\mathrm{,}{y}_j\mathrm{,}{t}_n\right)+\frac{{\partial }^{2}v}{\partial y^2}\left(x_i\mathrm{,}{y}_j\mathrm{,}{t}_{n-1}\right)\right]\\ -\frac{{\tau }^2}{4}\frac{{\partial }^{4}v}{\partial y^2\partial t^2}\left(x_i\mathrm{,}{y}_j\mathrm{,}{t}_n\right)+O\left({\tau }^4\right)\mathrm{.}\end{array}$

Substituting the above four formulas into (2.2) and replacing $u_{ij}^n$ with $\frac{1}{2}\left(u_{ij}^{n+1}+u_{ij}^{n-1}\right)$, we get

$\begin{array}{l}\mu {\delta }_{\hat{t}}{u}_{ij}^n+{\delta }_t^2{u}_{ij}^n-\frac{1}{4}\left[\frac{{\partial }^{2}v}{\partial x^2}\left(x_i\mathrm{,}{y}_j\mathrm{,}{t}_{n+1}\right)+2\frac{{\partial }^{2}v}{\partial x^2}\left(x_i\mathrm{,}{y}_j\mathrm{,}{t}_n\right)+\frac{{\partial }^{2}v}{\partial x^2}\left(x_i\mathrm{,}{y}_j\mathrm{,}{t}_{n-1}\right)\right]\\ -\frac{1}{4}\left[\frac{{\partial }^{2}v}{\partial y^2}\left(x_i\mathrm{,}{y}_j\mathrm{,}{t}_{n+1}\right)+2\frac{{\partial }^{2}v}{\partial y^2}\left(x_i\mathrm{,}{y}_j\mathrm{,}{t}_n\right)+\frac{{\partial }^{2}v}{\partial y^2}\left(x_i\mathrm{,}{y}_j\mathrm{,}{t}_{n-1}\right)\right]+\frac{k}{2}\left(u_{ij}^{n+1}+u_{ij}^{n-1}\right)\\ =f\left(x_i\mathrm{,}{y}_j\mathrm{,}{t}_n\right)+{\tau }^{2}g\left(x_i\mathrm{,}{y}_j\mathrm{,}{t}_n\right)+O\left({\tau }^4\right)\mathrm{,}\end{array}$ (2.3)

where

$\begin{array}{c}g\left(x_i\mathrm{,}{y}_j\mathrm{,}{t}_n\right)=\frac{1}{6}\frac{{\partial }^{3}u}{\partial t^3}\left(x_i\mathrm{,}{y}_j\mathrm{,}{t}_n\right)+\frac{1}{12}\frac{{\partial }^{4}u}{\partial t^4}\left(x_i\mathrm{,}{y}_j\mathrm{,}{t}_n\right)\\ +\frac{1}{4}\frac{{\partial }^{4}v}{\partial x^2\partial t^2}\left(x_i\mathrm{,}{y}_j\mathrm{,}{t}_n\right)+\frac{1}{4}\frac{{\partial }^{4}v}{\partial y^2\partial t^2}\left(x_i\mathrm{,}{y}_j\mathrm{,}{t}_n\right)\mathrm{.}\end{array}$

Performing the compact operator $A_x$, $A_y$ on both sides of (2.3), we get

$\begin{array}{l}\mu A_h{\delta }_{\hat{t}}{u}_{ij}^n+A_h{\delta }_t^2{u}_{ij}^n-\frac{1}{4}{A}_h\left[\frac{{\partial }^{2}v}{\partial x^2}\left(x_i,y_j,t_{n+1}\right)+2\frac{{\partial }^{2}v}{\partial x^2}\left(x_i,y_j,t_n\right)+\frac{{\partial }^{2}v}{\partial x^2}\left(x_i,y_j,t_{n-1}\right)\right]\\ -\frac{1}{4}{A}_h\left[\frac{{\partial }^{2}v}{\partial y^2}\left(x_i,y_j,t_{n+1}\right)+2\frac{{\partial }^{2}v}{\partial y^2}\left(x_i,y_j,t_n\right)+\frac{{\partial }^{2}v}{\partial y^2}\left(x_i,y_j,t_{n-1}\right)\right]+\frac{k}{2}{A}_h\left(u_{ij}^{n+1}+u_{ij}^{n-1}\right)\\ =A_h{f}_{ij}^n+{\tau }^2{A}_h{g}_{ij}^n+O\left({\tau }^4\right).\end{array}$ (2.4)

According to Lemma 2.1, and using the commutativity of $A_x{A}_y$, we have

$\mu A_h{\delta }_{\hat{t}}{u}_{ij}^n+A_h{\delta }_t^2{u}_{ij}^n-\frac{1}{4}{B}_h\left(v_{ij}^{n+1}+2v_{ij}^n+v_{ij}^{n-1}\right)+\frac{k}{2}{A}_h\left(u_{ij}^{n+1}+u_{ij}^{n-1}\right)=A_h{f}_{ij}^n+R_{i\mathrm{,}j\mathrm{,}n}^{\left(1\right)}\mathrm{,}$ (2.5)

where

$\begin{array}{c}{R}_{i\mathrm{,}j\mathrm{,}n}^{\left(1\right)}={\tau }^2{A}_x{A}_y{g}_{ij}^n+\frac{h_x^4}{960}{A}_y\left[\frac{{\partial }^{6}v}{\partial x^6}\left(x_i\mathrm{,}{y}_j\mathrm{,}{t}_{n+1}\right)+2\frac{{\partial }^{6}v}{\partial x^6}\left(x_i\mathrm{,}{y}_j\mathrm{,}{t}_n\right)+\frac{{\partial }^{6}v}{\partial x^6}\left(x_i\mathrm{,}{y}_j\mathrm{,}{t}_{n-1}\right)\right]\\ +\frac{h_y^4}{960}{A}_x\left[\frac{{\partial }^{6}v}{\partial y^6}\left(x_i\mathrm{,}{y}_j\mathrm{,}{t}_{n+1}\right)+2\frac{{\partial }^{6}v}{\partial y^6}\left(x_i\mathrm{,}{y}_j\mathrm{,}{t}_n\right)+\frac{{\partial }^{6}v}{\partial y^6}\left(x_i\mathrm{,}{y}_j\mathrm{,}{t}_{n-1}\right)\right]\\ +O\left({\tau }^4+h_x^6+h_y^6\right)\mathrm{.}\end{array}$

Assume u and v are sufficiently smooth, then there exists a positive constant C such that

$\left|R_{i\mathrm{,}j\mathrm{,}n}^{\left(1\right)}\right|\le C\left({\tau }^2+h_x^4+h_y^4\right)\mathrm{.}$

Omitting the small terms $R_{i\mathrm{,}j\mathrm{,}n}^{\left(1\right)}$ from (2.5) and letting $U_{ij}^n$, $V_{ij}^n$ be the approximate solution to $u_{ij}^n$, $v_{ij}^n$, we get difference discrete scheme of (2.1a)

$\mu A_h{\delta }_{\hat{t}}{U}_{ij}^n+A_h{\delta }_t^2{U}_{ij}^n-\frac{1}{4}{B}_h\left(V_{ij}^{n+1}+2V_{ij}^n+V_{ij}^{n-1}\right)+\frac{k}{2}{A}_h\left(U_{ij}^{n+1}+U_{ij}^{n-1}\right)=A_h{f}_{ij}^n\mathrm{.}$ (2.6)

Considering the Equation (2.1b) at node $\left(x_i\mathrm{,}{y}_j\mathrm{,}{t}_n\right)$, for $1\le i\le M_x-1$, $1\le j\le M_y-1$, $1\le n\le N-1$, we have

$a\frac{{\partial }^{2}u}{\partial x^2}\left(x_i\mathrm{,}{y}_j\mathrm{,}{t}_n\right)+a\frac{{\partial }^{2}u}{\partial y^2}\left(x_i\mathrm{,}{y}_j\mathrm{,}{t}_n\right)+v\left(x_i\mathrm{,}{y}_j\mathrm{,}{t}_n\right)=0.$ (2.7)

Similar to (2.5), (2.7) is rewritten in the following form

$a{B}_h{u}_{ij}^n+A_h{v}_{ij}^n=R_{i\mathrm{,}j\mathrm{,}n}^{\left(2\right)}\mathrm{,}$ (2.8)

where

$R_{i,j,n}^{\left(2\right)}=-\frac{h_x^4}{240}{A}_y\frac{{\partial }^{6}u}{\partial x^6}\left(x_i,y_j,t_n\right)-\frac{h_y^4}{240}{A}_x\frac{{\partial }^{6}u}{\partial y^6}\left(x_i,y_j,t_n\right)+O\left(h_x^6+h_y^6\right).$

Assume u and v are sufficiently smooth, then there exists a positive constant C such that

$\left|R_{i\mathrm{,}j\mathrm{,}n}^{\left(2\right)}\right|\le C\left(h_x^4+h_y^4\right)\mathrm{.}$

Omitting the small terms $R_{i\mathrm{,}j\mathrm{,}n}^{\left(2\right)}$ from (2.8) and letting $U_{ij}^n$, $V_{ij}^n$ be the approximate solution to $u_{ij}^n$, $v_{ij}^n$, we get difference discrete scheme of (2.1b)

$a{B}_h{U}_{ij}^n+A_h{V}_{ij}^n=0.$ (2.9)

Considering the initial condition ${u|}_{t=0}={\phi }_1\left(x,y\right)$, we have

$U_{ij}^0=u_{ij}^0={\phi }_{1,ij},V_{ij}^0=v_{ij}^0=-a\Delta {\phi }_{1,ij}.$ (2.10)

In order to get the computation formats for the first time layer, using Taylor expansion and combining (2.1) and (2.10), we obtain

$\begin{array}{l}\begin{array}{l}{u}_{ij}^1\hfill \end{array}=u_{ij}^0+\tau {\left(\frac{\partial u}{\partial t}\right)}_{ij}^0+\frac{{\tau }^2}{2}{\left(\frac{{\partial }^{2}u}{\partial t^2}\right)}_{ij}^0+\frac{{\tau }^3}{6}{\left(\frac{{\partial }^{3}u}{\partial t^3}\right)}_{ij}^0+O\left({\tau }^4\right)\\ ={\phi }_{1,ij}+\tau {\phi }_{2,ij}+\frac{{\tau }^2}{2}\left(f_{ij}^0-a{\Delta }^2{\phi }_{1,ij}-\mu {\phi }_{2,ij}-k{\phi }_{1,ij}\right)\\ +r_{ij}^{\left(1\right)},1\le i\le M_x-1,1\le j\le M_y-1,\end{array}$ (2.11)

and

$\begin{array}{l}{v}_{ij}^1=v_{ij}^0+\tau {\left(\frac{\partial v}{\partial t}\right)}_{ij}^0+\frac{{\tau }^2}{2}{\left(\frac{{\partial }^{2}v}{\partial t^2}\right)}_{ij}^0+O\left({\tau }^3\right)\\ =-a\Delta {\phi }_{1,ij}-a\tau \Delta {\phi }_{2,ij}+r_{ij}^{\left(2\right)},1\le i\le M_x-1,1\le j\le M_y-1,\end{array}$ (2.12)

where

$r_{ij}^{\left(1\right)}=\frac{{\tau }^3}{6}{\left(\frac{{\partial }^{3}u}{\partial t^3}\right)}_{ij}^0+O\left({\tau }^4\right),\mathrm{<mtext>\; \hspace{0.17em}\; </mtext>\; <mtext>\; \hspace{0.17em}\; </mtext>\; <mtext>\; \hspace{0.17em}\; </mtext>\; <mtext>\; \hspace{0.17em}\; </mtext>}{r}_{ij}^{\left(2\right)}=\frac{{\tau }^2}{2}{\left(\frac{{\partial }^{2}v}{\partial t^2}\right)}_{ij}^0+O\left({\tau }^3\right).$

Omitting the small terms $r_{ij}^{\left(1\right)}$ and $r_{ij}^{\left(2\right)}$ from (2.11) and (2.12) respectively, and letting $U_{ij}^n$, $V_{ij}^n$ be the approximate solution to $u_{ij}^n$, $v_{ij}^n$, we get

$\begin{array}{l}{U}_{ij}^1={\phi }_{1,ij}+\tau {\phi }_{2,ij}+\frac{{\tau }^2}{2}\left(f_{ij}^0-a{\Delta }^2{\phi }_{1,ij}-\mu {\phi }_{2,ij}-k{\phi }_{1,ij}\right),\\ V_{ij}^1=-a\Delta {\phi }_{1,ij}-a\tau \Delta {\phi }_{2,ij}.\end{array}$ (2.13)

Considering the boundary value conditions (2.1d), (2.1e), we have

$\begin{array}{l}{U}_{0,j}^n={\psi }_{1,j}^n,\mathrm{<mtext>\; \hspace{0.17em}\; </mtext>\; <mtext>\; \hspace{0.17em}\; </mtext>}{U}_{M_x,j}^n={\psi }_{2,j}^n,\mathrm{<mtext>\; \hspace{0.17em}\; </mtext>\; <mtext>\; \hspace{0.17em}\; </mtext>}{U}_{i,0}^n={\psi }_{3,i}^n,\mathrm{<mtext>\; \hspace{0.17em}\; </mtext>\; <mtext>\; \hspace{0.17em}\; </mtext>}{U}_{i,M_y}^n={\psi }_{4,i}^n,\\ V_{0,j}^n=-a{g}_{1,j}^n,\mathrm{<mtext>\; \hspace{0.17em}\; </mtext>\; <mtext>\; \hspace{0.17em}\; </mtext>}{V}_{M_x,j}^n=-a{g}_{2,j}^n,\mathrm{<mtext>\; \hspace{0.17em}\; </mtext>\; <mtext>\; \hspace{0.17em}\; </mtext>\; <mtext>\; \hspace{0.17em}\; </mtext>}{V}_{i,0}^n=-a{g}_{3,i}^n,\\ V_{i,M_y}^n=-a{g}_{4,i}^n,\mathrm{<mtext>\; \hspace{0.17em}\; </mtext>\; 1}\le i\le M_x-1,1\le j\le M_y-1,0\le n\le N.\end{array}$ (2.14)

Combining (2.6), (2.9), (2.10), (2.13), (2.14), we construct the following compact implicit difference scheme of (2.1)

$\{\begin{array}{l}(\text{a})\mu A_h{\delta }_{\hat{t}}{U}_{ij}^n+A_h{\delta }_t^2{U}_{ij}^n-\frac{1}{4}{B}_h\left(V_{ij}^{n+1}+2V_{ij}^n+V_{ij}^{n-1}\right)+\frac{k}{2}{A}_h\left(U_{ij}^{n+1}+U_{ij}^{n-1}\right)=A_h{f}_{ij}^n,\hfill \\ (\text{b})a{B}_h{U}_{ij}^n+A_h{V}_{ij}^n=0,\hfill \\ (\text{c})U_{ij}^0={\phi }_{1,ij},\mathrm{<mtext>\; \hspace{0.17em}\; </mtext>\; <mtext>\; \hspace{0.17em}\; </mtext>\; <mtext>\; \hspace{0.17em}\; </mtext>\; <mtext>\; \hspace{0.17em}\; </mtext>}{V}_{ij}^0=-a\Delta {\phi }_{1,ij},\hfill \\ (\text{d})U_{ij}^1={\phi }_{1,ij}+\tau {\phi }_{2,ij}+\frac{{\tau }^2}{2}\left(f_{ij}^0-a{\Delta }^2{\phi }_{1,ij}-\mu {\phi }_{2,ij}-k{\phi }_{1,ij}\right),\hfill \\ \mathrm{<mtext>\; \hspace{0.17em}\; </mtext>\; <mtext>\; \hspace{0.17em}\; </mtext>}\mathrm{<mtext>\; \hspace{0.17em}\; </mtext>\; <mtext>\; \hspace{0.17em}\; </mtext>}{V}_{ij}^1=-a\Delta {\phi }_{1,ij}-a\tau \Delta {\phi }_{2,ij},\hfill \\ (\text{e})U_{0,j}^n={\psi }_{1,j}^n,\mathrm{<mtext>\; \hspace{0.17em}\; </mtext>\; <mtext>\; \hspace{0.17em}\; </mtext>}{U}_{M_x,j}^n={\psi }_{2,j}^n,\mathrm{<mtext>\; \hspace{0.17em}\; </mtext>\; <mtext>\; \hspace{0.17em}\; </mtext>}{U}_{i,0}^n={\psi }_{3,i}^n,\mathrm{<mtext>\; \hspace{0.17em}\; </mtext>\; <mtext>\; \hspace{0.17em}\; </mtext>}{U}_{i,M_y}^n={\psi }_{4,i}^n,\hfill \\ \mathrm{<mtext>\; \hspace{0.17em}\; </mtext>}\mathrm{<mtext>\; \hspace{0.17em}\; </mtext>\; <mtext>\; \hspace{0.17em}\; </mtext>\; <mtext>\; \hspace{0.17em}\; </mtext>}{V}_{0,j}^n=-a{g}_{1,j}^n,\mathrm{<mtext>\; \hspace{0.17em}\; </mtext>\; <mtext>\; \hspace{0.17em}\; </mtext>}{V}_{M_x,j}^n=-a{g}_{2,j}^n,\mathrm{<mtext>\; \hspace{0.17em}\; </mtext>\; <mtext>\; \hspace{0.17em}\; </mtext>}{V}_{i,0}^n=-a{g}_{3,i}^n,\mathrm{<mtext>\; \hspace{0.17em}\; </mtext>\; <mtext>\; \hspace{0.17em}\; </mtext>}{V}_{i,M_y}^n=-a{g}_{4,i}^n.\hfill \end{array}$ (2.15)

According to the preceding derivation, the local truncation error of the difference scheme (2.15) is $O\left({\tau }^2+h_x^4+h_y^4\right)$.

3. Stability and Convergence Analysis

In this part, we prove the stability and convergence of difference schemes (2.15) by energy method. We introduce the space $S_h=\left\{u\mathrm{<mo>\; |\; </mo>}{u}_{ij}\in R^{\left(M_x+1\right)\times \left(M_y+1\right)}\mathrm{,0}\le i\le M_x\mathrm{,0}\le j\le M_y\right\}$, $S_{h0}=\left\{u\in S_h,u_{0,j}=u_{M_x,j}=u_{i,0}=u_{i,M_y}=0,0\le i\le M_x,0\le j\le M_y\right\}$.

$\forall u\mathrm{,}v\in S_{h0}$, we define inner products and norms on the space $S_{h0}$

$\left(u,v\right)={\displaystyle \underset{i=1}{\overset{M_x-1}{\sum }}}{\displaystyle \underset{j=1}{\overset{M_y-1}{\sum }}}{u}_{ij}{v}_{ij}{h}_x{h}_y,\left({\delta }_{x}u,{\delta }_{x}v\right)={\displaystyle \underset{i=1}{\overset{M_x}{\sum }}}{\displaystyle \underset{j=1}{\overset{M_y-1}{\sum }}}{\delta }_x{u}_{i-\frac{1}{2},j}{\delta }_x{v}_{i-\frac{1}{2},j}{h}_x{h}_y,$

$\begin{array}{l}\left({\delta }_x^{2}u,{\delta }_x^{2}v\right)={\displaystyle \underset{i=1}{\overset{M_x-1}{\sum }}}{\displaystyle \underset{j=1}{\overset{M_y-1}{\sum }}}{\delta }_x^2{u}_{ij}{\delta }_x^2{v}_{ij}{h}_x{h}_y,\\ \left({\delta }_x{\delta }_{y}u,{\delta }_x{\delta }_{y}v\right)={\displaystyle \underset{i=1}{\overset{M_x}{\sum }}}{\displaystyle \underset{j=1}{\overset{M_y}{\sum }}}{\delta }_x{\delta }_y{u}_{i-\frac{1}{2},j-\frac{1}{2}}{\delta }_x{\delta }_y{v}_{i-\frac{1}{2},j-\frac{1}{2}}{h}_x{h}_y,\end{array}$

${\Vert u\Vert }^2=\left(u,u\right),{\Vert {\delta }_{x}u\Vert }^2=\left({\delta }_{x}u,{\delta }_{x}u\right),{\Vert {\delta }_x^{2}u\Vert }^2=\left({\delta }_x^{2}u,{\delta }_x^{2}u\right),{\left|{\delta }_x{\delta }_{y}u\right|}^2=\left({\delta }_x{\delta }_{y}u,{\delta }_x{\delta }_{y}u\right).$

In order to give stability and convergence analysis of the difference scheme (2.15), the following lemmas are given.

Lemma 3.1 [13] For any grid function $u\mathrm{,}v\in S_{h0}$ , we have

(a) $\left(A_{x}u,v\right)=\left(u,A_{x}v\right)$, $\left(A_{y}u,v\right)=\left(u,A_{y}v\right)$,

(b) $\left(A_{h}u,v\right)=\left(u,A_{h}v\right)$, $\left(B_{h}u,v\right)=\left(u,B_{h}v\right)$,

(c) $\left({\delta }_x^{2}u,v\right)=\left(u,{\delta }_x^{2}v\right)$, $\left({\delta }_y^{2}u,v\right)=\left(u,{\delta }_y^{2}v\right)$.

Lemma 3.2 [13] [17] For any grid function $u\mathrm{,}v\in S_{h0}$ , we have

(a) $\frac{2}{3}{\Vert u\Vert }^2\le \left(A_{x}u\mathrm{,}u\right)\le {\Vert u\Vert }^2$, $\frac{2}{3}{\Vert u\Vert }^2\le \left(A_{y}u\mathrm{,}u\right)\le {\Vert u\Vert }^2$,

(b) $\frac{4}{9}{\Vert u\Vert }^2\le \left(A_{h}u\mathrm{,}u\right)\le {\Vert u\Vert }^2$,

(c) $h_x^2\left(A_y{\delta }_{x}u\mathrm{,}{\delta }_{x}u\right)\le 4\left(A_{y}u\mathrm{,}u\right)$, $h_y^2\left(A_x{\delta }_{y}u\mathrm{,}{\delta }_{y}u\right)\le 4\left(A_{x}u\mathrm{,}u\right)$.

For convenience, we set

$\begin{array}{l}{u}^n={\left(u_{\mathrm{1,1}}^n\mathrm{,}{u}_{\mathrm{2,1}}^n\mathrm{,}\cdots \mathrm{,}{u}_{M_x-\mathrm{1,1}}^n\mathrm{,}{u}_{\mathrm{1,2}}^n\mathrm{,}{u}_{\mathrm{2,2}}^n\mathrm{,}\cdots \mathrm{,}{u}_{M_x-\mathrm{1,2}}^n\mathrm{,}\cdots \mathrm{,}{u}_{\mathrm{1,}{M}_y-1}^n\mathrm{,}{u}_{\mathrm{2,}{M}_y-1}^n\mathrm{,}\cdots \mathrm{,}{u}_{M_x-\mathrm{1,}{M}_y-1}^n\right)}^{\text{T}}\mathrm{,}\\ v^n={\left(v_{\mathrm{1,1}}^n\mathrm{,}{v}_{\mathrm{2,1}}^n\mathrm{,}\cdots \mathrm{,}{v}_{M_x-\mathrm{1,1}}^n\mathrm{,}{v}_{\mathrm{1,2}}^n\mathrm{,}{v}_{\mathrm{2,2}}^n\mathrm{,}\cdots \mathrm{,}{v}_{M_x-\mathrm{1,2}}^n\mathrm{,}\cdots \mathrm{,}{v}_{\mathrm{1,}{M}_y-1}^n\mathrm{,}{v}_{\mathrm{2,}{M}_y-1}^n\mathrm{,}\cdots \mathrm{,}{v}_{M_x-\mathrm{1,}{M}_y-1}^n\right)}^{\text{T}}\mathrm{.}\end{array}$

Theorem 3.1. The difference scheme (2.15) is stable. Let ${\psi }_i=g_i=\mathrm{0\; <mtext>\; \hspace{0.17em}\; </mtext>}\left(i=1,2,3,4\right)$ and $\left\{u^n\mathrm{,}{v}^n\right\}$ refer to the solutions of the difference scheme

$\mu A_h{\delta }_{\hat{t}}{u}_{ij}^n+A_h{\delta }_t^2{u}_{ij}^n-\frac{1}{4}{B}_h\left(v_{ij}^{n+1}+2v_{ij}^n+v_{ij}^{n-1}\right)+\frac{k}{2}{A}_h\left(u_{ij}^{n+1}+u_{ij}^{n-1}\right)=P_{ij}^n\mathrm{,}$ (3.1)

$a{B}_h{u}_{ij}^n+A_h{v}_{ij}^n=Q_{ij}^n,$ (3.2)

$\begin{array}{l}{u}_{ij}^0={\phi }_{1,ij},\mathrm{<mtext>\; \hspace{0.17em}\; </mtext>\; <mtext>\; \hspace{0.17em}\; </mtext>\; <mtext>\; \hspace{0.17em}\; </mtext>\; <mtext>\; \hspace{0.17em}\; </mtext>}{v}_{ij}^0=-a\Delta {\phi }_{1,ij},\\ u_{ij}^1={\phi }_{1,ij}+\tau {\phi }_{2,ij}+\frac{{\tau }^2}{2}\left(f_{ij}^0-a{\Delta }^2{\phi }_{1,ij}-\mu {\phi }_{2,ij}-k{\phi }_{1,ij}\right),\\ v_{ij}^1=-a\Delta {\phi }_{1,ij}-a\tau \Delta {\phi }_{2,ij},\\ u_{0,j}^n=u_{M_x,j}^n=u_{i,0}^n=u_{i,M_y}^n=0,\mathrm{<mtext>\; \hspace{0.17em}\; </mtext>\; <mtext>\; \hspace{0.17em}\; </mtext>}{v}_{0,j}^n=v_{M_x,j}^n=v_{i,0}^n=v_{i,M_y}^n=0.\end{array}$ (3.3)

There is the stability estimate

$\begin{array}{l}{\Vert {\delta }_t{u}^{n+\frac{1}{2}}\Vert }^2+{\Vert u^{n+1}\Vert }^2+{\Vert u^n\Vert }^2+{\Vert v^{n+1}+v^n\Vert }^2\\ \le C\left({\Vert {\delta }_t{u}^{\frac{1}{2}}\Vert }^2+{\Vert v^1+v^0\Vert }^2+{\Vert u^0\Vert }^2+{\Vert u^1\Vert }^2+{\displaystyle {\sum }_{m=1}^n{\Vert P^m\Vert }^2}+{\displaystyle {\sum }_{m=1}^n{\Vert {\delta }_{\hat{t}}{Q}^m\Vert }^2}\right)\mathrm{.}\end{array}$

Proof. Taking the inner product on both sides of (3.1) with $2{\delta }_{\hat{t}}{u}^n$, we obtain

$\begin{array}{l}2\mu \left(A_h{\delta }_{\hat{t}}{u}^n\mathrm{,}{\delta }_{\hat{t}}{u}^n\right)+2\left(A_h{\delta }_t^2{u}^n\mathrm{,}{\delta }_{\hat{t}}{u}^n\right)-\frac{1}{2}\left(B_h\left(v^{n+1}+2v^n+v^{n-1}\right)\mathrm{,}{\delta }_{\hat{t}}{u}^n\right)\\ +k\left(A_h\left(u^{n+1}+u^{n-1}\right)\mathrm{,}{\delta }_{\hat{t}}{u}^n\right)=2\left(P^n\mathrm{,}{\delta }_{\hat{t}}{u}^n\right)\mathrm{.}\end{array}$ (3.4)

From (3.2), we have

$a{B}_h{\delta }_{\hat{t}}{u}^n+A_h{\delta }_{\hat{t}}{v}^n={\delta }_{\hat{t}}{Q}^n.$ (3.5)

Taking the inner product on both sides of (3.5) with $v^{n+1}+2v^n+v^{n-1}$, we obtain

$\begin{array}{l}a\left(B_h{\delta }_{\hat{t}}{u}^n\mathrm{,}{v}^{n+1}+2v^n+v^{n-1}\right)+\left(A_h{\delta }_{\hat{t}}{v}^n\mathrm{,}{v}^{n+1}+2v^n+v^{n-1}\right)\\ =\left({\delta }_{\hat{t}}{Q}^n\mathrm{,}{v}^{n+1}+2v^n+v^{n-1}\right)\mathrm{.}\end{array}$ (3.6)

Multiplying the both sides of (3.4) by a and multiplying the both sides of (3.6) by $\frac{1}{2}$, then summing the results and applying Lemma 3.1, we obtain

$\begin{array}{l}2a\mu \tau \left(A_h{\delta }_{\hat{t}}{u}^n\mathrm{,}{\delta }_{\hat{t}}{u}^n\right)+a\left(A_h{\delta }_t{u}^{n+\frac{1}{2}}\mathrm{,}{\delta }_t{u}^{n+\frac{1}{2}}\right)-a\left(A_h{\delta }_t{u}^{n-\frac{1}{2}}\mathrm{,}{\delta }_t{u}^{n-\frac{1}{2}}\right)\\ +\frac{ka}{2}\left(A_h{u}^{n+1}\mathrm{,}{u}^{n+1}\right)-\frac{ka}{2}\left(A_h{u}^{n-1}\mathrm{,}{u}^{n-1}\right)+\frac{1}{4}\left(A_h\left(v^{n+1}+v^n\right)\mathrm{,}{v}^{n+1}+v^n\right)\\ -\frac{1}{4}\left(A_h\left(v^n+v^{n-1}\right)\mathrm{,}{v}^n+v^{n-1}\right)\\ =2a\tau \left(P^n\mathrm{,}{\delta }_{\hat{t}}{u}^n\right)+\frac{\tau }{2}\left({\delta }_{\hat{t}}{Q}^n\mathrm{,}{v}^{n+1}+2v^n+v^{n-1}\right)\mathrm{.}\end{array}$

Applying Lemma 3.2, we obtain

$\begin{array}{l}8a\mu \tau {\Vert {\delta }_{\hat{t}}{u}^n\Vert }^2+4a{\Vert {\delta }_t{u}^{n+\frac{1}{2}}\Vert }^2+2ka{\Vert u^{n+1}\Vert }^2+{\Vert v^{n+1}+v^n\Vert }^2\\ \le 4a{\Vert {\delta }_t{u}^{n-\frac{1}{2}}\Vert }^2+2ka{\Vert u^{n-1}\Vert }^2+{\Vert v^n+v^{n-1}\Vert }^2+18a\tau \left(P^n\mathrm{,}{\delta }_{\hat{t}}{u}^n\right)\\ \mathrm{<mtext>\; \hspace{0.17em}\; </mtext>\; <mtext>\; \hspace{0.17em}\; </mtext>\; <mtext>\; \hspace{0.17em}\; </mtext>}+\frac{9\tau }{2}\left({\delta }_{\hat{t}}{Q}^n\mathrm{,}{v}^{n+1}+2v^n+v^{n-1}\right)\mathrm{.}\end{array}$ (3.7)

We set

$E^n=4a{\Vert {\delta }_t{u}^{n+\frac{1}{2}}\Vert }^2+2ka\left({\Vert u^{n+1}\Vert }^2+{\Vert u^n\Vert }^2\right)+{\Vert v^{n+1}+v^n\Vert }^2\mathrm{,}$

then

$E^0=4a{\Vert {\delta }_t{u}^{\frac{1}{2}}\Vert }^2+2ka\left({\Vert u^1\Vert }^2+{\Vert u^0\Vert }^2\right)+{\Vert v^1+v^0\Vert }^2\mathrm{.}$

We write (3.7) in the following form

$8a\mu \tau {\Vert {\delta }_{\hat{t}}{u}^n\Vert }^2+E^n\le E^{n-1}+18a\tau \left(P^n\mathrm{,}{\delta }_{\hat{t}}{u}^n\right)+\frac{9\tau }{2}\left({\delta }_{\hat{t}}{Q}^n\mathrm{,}{v}^{n+1}+2v^n+v^{n-1}\right)\mathrm{.}$ (3.8)

Noticing $\mu >0$, $a>0$, summing over n from 1 to m on both sides of (3.8) and applying Cauchy-Schwarz inequality, we get

$\begin{array}{c}{E}^m\le E^0+9a\tau {\displaystyle \underset{n=1}{\overset{m}{\sum }}}{\Vert P^n\Vert }^2+\frac{9\tau }{4}{\displaystyle \underset{n=1}{\overset{m}{\sum }}}{\Vert {\delta }_{\hat{t}}{Q}^n\Vert }^2+18a\tau {\displaystyle \underset{n=0}{\overset{m}{\sum }}}{\Vert {\delta }_t{u}^{n+\frac{1}{2}}\Vert }^2+9\tau {\displaystyle \underset{n=0}{\overset{m}{\sum }}}{\Vert v^{n+1}+v^n\Vert }^2\\ \le E^0+9\tau {\displaystyle \underset{n=0}{\overset{m}{\sum }}}{E}^n+9a\tau {\displaystyle \underset{n=1}{\overset{m}{\sum }}}{\Vert P^n\Vert }^2+\frac{9\tau }{4}{\displaystyle \underset{n=1}{\overset{m}{\sum }}}{\Vert {\delta }_{\hat{t}}{Q}^n\Vert }^2.\end{array}$ (3.9)

Applying discrete Gronwall inequality, we obtain

$E^m\le \exp \left(9n\tau \right)\left(E^0+9a\tau {\displaystyle \underset{n=1}{\overset{m}{\sum }}}{\Vert P^n\Vert }^2+\frac{9\tau }{4}{\displaystyle \underset{n=1}{\overset{m}{\sum }}}{\Vert {\delta }_{\hat{t}}{Q}^n\Vert }^2\right).$

The theorem has been proved.

We present the error estimates for the difference scheme (2.15) in Theorem 3.2.

Setting ${\xi }_{ij}^n=u_{ij}^n-U_{ij}^n$, ${\eta }_{ij}^n=v_{ij}^n-V_{ij}^n$, then subtracting (2.5) and (2.8) from (2.6) and (2.9) at grid point $\left(x_i\mathrm{,}{y}_j\mathrm{,}{t}_n\right)$, respectively, we get the error equations

$\begin{array}{l}\frac{\mu \tau }{2}{A}_h\left({\delta }_t{\xi }^{n+\frac{1}{2}}+{\delta }_t{\xi }^{n-\frac{1}{2}}\right)+A_h\left({\delta }_t{\xi }^{n+\frac{1}{2}}-{\delta }_t{\xi }^{n-\frac{1}{2}}\right)-\frac{\tau }{4}{B}_h\left({\eta }^{n+1}+2{\eta }^n+{\eta }^{n-1}\right)\\ +\frac{k\tau }{2}{A}_h\left({\xi }^{n+1}+{\xi }^{n-1}\right)=\tau R_n^{\left(1\right)},\mathrm{<mtext>\; \hspace{0.17em}\; </mtext>\; <mtext>\; \hspace{0.17em}\; </mtext>\; <mtext>\; \hspace{0.17em}\; </mtext>\; <mtext>\; \hspace{0.17em}\; </mtext>\; <mtext>\; \hspace{0.17em}\; </mtext>\; <mtext>\; \hspace{0.17em}\; </mtext>\; <mtext>\; \hspace{0.17em}\; </mtext>\; <mtext>\; \hspace{0.17em}\; </mtext>\; 1}\le n\le N-1,\end{array}$ (3.10)

and

$a{B}_h{\xi }^n+A_h{\eta }^n=R_n^{\left(2\right)},\mathrm{<mtext>\; \hspace{0.17em}\; </mtext>\; <mtext>\; \hspace{0.17em}\; </mtext>\; <mtext>\; \hspace{0.17em}\; </mtext>\; <mtext>\; \hspace{0.17em}\; </mtext>\; <mtext>\; \hspace{0.17em}\; </mtext>\; <mtext>\; \hspace{0.17em}\; </mtext>\; <mtext>\; \hspace{0.17em}\; </mtext>\; <mtext>\; \hspace{0.17em}\; </mtext>\; <mtext>\; \hspace{0.17em}\; </mtext>\; <mtext>\; \hspace{0.17em}\; </mtext>\; 0}\le n\le N.$ (3.11)

Theorem 3.2. Assuming that $\left\{u^n\mathrm{,}{v}^n\right\}$ and $\left\{U^n\mathrm{,}{V}^n\right\}$ are, respectively, the solutions of the initial-boundary problem (2.1) and the solutions of the difference scheme (2.15), we have

$\Vert {\delta }_t{\xi }^{n+\frac{1}{2}}\Vert +\Vert {\xi }^n\Vert +\Vert {\xi }^{n+1}\Vert +\Vert {\eta }^{n+1}+{\eta }^n\Vert \le C\left({\tau }^2+h^4\right)\mathrm{.}$

Proof. Similar to the proof of Theorem 3.1, we obtain

$\begin{array}{l}4a{\Vert {\delta }_t{\xi }^{n+\frac{1}{2}}\Vert }^2+{\Vert {\eta }^{n+1}+{\eta }^n\Vert }^2+2ka\left({\Vert {\xi }^k\Vert }^2+{\Vert {\xi }^{k+1}\Vert }^2\right)\\ \le \exp \left(9T\right)[4a{\Vert {\delta }_t{\xi }^{\frac{1}{2}}\Vert }^2+{\Vert {\eta }^1+{\eta }^0\Vert }^2+2ka\left({\Vert {\xi }^0\Vert }^2+{\Vert {\xi }^1\Vert }^2\right)\\ \begin{array}{c}\stackrel{}{}\\ \end{array}+\frac{9\tau }{4}{\displaystyle \underset{n=1}{\overset{k}{\sum }}}\left(4a{\Vert R_n^{\left(1\right)}\Vert }^2+{\Vert G^n\Vert }^2\right)],\end{array}$ (3.12)

where $G^n=\frac{R_{n+1}^{\left(2\right)}-R_{n-1}^{\left(2\right)}}{2\tau }=\frac{\partial }{\partial t}{R}_n^{\left(2\right)}+O\left({\tau }^2\right)$ which is bounded by $\Vert G^n\Vert \le C\left({\tau }^2+h^4\right)$. Then, we estimate the right side of (3.12).

Noticing the initial conditions, we know ${\xi }^0=0$, ${\xi }^1=r^{\left(1\right)}\le C\left({\tau }^3\right)$, $\left|{\delta }_t{\xi }^{\frac{1}{2}}\right|=\left|\frac{{\xi }^1-{\xi }^0}{\tau }\right|\le C{\tau }^2$, then

${\Vert {\delta }_t{\xi }^{\frac{1}{2}}\Vert }^2+{\Vert {\xi }^0\Vert }^2+{\Vert {\xi }^1\Vert }^2\le C{\left({\tau }^2\right)}^2\mathrm{.}$ (3.13)

Next, ${\eta }^0=0,{\eta }^1=r^{\left(2\right)}\le C\left({\tau }^2\right)$, then

${\Vert {\eta }^1+{\eta }^0\Vert }^2\le C{\left({\tau }^2\right)}^2\mathrm{.}$ (3.14)

Finally, noticing $\left|R_n^{\left(1\right)}\right|\le C\left({\tau }^2+h^4\right)\mathrm{,}\Vert G^n\Vert \le C\left({\tau }^2+h^4\right)$, it holds that

$9a\tau {\displaystyle \underset{n=1}{\overset{k}{\sum }}}{\Vert R_n^{\left(1\right)}\Vert }^2+\frac{9\tau }{4}{\displaystyle \underset{n=1}{\overset{k}{\sum }}}{\Vert G^n\Vert }^2\le C{\left({\tau }^2+h^4\right)}^2.$ (3.15)

Substituting (3.13), (3.14) and (3.15) into (3.12), we obtain

$4a{\Vert {\delta }_t{\xi }^{n+\frac{1}{2}}\Vert }^2+{\Vert {\eta }^{n+1}+{\eta }^n\Vert }^2+2ka\left({\Vert {\xi }^n\Vert }^2+{\Vert {\xi }^{n+1}\Vert }^2\right)\le C{\left({\tau }^2+h^4\right)}^2\mathrm{.}$

The theorem has been proved.

4. Numerical Experiment

In this section, we present two numerical examples. In addition to showing convergence orders and the effectiveness of the presented difference scheme (2.15), we apply the difference scheme to a practical problem and the results show the effect of the viscosity coefficient on the plate vibration.

Example 1. In this example, let $\mu =2$, $a=1$, $k=2$, $\Omega =\left(0,1\right)\times \left(0,1\right)$, $0<t\le 1$, the analytic solution is chosen to be $u\left(x,y,t\right)=\cos \pi t\sin \pi x\sin \pi y$, then $v\left(x,y,t\right)=2{\pi }^2\cos \pi t\sin \pi x\sin \pi y$, and the source term is $f\left(x,y,t\right)$,

$f\left(x,y,t\right)=\left(-2\pi \sin \pi t-{\pi }^2\cos \pi t+4{\pi }^4\cos \pi t+2\cos \pi t\right)\sin \pi x\sin \pi y.$

Taking the space step $h=\frac{1}{M_x}=\frac{1}{M_y}$, and the time step $\tau =\frac{1}{N}$, we solve the

problem by using the presented compact implicit difference scheme (2.5). Table 1 shows the error and spatial convergence order of the numerical solution U. Table 2 shows the error and time convergence order of the numerical solution U. Figure 1 and Figure 2 show the images of the numerical solution U and the analytic solution u when the mesh is divided into $M_x=M_y=32$, $N=M_x^2$. It can be seen from the tables that in the sense of the maximum norm and the L2 norm, the time convergence order has reached the second order, and the space convergence order has reached the fourth order, which shows the validity of the compact difference scheme established.

Given below is the images of the numerical solution and the analytic solution at $T=1$ with $\tau =h^2=\frac{1}{{32}^2}$.

Example 2. In this example, let $\Omega =\left(0,1\right)\times \left(0,1\right)$, $0<t\le 3$, $a=2.67$,

Table 1. The computational errors and convergence orders of U in space.

Table 2. The computational errors and convergence orders of U in time.

Figure 1. The numerical solution U when $h=\frac{1}{32}$.

${\phi }_1\left(x,y\right)=\sin \pi x\sin \pi y$, $f\left(x,y,t\right)={\phi }_2\left(x,y\right)=0$. We consider the problem which is more realistic in sense.

We observe vibration at the center point of the plate when the viscosity coefficient $\mu$ takes different values. Setting $k=1$, Figure 3 and Figure 4 show the images of displacement u when $\mu$ is set as 1 and 10 respectively. It can be seen from the figures that using the compact difference method to solve the problem (1.1) well shows the process of damped vibration. And the figures depict the effect of damping caused by the viscous term on vibration. We can see the values of

Figure 2. The analytic solution u when $h=\frac{1}{32}$.

Figure 3. The displacement u when $\mu =1$.

Figure 4. The displacement u when $\mu =10$.

the amplitude of the center point of the plate in different time periods. The amplitude becomes smaller and smaller as time goes on, and speed of decreasing amplitude becomes faster and faster with increase of the viscosity coefficient.

5. Conclusion

In this article, we develop a compact difference scheme for the viscoelastic plate vibration equation. The fourth-order differential equation is transformed into a second-order differential equation system by introducing intermediate variables, then the spatial derivative and time derivative are discretized by fourth-order compact difference method and the central difference method respectively. The stability and convergence of the scheme are proved by using the energy method. The difference scheme is stable and convergence rate is $O\left({\tau }^2+h^4\right)$. The numerical results verify the validity and accuracy of the scheme. In future work, we will consider the viscoelastic plate vibration equation with fractional derivative.

Funding

The work is supported by Shandong Provincial Natural Science Foundation, China (ZR202102250267).

Acknowledgements

The authors would like to thank editor and referees for their valuable advice for the improvement of this article.

Conflicts of Interest

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

References

[1]	Cheng, H. and Takiguchi, T. (2019) Proposal to Study Uniqueness Problem of Rayleigh Wave in Half-Space Kelvin Viscoelastic Media. Journal of Applied Mathematics and Physics, 7, 2074-2088. https://doi.org/10.4236/jamp.2019.79142
[2]	Yao, N. and Luo, Z. (2019) Existence of Solutions for Boundary Value Problems of Vibration Equation with Fractional Derivative. Journal of Applied Mathematics and Physics, 7, 1067-1076. https://doi.org/10.4236/jamp.2019.75072
[3]	Cao, Z.Y. (1989) Vibration Theory of Plates and Shells. China Railway Publishing House, Beijing.
[4]	Wang, Y.D. (2014) Finite Time Blow-Up and Gloal Solutions for Fourth Order Damped Wave Equations. Journal of Mathematical Analysis and Applications, 418, 713-733. https://doi.org/10.1016/j.jmaa.2014.04.015
[5]	Ciment, M. and Leventhal, S.H. (1975) Higher Order Compact Implicit Schemes for the Wave Equation. Mathematics of Computational, 29, 985-994. https://doi.org/10.1090/S0025-5718-1975-0416049-2
[6]	Lu, J.Y. and Ge, Y.B. (2020) A Fourth-Order Compact Difference Method for Solving the 1D Wave Equation. Journal of Ningxia University (Natural Science Edition), 41, 17-22
[7]	Ma, Y.Z., Li, X.G. and Ge, Y.B. (2010) A High-Accuracy Alternating Direction Implicit Method for Solving the Two-Dimensional Wave Equation. Journal of Sichuan Normal University (Natural Science), 33, 179-183.
[8]	Jiang, Y.Z. and Ge, Y.B. (2020) An Explicit Fourth-Order Compact Difference Scheme for Solving the 2D Wave Equation. Advances in Difference Equations, 2020, Article No. 415. https://doi.org/10.1186/s13662-020-02870-z
[9]	Deng, D.W. and Zhang, C.J. (2013) Application of a Fourth-Order Compact ADI Method to Solve a Two Dimensional Linear Hyperbolic Equation. International Journal of Computer Mathematics, 90, 273-291. https://doi.org/10.1080/00207160.2012.713475
[10]	Cheng, T. and Liu, L.B. (2010) Compact Difference Scheme with High Accuracy for Solving Four Order Parabolic Equation. College Mathematics, 26, 71-75.
[11]	Li, Q. and Yang, Q. (2019) Compact Difference Scheme for Two-Dimensional Fourth-Order Hyperbolic Equation. Advances in Difference Equations, 2019, Article No. 328. https://doi.org/10.1186/s13662-019-2094-4
[12]	Li, Q., Yang, Q. and Chen, H.Z. (2020) Compact Difference Scheme for Two-Dimensional Fourth-Order Nonlinear Hyperbolic Equation. Journal of Experimental Agriculture International, 36, 1938-1961. https://doi.org/10.1002/num.22511
[13]	Cui, M.R. (2010) High Order Compact Alternating Direction Implicit Method for the Generalized Sine-Gordon Equation. Journal of Computational and Applied Mathematics, 235, 837-849. https://doi.org/10.1016/j.cam.2010.07.016
[14]	Wang, Y.M. and Guo, B.Y. (2007) Fourth-Order Compact Finite Difference Method for Fourth-Order Nonlinear Elliptic Boundary Value Problems. Journal of Computational and Applied Mathematics, 221, 76-97. https://doi.org/10.1016/j.cam.2007.10.007
[15]	Shi, F.Y. and Li, S.C. (2018) A Three-Points Stable and Compact Difference Scheme for the Vibration Equation of Beams. Journal of Hebei Normal University (Natural Science Edition), 42, 19-23.
[16]	Sun, Z.Z. (2012) Numerical Solution of Partial Differential Equations. Science Press, Beijing.
[17]	Deng, D.W. and Zhang, C.J. (2013) A New Fourth-Order Numerical Algorithm for a Class of Three-Dimensional Nonlinear Evolution Equations. Numerical Methods for Partial Differential Equations, 29, 102-130. https://doi.org/10.1002/num.21701