A numerical study of fractional relaxation–oscillation equations involving Ã -Caputo fractional derivative

We provide a numerical method to solve a certain class of fractional differential equations involving ψ -Caputo fractional derivative. The considered class includes as particular case fractional relaxation–oscillation equations. Our approach is based on operational matrix of fractional integration of a new type of orthogonal polynomials. More precisely, we introduce ψ -shifted Legendre polynomial basis, and we derive an explicit formula for the ψ -fractional integral of ψ -shifted Legendre polynomials. Next, via an orthogonal projection on this polynomial basis, the problem is reduced to an algebraic equation that can be easily solved. The convergence of the method is justiﬁed rigorously and conﬁrmed by some numerical experiments.


Introduction
A relaxation oscillator is an oscillator based upon the performance of a physical system's resending to equilibrium after being disturbed. The relaxation-oscillation equation is the primary equation of relaxation and oscillation processes. The standard form of a relaxation equation is given by where λ ∈ R is a constant and f is a given function. Equation (1) models several physical phenomena, such as the Maxwell model, which describes the behavior of a viscoelastic material using a spring and a dashpot in series. In this case, λ = E η , where E is the elastic modulus, η is the viscosity coefficient, and f (t) denotes E multiplying the strain rate. The standard oscillation equation describes a simple physical process with a controlled phase shift. In the simplest linear case, the equation describes oscillation y of a system responding to an external forcing f : where λ = ω 2 0 with ω 0 is the natural (resonant) frequency of the oscillator. Note that in (2), it is supposed that the damping coefficient is zero.
In order to represent slow relaxation and damped oscillation, fractional derivatives are employed in the relaxation and oscillation models (1) and (2) (see [19,20]). Fractional relaxation-oscillation equation has the following form under the initial conditions y(0) = y 0 , if 0 < α < 1 ( 4 ) or where D α is a certain fractional derivative operator of order α. In the case 0 < α < 1, (3)-(4) describe the relaxation with the power law attenuation. In the case 1 < α < 2, (3)- (5) represent the damped oscillation with viscoelastic intrinsic damping of oscillator (see [5,26]). Recently, the numerical study of fractional relaxation-oscillation equations has attracted much attention. In [5], the authors studied the numerical solution of (3) (with f (t) = 0) by considering positive fractional derivative and fractal derivative. In [11], the authors used a Taylor matrix method in order to obtain the numerical solution of (3) by considering Caputo fractional derivative. This method is based on a fractional version of Taylor's formula established in [21]. In [12], the numerical solution of (3) in which the fractional derivative is given in the Caputo sense, is obtained by the optimal homotopy asymptotic method. In [7], the numerical solution of (3) with Caputo fractional derivative, is computed using a trapezoidal approximation of the fractional integral. In [25], a generalized wavelet collocation operational matrix method based on Haar wavelets is proposed to solve (3) in which the fractional derivative is given in the Caputo sense.
In this paper, motivated by the above cited works, we are concerned with the numerical solution of the fractional differential equation under the initial conditions where λ ∈ R is a constant, max m − 1, 1 2 < α < m, m ∈ N, ψ : I = [a, b] → [0, 1] is an increasing function that belongs to C m (I ) satisfying ψ (t) > 0, t ∈ I , and ψ(I ) = [0, 1], f : I → R is a given function, D α,ψ a is the ψ-Caputo fractional derivative of order α (see Definition 2.3), and for i = 0, 1, . . . , m − 1, Note that in the case ψ(t) = t, (6)-(7) includes as particular cases (3)-(4) for m = 1, and (3)-(5) for m = 2, where D α is the Caputo fractional derivative of order α. Our approach is based on operational matrix of fractional integration of a new type of orthogonal polynomials. More precisely, we introduce ψ-shifted Legendre polynomial basis, and we derive an explicit formula for the ψ-fractional integral of ψ-shifted Legendre polynomials. Next, by projecting the problem on this polynomial basis, we obtain an algebraic equation that can be easily solved. We present a rigorous proof of the convergence of the proposed method. Moreover, some numerical experiments are given to show the convergence of such procedure, comparing the exact solution with numerical approximation.
The operational matrix of integer integration has been used with different types of orthogonal polynomials such as Chebyshev polynomials [22], Legendre polynomials [23], Laguerre and Hermite polynomials [8], etc. Next, it was extended to the fractional case by many authors, see for example [3,4,9,10,13,15,24], and the references therein. Note that in all these cited works, only Caputo or Riemann-Liouville fractional derivatives were considered.
The paper is organized as follows. In Sect. 2, we present some preliminaries on fractional calculus and Hilbertian analysis, and we fix some notations. In Sect. 3, we introduce ψ-shifted Legendre polynomials, and study their properties. In Sect. 4, an explicit formula for the ψfractional integral of ψ-shifted Legendre polynomials is derived. The numerical scheme and the convergence of the method are discussed in Sect. 5. Finally, in Sect. 6, some numerical experiments are presented to demonstrate the efficiency of the proposed approach.
First, let us define some functional spaces that will be used later. Given a finite interval J ⊂ R, let with respect to the norm Let L 2 ψ (I ; R) be the functional space defined by It can be easily seen that L 2 ψ (I ; R) is a Hilbert space with respect to the scalar product We denote by · L 2 ψ (I ) the norm in L 2 ψ (I ; R) induced by the scalar product (·, ·) L 2 ψ (I ) , i.e., Further, given a function f : I → R, we define the function f : Using the change of variable s = ψ(t), we obtain Again, using the change of variable s = ψ(t), we obtain Definition 2.1 (see [16]) Let f ∈ L 1 ψ (I ; R). The ψ-fractional integral of order α > 0 of the function f is given by where is the Gamma function.
Note that in the particular case ψ(t) = t, I α,ψ a reduces to Riemann-Liouville fractional integral of order α. In the case ψ(t) = ln t (a > 0), I α,ψ a reduces to Hadamard fractional integral of order α (see [16] for more details).
is a linear and continuous operator. Moreover, we have .
Using Hölder's inequality, for all t ∈ I , we have which yields the desired result.
Definition 2.1 can be extended to a vector function as follows.

