Energy estimates for two-dimensional space-Riesz fractional wave equation

TL;DR

This paper develops an energy method to analyze the stability and convergence of the two-dimensional space-Riesz fractional wave equation, providing theoretical proofs and numerical verification of its accuracy and stability.

Contribution

It introduces a novel energy-based approach for error estimation of the 2D space-Riesz fractional wave equation with variable coefficients, including stability and convergence analysis.

Findings

Unconditional stability of the scheme is proved.

Convergence order is d7b5b2 in space and time.

Numerical results confirm theoretical error estimates.

Abstract

The fractional wave equation governs the propagation of mechanical diffusive waves in viscoelastic media which exhibits a power-law creep, and consequently provided a physical interpretation of this equation in the framework of dynamic viscoelasticity. In this paper, we first develop the energy method to estimate the one-dimensional space-Riesz fractional wave equation. For two-dimensional cases with the variable coefficients, the discretized matrices are proved to be commutative, which ensures to carry out of the priori error estimates. The unconditional stability and convergence with the global truncation error $\mathcal{O}(\tau^2+h^2)$ are theoretically proved and numerically verified. In particulary, the framework of the priori error estimates and convergence analysis are still valid for the compact finite difference scheme and the nonlocal wave equation.

Tables

Table 1. Table 1 : The maximum errors and convergent orders for ( 2.14 ) with τ = h 𝜏 ℎ \tau=h .

$\tau$	$\alpha =1.3,\theta =0.25$	Rate	$\alpha =1.6,\theta =0.5$	Rate	$\alpha =1.9,\theta =1$	Rate
1/40	8.8516e-05		9.0532e-05		7.5678e-05
1/80	2.2156e-05	1.9983	2.2329e-05	2.0195	1.9487e-05	1.9574
1/160	5.5242e-06	2.0038	5.5161e-06	2.0172	4.7477e-06	2.0372
1/320	1.3761e-06	2.0052	1.3636e-06	2.0163	1.1581e-06	2.0355

Table 2. Table 2 : The maximum errors and convergent orders for ( 2.27 ) with τ = h x = h y 𝜏 subscript ℎ 𝑥 subscript ℎ 𝑦 \tau=h_{x}=h_{y} and θ = 0.75 𝜃 0.75 \theta=0.75 .

$\tau$	$\alpha =1.3,\beta =1.7$	Rate	$\alpha =1.5,\beta =1.5$	Rate	$\alpha =1.7,\beta =1.3$	Rate
1/20	1.4066e-04		1.4066e-04		1.4500e-04
1/40	3.8290e-05	1.8772	3.7449e-05	1.9093	3.7041e-05	1.9688
1/80	9.4992e-06	2.0111	9.4992e-06	1.9790	9.4992e-06	1.9632
1/160	2.4049e-06	1.9818	2.4049e-06	1.9818	2.4049e-06	1.9818

Equations

$$\frac{{\partial }^{2}u(x,y,t)}{\partial t^2}=a(x,y)\frac{{\partial }^{\alpha }u(x,y,t)}{\partial |x{|}^{\alpha }}+b(x,y)\frac{{\partial }^{\beta }u(x,y,t)}{\partial |y{|}^{\beta }}+f(x,y,t).$$

$$\begin{array}[]{ll}&u(x,y,0)=\varphi(x,y)~{}~{}~{}{\rm for}~{}~{}~{}(x,y)\in\Omega,\\ &u_{t}(x,y,0)=\psi(x,y)~{}~{}~{}{\rm for}~{}~{}~{}(x,y)\in\Omega,\end{array}$$

$$u(x,y,t)=0\mathrm{for}(x,y)\in \partial \Omega$$

$$\frac{{\partial }^{\alpha }u(x,y,t)}{\partial |x{|}^{\alpha }}=-{\kappa }_{\alpha }\left({}_0{D}_x^{\alpha }{+}_x{D}_{x_r}^{\alpha }\right)u(x,y,t)\mathrm{with}{\kappa }_{\alpha }=\frac{1}{2\cos (\alpha \pi /2)},$$

$$\begin{array}{rl}{}_0{D}_x^{\alpha }u(x,y,t)& =\frac{1}{\Gamma (n-\alpha )}\frac{{\partial }^n}{\partial x^n}{\int }_0^x{\left(x-\xi \right)}^{n-\alpha -1}u(\xi ,y,t)d\xi ,\\ {}_x{D}_{x_r}^{\alpha }u(x,y,t)& =\frac{(-1{)}^n}{\Gamma (n-\alpha )}\frac{{\partial }^n}{\partial x^n}{\int }_x^{x_r}{\left(\xi -x\right)}^{n-\alpha -1}u(\xi ,y,t)d\xi .\end{array}$$

