Numerical algorithm for two-dimensional time-fractional wave equation of distributed-order with a nonlinear source term

TL;DR

This paper develops an ADI difference scheme for solving a complex two-dimensional time-fractional wave equation with distributed order and nonlinear source, ensuring stability, convergence, and verified accuracy through numerical tests.

Contribution

It introduces a novel ADI difference scheme for a complex fractional wave equation with distributed order and nonlinear source, analyzing its stability and convergence.

Findings

The scheme is unconditionally stable.

The method achieves high accuracy in numerical experiments.

Numerical results confirm the effectiveness of the proposed algorithm.

Abstract

In this paper, an alternating direction implicit (ADI) difference scheme for two-dimensional time-fractional wave equation of distributed-order with a nonlinear source term is presented. The unique solvability of the difference solution is discussed, and the unconditional stability and convergence order of the numerical scheme are analysed. Finally, numerical experiments are carried out to verify the effectiveness and accuracy of the algorithm.

Tables

Table 1. Table 1: Errors and convergence orders for Example  4.1 in temporal direction with h 1 = h 2 = π 500 subscript ℎ 1 subscript ℎ 2 𝜋 500 h_{1}=h_{2}=\frac{\pi}{500} and Δ ​ β = 1 160 Δ 𝛽 1 160 \Delta\beta=\frac{1}{160} .

$\tau$	${\Vert e^N\Vert }_{L^{\infty }}$	Order	${\Vert e^N\Vert }_{L^2}$	Order
1/10	0.0839	-	0.1225	-
1/20	0.0439	0.9344	0.0634	0.9502
1/40	0.0227	0.9515	0.0326	0.9596
1/80	0.0117	0.9526	0.0167	0.9650
1/160	0.0059	0.9877	0.0085	0.9743

Table 2. Table 2: Errors and convergence orders for Example  4.1 with an optimal step size ratio for τ 𝜏 \tau and Δ ​ β Δ 𝛽 \Delta\beta , and h 1 = h 2 = π 500 subscript ℎ 1 subscript ℎ 2 𝜋 500 h_{1}=h_{2}=\frac{\pi}{500} .

$\tau$	$\Delta \beta$	${\Vert e^N\Vert }_{L^{\infty }}$	Order	${\Vert e^N\Vert }_{L^2}$	Order
1/100	1/10	0.0093	-	0.0133	-
1/400	1/20	0.0024	1.9542	0.0034	1.9678
1/1600	1/40	6.0481e-04	1.9885	8.6411e-04	1.9762
1/6400	1/80	1.4751e-04	2.0357	2.1076e-04	2.0365

Table 3. Table 3: Errors and convergence orders for Example  4.1 with an optimal step size ratio for τ 𝜏 \tau , h 1 subscript ℎ 1 h_{1} , h 2 subscript ℎ 2 h_{2} , and Δ ​ β Δ 𝛽 \Delta\beta .

$\tau$	$h_1=h_2$	$\Delta \beta$	${\Vert e^N\Vert }_{L^{\infty }}$	Order	${\Vert e^N\Vert }_{L^2}$	Order
1/64	$\pi /2$	1/8	0.4602	-	0.7230	-
1/256	$\pi /4$	1/16	0.1195	1.9453	0.1689	2.0978
1/1024	$\pi /8$	1/32	0.0301	1.9892	0.0426	1.9872
1/4096	$\pi /16$	1/64	0.0075	2.0048	0.0107	1.9932
1/16384	$\pi /32$	1/128	0.0019	1.9809	0.0027	1.9866
1/65536	$\pi /64$	1/256	4.7098e-04	2.0123	6.6801e-04	2.0150

Equations

$$\displaystyle\int_{1}^{2}p(\beta){}_{0}^{C}D_{t}^{\beta}u(x,y,t)d\beta=\frac{\partial^{2}u(x,y,t)}{\partial x^{2}}+\frac{\partial^{2}u(x,y,t)}{\partial x^{2}}+f\big{(}x,y,t,u(x,y,t)\big{)},$$

