Finite element approximation of fractional hyperbolic integro-di ff erential equation

: In this article, we propose a Galerkin finite element method for numerically solving a type of fractional hyperbolic integro-di ff erential equation, which can be considered as the generalization of the classical hyperbolic Volterra integro-di ff erential equation. Along with Galerkin finite element method in spatial direction, we apply a second order symmetric di ff erence method in time. Next we discuss the regularity analysis of the weak solution and convergence analysis of the semi-discrete scheme. Then we further study the stability analysis and the error estimation of the fully discrete problems, according to the properties of fractional Ritz-Volterra projection, Ritz projection and so on. Numerical examples with comparisons among the proposed schemes verify our theoretical analyses.


Introduction
The field of fractional calculus, which deals with mathematical analysis and physical applications of derivatives and integrals of arbitrary order, has become a hot topic nowadays. Fractional derivatives and integral can be considered as the generalization of integer order derivative. However, there are quite different between the fractional derivative and the integer order derivative. One of the reasons is that fractional operators, such as Riemann-Liouville derivative or Caputo derivative, has the nonlocal characteristics and weakly singular properties. But just because of above characteristics, the fractional calculus performs more perfectly than the classical calculus, especially in the field of anomalous diffusion problems. Many published papers reveal that fractional models, such as the fractional differential equations and the fractional integro-differential equations show more realistic dynamic behavior than the classical differential equations and the classical integro-differential ones.
Fractional integro-differential equations often describe the anomalous diffusion phenomena which come from the dynamic behaviors of viscoelastic materials, heat conduction with memory, and so on [1][2][3][4]. They can be divided into parabolic integro-differential equations and hyperbolic integro-differential equations along the time axis. The fractional hyperbolic integro-differential equations can be modelled by the generalized constitutive relations between stress σ and strain ϵ of the linear viscoelasticity [5]. If the corresponding generalized constitutive relation satisfies where a(t) is the stress relaxation modulus and E 0 is the Young's modulus. Then substituting (1.1) into its motion equation ρu tt = divσ + f , we can obtain a type of generalized hyperbolic integro-differential equations where u is the displacement, ρ is the density, f is the external force. If the kernel K belongs to the form of power law (e.g. the form t β Γ(1+β) , where Γ is the Gamma function), (1.2) is retained to the temporal fractional integro-differential equations, in which the power law widely exists in complex systems [6]. Meanwhile, if the displacement u of continuous media satisfies the 2α (1/2 < α < 1) order Lévy stable distribution in spatial direction, by applying the power law approximation form of its Fourier transform and inverse Fourier transform, Eq (1.2) turns to be the following fractional hyperbolic integro-differential equation where ∂ 2α ∂|x| 2α is the 2α order Riesz fractional derivative, often describes the 2α order Lévy flights [7]. The existence and uniqueness of the analytic solution of the above fractional hyperbolic integro-differential equation can be proved by using the Fourier transform and the Laplace transform. It is omitted here, because the proof is very similar to the corresponding parabolic problem. One can refer to [8,9].
There are several methods to study the fractional hyperbolic equations. Dassios and Font [10] studied the analytical solution of the time-fractional hyperbolic heat equation, in which the fractional derivatives contain three kinds of definitions. Kumar and Rai [11] presented a fractional hyperbolic bioheat transfer model and used a hybrid numerical scheme based on fractional Legendre wavelet method and finite difference scheme to study the numerical solution. Ashyralyev, Dal and Pinar [12] studied an initial boundary value problem for the fractional hyperbolic equation by difference scheme and discussed the stability details. Qiu et al. [13] constructed a formally second-order BDF finite difference scheme for a integro-differential equations with the multi-term kernels. Qiu et al. [14] presented a formally second-order backward differentiation formula for the Volterra integro-differential equation with a weakly singular kernel. As we all known, finite element method, finite difference method and spectral method are the classical numerical methods. They have been widely applied not only in the classical hyperbolic equations, but also in the fractional parabolic differential equations, e.g. [9,[15][16][17][18][19][20][21][22][23][24][25][26][27]. However, for the fractional hyperbolic integro-differental equation, there are relatively few. Moreover, it is quite difficult to get the high accuracy and the high convergence order methods, because the fractional integral and the fractional derivative are mixed into one term.
In this paper, we consider a Galerkin finite element method to solve the initial-boundary value problem of the fractional hyperbolic integro-differential equation together with the homogeneous Dirichlet boundary conditions and the initial conditions where J 1+β 0,t (0 < β < 1) represents the temporal Riemann-Liouville integral, κ is assumed to be a nonnegative constant, which represents the diffusion rate of particles. Here we use the time stepping method based on the symmetric difference approach (U n−1/2 + U n+1/2 )/2 to approximate u(t n ) of the Riesz derivative part in (1.4), and center difference approach (U n−1 − 2U n + U n+1 )/2 to approximate u tt (t n ), combining with the high order quadrature schemes based on the product trapezoidal formula to treat the Riemann-Liouville integral term. Meanwhile, we use the Galerkin finite element method with r − 1 order piecewise polynomial as the shape function in space. The expected goal of our convergence order is O(h r + k 2 ). Theoretical analyses and numerical experiments of the designed algorithm will be presented in the following paper.
We organize the following sections. In Section 2, we introduce the preliminary definitions and properties of fractional integral and fractional derivatives. In Section 3 and Section 4, we construct a Galerkin finite element scheme for the fractional hyperbolic integro-differential equation, and then present the regularity analysis of the weak solution and convergence analysis of the semi-discrete scheme separately. In Section 5, we derive the fully discrete scheme based on the symmetric difference method in time direction. Then we further discuss the stability analysis and error estimate of the fully discrete scheme. In Section 6, we present some numerical experiments to illustrate the efficiency of the theoretical analyses.

Preliminaries
We first introduce some definitions and notations of fractional calculus [28,29], and some properties of fractional derivative space as well.
Definition 2.1 The Riemann-Liouville integral of order α is defined as where α > 0.

Definition 2.2
The Riesz derivatives of order α is defined as where the left and right Riemann-Liouville derivatives are defined as The fractional Sobolev space needs to be denoted. On the fractional Sobolev space H α (Ω), we denote the semi-norm | · | α by and or One can prove that they are equivalent to each other [15,17,18,20,29]. Then we define the norm ∥ · ∥ α by ∥u∥ α = (∥u∥ 2 + |u| 2 α ) 1/2 .

The weak formulation
According to the above Lemma 2.2, the original problem turns to be the following weak form with the homogeneous Dirichlet boundary conditions and the initial conditions Grönwall's inequality is important to our regularity analysis and error estimate.
where x(t) ≥ 0 and f (t) are absolutely integrable, then Now we consider the regularity of the weak solution. We always denote C as a generic positive constant which may be changed at different situations from now. Theorem 3.1 The weak solution u(x, t) in (3.1) satisfies the following energy estimate i.e., By integration in t and using Lemma 2.2, we get where the initial values u 0 and u 1 are defined in (1.5). Then integrating by parts obtains For the first term of the right hand side of (3.6), we have And for the second one, we get For the discussion of the right hand side of (3.7), we can prove it similar to (3.8) and (3.9), which is omitted here. For the last term of the right hand side of (3.5), we get Hence, we obtain . Then employing Grönwall's inequality (Lemma 3.1) gets that Thus we finish the proof of Theorem 3.1. Remark 1. The existence and uniqueness of the weak solution for (3.1) can be derived from Theorem 3.1. The higher order regularities of the weak solution can be guaranteed only by higher differentiabilities of the data and compatibilities, which are omitted here.

Error estimation of the semi-discretization
In this section, we will give the error estimate for the following semi-discretization for v ∈ S h r , and S h r = {v ∈ H α 0 (Ω) ∩ C(Ω), ∀v| r ∈ P r−1 (K)} is the finite element subspace, in which h is the stepsize in space variable x, and P r−1 (K) is the r − 1 degree polynomial space on the subinterval K of Ω.
In order to carry out the work of error estimation, we first define the projection P h u ∈ S h r satisfying Obviously, we have that ∥P h u − u∥ ≤ C∥u∥.
Then we define R h u as the Ritz projection of function u ∈ H α 0 (Ω) ∩ H r (Ω) satisfying The optimal error estimates about the P h u and R h u are very useful for our later discussion.
Remark 2. The initial value u 0h ∈ S h r is an approximation of u 0 , which can take P h u 0 or R h u 0 , but here it satisfies Then we define the fractional Ritz-Volterra projection (4.7) Here we will use the fractional Ritz-Volterra projection V h and give some useful lemmas to study the error estimate.
The following energy estimation of ρ tt is novel and crucial to our convergence analysis.
Proof. By using the convolution property of Riemann-Liouville integral, we can rewrite (4.7) as follows By differentiating with respect to t of both sides of above equation, it becomes Then by differentiating with respect to t again, (4.10) Taking v = V h u tt − R h u tt and using the definition of Ritz projection (4.3) obtain that Hence, Therefore Then by using the Grönwall's inequality, one has and finally t 0 ∥ρ tt (s)∥ α ds ≤ Ch r−α ( t β Γ(1+β) ∥u 0 ∥ r + t β+1 Γ(2+β) ∥u 1 ∥ r + t 0 ∥u tt ∥ r ds). Consider now the L 2 −estimate. According to the characteristics of Riesz derivatives, about the left and right Riemann-Liouville derivatives, and also using Eq (4.10), we get that where ω satisfies ∂ 2α ∂|x| 2α ω = y in H α 0 (Ω), and ∥ω∥ 2α ≤ C.
in which ∥R h ω − ω∥ α ≤ h α ∥ω∥ 2α is used. By using the Grönwall's inequality again (Lemma 3.1), one has After integration, we get Thus the proof is completed.
∥u tt (s)∥ r ds). Proof. By using V h u(t) ∈ S h r as an intermediate function, the error can be defined as ε = u h − u = (u h − V h u) + (V h u − u) = θ + ρ. Therefore, the error equation can be rewritten as where Lemma 2.1 is used for the computation of the Riemann-Liouville integral. By integration in t and using Lemma 4.4, we get Hence, we obtain ∥θ t ∥ + ∥θ∥ α ≤ C(∥θ t (0)∥ + J 1+β 0,t ∥θ(t)∥ α + t 0 ∥ρ tt ∥ds). Then we get in which the Grönwall's inequality is used. By using Lemma 4.3, we have is used. This finishes the error estimates of θ. Then combining Lemma 4.1 for error estimates of ρ, the proof of Theorem 4.1 is completed.