$$\begin{array}{rl}{}_0{D}_x^{\alpha }u(x_i)& ={\delta }_{x,+}^{\alpha }u(x_i)+\mathcal{O}(h^2)\mathrm{with}{\delta }_{x,+}^{\alpha }u(x_i)=\frac{1}{h^{\alpha }}\sum _{m=0}^{i+1}{\phi }_m^{\alpha }u(x_{i-m+1}),\\ {}_x{D}_{x_r}^{\alpha }u(x_i)& ={\delta }_{x,-}^{\alpha }u(x_i)+\mathcal{O}(h^2)\mathrm{with}{\delta }_{x,-}^{\alpha }u(x_i)=\frac{1}{h^{\alpha }}\sum _{m=0}^{N_x-i+1}{\phi }_m^{\alpha }u(x_{i+m-1}),\end{array}$$

$$\begin{array}{r}{\phi }_0^{\alpha }=\frac{\alpha }{2}{g}_0^{\alpha },{\phi }_m^{\alpha }=\frac{\alpha }{2}{g}_m^{\alpha }+\frac{2-\alpha }{2}{g}_{m-1}^{\alpha },m\ge 1,\end{array}$$

$$g_m^{\alpha }=(-1{)}^m\left(\begin{array}{c}\alpha \\ m\end{array}\right),i.e.,g_0^{\alpha }=1,g_m^{\alpha }=\left(1-\frac{\alpha +1}{m}\right)g_{m-1}^{\alpha },m\ge 1.$$

$$\begin{array}{r}\frac{{\partial }^{\alpha }u(x_i)}{\partial |x{|}^{\alpha }}={\nabla }_h^{\alpha }u(x_i)+\mathcal{O}(h^2)\end{array}$$

$${\nabla }_h^{\alpha }u(x_i)=-{\kappa }_{\alpha }\left({\delta }_{x,+}^{\alpha }+{\delta }_{x,-}^{\alpha }\right)u(x_i)=-\frac{{\kappa }_{\alpha }}{h^{\alpha }}\sum _{l=0}^{N_x}{\phi }_{i,l}^{\alpha }u(x_l),$$