Definition 2.2
Let F : I → R N , N ∈ N, be a vector function given by The ψ-fractional integral of order α > 0 of the vector function F is given by The ψ-fractional integral operator has the following properties (see [16]): and For n ∈ N ∪ {0}, we define the differential operator (δ ψ ) n by Recently, Almeida et al. [1,2] introduced the concept of ψ-Caputo fractional derivative as follows.
For some special cases of ψ, we obtain from (9) different Caputo-type fractional derivatives. For example, if ψ(t) = t, then (9) reduces to the Caputo fractional derivative of order α (see [16]). In the case ψ(t) = ln t, (9) reduces to the Caputo-Hadamard fractional derivative (see [14]). When function f is of class C n , the ψ-Caputo fractional derivative of f can be represented by the expression (cf. [1, Theorem 3]) For example, (see [1]) Lemma 2.3 (see [2]) Let n − 1 < α < n, n ∈ N. Then: 2. If f ∈ C n−1 (I ; R), then Further, let H be a Hilbert space with respect to a certain scalar product (·, ·) H . We denote by · H the norm on H induced by (·, ·) H . Let (e i ) i∈N∪{0} be a Hilbertian basis of H . For K ∈ N, we denote by We recall the following standard result from functional analysis (see, for example [6]).

Lemma 2.4 For K ∈ N, K is a linear and continuous operator on H , satisfying
Proof For all K ∈ N, we have Passing to the limit as K → ∞, we obtain Then

Ã-Shifted Legendre polynomials
The analytic form of the normalized shifted Legendre polynomials L i (s) of degree i ∈ N∪{0}, defined on [0, 1], is given by (see, for example [18]): We introduce ψ-shifted Legendre polynomials as follows. The ψ-shifted Legendre polynomials P i,ψ (t) of degree i ∈ N ∪ {0}, defined on I , are given by Therefore, by Lemma 3.1, P i,ψ : i ∈ N ∪ {0} is an orthonormal basis of L 2 ψ (I ; R). Next, we have to prove that for any f ∈ L 2 ψ (I ; R), we have Let f ∈ L 2 ψ (I ; R). Then, by Lemma 2.1, we know that f ∈ L 2 ([0, 1]; R). Hence, by Lemma 3.1, we have Then, for all t ∈ I , we have Again, using Lemma 2.1, we obtain which proves (11).
Further, given a function f ∈ L 2 ψ (I ; R) and K ∈ N, let K ( f ) be the orthogonal projection of f on By Lemmas 2.4 and 3.2, we deduce immediately the following fact.

