The spectral collocation method for solving a fractional integro-differential equation

1.   Introduction

All kinds of natural phenomena in nature can be explained by fractional calculus. For example, fractional integral and derivative theory are used to solve problems in nonlinear vibration of viscoelastic damping, control signal and processing, thermoelectric viscoelastic fluid, continuum and statistical mechanics, and other fields ^[1,2,3,4,5]. It is well known that generalized fractional operators are obtained by the extension of Caputo and Riemann-Liouville fractional derivatives. In convolution, because of the complexity of weight and kernel function, it is more difficult to design the higher-order numerical scheme of calculus equations with generalized fractional operators.

In the past decades, many numerical methods have been proposed for the differential and integral equations of different type problems. Such as, collocation method based on the Muntz- Legendre polynomials ^[6], Wavelet-Galerkin method and homotopy perturbation method ^[7], Legendre wavelets method ^[8], the rational second kind Chebyshev pseudospectral method ^[9], An hp-version spectral collocation method ^[10], and so on ^{[11,12,13,14]}. As one of the efficient numerical computing technique, the spectral collocation method provides a powerful tool for solving fractional calculus equations in recent years. Therefore, we also hope to provide a useful numerical method for fractional calculus equations.

In this paper, we considered the following fractional integro-differential equation:

with an initial condition

where $0 < \alpha < 1, \ 0 < \mu < 1, \ t\in I = [0, 1],$ $\ D_t^\alpha$ is the classical fractional derivative of order $\alpha$, $f(t)$ and $K(s, t)$ are known functions, and $y(t)$ is unknown function. This type of fractional integral-differential equation appear in several other places as in the modelling of particle motion in physics and mechanics in ^[15].

After careful planning, this article is developed as follows: the first section introduces the research background of fractional calculus equation model. Then, some results of fractional integro-differential and the properties of the polynomials in Section 2. In Section 3, based on the polynomials, a new method of the spectral collocation and Gauss-Legendre quadrature rule is proposed for solving fractional integro-differential problems. In Section 4, we analyze the error truncation of the numerical method. It considers some examples to illustrate the high accuracy of the proposed approach in Section 5. Finally, the conclusion is given.

2.   Properties of the shifted Legendre polynomials

2.1.   Fractional calculus

Several definitions of fractional calculus has been developed over time, it can be seen in ^[16,17].

Definition 1. ^[16] The Riemann-Liouville fractional integral operator of order $\alpha\geq 0$ is defined

Definition 2. ^[16] The fractional derivative of order $\alpha$ in the Caputo sense is defined

where $m-1 < \alpha\leq m,$ $m\in \mathbb{N},$ $x > 0$.

The properties of the operator $I^{\alpha}$ and $D^{\alpha}$ as follows

2.2.   The shifted Legendre polynomials

According to the definition in ^[16,18], we make it easier to solve Eq $(1.1)$, we first introduce its definition about the shifted Legendre polynomials of degree $i$ over the interval $[0, 1]$, as follows

and it is orthogonal for the legendre polynomial $L_i$ on $[-1, 1]$.

The shifted Legendre polynomial $\varphi_i(t)$ satisies the following form

and the complete basis is formed on $[0, 1]$. Then, the analytic form of $\varphi_i(t)$ may be expressed with $t$ as

where $b_{s, i} = \frac{(- 1)^{s + i}(i + s)!}{(i - s)!(s!)^2}$, and $\varphi _i(0) = (-1)^i$, $\varphi_i(1) = 1$.

It is orthogonal under the constraint $w = 1$ for shifted Legendre polynomials on $t \in [0, 1]$, that is $\int_0^1 \varphi_i(s)\varphi_j(s)ds = \frac{1}{i+j + 1}\delta_{i, j}$, here $\delta_{i, j}$ is Kronecker symbol.

2.3.   Function approximation

In reference [18], the following symbols are defined in $L^2(I)$ for weighted inner product and its norm as

Property 1. ^[18] Suppose $u \in L^2(I)$ is arbitrary, there exists a unique $q_m^ * \in {P_m},$ the following relation holds

where$m \in \mathbb{N},$ ${P_m}$ is the space of polynomials of order $m,$$q_m^ * (x) = \sum\limits_{k = 0}^m \hat u_k \phi _k(x), \ \hat u_k = \frac{(u, \phi _k)}{\left\| \phi_k \right\|_0^2},$and $\left\{ \phi _k \right\}_{k = 0}^m$\ forms an $L^2$-orthogonal basis for $P_m$.

By Property $1$, for any $y \in L^2(I)$, there is a unique best approximation $y_N$, such that