$$\varphi_{i,l}^{\alpha}=\left\{\begin{array}[]{l@{\quad} l}\varphi_{i-l+1}^{\alpha},&l<i-1,\\ \varphi_{0}^{\alpha}+\varphi_{2}^{\alpha},&l=i-1,\\ 2\varphi_{1}^{\alpha},&l=i,\\ \varphi_{0}^{\alpha}+\varphi_{2}^{\alpha},&l=i+1,\\ \varphi_{l-i+1}^{\alpha},&l>i+1.\end{array}\right.$$

$$\begin{array}{rl}& \frac{1}{h^{\alpha }}{\left[\sum _{l=0}^{N_x}{\phi }_{1,l}^{\alpha }u(x_l),\sum _{l=0}^{N_x}{\phi }_{2,l}^{\alpha }u(x_l),\dots ,\sum _{l=0}^{N_x}{\phi }_{N_x-1,l}^{\alpha }u(x_l)\right]}^T=\left({\delta }_{x,+}^{\alpha }+{\delta }_{x,-}^{\alpha }\right)u=\frac{1}{h^{\alpha }}{A}_{\alpha }u,\end{array}$$

$$\begin{array}{rl}& {\nabla }_h^{\alpha }u=-{\kappa }_{\alpha }\left({\delta }_{x,+}^{\alpha }+{\delta }_{x,-}^{\alpha }\right)u=\frac{-{\kappa }_{\alpha }}{h^{\alpha }}{A}_{\alpha }u,\end{array}$$

$$A_{\alpha }=B_{\alpha }+B_{\alpha }^T\mathrm{with}{B}_{\alpha }=\left[\begin{array}{cccccc}{\phi }_1^{\alpha }& {\phi }_2^{\alpha }& {\phi }_3^{\alpha }& \cdots & {\phi }_{N_x-2}^{\alpha }& {\phi }_{N_x-1}^{\alpha }\\ {\phi }_0^{\alpha }& {\phi }_1^{\alpha }& {\phi }_2^{\alpha }& {\phi }_3^{\alpha }& \cdots & {\phi }_{N_x-2}^{\alpha }\\ & {\phi }_0^{\alpha }& {\phi }_1^{\alpha }& {\phi }_2^{\alpha }& \ddots & ⋮\\ & & \ddots & \ddots & \ddots & {\phi }_3^{\alpha }\\ & & & \ddots & {\phi }_1^{\alpha }& {\phi }_2^{\alpha }\\ & & & & {\phi }_0^{\alpha }& {\phi }_1^{\alpha }\end{array}\right].$$

$$\frac{{\partial }^{2}u(x,t)}{\partial t^2}=a(x)\frac{{\partial }^{\alpha }u(x,t)}{\partial |x{|}^{\alpha }}+f(x,t)$$

$$\begin{array}[]{ll}&u(x,0)=\varphi(x)~{}~{}~{}{\rm for}~{}~{}~{}x\in\Omega,\\ &u_{t}(x,0)=\psi(x)~{}~{}~{}{\rm for}~{}~{}~{}x\in\Omega.\end{array}$$

$$\frac{{\partial }^{2}u(x,t)}{\partial t^2}={\delta }_t^{2}u(x_i,t_k)+\mathcal{O}({\tau }^2)\mathrm{with}{\delta }_t^{2}u(x_i,t_k)=\frac{u_i^{k+1}-2u_i^k+u_i^{k-1}}{{\tau }^2}.$$

$$\theta u_i^{k+1}+(1-2\theta )u_i^k+\theta u_i^{k-1},\frac{1}{4}\le \theta \le 1,$$

$$\begin{array}{rl}& \frac{u(x_i,t_{k+1})-2u(x_i,t_k)+u(x_i,t_{k-1})}{{\tau }^2}\\ & =a(x_i){\nabla }_h^{\alpha }\left[\theta u(x_i,t_{k+1})+(1-2\theta )u(x_i,t_k)+\theta u(x_i,t_{k-1})\right]+f(x_i,t_k)+R_i^k\end{array}$$

$$R_i^k\le C_{u,\alpha }({\tau }^2+h^2),$$

$${\delta }_t^2{u}_i^k=\theta a_i{\nabla }_h^{\alpha }{u}_i^{k+1}+(1-2\theta )a_i{\nabla }_h^{\alpha }{u}_i^k+\theta a_i{\nabla }_h^{\alpha }{u}_i^{k-1}+f_i^k,$$

$$\begin{array}{rl}& u_i^{k+1}+\theta \frac{{\tau }^2}{h^{\alpha }}{\kappa }_{\alpha }{a}_i\left[\sum _{m=0}^{i+1}{\phi }_m^{\alpha }{u}_{i-m+1}^{k+1}+\sum _{m=0}^{N_x-i+1}{\phi }_m^{\alpha }{u}_{i+m-1}^{k+1}\right]\\ & =2u_i^k-(1-2\theta )\frac{{\tau }^2}{h^{\alpha }}{\kappa }_{\alpha }{a}_i\left[\sum _{m=0}^{i+1}{\phi }_m^{\alpha }{u}_{i-m+1}^k+\sum _{m=0}^{N_x-i+1}{\phi }_m^{\alpha }{u}_{i+m-1}^k\right]\\ & -u_i^{k-1}-\theta \frac{{\tau }^2}{h^{\alpha }}{\kappa }_{\alpha }{a}_i\left[\sum _{m=0}^{i+1}{\phi }_m^{\alpha }{u}_{i-m+1}^{k-1}+\sum _{m=0}^{N_x-i+1}{\phi }_m^{\alpha }{u}_{i+m-1}^{k-1}\right]+{\tau }^2{f}_i^k.\end{array}$$

$$\begin{array}{rl}u(x_i,\tau )& =u(x_i,0)+\tau \frac{\partial u(x_i,0)}{\partial t}+\frac{{\tau }^2}{2}\frac{{\partial }^{2}u(x_i,0)}{\partial t^2}+\frac{1}{2}{\int }_0^{\tau }(\tau -t{)}^2\frac{{\partial }^{3}u(x_i,t)}{\partial t^3}dt\\ & =\phi (x_i)+\tau \psi (x_i)+\frac{{\tau }^2}{2}\left[a(x_i)\frac{{\partial }^{\alpha }u(x_i,0)}{\partial |x{|}^{\alpha }}+f(x_i,0)\right]+\frac{1}{2}{\int }_0^{\tau }(\tau -t{)}^2\frac{{\partial }^{3}u(x_i,t)}{\partial t^3}dt.\end{array}$$

$$\begin{array}{rl}{u}_i^1& =\phi (x_i)+\tau \psi (x_i)+\frac{{\tau }^2}{2}\left[a(x_i){\nabla }_h^{\alpha }u(x_i,0)+f(x_i,0)\right]\end{array}$$

$$U^k={\left[u_1^k,u_2^k,\dots ,u_{N_x-1}^k\right]}^{\mathrm{T}},F^k={\left[f_1^k,f_2^k,\dots ,f_{N_x-1}^k\right]}^{\mathrm{T}}.$$

$$\begin{array}{rl}& \left[I+\theta \frac{{\tau }^2}{h^{\alpha }}{\kappa }_{\alpha }D{A}_{\alpha }\right]U^{k+1}=\left[2I-\left(1-2\theta \right)\frac{{\tau }^2}{h^{\alpha }}{\kappa }_{\alpha }D{A}_{\alpha }\right]U^k-\left[I+\theta \frac{{\tau }^2}{h^{\alpha }}{\kappa }_{\alpha }D{A}_{\alpha }\right]U^{k-1}+\tau F^k,\end{array}$$

$$D=\left[\begin{array}{cccc}{a}_1& & & \\ & a_2& & \\ & & \ddots & \\ & & & a_{N_x-1}\end{array}\right].$$

$$\frac{{\partial }^{2}u(x,y,t)}{\partial t^2}=\frac{u_{i,j}^{k+1}-2u_{i,j}^k+u_{i,j}^{k-1}}{{\tau }^2}+\mathcal{O}({\tau }^2),$$

$$\begin{split}&\frac{u(x_{i},y_{j},t_{k+1})-2u(x_{i},y_{j},t_{k})+u(x_{i},y_{j},t_{k-1})}{\tau^{2}}\\ &=a(x_{i},y_{j})\nabla^{\alpha}_{h_{x}}\Big{(}\theta u(x_{i},y_{j},t_{k+1})+(1-2\theta)u(x_{i},y_{j},t_{k})+\theta u(x_{i},y_{j},t_{k-1})\Big{)}\\ &\quad+b(x_{i},y_{j})\nabla^{\beta}_{h_{y}}\Big{(}\theta u(x_{i},y_{j},t_{k+1})+(1-2\theta)u(x_{i},y_{j},t_{k})+\theta u(x_{i},y_{j},t_{k-1})\Big{)}+f(x_{i},y_{j},t_{k})+R_{i,j}^{k},\end{split}$$

$$R_{i,j}^k\le C_{u,\alpha ,\beta }({\tau }^2+h_x^2+h_y^2).$$

Full text

Keywords: Riesz fractional wave equation; Nonlocal wave equation; Priori error estimates; Energy method; Numerical stability and convergence

1 Introduction

The fractional wave equation is obtained from the classical wave equation by replacing the second-order derivative with a fractional derivative of order $\alpha$, $1<\alpha\leq 2$. Mainardi [19] pointed out that the fractional wave equation governs the propagation of mechanical diffusive waves in viscoelastic media which exhibits a power-law creep, and consequently provided a physical interpretation of this equation in the framework of dynamic viscoelasticity. In this paper, we study a second-order accurate numerical method in both space and time for the two-dimensional space-Riesz fractional wave equation with the variable coefficients whose prototype is, for $1<\alpha,\beta\leq 2$,

[TABLE]

The initial conditions are

[TABLE]

and the Dirichlet boundary condition

[TABLE]

with $\Omega=(0,x_{r})\times(0,y_{r})$. The function $f(x,y,t)$ is a source term and all the coefficients are positive, i.e., $0<a_{0}\leq a(x,y)\leq a_{1}$ and $0<b_{0}\leq b(x,y)\leq b_{1}$.

The space-Riesz fractional derivative appears in the continuous limit of lattice models with long-range interactions [30], for $n\in\mathbb{N}$, $n-1\leq\alpha<n$, which is defined as [25]

[TABLE]

where

[TABLE]

For the Caputo-Riesz time-space fractional wave equation ${{}^{c}}\!{D}_{t}^{\gamma}u(x,t)=\frac{\partial^{\alpha}u(x,t)}{\partial|x|^{\alpha}}$ with $1<\alpha,\gamma\leq 2$, Mainardi (2001) et al. obtained the fundamental solution of the space-time fractional diffusion equation [20]. Metzler and Nonnenmacher (2002) investigated the physical backgrounds and implications of a space-and time-fractional diffusion and wave equation [22]. The numerical solution of space-time fractional diffusion-wave equations are discussed in [2, 12], but they are lack of the stability and convergence analysis. To rewrite the fractional diffusion-wave equation as the the Volterra type integro-differential equations, the stability and convergence analysis are given with the zero initial conditions [6]. For $1<\gamma<2$ and $\alpha=2$, it has been proposed by various authors [7, 8, 17, 21, 23, 31, 34, 35, 36]. For example, based on the second-order fractional Lubich’s methods [18], Cuesta (2006) et al. derived the second-order error bounds of the time discretization in a Banach space with the $\frac{\partial^{2}u}{\partial x^{2}}$ a sectorial operator [8] and Yang (2014) et al. obtained the second-order convergence schemes with $1\leq\gamma\leq 1.71832$ [34]. For $\gamma=2$ and $1<\alpha<2$, it seems that achieving a second-order accurate scheme for (1.1) is not an easy task with the nonzero initial conditions. This paper focuses on providing the weighted numerical scheme to solve the space-Riesz fractional wave equation with the nonzero initial conditions and the variable coefficients in one-dimensional and two-dimensional cases. The unconditional stability and convergence with the global truncation error $\mathcal{O}(\tau^{2}+h^{2})$ are theoretically proved and numerically verified by the energy method, which can be easily extended to the nonlocal wave equation [11].

The rest of the paper is organized as follows. The next section proposes the second-order accurate scheme for (1.1). In Section 3, we carry out a detailed stability and convergence analysis with the second order accuracy in both time and space directions for the derived schemes. To show the effectiveness of the schemes, we perform the numerical experiments to verify the theoretical results in Section 4. The paper is concluded with some remarks in the last section.

2 Discretization Schemes

Let the mesh points $x_{i}=ih$, $i=0,1,\ldots,N_{x}$, and $t_{k}=k\tau$, $0\leq k\leq{N_{t}}$ with $h=x_{r}/{N_{x}}$, $\tau=T/{N_{t}}$, i.e., $h$ is the uniform space stepsize and $\tau$ the time stepsize. And ${u_{i}^{k}}$ denotes the approximated value of $u(x_{i},t_{k})$, $a_{i}=a(x_{i})$, $f_{i}^{k}=f(x_{i},t_{k})$.

Nowadays, there are already many types of high order discretization schemes for the Riemann-Liouville space fractional derivatives [3, 13, 15, 24, 27, 29, 32]. Here, we take the following schemes to approach (1.3), see in [4, 32]

[TABLE]

where

[TABLE]

and

[TABLE]

Using (1.3) and (2.1), we obtain the approximation operator of the space-Riesz fractional derivative

[TABLE]

with

[TABLE]

where $i=1,\ldots,N_{x}-1$ (together with the zero Dirichlet boundary conditions) and

[TABLE]

Taking $u=[u({x_{1}}),u({x_{2}}),\cdots,u({x_{N_{x}-1}})]^{\rm T}$, and using (2.1), (2.2), there exists

[TABLE]

it yields

[TABLE]

where the matrix

[TABLE]

2.1 Numerical scheme for one-dimensional space-Riesz fractional wave equation

We now examine the full discretization scheme to the one-dimensional space-Riesz fractional wave equation, i.e,

[TABLE]

with $0<a_{0}\leq a(x)\leq a_{1}$ and the zero Dirichlet boundary condition. The initial conditions are

[TABLE]

In the time direction derivative, we use the following center difference scheme

[TABLE]

In order to achieve an unconditional stable algorithm, we use the weighted algorithm for the space-Riesz fractional derivative, i.e.,

[TABLE]

to approximate $u(x_{i},t_{k})$. From (2.2) and the above equations, we can rewrite (2.5) as

[TABLE]

with the local truncation error

[TABLE]

where the constant $C_{u,\alpha}$ is independent of $h$ and $\tau$. Therefore, the full discretization of (2.5) has the following form

[TABLE]

i.e.,

[TABLE]

Using (2.5), (2.6) and Taylor expansion with integral form of the remainder, there exists

[TABLE]

Then we can obtain $u_{i}^{1}$, i.e.,

[TABLE]

with the local truncation error $\mathcal{O}(\tau^{3}+\tau^{2}h^{2})$, see Section 3.

For the convenience of implementation, we use the matrix form of the grid functions

[TABLE]

Hence, the finite difference scheme (2.11) can be recast as

[TABLE]

where $A_{\alpha}$ is defined by (2.4) and the diagonal matrix

[TABLE]

2.2 Numerical scheme for two-dimensional space-Riesz fractional wave equation

Let the mesh points $x_{i}=ih_{x}$, $i=0,1,\ldots,N_{x}$ and $y_{j}=jh_{y}$, $j=0,1,\ldots,N_{y}$ and $t_{k}=k\tau$, $0\leq k\leq{N_{t}}$ with $h_{x}=x_{r}/{N_{x}}$, $h_{y}=y_{r}/N_{y}$, $\tau=T/N_{t}$. Similarly, we take $u_{i,j}^{k}$ as the approximated value of $u(x_{i},y_{j},t_{k})$, $a_{i,j}=a(x_{i},y_{j})$, $b_{i,j}=b(x_{i},y_{j})$, $f_{i,j}^{k}=f(x_{i},y_{j},t_{k})$. We use the center difference scheme to do the discretization in time direction derivative,

[TABLE]

and the weighted schemes for the space-Riesz fractional derivative, i.e., $\theta u_{i,j}^{k+1}+(1-2\theta)u_{i,j}^{k}+\theta u_{i,j}^{k-1}$ to approximate $u(x_{i},y_{j},t_{k})$. Therefore (1.1) can be rewritten as

[TABLE]

where the local truncation error is

[TABLE]

Similarly, we denote

[TABLE]

Therefore, the resulting discretization of (1.1) has the following form

[TABLE]

i.e.,

[TABLE]

Using (2.12) and (2.13), we can obtain

[TABLE]

with the local truncation error $\mathcal{O}(\tau^{3}+\tau^{2}h_{x}^{2}+\tau^{2}h_{y}^{2})$, see Section 3.

For the two-dimensional space-Riesz fractional wave equation (1.1), the relevant perturbation equation of (2.20) is of the form

[TABLE]

Comparing (2.22) with (2.20), the splitting term is given by

[TABLE]

since \big{(}u_{i,j}^{k+1}-2u_{i,j}^{k}+u_{i,j}^{k-1}\big{)} is an $\mathcal{O}(\tau^{2})$ term, it implies that the perturbation contributes an $\mathcal{O}(\tau^{6})$ error component to the truncation error of (2.20). Thus we can rewrite (1.1) as

[TABLE]

where

[TABLE]

Hence, the system (2.22) can be solved by the alternating direction implicit method (D-ADI) [9, 10]:

[TABLE]

where $u_{i,j}^{*}$ is an intermediate solution. Take

[TABLE]

and denote

[TABLE]

where $I$ denotes the unit matrix and the symbol $\otimes$ the Kronecker product [16], and $A_{\alpha}$, $A_{\beta}$ are defined by (2.4). Therefore, we can rewrite (2.25) as the following form

[TABLE]

where

[TABLE]

and

[TABLE]

3 Convergence and Stability Analysis

To rewrite the fractional diffusion-wave equation as the the Volterra type integro-differential equations, the stability and convergence analysis are given with the zero initial conditions [6]. Here, we first develop the energy method to estimate the space-Riesz fractional wave equation with the nonzero initial conditions. For two-dimensional cases with the variable coefficients, the discretized matrices are proved to be commutative, which ensures to carry out of the priori error estimates.

Lemma 3.1**.**

[33]** Let $\nabla^{\alpha}_{h}$ be given in (2.3) and $1<\alpha\leq 2$. Then there exists an symmetric positive definite matrix $\Lambda_{h}^{\alpha}$ such that

[TABLE]

Lemma 3.2**.**

(Discrete Gronwall Lemma [26]) Assume that $\{a_{k}\}$ and $\{b_{k}\}$ is a nonnegative sequence, and the sequence $\varphi^{k}$ satisfies

[TABLE]

where $c_{0}\geq 0.$ Then the sequence $\{\varphi^{k}\}$ satisfies

[TABLE]

Lemma 3.3**.**

[16, p. 141]** Let $A\in\mathbb{R}^{n\times n}$ have eigenvalues $\{\lambda_{i}\}_{i=1}^{n}$ and $B\in\mathbb{R}^{m\times m}$ have eigenvalues $\{\mu_{j}\}_{j=1}^{m}$. Then the $mn$ eigenvalues of $A\otimes B$ are

[TABLE]

Lemma 3.4**.**

[16, p. 140]** Let $A\in\mathbb{R}^{m\times n}$, $B\in\mathbb{R}^{r\times s}$, $C\in\mathbb{R}^{n\times p}$, and $D\in\mathbb{R}^{s\times t}$. Then

[TABLE]

Moreover, for all $A$ and $B$, $(A\otimes B)^{T}=A^{T}\otimes B^{T}$.

Lemma 3.5**.**

Let $\mathcal{A}_{x}=I\otimes A_{\alpha}$ and $\mathcal{A}_{y}=A_{\beta}\otimes I$ be defined by (2.26). Then

[TABLE]

where we denote $\Lambda_{x}:=I\otimes\Lambda_{\alpha}$ and $\Lambda_{y}:=\Lambda_{\beta}\otimes I$.

Proof.

From [33] or Lemma 3.1, there exists $-A_{\alpha}=\Lambda_{\alpha}\cdot\Lambda_{\alpha}$ and $-A_{\beta}=\Lambda_{\beta}\cdot\Lambda_{\beta}$, since $-A_{\alpha}$ and $-A_{\beta}$ are the symmetric positive definite matrices. Taking $\Lambda_{x}:=I\otimes\Lambda_{\alpha}$ and $\Lambda_{y}:=\Lambda_{\beta}\otimes I$ and using Lemma 3.4, the results are obtained. ∎

Lemma 3.6**.**

Let $\nabla^{\alpha}_{h_{x}}$ and $\nabla^{\beta}_{h_{y}}$ be given in (2.18) with $1<\alpha,\beta\leq 2$. Then there exist the symmetric positive definite matrices $\Lambda_{h_{x}}^{\alpha}$ and $\Lambda_{h_{y}}^{\beta}$, respectively, such that

[TABLE]

and

[TABLE]

Proof.

According to (2.18) and (2.26), it implied that

[TABLE]

From Lemmas 3.3 and 3.5, we know that $\mathcal{A}_{x}=I\otimes A_{\alpha}$ is a symmetric negative definite, which leads to $-\nabla^{\alpha}_{h_{x}}$ (or $-\nabla^{\beta}_{h_{y}}$ ) is the symmetric positive definite. The proof is completed. ∎

3.1 Convergence and stability for one-dimensional space-Riesz fractional wave equation

First, we introduce some relevant notations and properties of discretized inner product given in [14, 28]. Denote $u^{k}=\{u_{i}^{k}|0\leq i\leq N_{x},0\leq k\leq N_{t}\}$ and $v^{k}=\{v_{i}^{k}|0\leq i\leq N_{x},0\leq k\leq N_{t}\}$, which are grid functions. And

[TABLE]

Lemma 3.7**.**

Let $\frac{1}{4}\leq\theta\leq 1$, $1<\alpha\leq 2$ and $\{u_{i}^{k}\}$ be the solution of the difference scheme

[TABLE]

with the initial conditions and the Dirichlet boundary conditions

[TABLE]

Then

[TABLE]

where the energy norm is defined by

[TABLE]

Proof.

Multiplying (2.10) by $h(u^{k+1}_{i}-u^{k-1}_{i})$, respectively, it yields

[TABLE]

and

[TABLE]

Then summing up for $i$ from 1 to $N_{x}-1$ for the above equations, respectively, there exists

[TABLE]

and

[TABLE]

where

[TABLE]

According to Lemma 3.1, which leads to

[TABLE]

and

[TABLE]

Combine (2.10), (3.2) and (3.3), we obtain

[TABLE]

i.e.,

[TABLE]

Adding $\theta\parallel\sqrt{a}\Lambda_{h}^{\alpha}u^{k}\parallel^{2}$ on both sides of the above equation, there exists

[TABLE]

Denoting

[TABLE]

i.e.,

[TABLE]

where we use

[TABLE]

From

[TABLE]

and (3.5), (3.4), we obtain

[TABLE]

i.e,

[TABLE]

Therefore, for $\tau\leq 2/3$, it yields

[TABLE]

Using the discrete Gronwall inequality (see Lemma 3.2), we have

[TABLE]

The proof is completed. ∎

Theorem 3.1**.**

Let $u(x_{i},t_{k})$ be the exact solution of (2.5) with $1<\alpha\leq 2$, $\frac{1}{4}\leq\theta\leq 1$; $u^{k}_{i}$ be the solution of the finite difference scheme (2.10) and $e_{i}^{k}=u(x_{i},t_{k})-u^{k}_{i}$. Then

[TABLE]

where the energy norm is defined by

[TABLE]

Proof.

Subtracting (2.10) from (2.8), it yields

[TABLE]

Using Lemma 3.7, we obtain

[TABLE]

where

[TABLE]

Next we estimate the local error truncation of $E^{0}$. Since $e_{i}^{0}=0$ and

[TABLE]

where $\xi_{i}\in(0,x_{r})$ and

[TABLE]

it implies that

[TABLE]

Here, the coefficients $C_{l,\alpha},1\leq l\leq 2$ are the constants independent of $h$ and $\tau$.

According to (2.2) and the above equations, there exists

[TABLE]

where $\xi_{i},\overline{\xi}_{i},\widetilde{\xi_{i}}\in(0,x_{r})$ and $C_{l,\alpha},1\leq l\leq 5$ are the constants independent of $h$ and $\tau$. Using (3.9), (3.10) and the above equation, we have

[TABLE]

with a constant $C_{\alpha}$. From (2.9), (3.8) and (3.11), it means that

[TABLE]

with $\widetilde{C}_{\alpha}=2\max\{C_{\alpha}^{2}x_{r},\frac{3}{2}C^{2}_{u,\alpha}T\}$. The proof is completed. ∎

Theorem 3.2**.**

The difference scheme (2.14) with $1<\alpha\leq 2$ and $\frac{1}{4}\leq\theta\leq 1$ is unconditionally stable.

Proof.

From Lemma 3.7, the proof is completed. ∎

3.2 Convergence and stability for two-dimensional space-Riesz fractional wave equation

Denote $u^{k}=\{u_{i,j}^{k}|0\leq i\leq N_{x},0\leq j\leq N_{y},0\leq k\leq N_{t}\}$ and $v^{k}=\{v_{i,j}^{k}|0\leq i\leq N_{x},0\leq j\leq N_{y},0\leq k\leq N_{t}\}$, which are grid functions. And

[TABLE]

Lemma 3.8**.**

Let $\frac{1}{4}\leq\theta\leq 1$, $1<\alpha,\beta\leq 2$ and $\{u_{ij}^{k}\}$ be the solution of the difference scheme

[TABLE]

with the initial conditions and the Dirichlet boundary conditions

[TABLE]

Then

[TABLE]

where the energy norm is defined by

[TABLE]

Proof.

Multiplying (3.13) by $h_{x}h_{y}(u^{k+1}_{i,j}-u^{k-1}_{i,j})$ and using Lemmas 3.5, 3.6, there exists

[TABLE]

and

[TABLE]

Then summing up for $i$ from 1 to $N_{x}-1$ and for $j$ from 1 to $N_{y}-1$, we have

[TABLE]

and

[TABLE]

where

[TABLE]

According to Lemma 3.6, we have

[TABLE]

and

[TABLE]

From (3.14) and (3.15), we obtain

[TABLE]

i.e.,

[TABLE]

Adding $\theta\parallel\sqrt{a}\Lambda_{h_{x}}^{\alpha}u^{k}\parallel^{2}+\theta\parallel\sqrt{b}\Lambda_{h_{y}}^{\beta}u^{k}\parallel^{2}$ on both sides of the above equation, we have

[TABLE]

Denoting

[TABLE]

we have

[TABLE]

We rewrite (3.16) as the following form

[TABLE]

where we use

[TABLE]

and

[TABLE]

According to

[TABLE]

and (3.18), (3.17), there exists

[TABLE]

i.e.,

[TABLE]

For $\tau\leq 2/3$, which leads to

[TABLE]

From Lemma 3.2, there exists

[TABLE]

The proof is completed. ∎

Theorem 3.3**.**

Let $u(x_{i},y_{j},t_{k})$ be the exact solution of (1.1) with $1<\alpha,\beta\leq 2$, $u^{k}_{i,j}$ be the solution of (2.22) and $e_{ij}^{k}=u(x_{i},y_{j},t_{k})-u_{ij}^{k}$. Then

[TABLE]

where the energy norm is defined by

[TABLE]

Proof.

Subtracting (2.22) from (2.23), it yields

[TABLE]

Using Lemma 3.8, there exists

[TABLE]

with the energy norm

[TABLE]

Next we estimate the local error truncation of $E^{0}$. Since $e_{i,j}^{0}=0$,   $\theta^{2}\tau^{6}\parallel\sqrt{ab}\Lambda_{h_{x}}^{\alpha}\Lambda_{h_{y}}^{\beta}e^{1}_{\overline{t}}\parallel^{2}=\mathcal{O}(\tau^{6})$ in (3.14) and

[TABLE]

Here the coefficients $C_{1,\alpha}$ and $C_{1,\beta}$ are the constants independent of $h$, $\tau$ and

[TABLE]

Then we obtain

[TABLE]

From (1.1) and the above equations, there exists

[TABLE]

where $\xi_{i},\widehat{\xi}_{i},\overline{\xi_{i}},\widetilde{\xi_{i}}\in(0,x_{r})$ and $C_{l,\alpha}$ and $C_{l,\alpha,\beta},1\leq l\leq 3$ are the constants. Similarly, we have

[TABLE]

with the constants $\widetilde{C}_{1,\alpha,\beta}$ and $\widetilde{C}_{3,\alpha,\beta}$.

According to (3.20), (3.21) and the above equations, we get

[TABLE]

where $C_{\alpha,\beta}$ is a constant. Hence, using (2.17), (3.19) and (3.22), there exists

[TABLE]

with $\widetilde{C}_{\alpha,\beta}=2\max\{C_{\alpha,\beta}^{2}x_{r}y_{r},\frac{3}{2}C^{2}_{u,\alpha,\beta}T\}$. The proof is completed. ∎

Theorem 3.4**.**

The difference scheme (2.27) with $1<\alpha,\beta\leq 2$ and $\frac{1}{4}\leq\theta\leq 1$ is unconditionally stable.

Proof.

From Lemma 3.7, the result is obtained. ∎

Remark 3.1**.**

The operator $\mathcal{L}$ appears in the nonlocal wave equation [11]

[TABLE]

From [5], we known that the approximation operator of $-\mathcal{L}_{\delta}$ is also the symmetric positive definite. Hence, the framework of the stability and convergence analysis are still valid for the nonlocal wave equation.

4 Numerical results

In this section, we numerically verify the above theoretical results and the $l_{\infty}$ norm is used to measure the numerical errors.

Example 4.1**.**

Consider the space-Riesz fractional wave equation (2.5), on a finite domain $0<x<1$, $0<t\leq 1$ with the coefficient $d(x)=x^{\alpha}$, the forcing function is

[TABLE]

with the initial conditions $u(x,0)=x^{2}(1-x)^{2}$, $\frac{\partial}{\partial t}u(x,0)=-x^{2}(1-x)^{2}$, and the boundary conditions $u(0,t)=u(1,t)=0$. The exact solution of the fractional PDEs is

[TABLE]

Table 1 shows that the scheme (2.14) is second order convergent in both space and time directions.

Example 4.2**.**

Consider the two-dimensional space-Riesz fractional wave equation (1.1), on a finite domain $0<x<1,\,0<y<1$, $0<t\leq 1/2$ with the variable coefficients

[TABLE]

and the initial conditions $u(x,y,0)=\sin(1)x^{2}(1-x)^{2}y^{2}(1-y)^{2}$, $\frac{\partial}{\partial t}u(x,y,0)=\cos(1)x^{2}(1-x)^{2}y^{2}(1-y)^{2}$ with the zero Dirichlet boundary conditions on the rectangle. The exact solution of the PDEs is

[TABLE]

Using the above conditions, it is easy to obtain the forcing function $f(x,y,t)$.

Table 2 shows that the scheme (2.27) is second order convergent in both space and time directions.

5 Conclusion

In this work we have developed the energy method to estimate the two-dimensional space-Riesz fractional wave equation with the variable coefficients. To the best of our knowledge, the convergence and stability are lack of study for the one-dimensional space-Riesz fractional wave equation with the nonzero conditions. In this paper, the priori error estimates have been established and the convergence analysis and stability of the proposed method have been proved. For two-dimensional cases with the variable coefficients, the discretized matrices are proved to be commutative, which ensures to carry out of the priori error estimates. Numerical results have been given to illustrate the robustness and efficiency of the presented method with the second order convergence. We remark that though this current paper focus on the space-Riesz fractional wave equation, the energy estimates is still valid for the compact finite difference schemes and the nonlocal wave equation [11].

Acknowledgments

This work was supported by NSFC 11601206, the Fundamental Research Funds for the Central Universities under Grant No. lzujbky-2016-105, and SIETP 201710730065.