A multiple-step adaptive pseudospectral method for solving multi-order fractional differential equations

2.1 Definitions of fractional derivatives

In this subsection, we start by recalling the essential definitions of the fractional calculus. There are various definitions related to fractional integration and differentiation, such as Grünwald–Letnikov’s definition, Riemann–Liouville’s definition and Caputo definition [1]. Suppose α > 0, n = [α] + 1, where [α] denotes the integer part of α, and the function f(x) has n continuous bounded derivatives in [0, T], f ∈ Cⁿ[0, T], for every T > 0.

Recall that for α ∈ ℕ (ℕ = {1, 2, …}), the Caputo differential operator coincides with the usual differential operator of an integer order. Further, it is seen that fractional differential operators are defined through convolution integrals. Therefore, unlike integer order derivatives, they are non-local operators. They depend on all function values from its lower limit t = 0 up to the evaluation point t = x. This characteristic makes their evaluation more complicated.

Recalling the fact that the Grünwald–Letnikov definition coincides with the Riemann–Liouville derivative, it is pointed out that the Riemann–Liouville derivative has certain disadvantages when trying to model real-world phenomena with fractional differential equations [41]. Therefore, we shall use a modified fractional differential operator D^α proposed by Caputo in his work on the theory of viscoelasticity [42]. The reason for adopting the Caputo definition is as follows: in order to produce a unique solution for differential equations, it is required to specify additional conditions. For the case of differential equations with the Caputo derivatives, these additional conditions are just the traditional conditions, which are akin to those of classical differential equations. In contrast, for differential equations with the Riemann-Liouville derivatives, we are forced to specify some fractional derivatives of the unknown solution at the initial point x = 0, which are functions of x. In practical applications, these initial values are frequently not available, and it may not even be clear what their physical meaning is (see [43]). For more details about the mathematical properties of fractional derivatives, see [1, 2, 3, 4].

2.2 Basic properties of Legendre and shifted Legendre polynomials

The well-known Legendre polynomials, P_i(z), i = 0, 1, …, can be determined by the following recurrence relation:

where P₀(z) = 1 and P₁(z) = z. If P_i(z) is normalized so that P_i(1) = 1, then for any i, the Legendre polynomials in terms of power of z are

where i2 denotes the integer part of i2. Also, they are orthogonal with respect to L² inner product on the interval [−1, 1] with the weight function w(z) = 1, that is

where δ_ij is the Kronecker delta. The Legendre–Gauss (LG) collocation points −1 < z₁ < z₂ < … < z_N−1 < 1 are the roots of P_N−1(z). There are no explicit formulas for the LG points; however, they can be computed numerically using existing subroutines. The LG points have the property that

is exact for polynomials of degree at most 2N − 3, where

are the LG quadrature weights. For more details about Legendre polynomials, see [35].

The shifted Legendre polynomials on the interval x ∈ [a, b] are defined by

which are obtained by an affine transformation from the Legendre polynomials. The set of shifted Legendre polynomials is a complete L²(a, b)-orthogonal system with the weight function w(x) = 1. Thus, any function f ∈ L²(a, b) can be expanded in terms of shifted Legendre polynomials. The shifted Legendre–Gauss (ShLG) collocation points a < x₁ < x₂ < … < x_N−1 < b on the interval [a, b], are obtained by shifting the LG points, z_j, using the transformation

Thanks to the property of the standard LG quadrature, it follows that for any polynomial p of degree at most 2N − 3 on (a, b),

where w^j=b−a2wj, 1 ≤ j ≤ N − 1 are the ShLG quadrature weights.

4.1 Problem replacement

Let h > 0 be the step-size. Consider the subintervals I_k = [(k − 1)h, kh], k = 1, 2, … and let y_k(x) be the solution of the FIVP in the subinterval I_k. The FIVP (3.1)–(3.2) can be replaced with the following sequence of FIVPs on subintervals I_k:

where the notation y_(k) means that the function y is considered up to the k^th subinterval, i.e.,

Note that the reason for introducing this notation is the non-locality of the Caputo differential operator. It is important to mention that, according to (4.2), the initial conditions for step k (except for the first step where the initial conditions are available in (3.2)) are determined based on the solution obtained earlier at step k − 1. We choose to solve this sequence of problems via Legendre–Gauss adaptive pseudospectral method because of its high accuracy in large domain calculations.

4.2 Piecewise polynomial interpolation

Consider the ShLG collocation points (k − 1)h < x_k1 < ⋯ < x_k,N−1 < kh on the k^th subinterval I_k, k = 1, 2, …, obtained using the relation (2.5). Obviously,

Also, consider two additional noncollocated points x_k0 = (k − 1)h and x_kN = kh. Let us consider the space of polynomials of degree at most N on the subinterval I_k as

where

is a basis of Lagrange polynomials on subinterval I_k that satisfy L_ki(x_kj) = δ_ij where δ_ij is the Kronecker delta function. The L²(I_k)-orthogonal projection 𝓘_N : L²(I_k) → 𝓟_N(I_k) is a mapping in a way that for any y_k ∈ L²(I_k) we have

or equivalently

where y_ki = y_k(x_ki). Here, it can be easily seen that for i = 0, 1, …, N and k = 1, 2, …, we have

Thus, by utilizing Eq. (4.5) for Eq. (4.4), the approximation of y_k(x) within each subinterval I_k can be restated as