Fully discrete scheme based on a symmetric difference approximation
We now turn to discuss the fully discrete scheme based on a symmetric difference approximation. Let f n = f (t n ), t n = nk, n = 1, 2, · · · , N, where k = T/N is the steplength in time variable t. Denote by U n the approximation solution and by ∂U n = (U n+1 − U n )/k, ∂U n = (U n − U n−1 )/k the forward and backward difference quotient of U n respectively, then ∂∂U n = ∂(U n − U n−1 )/k = (U n+1 − 2U n + U n−1 )/k 2 is the center difference quotient of second order to approximate the second time derivative term u tt . We also denote the average There are several ways to approximate the fractional integral. Here we select the product trapezoidal technique to approximate the fractional integral J 1+β 0,t n g(t), 1 ≤ n ≤ N, under the condition g(t) ∈ C 2 ([0, T ]). The following is the truncate error estimation.
The discrete Grönwall's inequality is very important to our stability analysis and the convergence analysis. Lemma 5.2. (Discrete Grönwall's inequality) ( [33][34][35][36]) Assume that ω n ≥ 0, f n ≥ 0 and that for n = 0, 1, · · · , y n ≥ 0 satisfies Next we denote by I h u the piecewise polynomial interpolation operator of u in S h r satisfying I h u(t n ) = u(t n ), n = 0, 1, · · · , N. Now we give the fully discrete scheme with given initial values U 0 and U 1 ∈ S h r , and σ n (g) = k 1+β Γ(3+β) n j=0 b n− j g(t j ). In fact, the above scheme is equivalent to the following form where U −1/2 can be approximated by U 0 and U 1 , and the second-order accuracy should be guaranteed. Furthermore, we can move the terms U n+1 in the right hand side of (5.2) to the left hand side for explicit processing.
Then we discuss the stability analysis of the fully discrete scheme (5.1) in the following form. Theorem 5.1. The solution of (5.1) satisfies the following stability conclusion . Multiplying both sides of (5.4) by k and summing from n = 1 to N obtain For the term k| N n=1 I n 1 |, we have in which b n− j = (n − j + 1) 2+β − 2(n − j) 2+β + (n − j − 1) 2+β is monotonically increasing for n = 1, · · · , N. And for the term k| N n=1 I n 2 |, we have After some adjustments and applying the discrete Grönwall's inequality (Lemma 5.2) obtain Then we finish the proof of the theorem. Next we discuss the error estimate of the fully discrete scheme. Denote by the error U n − u n = (U n − V h u n ) + (V h u n − u n ) = θ n + ρ n , we have where θ 0 = u 0 − U 0 , θ 1 = u 1 − U 1 , and e n = e n 1 + e n 2 + e n 3 + e n 4 , e n 1 = ∂∂ρ n , e n 2 = ∂∂u(t n ) − u tt (t n ), e n To construct the discrete initial values U 0 and U 1 , let u 2 = u tt (0) = f (0) − ∂ 2α ∂|x| 2α u 0 and define U 1 = V 0 + kV 1 + V 2 k 2 /2, where U 0 = V 0 , V j = P h u j , j = 0, 1, 2. According to Lemma 4.1, we get the following conclusions Then for θ n = U n − V h u n , there is In the following, we give the error estimation of the fully discrete scheme (5.1). Theorem 5.2. The solution U of (5.1) and the solution u of (1.4) at t n+1/2 satisfy the following conclusion Proof. The proof is similar to Theorem 5.1. Taking v = ∂θ n+1/2 in (5.6), we have where (∂∂θ n , ∂θ n+1/2 ) = (∂θ n − ∂θ n−1 , ∂θ n + ∂θ n−1 )/2k = 1 2 ∂∥∂θ n ∥ 2 , . For the first term of the right sides of (5.9), which is similar to (5.5), we obtain Summing (5.9) from n = 1 to N obtains that Then applying Lemma 5.2, we get ∥e n ∥. (5.10) By the Taylor expansion, we know that Combined with (5.7) and (5.10), we finish the proof of Theorem 5.2 Remark 3. The above stability analysis and the convergence analysis of the fully discrete schemes can be extended to high-dimensional cases without difficulty, which are omitted here.