The relationship between the function $y(t) \in L^2(I)$ and the shifted Legendre polynomials as follows

Then, the shifted Legendre polynomials with first (N+1)-terms are considered as

Similarly, for any function $K(s, t) \in L^2(I)$ with the variable $s$, it can be expanded by $\varphi_i(s)$ as

which is used the Gauss-Legendre quadrature rule to compute the coefficients in this paper.

2.4.   Fractional operational matrices of derivative and integral

An approximated of the derivative of the shifted Legendre polynomials is given as

where $\boldsymbol{D}^\gamma$ is the derivative operation matrix of the shifted Legendre polynomials and $\Phi(s) = (\varphi_0(s), \ldots, \varphi_N(s))^T$.

Lemma 1. ^[19]Suppose ${{\boldsymbol{D}}}^\gamma$ is an operation matrix of fractional derivatives of order $\gamma > 0$ for the generalized shifted Legendre functions of $m \times m$ in the Caputo sense, then the $d_{ij}^{(\gamma)}$ are given for ${\boldsymbol{D}}^\gamma$ as

where

and $b_{s, i} = \frac{(- 1)^{i + s}(i + s)!}{(i - s)!(s!)^2}.$

Then$,$ we integrate for $\Phi(s)$ as follows

so called ${\boldsymbol{P}}^\upsilon$ is integral operator matrice of the shifted Legendre functions.

Lemma 2. ^[16]Suppose ${{{\boldsymbol{P}}}^\upsilon }$ is the $m \times m$ the generalized shifted Legendre functions operation matrix of Riemann-Liouville fractional integral of order $\upsilon$ $,$ then ${P^{(\upsilon)}_{i, j}}$ can be obtained for ${{{\boldsymbol{P}}}^\upsilon }$ as

3.   Numerical algorithm

In this part, based on the shifted Legendre polynomials, we consider the idea of combination of Gauss-Legendre quadrature rule and spectral collocation method to solve the Eq $(1.1)$. Then, using this numerical algorithm, we first expand $y(t)$ and $K(s, t)$ by the shifted Legendre polynomials of first (N+1)-terms as

where the coefficients are obtained by $(2.1)$ and $(2.2)$, respectively.

Substituting $(3.1)$ and $(3.2)$ in Eq $(1.1)$, we obtain

based on the orthogonality of polynomials on certain interval, have

simplified as

where $H(t) = [h_0(t), \ldots, h_N(t)]^T, h_i = k_i(t)/(2i + 1), i = 0, 1, \ldots, N.$

To solving the coefficients $Z$, we substitute collocation points into Eq $(3.3)$ lead to

where the points of shifted Chebyshev-Gauss are $t_r = (1-\cos((2r + 1)\pi/(2N + 2)))/2, r = 0, \ \ldots, \ N - 1.$ Thus, a system is formed by the $N$ nonlinear equation. In addition, an algebraic equation with the initial condition is provided to obtain $Z$ with

The value of $\{z_i\}_{i = 0}^N$ can be obtained by solving Eqs $(3.4)$ and $(3.5)$, the approximate solution ${y_N}(t) = Z^T\Phi(t)$ can be obtained. That is, the partial terms of vector $Z$ in Eq (3.4) are deformed and $Z$ is proposed as

set

where $i = 0, \ 1, \ 2, \ \cdots, \ N, \ r = 0, \ 1, \ 2, \ \cdots, \ N - 1.$

Then, by Eqs $(3.5)$ and $(3.6)$ we have

let $\hat{A} = [A, \Phi(0)], \ $$\hat{F} = [F, y_0],$ we can be obtained $Z$ by following

4.   Error analysis

The errors of the proposed method are analyzed, which is the work of this section. Then we introduce the following approximation operator, we can see in ^[20]. In the time direction, we refer to $\Pi_N$ as the traditional orthogonal projection operator of $L^2(I)$ and $P_N(I)$ is a space of polynomials of degree up to $N$ with the time variable $t$. $\Pi_N:L^2(I)\to P_N(I)$, that is, for any $y \in L^2(I)$, $\Pi_N y \in P_N(I)$, it satisfies

In the following lemma, we give an error estimate for approximating the operator $\Pi_N$.

Lemma 3. ^[20] For $0 < \alpha < 1,$ $\gamma > 1,$ if $y \in H^\alpha(I) \cap H^\gamma(I)$, then we have

where $C$ is a constant that is different from $N$.

Lemma 4. ^[20]Set the function $y\left(t \right)$ be expanded as ${y_N}(t) = {Z^T}\Phi(t)$ by shifted Legendre polynomials in $[0, 1],$ where

and

