Finite element approximation of fractional hyperbolic integro-differential equation

1.   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 $\sigma$ and strain $\epsilon$ 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 $\rho u_{tt} = {\mbox{div}} \sigma+ f$, we can obtain a type of generalized hyperbolic integro-differential equations

where $u$ is the displacement, $\rho$ is the density, $f$ is the external force. If the kernel $K$ belongs to the form of power law (e.g. the form $\frac{t^{\beta}}{\Gamma(1+\beta)}$, where $\Gamma$ 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\alpha \; (1/2 < \alpha < 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 $\frac{\partial^{2\alpha} }{\partial |x|^{2\alpha}}$ is the $2\alpha$ order Riesz fractional derivative, often describes the $2\alpha$ 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_{0, t}^{1+\beta}(0 < \beta < 1)$ represents the temporal Riemann-Liouville integral, $\kappa$ 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.

2.   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 $\alpha$ is defined as

where $\alpha > 0$.

Definition 2.2 The Riesz derivatives of order $\alpha$ is defined as

where the left and right Riemann-Liouville derivatives are defined as

where $C_\alpha = -\frac{1}{2\cos(\pi \alpha/2)}, n-1 < \alpha < n.$

The fractional Sobolev space needs to be denoted. On the fractional Sobolev space $H^\alpha(\Omega)$, we denote the semi-norm $|\cdot|_{\alpha}$ by

or

and or

One can prove that they are equivalent to each other ^{[15,17,18,20,29]}. Then we define the norm $\|\cdot\|_{\alpha}$ by

Lemma 2.1. (^[29]) The Riemann-Liouville integral operator has the following properties

where $p, q > 0$.

Lemma 2.2. (^[15]) There exist constants $C_1, C_2, C_3, C_4 > 0$ such that for any $u\in H_0^{2\alpha}(\Omega)$, and $v\in H_0^{\alpha}(\Omega)$,

3.   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.

Lemma 3.1. (Grönwall's inequality) (^[5]) Assume that the nonnegative function $y(t)$ satisfies the integral inequality

where $x(t)\geq 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

Proof. Setting $v = u_t$ in (3.1), it becomes

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

where

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.

4.   Error estimation of the semi-discretization

In this section, we will give the error estimate for the following semi-discretization

for $v\in S_r^h$, and $S_r^h = \{v\in H_0^{\alpha}(\Omega)\cap C(\Omega), \forall v|_r\in 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 $\Omega$.

In order to carry out the work of error estimation, we first define the projection $P_hu\in S_r^h$ satisfying

Obviously, we have that

Then we define $R_h u$ as the Ritz projection of function $u\in H_0^{\alpha}(\Omega)\cap H^{r}(\Omega)$ satisfying

The optimal error estimates about the $P_h u$ and $R_h u$ are very useful for our later discussion.

Lemma 4.1. (^[19]) Assume $u\in H_0^{\alpha}(\Omega)\cap H^{r}(\Omega)$, then we have that

Remark 2. The initial value $u_{0h}\in S_r^h$ is an approximation of $u_0$, which can take $P_hu_0$ or $R_hu_0$, but here it satisfies

Then we define the fractional Ritz-Volterra projection $V_hu\in S_r^h$ satisfying

Here we will use the fractional Ritz-Volterra projection $V_h$ and give some useful lemmas to study the error estimate.

Lemma 4.2. (^[9]) Denote by $\rho = V_h u-u$, then we have

The following energy estimation of $\rho_{tt}$ is novel and crucial to our convergence analysis.

Lemma 4.3. For the operator $\rho_{tt}$, we have

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,

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

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 $\omega$ satisfies $\frac{\partial^{2\alpha}}{\partial |x|^{2\alpha}}\omega = y$ in $H_0^{\alpha}(\Omega)$, and $\|\omega\|_{2\alpha}\leq C$. Choosing $v = R_h \omega$, it becomes

in which $\|R_h \omega- \omega\|_{\alpha}\leq h^{\alpha}\|\omega\|_{2\alpha}$ is used. By using the Grönwall's inequality again (Lemma 3.1), one has

After integration, we get

Therefore,

Thus the proof is completed.

Lemma 4.4. (^[5]) For each $\epsilon$ there is a constant $C_\epsilon = C_\epsilon(t)$ such that

Theorem 4.1. Assume that the initial values $u_{0h} = R_h u_0$ and $\|u_{1h}-u_1\|\leq Ch^r \|u_1\|_{r}$. Let $u(t)\in H_0^{\alpha}(\Omega)\cap H^{r}(\Omega)$ solve (3.1), and $u_h(t)\in S_r^h$ solve (4.1), then there exists a constant $C$ satisfies

Proof. By using $V_hu(t)\in S_r^h$ as an intermediate function, the error can be defined as $\varepsilon = u_h-u = (u_h-V_hu)+(V_hu-u) = \theta+\rho.$ Therefore, the error equation can be rewritten as

where $\theta(0) = 0$, $\theta_{t}(0) = (V_hu)_t(0)-u_{1h}$.

Taking $v = \theta_t$, it becomes

i.e.

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

where

Hence, we obtain

Then we get

in which the Grönwall's inequality is used. By using Lemma 4.3, we have

where

is used. This finishes the error estimates of $\theta$. Then combining Lemma 4.1 for error estimates of $\rho$, the proof of Theorem 4.1 is completed.

5.   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, \cdots, N$, where $k = T/N$ is the steplength in time variable $t$. Denote by $U^n$ the approximation solution and by ${\partial} U^n = (U^{n+1}-U^{n})/k, \overline{\partial} U^n = (U^{n}-U^{n-1})/k$ the forward and backward difference quotient of $U^n$ respectively, then

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_{0, t_{n}}^{1+\beta}g(t), 1\leq n\leq N,$ under the condition $g(t)\in C^2([0, T]).$ The following is the truncate error estimation.

Lemma 5.1. (^[30,31,32]) Suppose $u(t)\in C^2([0, T])$, then we have

where

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 $\omega_n\geq 0, f_n\geq 0$ and that for $n = 0, 1, \cdots, y_n\geq 0$ satisfies

then we have

Next we denote by $I_h u$ the piecewise polynomial interpolation operator of $u$ in $S_r^h$ satisfying

Lemma 5.3. (^[19]) Assume $u\in H_0^\alpha(\Omega)\cap H^{r}(\Omega)$, then we have that

Now we give the fully discrete scheme

with given initial values $U^0$ and $U^1\in S_r^h$, and $\sigma^n (g) = \frac{k^{1+\beta}}{\Gamma(3+\beta)}\sum_{j = 0}^{n}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

for $n\geq 1, t_{n+1}\leq T.$

Proof. Choosing $v = \overline{\partial} U^{n+1/2}$ in (5.1) to obtain

Note that

Multiplying both sides of (5.4) by $k$ and summing from $n = 1$ to $N$ obtain

For the term $k|\sum_{n = 1}^N I_1^n |$, we have

in which $b_{n-j} = (n-j+1)^{2+\beta}-2(n-j)^{2+\beta}+(n-j-1)^{2+\beta}$ is monotonically increasing for $n = 1, \cdots, N$. And for the term $k|\sum_{n = 1}^NI_2^n|$, 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) = \theta^n+\rho^n,$ we have

where $\theta^0 = u_0-U^0, \theta^1 = u_1-U^1,$ and

To construct the discrete initial values $U^0$ and $U^1$, let $u_2 = u_{tt}(0) = f(0)-\frac{\partial^{2\alpha}}{\partial |x|^{2\alpha}}u_0$ and define $U^1 = V_0+kV_1+V_2k^2/2,$ where $U^0 = V_0, V_j = P_hu_j, j = 0, 1, 2.$ According to Lemma 4.1, we get the following conclusions

Then for $\theta^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 = \overline{\partial}\theta^{n+1/2}$ in (5.6), we have

where

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

By the Taylor expansion, we know that

Therefore

and

Thus

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.

6.   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

with homogeneous boundary conditions. And the corresponding parameters satisfy $\Omega = [-1, 1], T = 1, \alpha = 0.9, \beta = 0.1$.

Case Ⅰ. We choose the source term $f = 2(1-x)(1+x)-(t^2+\frac{2t^{3+\beta}}{\Gamma(4+\beta)})\frac{\partial^{2\alpha}}{\partial |x|^{2\alpha}}(1-x)(1+x)$, 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 Ⅰ, which supports the predicted rates of the convergence. Figure 1 shows the exact solution and the numerical solution of case Ⅰ with $\alpha = 0.9, \beta = 0.1$, $h = 1/64, k = 0.01h$ at $t = 0.5, x\in[-1, 1]$. Figure 2 shows the error between the exact solution and the numerical solution of case Ⅰ with $\alpha = 0.9, \beta = 0.1$, $h = 1/64, k = 0.01h$, at $x = 0, t\in[0, 1]$. And Figure 3 shows the numerical solution of case Ⅰ under the same conditions.

Case Ⅱ. We choose the exact solution $u(x, t) = \sin(\pi 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 = \pi(1-x)(1+x).$

Table 2 shows the numerical results and convergence rates of case Ⅱ, which support the predicted rates of the convergence. Figure 4 shows the exact solution and the numerical solution of case Ⅱ with $\alpha = 0.9, \beta = 0.1$, $h = 1/64, k = 0.01h$, at $t = 0.5, x\in[-1, 1]$. Figures 5 and 6 show the error between the exact solution and the numerical solution, the numerical solution of case Ⅱ with $\alpha = 0.9, \beta = 0.1$, $h = 1/64, k = 0.01h$ for $x\in [-1, 1], t\in [0, 1]$ separately.

Example 6.2. In this example, we consider the following fractional hyperbolic equation

with homogeneous Dirichlet boundary conditions in $\Omega = [0, 1], T = 1$.

We choose the source term $f = 2(1-x)x^{\alpha}-(t^2+\frac{2t^{3+\beta}}{\Gamma(4+\beta)})D_{a, x}^{2\alpha} (1-x)x^{\alpha}$, then the exact solution is $u(x, t) = t^{2}(1-x)x^\alpha$, which has a weak singularity at the boundary point $x = 0$ if $0.5 < \alpha < 1$.

Table 3 shows the errors and convergence rates with parameters $\alpha = 0.6, \beta = 0.1$, $h = 1/64, k = 0.01h$ for $x\in[0, 1], t\in [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 $\alpha = 0.9, \beta = 0.1$, $h = 1/64$ for $x\in[0, 1], t\in [0, 1]$ separately. And Figure 9 shows the numerical solution with different values of $\alpha$ 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.

7.   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.

Acknowledgments

This work was sponsored by Natural Science Foundation of Shanghai under grant No. 19ZR1422000, Natural Science Foundation of Anhui Province under grant No. 2008085MA11, and Natural Science Foundation of Education Department of Anhui Province under grant No. KJ2018A0385.

The authors wish to thank the anonymous reviewers for their comments and proposals, which greatly improve the quality of the paper. The authors also wish to thank Prof. Changpin Li for his helpful discussions and suggestions.

Conflict of interest

All authors declare no conflicts of interest in this paper.