Numerical experiments
In order to test the effectiveness of the designed numerical algorithm, we present the following numerical experiments in this section.
In the Galerkin finite element approximation, we select the hat function as the shape function, followed by the symmetric difference scheme and fractional trapezoidal formula for the time stepping. Through the theoretical analysis of the previous sections, the expected goal of the convergence order with L 2 norm should be O(k 2 + h 2 ). Example 6.1. In this example, we study the following fractional hyperbolic equation , then the exact solution is u(x, t) = t 2 (1 − x)(1 + x), which determines that the initial values u 0 = u 1 = 0. Table 1 shows the numerical results and convergence rates of case I, which supports the predicted rates of the convergence. Figure 1 shows the exact solution and the numerical solution of case I with α = 0.9, β = 0.1, h = 1/64, k = 0.01h at t = 0.5, x ∈ [−1, 1]. Figure 2 shows the error between the exact solution and the numerical solution of case I with α = 0.9, β = 0.1, h = 1/64, k = 0.01h, at x = 0, t ∈ [0, 1]. And Figure 3 shows the numerical solution of case I under the same conditions.    Case II. We choose the exact solution u(x, t) = sin(πt)(1 − x)(1 + x). And the source term f is obtained numerically by using the fractional trapezoidal formula. Then the initial values satisfy u 0 = 0, u 1 = π(1 − x)(1 + x). Table 2 shows the numerical results and convergence rates of case II, which support the predicted rates of the convergence. Figure 4 shows the exact solution and the numerical solution of case II with α = 0.9, β = 0.1, h = 1/64, k = 0.01h, at t = 0.5, x ∈ [−1, 1]. Figures 5 and 6 show the error between the exact solution and the numerical solution, the numerical solution of case II with α = 0.9, β = 0.1, h = 1/64, k = 0.01h for x ∈ [−1, 1], t ∈ [0, 1] separately.     Example 6.2. In this example, we consider the following fractional hyperbolic equation with homogeneous Dirichlet boundary conditions in Ω = [0, 1], T = 1. We choose the source term f = 2(1 − x)x α − (t 2 + 2t 3+β Γ(4+β) )D 2α a,x (1 − x)x α , then the exact solution is u(x, t) = t 2 (1 − x)x α , which has a weak singularity at the boundary point x = 0 if 0.5 < α < 1. Table 3 shows the errors and convergence rates with parameters α = 0.6, β = 0.1, h = 1/64, k = 0.01h for x ∈ [0, 1], t ∈ [0, 1]. Figures 7 and 8 show the numerical solution and the absolute error between the exact solution and the numerical solution of Example 6.2 with α = 0.9, β = 0.1, h = 1/64 for x ∈ [0, 1], t ∈ [0, 1] separately. And Figure 9 shows the numerical solution with different values of α at time t = 0.5. From Figures 8 and 9, we can see that the numerical solution is basically coincided with the exact solution. Note that the selected exact solution has a weak singularity at the boundary point x = 0, therefore the scheme does not work very well near zero.

Conclusions
In this paper, we use the Galerkin finite element method and the symmetric difference method to solve the fractional hyperbolic integro-differential equation, where the space fractional derivative is in Riesz sense and the integro-differential term is compounded of the Riesz space fractional derivative and the Riemann-Liouville time fractional integral. We apply the fractional trapezoidal formula to treat the fractional integral and employ enough points to ensure the convergence order. Numerical examples are presented to test the effectiveness of the convergence analysis. From the numerical results, we can see that the designed numerical algorithm performs well and the convergence orders conform to the convergence analysis.
As is known to all, fractional calculus has weak singularity and nonlocality from its origin [37]. Not only the fractional differential equation, but also the fractional integro-differential equation, their solutions both behave the weak singularities. In this paper, we design a solution with a weak singularity at the boundary point x = 0, which is verified by numerical experiments. Meanwhile, because of its nonlocality, although the above theoretical analyses can be extended to the high-dimensional cases without difficulty, the capacities of computation and memory will become large. So how to reduce the computationally expensive and the storage requirement comes into being the main problem. Maybe the fast algorithm is a good choice. In future, we will continue to study these problems.