then$,$ such that

where $C$ is a constant independent of $N$.

Similar to the Lemma 3.6 in ^[21], it is easy to proof the following Lemma.

Lemma 5. For $0 < \mu < 1$ and $\mu$ is arbitrary, if $0 < \sigma < min(1/2, \mu)$, then

where $S y(s) = \frac{1}{\Gamma(\mu)}{\int_0^t {(t - s)}^{\mu - 1}}y(s)ds$, $S y_N(s) = \frac{1}{\Gamma(\mu)}{\int_0^t {(t - s)}^{\mu - 1}}y_N(s)ds$ and $C$ is a constant independent of $N$.

Next, replace $y(t)$ and $K(s, t)$ with $y_N(t)$ and $K_N(s, t)$ in (1.1), rewriting the equation as follows

we define the following operators

In the following section, we estimate the error for $L(y(t)-y_N(t))$.

Theorem 1. Let $y$ and ${y_N}$ be the analytic solution and numerical solutions of (1.1) and (1.2). Also assume that $\| y \|_{\gamma } < \infty, \ \gamma > 1,$ $||||K(s, t)||_m||_0\le K, m\geq 0,$ then

where $K$ is a real numbers and $C$ is a constant independent of $N$.

Proof. By Eqs (1.1) and (4.1) and triangle inequality, we obtain

where $R_1, R_2, R_3$ are defined by the right term of $(4.3)$, respectively.

For $R_1$, it is easy to obtain based on the Lemma 3 as

For $R_2$, based on the Lemma 5 it is directly to obtain that

where $0 < \sigma < min(1/2, \mu)$.

Next, we will estimate the $R_3$,

When $N \to \infty$, we get ${y_N}(s) \to y(s)$, then we have

where $C$ is a constant that is different from $N$.

Using (4.4)–(4.7), we have

At this point, the above theorem has been proved.

5.   Numerical examples

In this section, it is verified about the validity and applicability of the method by the following examples.

Example 1. Consider the fractional integro-differential equation as follows

where

its exact solution is $y(t) = e^t- 1,$ where $E_{1, 1 - \alpha}$ and $E_{1, 1 + \mu }$ is the classical Mittag-Leffler function.

Example 2. The fractional integro-differential equation is considered as

where

its exact solution is $y(t) = e^t - t - 1.$

In Eqs $(5.1)$ and $(5.2)$ with the initial condition $y(0) = 0$. In the numerical implementation, we provide a value of $\mu$ and the nodes and weights of the Gauss-Legendre quadrature rule. With the increasing value of $N$, the error result of different $\alpha$ are shown in Figure 1. A logarithmic scale has been used for error-axis in figures. Clearly, the errors show an exponential decay, since in these semi-log representations one observes that the error variations are essentially linear versus the polynomial degrees for $\alpha$. The results show that the method with a high accuracy to solving fractional integro-differential equations, and the error of $\|y(t)-y_N(t)\|_0$ converge exponentially with the increase of $N$ and reach machine precision. Meanwhile, the $L_2$ error has the advantage that it can better reflect the actual situation error of $\|y(t)-y_N(t)\|_0$ in Tables 1 and 2.

In order to reflect the effectiveness of the proposed method, the results be compared with euler method in following example.

Example 3. Consider the FIDE

with the initial condition $y(0) = 0.$ Its exact solution is $y(t) = t^3.$

From Table 3, it is clear that the error of absolute with the present method by polynomials of degree $N$ approximation. Clearly, in Example $3$ the method gives more accurate results by using lower order polynomials and fewer points than the euler method, and our method is easier to implement and understand.

6.   Conclusions

Based on the shifted Legendre polynomials, the idea of combining Gauss-Legendre quadrature rule and spectral collocation method, a new spectral collocation method is proposed for solving a fractional integro-differential equation in this paper. The error of this method is analyzed by Theorem $1$. It is found that the error decreases exponentially from Figure 1, Tables 1 and 2. In Table 3, the method gives more accurate results than the euler method. The results show that this method is high accuracy and easy to be implemented. In the future, based on the idea of ^[22,23], we will use the shifted Legendre polynomials for stochastic fractional integro-differential problems and shape optimization.

Acknowledgments

This research was supported by National Natural Science Foundation of China (Grant No. 11961009, 11901135), Foundation of Guizhou Science and Technology Department (Grant No. [2020]1Y015), Foundation for graduate students of Guizhou Provincial Department of Education (Grant No. YJSCXJH[2020]136). The authors thank the anonymous referees for their valuable suggestions to improve the quality of this work significantly.

Conflict of interest

The authors declare that they have no competing interests.