Operational matrices of integrations
Let F : I → R K , K ∈ N, be a vector function given by Suppose that F ∈ L 2 ψ (I ; R K ), i.e., F i ∈ L 2 ψ (I ; R), i = 0, 1, . . . , K − 1. We shall use the notation where K is the orthogonal projection operator defined by (12). We define the binary relation on L 2 ψ (I ; R K ) by We note that is not symmetric.
Let M α K ×K be the square matrix of size K , given by We have the following result.

Numerical scheme and convergence
Let us consider the fractional boundary value problem (6)- (7). We suppose that f ∈ L 2 ψ (I ; R) and (6)-(7) admits a unique solution y ∈ L 2 ψ (I ; R) (for existence and uniqueness of solution results, we suggest the reader the work [2]). Under the considered assumptions, from (6), we have D α,ψ a y ∈ L 2 ψ (I ; R). For K ∈ N (K is supposed to be large enough), we have where K is the orthogonal projection operator given by (12). Let H K ,α ∈ R K be the vector defined by Then, we have K D α,ψ a y = H T K ,α φ K ,ψ (t), t ∈ I , where φ K ,ψ is the vector function defined by (14). Using Lemmas 2.4 and 3.2, we deduce immediately the following convergence result.

Lemma 5.1 We have
On the other hand, since α > 1 2 , by Lemma 2.2, we know that is a linear and continuous operator. Therefore, by Lemma 5.1, we deduce the following result.

Lemma 5.2 We have
Next, we shall prove the following convergence result.

Lemma 5.3 We have
where K is the operator defined by (13).

Proof We have
On the other hand, observe that Hence, we obtain Next, using Lemma 2.4, we obtain Therefore, by Lemma 5.2, we deduce that Again, by Lemma 2.4, we have Finally, combining (18), (19) and (20), the desired result follows. Now, using Theorem 4.1 and Lemma 5.3, we deduce the following result.

Lemma 5.4 We have
Next, using the property (10) and the initial conditions (7), we obtain Let Z K ∈ R K be the (known) vector satisfying Let {y K } ⊂ L 2 ψ (I ; R) be the sequence defined by Theorem 5. 1 We have .
Using (21) and (22), we obtain Using Lemma 5.4 and passing to the limit as K → ∞, we obtain the desired result.
From Theorem 5.1, the solution y to (6)- (7) can be approximated by the sequence {y K } defined by (23). However, this sequence depends on the unknown vector H K ,α ∈ R K . Therefore, this vector should be computed before using the approximation given by Theorem 5.1.

Lemma 5.5
For all K ∈ N, we have Proof The result follows immediately by Lemma 2.5 and Theorem 5.1.
Let Q K ∈ R K be the (known) vector satisfying Remark 5.1 Note that it is supposed in (24) that the matrix A K is invertible. If it is not the case, then one may increase, iteratively, the number of the ψ-shifted Legendre coefficients by one, until A K becomes invertible.
Finally, after solving (24), the desired solution can be approximated via the sequence {y K } given by (23).

Numerical results
In this section, we illustrate the proposed method with some numerical experiments.
Let us consider the fractional oscillator equation under the initial conditions y(a) = y (a) = 0, where I = [a, b] and ψ : I → [0, 1] is an increasing function that belongs to C 2 (I ) such that ψ (t) > 0, t ∈ I and ψ(I ) = [0, 1]. We remark that problem (25)-(26) has a unique solution (cf [2]). We denote by y * the exact solution of (25)- (26). It can be easily seen that The numerical solution of (25)-(26) is denoted by y. We denote by E(t) the absolute error at the point t ∈ I , that is, Different choices of the function ψ are considered in this example.
The approximate solution of (31)-(28) and the absolute error at different points t ∈ [0, 1], in the case K = 6, are shown in Table 4.
Comparison of exact solutions and numerical solutions of (25)- (26) for all the considered cases are shown in Fig. 1, in the case K = 6.

Conclusion
A numerical approach based on operational matrix of fractional integration of a new type of orthogonal polynomials is introduced in this paper for solving a certain class of fractional differential equations involving ψ-Caputo fractional derivative. The convergence of the method is proved and the numerical experiments presented in Sect. 6 confirm the efficiency of this approach.