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

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.


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 singleorder 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.  [10]. In [11], for the one-dimensional distributed-order diffusion-wave equation, R. Gorenflo  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 meth-ods 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 integrat-  [16]. By using the Grünwald-Letnikov formula, Gao et al. proposed two difference schemes to solve the one-dimensional distributedorder 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 where Ω = (0, L 1 ) × (0, L 2 ), and ∂Ω is the boundary of Ω. The fractional and the function p(β) is served as weight for the order of differentiation such that p(β) > 0 and 2 1 p(β)dβ = c 0 > 0. We assume that p(β), φ(x, y, t), ψ 1 (x, y), ψ 2 (x, y) and f (x, y, t, u) are continuous, and the nonlinear source term f satisfies a Lipschitz condition of the form where L f is a positive constant.
The main procedure of developing numerical scheme for solving problem (1.1)−(1.3) is as follows. Firstly a suitable numerical quadrature formula is adopted to discrete the integral in (1.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.

The derivation of the ADI scheme
This section focuses on deriving the ADI scheme for the problems (1.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 τ = 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, · · · , M 1 , y j = jh 2 for j = 0, 1, · · · , M 2 and t n = nτ for n = 0, 1, · · · , N . DenoteΩ h = We introduce the following notations: Consider Eq. (1.1) at the point (x i , y j , t n ), and we write it as Take an average of Eq. (2.1) on time level t = t n and t = t n−1 , then we have 2]. Let K be a positive integer, and ∆β = 1/K be the uniform step size. Take β l = 1 + 2l−1 2 ∆β, 1 ≤ l ≤ K, then the mid-point quadrature rule is used for approximating the Next, we solve the multi-term time fractional wave equation ( In the meantime, using the second order finite difference to approximate the second order derivatives in (2.4), it is obtained . Subsequently, the nonlinear source term is dealt with in the following manner to avoid a system of nonlinear equations when computing: From (2.6), we can deduce that there exists a positive constant C 1 such that where C 2 is a positive constant. Thus there exists a positive constant C 3 such Besides, it is obvious that it can be concluded that .
In addition, | ln τ | ≤ Cτ −ε for any positive and small ε when τ is sufficiently small, thus the term O(τ 2 | ln τ |) is almost the same as O(τ 2 ) when τ is sufficiently small. Adding the high order term τ on both sides of (2.9), we derive and it is clear that Also, for the initial and boundary value conditions, we have Let u n ij be the numerical approximation to u(x i , y j , t n ). Neglecting the small term R Notice a (β l ) 0 = 1, then Eq. (2.13) can be rewritten as: where I denotes the identity operator. Let 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 ≤ n ≤ N ), firstly, for all fixed y = y j (1 ≤ j ≤ M 2 −1), solving a set of M 1 −1 equations at the mesh points afterwards, for all fixed x = x i (1 ≤ i ≤ M 1 − 1), by computing a set of M 2 − 1 equations at the mesh points y j (1 ≤ j ≤ M 2 − 1), the solution u n ij can be obtained:

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Ω
For any grid function v ∈ V h , the following discrete norms and Sobolev seminorm The discrete Gronwall's inequality is also introduced below since it is necessary to prove the stability and convergence of the proposed method. Assume that u n ij is the approximate solution of u n ij , which is the exact solution of the scheme (2.13)−(2.15). Denote Whereafter using the discrete Green formula, we get and (3.5) Analogous to (3.5), it is also obtained On the basis of (1.4), there holds that (3.7) From Equations (3.3)−(3.7), the inequality below is derived (3.8) According to Lemma 3.1 and Lemma 3.2, we deduce from (3.8) that Finally, taking Lemma 3.4, it follows that This completes the proof.
In the following we consider the convergence of the difference approximation.
Noticing that U n ij is the exact solution of the system (1.1)−(1.3) and u n ij is the numerical solution of the difference scheme (2.13)−(2.15), we denote the error Subscribing (2.13)−(2.15) from (2.10)−(2.12), we get the error equations and (3.14) As for the remainder, it is deduced that (3.15) From (3.10)−(3.15) it follows that i.e., According to Lemma 3.1, we obtain Therefore, where Lemma 3.4 is applied. This completes the proof.

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 ∞ norms are both taken to measure the numerical errors. Denote , and e N L ∞ := max Example 4.1.
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 ∆β 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 ∞ 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 ∆β and τ 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.

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.