$$(x,y)\in \Omega ,t\in (0,T],$$

$$u(x,y,t)=\varphi (x,y,t),(x,y)\in \partial \Omega ,0\le t<T,$$

$$u(x,y,0)={\psi }_1(x,y),u_t(x,y,0)={\psi }_2(x,y),(x,y)\in \Omega ,$$

$${}_0^C{D}_t^{\beta }v(t)=\{\begin{array}{rl}& \frac{\partial v(t)}{\partial t}-\frac{\partial v(0)}{\partial t},\beta =1,\\ & \frac{1}{\Gamma (2-\beta )}{\int }_0^t(t-\xi {)}^{1-\beta }\frac{{\partial }^{2}v(\xi )}{\partial {\xi }^2}d\xi ,1<\beta <2,\\ & \frac{{\partial }^{2}v(t)}{\partial t^2},\beta =2,\end{array}$$

$$|f(x,y,t,u_1)-f(x,y,t,u_2)|\le L_f|u_1-u_2|,$$

$$u_{ij}^{n-\frac{1}{2}}=\frac{1}{2}(u_{ij}^n+u_{ij}^{n-1}),{\delta }_t{u}_{ij}^{n-\frac{1}{2}}=\frac{1}{\tau }(u_{ij}^n-u_{ij}^{n-1}),$$

$${\delta }_x{u}_{i-\frac{1}{2},j}^n=\frac{1}{h_1}(u_{ij}^n-u_{i-1,j}^n),{\delta }_x^2{u}_{ij}^n=\frac{1}{h_1}({\delta }_x{u}_{i+\frac{1}{2},j}^n-{\delta }_x{u}_{i-\frac{1}{2},j}^n),$$

$${\delta }_y{u}_{i,j-\frac{1}{2}}^n=\frac{1}{h_2}(u_{ij}^n-u_{i,j-1}^n),{\delta }_y^2{u}_{ij}^n=\frac{1}{h_2}({\delta }_x{u}_{i,j+\frac{1}{2}}^n-{\delta }_x{u}_{i,j-\frac{1}{2}}^n),$$

$${\Delta }_h{u}_{ij}={\delta }_x^2{u}_{ij}+{\delta }_y^2{u}_{ij}.$$

$${\int }_1^{2}p(\beta ){}_0^C{D}_t^{\beta }u(x_i,y_j,t_n)d\beta$$

$$=$$

$$\displaystyle\frac{1}{2}\bigg{(}\int_{1}^{2}p(\beta){}_{0}^{C}D_{t}^{\beta}u(x_{i},y_{j},t_{n})d\beta+\int_{1}^{2}p(\beta){}_{0}^{C}D_{t}^{\beta}u(x_{i},y_{j},t_{n-1})d\beta\bigg{)}$$

$$=$$

$$\displaystyle+\frac{1}{2}\bigg{[}f\big{(}x_{i},y_{j},t_{n},u(x_{i},y_{j},t_{n})\big{)}+f\big{(}x_{i},y_{j},t_{n-1},u(x_{i},y_{j},t_{n-1})\big{)}\bigg{]}.$$

$${\int }_1^{2}p(\beta ){}_0^C{D}_t^{\beta }{U}_{ij}^{n-\frac{1}{2}}d\beta =$$

$$\displaystyle+\frac{1}{2}\bigg{[}f\big{(}x_{i},y_{j},t_{n},U_{ij}^{n}\big{)}+f\big{(}x_{i},y_{j},t_{n-1},U_{ij}^{n-1}\big{)}\bigg{]}$$

$$\Delta \beta \sum _{l=1}^{K}p({\beta }_l){}_0^C{D}_t^{{\beta }_l}{U}_{ij}^{n-\frac{1}{2}}+R_1=\frac{{\partial }^2}{\partial x^2}{U}_{ij}^{n-\frac{1}{2}}+\frac{{\partial }^2}{\partial y^2}{U}_{ij}^{n-\frac{1}{2}}$$

$$\displaystyle+\frac{1}{2}\bigg{[}f\big{(}x_{i},y_{j},t_{n},U(x_{i},y_{j},t_{n})\big{)}+f\big{(}x_{i},y_{j},t_{n-1},U(x_{i},y_{j},t_{n-1})\big{)}\bigg{]},$$

$${}_0^C{D}_t^{{\beta }_l}{U}_{ij}^{n-\frac{1}{2}}$$

$$=$$

$$a_k^{({\beta }_l)}=(k+1{)}^{2-{\beta }_l}-k^{2-{\beta }_l},k=0,1,2,\cdots ,$$

$$\mid R_2^l\mid \le$$

$$\underset{0\le t\le t_n}{\max }\mid \frac{{\partial }^{3}u(x_i,y_j,t)}{\partial t^3}\mid {\tau }^{3-{\beta }_l},l=1,2,\cdots ,K.$$

$$\frac{{\partial }^{2}g(x_i)}{\partial x^2}=\frac{g(x_{i+1})-2g(x_i)+g(x_{i-1})}{(\Delta x{)}^2}-\frac{(\Delta x{)}^2}{12}\frac{{\partial }^{4}g({\xi }_i)}{\partial x^4},{\xi }_i\in (x_{i-1},x_{i+1})$$

$$\displaystyle\Delta\beta\sum_{l=1}^{K}p(\beta_{l})\frac{\tau^{1-\beta_{l}}}{\Gamma(3-\beta_{l})}\bigg{[}a_{0}^{(\beta_{l})}\delta_{t}U_{ij}^{n-\frac{1}{2}}-\sum_{k=1}^{n-1}\big{(}a_{n-k-1}^{(\beta_{l})}-a_{n-k}^{(\beta_{l})}\big{)}\delta_{t}U_{ij}^{k-\frac{1}{2}}$$

$$\displaystyle-a_{n-1}^{(\beta_{l})}\psi_{2}(x_{i},y_{j})\bigg{]}+\sum_{l=1}^{K}\Delta\beta p(\beta_{l})R_{2}^{l}+R_{1}$$

$$=$$

$$f(x_i,y_j,t_n,U_{ij}^n)=f(x_i,y_j,t_{n-1},U_{ij}^{n-1})+\mathcal{O}(\tau ).$$

$$\displaystyle\Delta\beta\sum_{l=1}^{K}p(\beta_{l})\frac{\tau^{1-\beta_{l}}}{\Gamma(3-\beta_{l})}\bigg{[}a_{0}^{(\beta_{l})}\delta_{t}U_{ij}^{n-\frac{1}{2}}-\sum_{k=1}^{n-1}\big{(}a_{n-k-1}^{(\beta_{l})}-a_{n-k}^{(\beta_{l})}\big{)}\delta_{t}U_{ij}^{k-\frac{1}{2}}-a_{n-1}^{(\beta_{l})}\psi_{2}(x_{i},y_{j})\bigg{]}$$

$$\displaystyle=\delta_{x}^{2}U_{ij}^{n-\frac{1}{2}}+\delta_{y}^{2}U_{ij}^{n-\frac{1}{2}}+f\big{(}x_{i},y_{j},t_{n-1},U_{ij}^{n-1}\big{)}+R_{ij}^{n-\frac{1}{2}}+\widetilde{R}_{ij}^{n-\frac{1}{2}},$$

$$R_{ij}^{n-\frac{1}{2}}=-\sum _{l=1}^K\Delta \beta p({\beta }_l)R_2^l+\mathcal{O}(h_1^2+h_2^2)+\mathcal{O}(\Delta {\beta }^2)$$

$${\tilde{R}}_{ij}^{n-\frac{1}{2}}=\mathcal{O}(\tau ).$$

$$\bigg{|}-\sum_{l=1}^{K}\Delta\beta p(\beta_{l})R_{2}^{l}\bigg{|}\leq C_{1}\tau^{1+\frac{1}{2}\Delta\beta}\sum_{l=1}^{K}\Delta\beta p(\beta_{l}).$$

$$\sum _{l=1}^K\Delta \beta p({\beta }_l)\sim {\int }_1^{2}p(\beta )d\beta =c_0,$$

$$\sum _{l=1}^K\Delta \beta p({\beta }_l)\le C_2,$$

$$\Big{|}R_{ij}^{n-\frac{1}{2}}\Big{|}\leq C_{3}\left(\tau^{1+\frac{1}{2}\Delta\beta}+h_{1}^{2}+h_{2}^{2}+\Delta\beta^{2}\right).$$

Full text

1 Introduction

The idea of distributed-order differential equation was first introduced by Caputo in his work for modeling the stress-strain behavior of an anelastic medium in 1960s [1]. Being different from the differential equations with the single-order fractional derivative and the ones with sums of fractional derivatives, i.e., multi-term fractional differential equations (FDEs), the distributed-order differential equations are derived by integrating the order of differentiation over a certain range. It can be regarded as a generalization of the aforementioned two classes of FDEs. A typical application of this kind of FDEs is in the retarding sub-diffusion process, where a plume of particles spreads at a logarithmic rate, which leads to ultraslow diffusion (see [2][3][4]). Another example is the fractional Langevin equation of distributed-order, which was proposed to model the kinetics of retarding sub-diffusion whose scaling exponent decreases in time, and then was applied to simulate the strongly anomalous ultraslow diffusion with the mean square displacement growing as a power of logarithm of time [5]. The distributed-order FDEs were also found playing important role in other various research fields, such as control and signal processing [6], modelling dielectric induction and diffusion [7], identification of systems [8], and so on.

Till now, there have been many important progresses for the research on analytical solutions of distributed-order FDEs. For the kinetic description of anomalous diffusion and relaxation phenomena, A. V. Chechkin et al. presented the diffusion-like equation with time fractional derivative of distributed-order in [9], where the positivity of the solutions of the proposed equation was proved and the relation to the continuous-time random walk theory was established. T. M. Atanackovic et al. analysed a Cauchy problem for a time distributed-order diffusion-wave equation by means of the theory of an abstract Volterra equation [10]. In [11], for the one-dimensional distributed-order diffusion-wave equation, R. Gorenflo et al. gave the interpretation of the fundamental solution of the Cauchy problem as a probability density function of the space variable $x$ evolving in time $t$ in the transform domain by employing the technique of the Fourier and Laplace transforms. Using the Laplace transform method, Z. Li et al investigated the asymptotic behavior of solutions to the initial-boundary-value problem for the distributed-order time-fractional diffusion equations [12].

In most instances, the analytical solutions of distributed-order differential equations are not easy to available, thus it stimulates researchers to develop numerical algorithms for approximate solutions. To our knowledge, the research on numerically solving the distributed-order differential equations are still in its infancy. The literatures [13][14][15] concerned on developing numerical methods for solving distributed-order ordinary differential equations. In terms of the distributed-order partial differential equations, most of the work are about the one-dimensional time distributed-order differential equations, and the integrating range of the order of time derivative is the interval $[0,1]$, which is named as time distributed-order diffusion equation. N. J. Ford et al. developed an implicit finite difference method for the solution of the diffusion equation with distributed order in time [16]. By using the Grünwald-Letnikov formula, Gao et al. proposed two difference schemes to solve the one-dimensional distributed-order differential equations, and the extrapolation method was applied to improve the approximate accuracy [17]. In [18], the authors handled the same distributed-order differential equations by employing a weighted and shifted Grünwald-Letnikov formula to derive several second-order convergent difference schemes. When the order of the time derivative is distributed over the interval $[1,2]$, it is called the time distributed-order wave equation. The study of the numerical solution of this kind of equation is rather more limited. Ye et al. derived and analysed a compact difference scheme for a distributed-order time-fractional wave equation in [19].

When considering the high-dimensional models, Gao et al. investigated ADI schemes for two-dimensional distributed-order diffusion equations [20][21], and they also developed two ADI difference schemes for solving the two-dimensional time distributed-order wave equations [22]. Due to the widespread use of the nonlinear models [23][24], M. L. Morgado et al. developed an implicit difference scheme for one-dimensional time distributed-order diffusion equation with a nonlinear source term [25]. For further discussion on the numerical approaches for solving the high-dimensional distributed-order partial differential equations, this paper is devoted to develop effective numerical algorithm for two-dimensional time-fractional wave equation of distributed-order with a nonlinear source term

[TABLE]

where $\Omega=(0,L_{1})\times(0,L_{2})$, and $\partial\Omega$ is the boundary of $\Omega$. The fractional derivative ${}_{0}^{C}D_{t}^{\beta}v(t)$ in (1) is given in the Caputo sense

[TABLE]

and the function $p(\beta)$ is served as weight for the order of differentiation such that $p(\beta)>0$ and $\int_{1}^{2}p(\beta)d\beta=c_{0}>0$. We assume that $p(\beta)$, $\phi(x,y,t)$, $\psi_{1}(x,y)$, $\psi_{2}(x,y)$ and $f(x,y,t,u)$ are continuous, and the nonlinear source term $f$ satisfies a Lipschitz condition of the form

[TABLE]

where $L_{f}$ is a positive constant.

The main procedure of developing numerical scheme for solving problem (1)$-$(1.3) is as follows. Firstly a suitable numerical quadrature formula is adopted to discrete the integral in (1), and a multi-term time fractional wave equation is left whereafter. Then we develop an ADI finite difference scheme which is uniquely solvable for the multi-term time fractional wave equation. By using the discrete energy method, we prove the derived numerical scheme is unconditionally stable and convergent.

The rest of this paper is organized in the following way. In Section 2, the ADI finite difference scheme is constructed and described detailedly. In Section 3, we give analysis on solvability, stability and convergence for the derived difference scheme. Numerical results are illustrated in Section 4 to confirm the effectiveness and accuracy of our method, and some conclusions are drawn in the last section.

2 The derivation of the ADI scheme

This section focuses on deriving the ADI scheme for the problems (1)$-$(1.3).

Let ${M_{1}}$, ${M_{2}}$ and $N$ be positive integers, and $h_{1}=L_{1}/M_{1}$, $h_{2}=L_{2}/M_{2}$ and $\tau=T/N$ be the uniform sizes of spatial grid and time step, respectively. Then a spatial and temporal partition can be defined as $x_{i}=ih_{1}$ for $i=0,1,\cdots,M_{1}$, $y_{j}=jh_{2}$ for $j=0,1,\cdots,M_{2}$ and $t_{n}=n\tau$ for $n=0,1,\cdots,N$. Denote $\bar{\Omega}_{h}=\{(x_{i},y_{j})\mid 0\leq i\leq M_{1},0\leq j\leq M_{2}\}$ and $\Omega_{\tau}=\{t_{n}\mid t_{n}=n\tau,0\leq n\leq N\}$, then the domain $\bar{\Omega}\times[0,T]$ is covered by $\bar{\Omega}_{h}\times\Omega_{\tau}$. Let $u=\{u_{ij}^{n}\mid 0\leq i\leq M_{1},0\leq j\leq M_{2},0\leq n\leq N\}$ be a grid function on $\bar{\Omega}_{h}\times\Omega_{\tau}$. We introduce the following notations:

[TABLE]

[TABLE]

[TABLE]

and

[TABLE]

Consider Eq. (1) at the point $(x_{i},y_{j},t_{n})$, and we write it as

[TABLE]

Take an average of Eq. (2.1) on time level $t=t_{n}$ and $t=t_{n-1}$, then we have

[TABLE]

Denote by $U_{ij}^{n}=u(x_{i},y_{j},t_{n})$ the grid functions on $\bar{\Omega}_{h}\times\Omega_{\tau}$ with $0\leq i\leq M_{1}$, $0\leq j\leq M_{2}$, $0\leq n\leq N$. Eq. (2.2) can be expressed as

[TABLE]

Firstly we discretize the integral term in (2.3). Suppose $p(\beta)\in C^{2}[1,2]$, ${}_{0}^{C}D_{t}^{\beta}u(x_{i},y_{j},t)|_{t=t_{n-1}}$ and ${}_{0}^{C}D_{t}^{\beta}u(x_{i},y_{j},t)|_{t=t_{n}}\in C^{2}[1,2]$. Let $K$ be a positive integer, and $\Delta\beta=1/K$ be the uniform step size. Take $\beta_{l}=1+\frac{2l-1}{2}\Delta\beta$, $1\leq l\leq K$, then the mid-point quadrature rule is used for approximating the integral in (2.3)

[TABLE]

where $R_{1}=\mathcal{O}(\Delta\beta^{2})$.

Next, we solve the multi-term time fractional wave equation (2.4) with the initial and boundary conditions (1.3) and (1.2). Suppose $u(x,y,t)\in C_{x,y,t}^{4,4,3}(\bar{\Omega}\times[0,T])$. According to Theorem 8.2.5 in [26], the Caputo derivative ${}_{0}^{C}D_{t}^{\beta_{l}}U_{ij}^{n-\frac{1}{2}}$, $1<\beta_{l}<2$ have the fully discrete difference scheme

[TABLE]

where

[TABLE]

and

[TABLE]

In the meantime, using the second order finite difference

[TABLE]

to approximate the second order derivatives in (2.4), it is obtained

[TABLE]

where $R_{3}=\mathcal{O}(h_{1}^{2}+h_{2}^{2})$. Subsequently, the nonlinear source term is dealt with in the following manner to avoid a system of nonlinear equations when computing:

[TABLE]

Substituting (2.8) in (2.7), we are left with

[TABLE]

where

[TABLE]

and

[TABLE]

From (2.6), we can deduce that there exists a positive constant $C_{1}$ such that

[TABLE]

Since

[TABLE]

we get

[TABLE]

where $C_{2}$ is a positive constant. Thus there exists a positive constant $C_{3}$ such that

[TABLE]

Besides, it is obvious that

[TABLE]

where $C_{4}$ is a positive constant.

Denote

[TABLE]

Since

[TABLE]

it can be concluded that

[TABLE]

In addition, $|\ln\tau|\leq C\tau^{-\varepsilon}$ for any positive and small $\varepsilon$ when $\tau$ is sufficiently small, thus the term $\mathcal{O}(\tau^{2}|\ln\tau|)$ is almost the same as $\mathcal{O}(\tau^{2})$ when $\tau$ is sufficiently small. Adding the high order term

[TABLE]

on both sides of (2.9), we derive

[TABLE]

where

[TABLE]

and it is clear that

[TABLE]

Also, for the initial and boundary value conditions, we have

[TABLE]

[TABLE]

Let $u_{ij}^{n}$ be the numerical approximation to $u(x_{i},y_{j},t_{n})$. Neglecting the small term $R_{ij}^{n-\frac{1}{2}}$, $\widetilde{R}_{ij}^{n-\frac{1}{2}}$ and $\widehat{R}_{ij}^{n-\frac{1}{2}}$ in (2.10), and using $u_{ij}^{n}$ instead of $U_{ij}^{n}$ in (2.10)$-$(2.12), we construct the difference scheme for (1)$-$(1.3) as follows:

[TABLE]

where

[TABLE]

and

[TABLE]

Notice $a_{0}^{(\beta_{l})}=1$, then Eq. (2) can be rewritten as:

[TABLE]

or

[TABLE]

where $I$ denotes the identity operator.

Let

[TABLE]

Together with (2.14) and (2.15) the ADI difference scheme is derived, and the procedure can be executed as follows:

On each time level $t=t_{n}$ $(1\leq n\leq N)$, firstly, for all fixed $y=y_{j}$ $(1\leq j\leq M_{2}-1)$, solving a set of $M_{1}-1$ equations at the mesh points $x_{i}$ $(1\leq i\leq M_{1}-1)$ to get the intermediate solution $u_{ij}^{\ast}$:

[TABLE]

afterwards, for all fixed $x=x_{i}$ $(1\leq i\leq M_{1}-1)$, by computing a set of $M_{2}-1$ equations at the mesh points $y_{j}$ $(1\leq j\leq M_{2}-1)$, the solution $u_{ij}^{n}$ can be obtained:

[TABLE]

3 Analysis of the ADI difference scheme

3.1 Solvability

It is clear that the ADI scheme (2.16)$-$(2.17) is a linear tridiagonal system in unknowns, and the coefficient matrices are strictly diagonally dominant. Thus the scheme (2.16)$-$(2.17) has a unique solution. This result can be written as following.

Theorem 3.1**.**

The ADI difference scheme (2.16)$-$(2.17) is uniquely solvable.

3.2 Stability

In this subsection we prove the unconditional stability and the convergence of the difference scheme (2.16)$-$(2.17). We start with some auxiliary definitions and useful results.

Denote the space of grid functions on $\bar{\Omega}_{h}$

[TABLE]

For any grid function $v\in\mathcal{V}_{h}$, the following discrete norms and Sobolev seminorm are introduced:

[TABLE]

[TABLE]

[TABLE]

Lemma 3.1**.**

[27]** For any grid function $v\in\mathcal{V}_{h}$, $\|v\|\leq\frac{1}{2\sqrt{3}}|v|_{1}$.

Lemma 3.2**.**

[26]** For any grid function $v\in\mathcal{V}_{h}$, $|v|_{1}\leq\frac{1}{2\sqrt{3}}\|\Delta_{h}v\|$.

Lemma 3.3**.**

[26]** For any $G=\{G_{1},G_{2},G_{3},\ldots\}$ and $q$, we have

[TABLE]

where

[TABLE]

The discrete Gronwall’s inequality is also introduced below since it is necessary to prove the stability and convergence of the proposed method.

Lemma 3.4**.**

[28]** Assume that $k_{n}$ and $p_{n}$ are nonnegative sequences, and the sequence $\Phi_{n}$ satisfies

[TABLE]

where $g_{0}\geq 0$. Then the sequence $\Phi_{n}$ satisfies

[TABLE]

Since the ADI difference scheme (2.16)$-$(2.17) is equivalent to (2)$-$(2.15) if the intermediate variable $u^{\ast}$ is eliminated, we analyze the stability and convergence by employing the difference scheme (2)$-$(2.15).

Assume that $\widetilde{u}_{ij}^{n}$ is the approximate solution of $u_{ij}^{n}$, which is the exact solution of the scheme (2)$-$(2.15). Denote $\varepsilon_{ij}^{n}=u_{ij}^{n}-\widetilde{u}_{ij}^{n},\ 0\leq i\leq M_{1},\ 0\leq j\leq M_{2},\ 0\leq n\leq N$, then we have the perturbation error equations

[TABLE]

where

[TABLE]

Theorem 3.2**.**

Assume that the condition (1.4) is satisfied, then the difference scheme (2.16)$-$(2.17) is unconditionally stable.

Proof.

Let

[TABLE]

then Eq. (3.2) is equivalent to

[TABLE]

Multiplying (3.2) by $h_{1}h_{2}\tau\delta_{t}\varepsilon_{ij}^{n-\frac{1}{2}}$, summing up for $i$ from $1$ to $M_{1}-1$, for $j$ from $1$ to $M_{2}-1$ and for $n$ from $1$ to $m$, we analyze each term in the derived equation. Firstly, by employing Lemma 3.3, we have

[TABLE]

where

[TABLE]

Whereafter using the discrete Green formula, we get

[TABLE]

and

[TABLE]

Analogous to (3.5), it is also obtained

[TABLE]

On the basis of (1.4), there holds that

[TABLE]

From Equations (3.3)$-$(3.7), the inequality below is derived

[TABLE]

According to Lemma 3.1 and Lemma 3.2, we deduce from (3.8) that

[TABLE]

Finally, taking Lemma 3.4, it follows that

[TABLE]

This completes the proof. ∎

In the following we consider the convergence of the difference approximation. Noticing that $U_{ij}^{n}$ is the exact solution of the system (1)$-$(1.3) and $u_{ij}^{n}$ is the numerical solution of the difference scheme (2)$-$(2.15), we denote the error

[TABLE]

Subscribing (2)$-$(2.15) from (2.10)$-$(2.12), we get the error equations

[TABLE]

Theorem 3.3**.**

Suppose that the continuous problem (1)$-$(1.3) has solution $u(x,y,t)\in C_{x,y,t}^{4,4,3}(\bar{\Omega}\times[0,T])$. Then there is a positive constant $C$ such that

[TABLE]

Proof.

The proof of convergence is similar to that of Theorem 3.2. Multiplying (3.2) by $h_{1}h_{2}\tau\delta_{t}e_{ij}^{n-\frac{1}{2}}$, summing up for $i$ from $1$ to $M_{1}-1$, for $j$ from $1$ to $M_{2}-1$ and for $n$ from $1$ to $m$, we estimate each term in the resulted equation.

By using analogous strategies as (3.3)$-$(3.7), we get (3.10)$-$(3.14) correspondingly.

[TABLE]

[TABLE]

[TABLE]

[TABLE]

and

[TABLE]

As for the remainder, it is deduced that

[TABLE]

From (3.10)$-$(3.15) it follows that

[TABLE]

i.e.,

[TABLE]

According to Lemma 3.1, we obtain

[TABLE]

Therefore,

[TABLE]

where Lemma 3.4 is applied. This completes the proof. ∎

4 Numerical results

In this section, a numerical example is tested to demonstrate the effectiveness of the proposed scheme, and verify the theoretical results including convergence orders and numerical stability. The discrete $L^{2}$ and $L^{\infty}$ norms are both taken to measure the numerical errors. Denote

[TABLE]

and

[TABLE]

Example 4.1**.**

[TABLE]

whose analytical solution is known and is given by

[TABLE]

In Figure 1 we illustrate the relative error, which verifies the convergence of the algorithm we proposed.

In Figure 2 we present a comparison of the exact and numerical solutions. It can be seen that the numerical solution is in good agreement with the exact solution.

In Table 1, the numerical accuracy of difference scheme (2.16)$-$(2.17) in time is recorded. Let the step sizes $h_{1}$, $h_{2}$, and $\Delta\beta$ be fixed and small enough such that the dominated error arise from the approximation of the time derivatives. Varying the step sizes in time, the numerical errors in discrete both $L^{\infty}$ and $L^{2}$ norms and the associated convergence orders are shown in this table respectively, which can be found in agreement with the theoretical analysis.

In Table 2, we take the fixed and small enough step sizes in space, and adopt an optimal step size ratio in time and distributed order. As $\Delta\beta$ and $\tau$ vary, we compute the errors and convergence orders listed in the table, which indicates that the convergence order in time and distributed order are about one and two, respectively.

Table 3 displays the computational results with an optimal step size ratio in time, space and distributed order. We can conclude from this table that the convergence orders with respect to time, space and distributed order are approximately one, two and two, respectively, which is in good agreement with our theoretical results analyzed in Section 3.

5 Conclusion

In this paper, we construct efficient numerical scheme for solving two-dimensional time-fractional wave equation of distributed-order with a nonlinear source term, and provide the theoretical analysis on stability and convergence by the discrete energy method. Numerical results are provided by figures and tables, which show the algorithm proposed in this work is effective and feasible. In the future work, the promotion of computational efficiency will be considered so that the more complicated problems can be handled.

Acknowledgements

This research was supported by National Natural Science Foundations of China (No.11471262). The authors would like to express their gratitude to the referees for their very helpful comments and suggestions on the manuscript.