where Y_k and L_k(x) are (N + 1) × 1 matrices given by Y_k = [y_k0, …, y_kN]^T and L_k(x) = [L₁₀(x − x_k0), …, L_1N(x − x_k0)]^T. It is observed that the series (4.6) includes the Lagrange polynomials associated with the noncollocated points x_k0 = (k − 1)h and x_kN = kh, i.e., the initial and terminal points of each subinterval. This then allows more accurate initial conditions to be obtained for the sequence of problems (4.1)–(4.2). Moreover, it is seen from (4.6) that in the present multiple-step scheme, it is only needed to produce the basis of Lagrange polynomials L_1i(x) at the first step, which reduces the number of arithmetic calculations and storage requirements especially when the number of steps is large. Now consider the interval [0, T] such that for a positive integer K, ⋃k=1KIk = [0, T]. For N ⩾ 1 we introduce the piecewise polynomials space

which is the space of the continuous functions over [0, T] whose restrictions on each subinterval I_k are polynomials of degree ≤ N. Then, for any continuous function y in [0, T], the piecewise interpolation polynomial, ψ_N(y), coincides on each subinterval I_k with the interpolating polynomial 𝓘_N(y) of y_k = y|_Ik at the ShLG points.

Next, in a typical collocation method the derivatives of the unknown function y_k are approximated by the analytic derivatives of the interpolating polynomial [35, 36]. The polynomials of degree N − l, (𝓘_N(y_k))^(l), are called Legendre–Gauss interpolation derivatives of y_k [35]. By differentiating (4.6), l times, the ShLG interpolation derivatives of the function y on each subinterval I_k, k = 1, 2, …, are obtained as

where Lk(l)(x)=dldxlLk(x).

4.3 Error estimates

In this subsection, we give some estimates for the errors of piecewise Legendre–Gauss interpolation and its integer-order and fractional-order derivatives. We recall that denoted by y^(l) the distributional derivative of function y of order l, the Sobolev space of integer order r ≥ 0 in the interval (a, b) is defined by [35]

The associated norm is

According to the error bounds (5.4.33) and (5.4.34) of [35], the following error estimates hold.

Similar results hold for the piecewise interpolation error y − ψ_N(y) in the norms of the Sobolev spaces (see Theorem 1 of [40]). In the next theorem, we shall give an estimate in the maximum norm.

The following theorem gives an estimate for the piecewise interpolation fractional derivatives error, D^αy − D^αψ_N(y), in the maximum norm. In what follows, similarly to the definition of y_(k) in Subsection 4.1, the notation ψNk means that the piecewise interpolation is considered up to the k^th subinterval I_k.

4.4 Numerical evaluation of Caputo fractional derivatives

As in a typical pseudospectral method [35], it is crucial to evaluate the values of D^αy(x) at the ShLG collocation points x_kj, k = 1, 2, …, j = 1, …, N − 1, using an accurate and stable numerical method. In this subsection, we present a stable scheme for approximating the values of D^αy(x_kj). Our technique is based on the direct usage of the Caputo definition. From the definition of the Caputo fractional derivative (2.2) and using Eqs. (4.6) and (4.7), we have

To calculate the integrals on the right-hand side of Eq. (4.15), the Gaussian quadrature rule (2.6) is used in each subinterval. Therefore,

Moreover, in order to use the Gaussian integration formula in subinterval I_k for the last integral in Eq. (4.15), the change of variables τ=hx1jt−xk0+xk0 should be applied, which yields

Substituting (4.16) and (4.17) into (4.15), results

Unfortunately, the numerical evaluation (4.18) can be problematic in finite arithmetic. In fact, it is not stable when α changes between two consecutive positive integers. This problem arises from the fact that the coefficients x1j−x1jx1ihn−α−1 in (4.17) become large when α increases. We remedy this deficiency by utilizing integration by parts. By integrating by parts, Eq. (4.17) reads

Finally, substituting (4.16) and (4.19) into (4.15) we arrive at the following approximation:

The stability and more accuracy of the approximation (4.20) compared to that of the approximation (4.18) is illustrated in Fig. 1, where we plotted the maximum absolute errors E_α = max{Dαy(xkj)−dk,jα,, k = 1, 2; j = 1, …, N − 1} obtained by (4.18) and (4.20) for the function y(x) = x^5/2, x ∈ [0, 1] for 0 < α < 2 and different values of N. As we can see, the approximate values obtained by using (4.18) become quite inaccurate when α goes to 1 or 2. In contrast, (4.20) leads to much more accurate approximations of D^αy(x_kj). In addition, Fig. 1 illustrates the convergence of the approximation (4.20). Indeed, the convergence of the approximation (4.20) is a consequence of Theorems 1–3 and the high accuracy of the Gaussian quadrature rule provided that the function y is sufficiently smooth.

Fig. 1

Maximum absolute errors E_α for the function y(x) = x^5/2 by using (4.18) (left) and by using (4.20) (right).

4.5 The multiple-step collocation method

Consider the problem replacement as explained in Subsection 4.1. Collocating the k^th FIVP in subinterval I_k, k = 1, 2, … given by Eq. (4.1) at the ShLG collocation points x_kj for j = 1, 2, …, N − l + 1 using Eqs. (4.6) and (4.20), we obtain the collocation condition

Moreover, the initial conditions of the k^th FIVP in subinterval I_k given by Eq. (4.2), are approximated using Eqs. (4.6) and (4.7) as

We note that the vector Y_k−1 is obtained earlier at step k − 1. Also, the initial conditions of the first step are given by (3.2). At step k, k = 1, 2, …, Eq. (4.21) together with Eq. (4.22) give a system of equations with N + 1 set of algebraic equations, which can be solved with one of the known methods to find the unknowns y_kj of the vector Y_k. Consequently, the unknown function y_k(x) in Eq. (4.6) can be calculated at the k^th step, which is indeed the approximate solution of the FIVP (3.1)–(3.2) on the subinterval [(k-1)h, kh).