An Expansion Formula with Higher-Order Derivatives for Fractional Operators of Variable Order

Abstract

We obtain approximation formulas for fractional integrals and derivatives of Riemann-Liouville and Marchaud types with a variable fractional order. The approximations involve integer-order derivatives only. An estimation for the error is given. The efficiency of the approximation method is illustrated with examples. As applications, we show how the obtained results are useful to solve differential equations, and problems of the calculus of variations that depend on fractional derivatives of Marchaud type.

1. Introduction

Fractional calculus is a natural extension of the integer-order calculus by considering derivatives and integrals of arbitrary real or complex order 𝛼∈𝕂, with 𝕂=ℝ or 𝕂=ℂ. The subject was born from a famous correspondence between L’Hopital and Leibniz in 1695 and then developed by many famous mathematicians, like Euler, Laplace, Abel, Liouville, and Riemann, just to mention a few names. Recently, fractional calculus has attracted the attention of a vast number of researchers, not only in mathematics, but also in physics and in engineering, and has proven to better describe certain complex phenomena in nature [1, 2].

Since the order 𝛼 of the integrals and derivatives may take any value, another interesting extension is to consider the order not as a constant during the process but as a variable 𝛼(𝑡) that depends on time. This provides an extension of the classical fractional calculus and it was introduced by Samko and Ross in 1993 [3] (see also [4]). The variable order fractional calculus is nowadays recognized as a useful tool, with successful applications in mechanics, in the modeling of linear and nonlinear viscoelasticity oscillators, and in other phenomena where the order of the derivative varies with time. For more on the subject, and its applications, we mention [5–11]. For a numerical approach see, for example, [12–14]. Results on differential equations and the calculus of variations with fractional operators of variable order can be found in [15, 16] and references therein. In this paper we show how fractional derivatives and integrals of variable order can be approximated by classical integer-order operators.

The outline of the paper is the following. In Section 2 we present the necessary definitions, namely, the fractional operators of Riemann-Liouville and Marchaud of variable order. Some properties of the operators are also given. The main core of the paper is Section 3, where we prove the expansion formulas for the considered fractional operators, with the size of the expansion being the derivative of order 𝑛∈ℕ. In Section 4 we show the accuracy of our method with some examples and how the approximations can be applied in different situations to solve problems involving variable order fractional operators.

2. Fractional Calculus of Variable Order

In the following, the order of the fractional operators is given by a function 𝛼∈𝐶¹([𝑎,𝑏],]0,1[); 𝑥(⋅) is assumed to ensure convergence for each of the involved integrals. For a complete and rigorous study of fractional calculus we refer to [17].

Definition 1. Let 𝑥(⋅) be a function with domain [𝑎,𝑏]. Then, for 𝑡∈[𝑎,𝑏], (i)the left Riemann-Liouville fractional integral of order 𝛼(⋅) is given by _𝑎𝐼_𝑡^𝛼(𝑡)𝑥(𝑡)=1Γ(𝛼(𝑡))^𝑡_𝑎(𝑡−𝜏)^{𝛼(𝑡)−1}𝑥(𝜏)𝑑𝜏, (1)(ii)the right Riemann-Liouville fractional integral of order 𝛼(⋅) is given by _𝑡𝐼_𝑏^𝛼(𝑡)𝑥(𝑡)=1Γ(𝛼(𝑡))^𝑏_𝑡(𝜏−𝑡)^{𝛼(𝑡)−1}𝑥(𝜏)𝑑𝜏, (2)(iii)the left Riemann-Liouville fractional derivative of order 𝛼(⋅) is given by _𝑎𝐷_𝑡^𝛼(𝑡)𝑥(𝑡)=1Γ((𝑡𝑑1−𝛼))𝑑𝑡^𝑡_𝑎(𝑡−𝜏)^{−𝛼(𝑡)}𝑥(𝜏)𝑑𝜏, (3)(iv)the right Riemann-Liouville fractional derivative of order 𝛼(⋅) is given by _𝑡𝐷_𝑏^𝛼(𝑡)𝑥(𝑡)=−1Γ((𝑡𝑑1−𝛼))𝑑𝑡^𝑏_𝑡(𝜏−𝑡)^{−𝛼(𝑡)}𝑥(𝜏)𝑑𝜏, (4)(v)the left Marchaud fractional derivative of order 𝛼(⋅) is given by _𝑎𝔻_𝑡^𝛼(𝑡)𝑥(𝑡)=𝑥(𝑡)Γ(1−𝛼(𝑡))(𝑡−𝑎)^𝛼(𝑡)+𝛼(𝑡)Γ(1−𝛼(𝑡))^𝑡_𝑎𝑥(𝑡)−𝑥(𝜏)(𝑡−𝜏)^1+𝛼(𝑡)𝑑𝜏, (5)(vi)the right Marchaud fractional derivative of order 𝛼(⋅) is given by _𝑡𝔻_𝑏^𝛼(𝑡)𝑥(𝑡)=𝑥(𝑡)Γ(1−𝛼(𝑡))(𝑏−𝑡)^𝛼(𝑡)+𝛼(𝑡)Γ(1−𝛼(𝑡))^𝑏_𝑡𝑥(𝑡)−𝑥(𝜏)(𝜏−𝑡)^1+𝛼(𝑡)𝑑𝜏. (6)

Remark 2. It follows from Definition 1 that _𝑎𝐷_𝑡^𝛼(𝑡)𝑥(𝑡)=𝑑𝑑𝑡^𝑎𝐼_𝑡^{1−𝛼(𝑡)}𝑥(𝑡),_𝑡𝐷_𝑏^𝛼(𝑡)𝑥(𝑡)𝑑=−𝑑𝑡^𝑡𝐼_𝑏^{1−𝛼(𝑡)}𝑥(𝑡). (7)

Example 3 (see [3]). Let 𝑥 be the power function 𝑥(𝑡)=(𝑡−𝑎)^𝛾. Then, for 𝛾>−1, we have _𝑎𝐼_𝑡^𝛼(𝑡)𝑥(𝑡)=Γ𝛾+1Γ(𝑡)()𝛾+𝛼+1𝑡−𝑎^{𝛾+𝛼(𝑡)},_𝑎𝔻_𝑡^𝛼(𝑡)𝑥(𝑡)=Γ𝛾+1Γ(𝑡)()𝛾−𝛼+1𝑡−𝑎^{𝛾−𝛼(𝑡)},_𝑎𝐷_𝑡^𝛼(𝑡)𝑥(𝑡)=Γ𝛾+1Γ(𝑡)()𝛾−𝛼+1𝑡−𝑎^{𝛾−𝛼(𝑡)}−𝛼⁽¹⁾(𝑡)Γ𝛾+1Γ(𝑡)()𝛾−𝛼+2𝑡−𝑎^{𝛾−𝛼(𝑡)+1}×()(𝑡)((𝑡,ln𝑡−𝑎−𝜓𝛾−𝛼+2+𝜓1−𝛼)) (8) where 𝜓 is the Psi function, that is, the derivative of the logarithm of the Gamma function: 𝜓(𝑡)=𝑑(Γ(𝑡=Γ𝑑𝑡ln))^′(𝑡)Γ(𝑡). (9)

From Example 3 we see that _𝑎𝐷_𝑡^𝛼(𝑡)𝑥(𝑡)≠_𝑎𝔻_𝑡^𝛼(𝑡)𝑥(𝑡). Also, the symmetry on power functions is violated, when we consider _𝑎𝐼_𝑡^𝛼(𝑡)𝑥(𝑡) and _𝑎𝐷_𝑡^𝛼(𝑡)𝑥(𝑡), but holds for _𝑎𝐼_𝑡^𝛼(𝑡)𝑥(𝑡) and _𝑎𝔻_𝑡^𝛼(𝑡)𝑥(𝑡). Later we explain this better, when we deduce the expansion formula for the Marchaud fractional derivative. In contrast with the constant fractional order case, the law of exponents fails for fractional integrals of variable order. However, a weak form holds (see [3]): if 𝛽(𝑡)≡𝛽, 𝑡∈]0,1[, then _𝑎𝐼_𝑡𝑎^𝛼(𝑡)𝐼^𝛽_𝑡𝑥(𝑡)=_𝑎𝐼_𝑡^{𝛼(𝑡)+𝛽}𝑥(𝑡).

3. Expansion Formulas with Higher-Order Derivatives

The main results of the paper provide approximations of the fractional derivatives of a given function 𝑥 by sums involving only integer derivatives of 𝑥. The approximations use the generalization of the binomial coefficient formula to real numbers: (𝑡)𝑘()−𝛼−1^𝑘=𝛼(𝑡)𝑘=+𝑘−1Γ(𝛼(𝑡)+𝑘)Γ(𝛼(𝑡))𝑘!. (10)

Theorem 4. Fix 𝑛∈ℕ and 𝑁≥𝑛+1, and let 𝑥(⋅)∈𝐶^𝑛+1([𝑎,𝑏],ℝ). Define the (left) moment of 𝑥 of order 𝑘 by 𝑉_𝑘(𝑡)=(𝑘−𝑛)^𝑡_𝑎(𝜏−𝑎)^{𝑘−𝑛−1}𝑥(𝜏)𝑑𝜏. (11) Then, _𝑎𝐷_𝑡^𝛼(𝑡)𝑥(𝑡)=𝑆₁(𝑡)−𝑆₂(𝑡)+𝐸_1,𝑁(𝑡)+𝐸_2,𝑁(𝑡) (12) with 𝑆₁(𝑡)=(𝑡−𝑎)^{−𝛼(𝑡)𝑛}_𝑘=0𝐴(𝛼(𝑡),𝑘)(𝑡−𝑎)^𝑘𝑥^(𝑘)(𝑡)+^𝑁_𝑘=𝑛+1𝐵(𝛼(𝑡),𝑘)(𝑡−𝑎)^𝑛−𝑘𝑉_𝑘(𝑡), (13) where 𝐴(𝛼(𝑡))=1,𝑘Γ((𝑡×𝑘+1−𝛼))1+^𝑁_{𝑝=𝑛+1−𝑘}Γ𝑝−𝑛+𝛼(𝑡)Γ(𝛼(𝑡))!,𝐵(𝛼(𝑡))=Γ((𝑡−𝑘𝑝−𝑛+𝑘𝑘=0,…,𝑛,,𝑘𝑘−𝑛+𝛼))Γ((𝑡Γ((𝑡)!,𝑆−𝛼))1+𝛼))(𝑘−𝑛₂(𝑡)=𝑥(𝑡)𝛼⁽¹⁾(𝑡)Γ((𝑡()1−𝛼))𝑡−𝑎^{1−𝛼(𝑡)}×ln(𝑡−𝑎)(𝑡)−11−𝛼(1−𝛼(𝑡))²−ln(𝑡−𝑎)^𝑁_𝑘=0(𝑡)𝑘−𝛼(−1)^𝑘+𝑘+1_𝑁𝑘=0(𝑡)𝑘()−𝛼−1^𝑘𝑁_𝑝=11𝑝+𝛼𝑘+𝑝+1₍₁₎(𝑡)Γ((𝑡()1−𝛼))𝑡−𝑎^{1−𝛼(𝑡)}×ln(𝑡−𝑎)^𝑁+𝑛+1_𝑘=𝑛+1(𝑡)×()−𝛼𝑘−𝑛−1−1^{𝑘−𝑛−1}()𝑘−𝑛𝑡−𝑎^𝑛−𝑘𝑉_𝑘(𝑡)−^𝑁+𝑛+1_𝑘=𝑛+1(𝑡)()−𝛼𝑘−𝑛−1−1^{𝑘−𝑛−1}×^𝑁_𝑝=11𝑝()𝑘+𝑝−𝑛𝑡−𝑎^{𝑛−𝑘−𝑝}𝑉_𝑘+𝑝(𝑡). (14) The error of the approximation _𝑎𝐷_𝑡^𝛼(𝑡)𝑥(𝑡)≈𝑆₁(𝑡)−𝑆₂(𝑡) is given by 𝐸_1,𝑁(𝑡)+𝐸_2,𝑁(𝑡), where 𝐸_1,𝑁(𝑡) and 𝐸_2,𝑁(𝑡) are bounded by 𝐸_1,𝑁(𝑡)≤𝐿_𝑛+1(𝑡)×exp(𝑛−𝛼(𝑡))²+𝑛−𝛼(𝑡)Γ(𝑛+1−𝛼(𝑡))(𝑛−𝛼(𝑡))𝑁^{𝑛−𝛼(𝑡)}×(𝑡−𝑎)^{𝑛+1−𝛼(𝑡)},𝐸_2,𝑁(𝑡≤𝐿)¹(𝑡)𝛼⁽¹⁾(𝑡)(𝑡−𝑎)^{2−𝛼(𝑡)}exp𝛼²(𝑡)−𝛼(𝑡)Γ(2−𝛼(𝑡))𝑁^{1−𝛼(𝑡)}×(+1|ln𝑡−𝑎)|𝑁 (15)(16) with 𝐿_𝑗(𝑡)=max_{[]𝜏∈𝑎,𝑡}𝑥^(𝑗)(𝜏),𝑗∈{1,𝑛+1}. (17)

Proof. Starting with equality _𝑎𝐷_𝑡^𝛼(𝑡)𝑥(𝑡)=1Γ((𝑡𝑑1−𝛼))𝑑𝑡^𝑡_𝑎(𝑡−𝜏)^{−𝛼(𝑡)}𝑥(𝜏)𝑑𝜏, (18) doing the change of variable 𝑡−𝜏=𝑢−𝑎 over the integral, and then differentiating it, we get _𝑎𝐷_𝑡^𝛼(𝑡)𝑥(𝑡)=1Γ((𝑡𝑑1−𝛼))𝑑𝑡^𝑡_𝑎(𝑢−𝑎)^{−𝛼(𝑡)}𝑥()=1𝑡−𝑢+𝑎𝑑𝑢Γ((𝑡𝑥(𝑎)1−𝛼))(𝑡−𝑎)^𝛼(𝑡)+^𝑡_𝑎𝑑)𝑑𝑡(𝑢−𝑎^{−𝛼(𝑡)}(=1×𝑥𝑡−𝑢+𝑎)𝑑𝑢Γ((𝑡×𝑥(𝑎)1−𝛼))(𝑡−𝑎)^𝛼(𝑡)+^𝑡_𝑎−𝛼⁽¹⁾(𝑡)(𝑢−𝑎)^{−𝛼(𝑡)}()𝑥()+()×ln𝑢−𝑎𝑡−𝑢+𝑎𝑢−𝑎^{−𝛼(𝑡)}×𝑥⁽¹⁾(𝑡−𝑢+𝑎)𝑑𝑢=𝑆₁(𝑡)−𝑆₂(𝑡) (19) with 𝑆₁(𝑡)=1Γ((𝑡𝑥(𝑎)1−𝛼))(𝑡−𝑎)^𝛼(𝑡)+^𝑡_𝑎(𝑡−𝜏)^{−𝛼(𝑡)}𝑥⁽¹⁾(𝜏),𝑆𝑑𝜏₂(𝑡)=𝛼⁽¹⁾(𝑡)Γ(1−𝛼(𝑡))^𝑡_𝑎(𝑡−𝜏)^{−𝛼(𝑡)}ln(𝑡−𝜏)𝑥(𝜏)𝑑𝜏. (20)(21)
The equivalence between (20) and (13) follows from the computations of [18]. To show the equivalence between (21) and (14) we start in the same way as done in [19], to get 𝑆₂(𝑡)=𝛼⁽¹⁾(𝑡)Γ((𝑡𝑥(𝑡)1−𝛼))^𝑡_𝑎(𝑡−𝑢)^{−𝛼(𝑡)}ln(𝑡−𝑢)𝑑𝑢−^𝑡_𝑎𝑥⁽¹⁾(𝜏)×^𝜏_𝑎(𝑡−𝑢)^{−𝛼(𝑡)}()=𝛼ln𝑡−𝑢𝑑𝑢𝑑𝜏₍₁₎(𝑡)Γ((𝑡𝑥(𝑡)1−𝛼)))(𝑡−𝑎^{1−𝛼(𝑡)}×ln(𝑡−𝑎)(𝑡)−11−𝛼(1−𝛼(𝑡))²−^𝑡_𝑎𝑥⁽¹⁾(𝜏)×^𝜏_𝑎(𝑡−𝑎)^{−𝛼(𝑡)}1−𝑢−𝑎𝑡−𝑎^{−𝛼(𝑡)}×ln(𝑡−𝑎)+ln1−𝑢−𝑎.𝑡−𝑎𝑑𝑢𝑑𝜏 (22) Now, applying Taylor’s expansion over 1−𝑢−𝑎𝑡−𝑎^{−𝛼(𝑡)} (23) and ln1−𝑢−𝑎,𝑡−𝑎 (24) we deduce that 𝑆₂(𝑡)=𝛼⁽¹⁾(𝑡)Γ((𝑡×𝑥(𝑡)1−𝛼)))(𝑡−𝑎^{1−𝛼(𝑡)}ln(𝑡−𝑎)(𝑡)−11−𝛼(1−𝛼(𝑡))²−^𝑡_𝑎𝑥⁽¹⁾(𝜏)^𝜏_𝑎(𝑡−𝑎)^{−𝛼(𝑡)}()×ln𝑡−𝑎_𝑁𝑘=0(𝑡)𝑘()−𝛼−1^𝑘(𝑢−𝑎)^𝑘(𝑡−𝑎)^𝑘−𝑑𝑢_𝜏𝑎(𝑡−𝑎)^{𝑁−𝛼(𝑡)}_𝑘=0(𝑡)𝑘()−𝛼−1^𝑘×(𝑢−𝑎)^𝑘(𝑡−𝑎)^𝑘𝑁_𝑝=11𝑝(𝑢−𝑎)^𝑝(𝑡−𝑎)^𝑝𝑑𝑢𝑑𝜏+𝐸_2,𝑁(𝑡)=𝛼⁽¹⁾(𝑡)Γ((𝑡𝑥(𝑡)1−𝛼)))(𝑡−𝑎^{1−𝛼(𝑡)}×ln(𝑡−𝑎)(𝑡)−11−𝛼(1−𝛼(𝑡))²−^𝑡_𝑎𝑥⁽¹⁾(𝜏)(𝑡−𝑎)^{−𝛼(𝑡)}()×ln𝑡−𝑎_𝑁𝑘=0(𝑡)𝑘−𝛼(−1)^𝑘(𝑡−𝑎)^𝑘×^𝜏_𝑎(𝑢−𝑎)^𝑘+𝑑𝑢𝑑𝜏_𝑡𝑎𝑥⁽¹⁾(𝜏)(𝑡−𝑎)^{−𝛼(𝑡)}×^𝑁_𝑘=0(𝑡)𝑘−𝛼(−1)^𝑘(𝑡−𝑎)^𝑘×^𝑁_𝑝=11𝑝(𝑡−𝑎)^𝑝×^𝜏_𝑎(𝑢−𝑎)^𝑘+𝑝𝑑𝑢𝑑𝜏+𝐸_2,𝑁(𝑡)=𝛼⁽¹⁾(𝑡)(𝑡−𝑎)^{−𝛼(𝑡)}Γ((𝑡×𝑥(𝑡()1−𝛼)))(𝑡−𝑎)ln𝑡−𝑎(𝑡)−11−𝛼(1−𝛼(𝑡))²−ln(𝑡−𝑎)^𝑁_𝑘=0(𝑡)𝑘−𝛼(−1)^𝑘(𝑡−𝑎)^𝑘()×𝑘+1_𝑡𝑎𝑥⁽¹⁾(𝜏)(𝜏−𝑎)^𝑘+1+𝑑𝜏_𝑁𝑘=0(𝑡)𝑘−𝛼(−1)^𝑘(𝑡−𝑎)^𝑘𝑁_𝑝=11𝑝(𝑡−𝑎)^𝑝×𝑘+𝑝+1_𝑡𝑎𝑥⁽¹⁾(𝜏)(𝜏−𝑎)^𝑘+𝑝+1𝑑𝜏+𝐸_2,𝑁(𝑡). (25) Integrating by parts, we conclude with the two following equalities: ^𝑡_𝑎𝑥⁽¹⁾(𝜏)(𝜏−𝑎)^𝑘+1𝑑𝜏=𝑥(𝑡)(𝑡−𝑎)^𝑘+1−𝑉_𝑘+𝑛+1(𝑡),_𝑡𝑎𝑥⁽¹⁾(𝜏)(𝜏−𝑎)^𝑘+𝑝+1𝑑𝜏=𝑥(𝑡)(𝑡−𝑎)^𝑘+𝑝+1−𝑉_{𝑘+𝑝+𝑛+1}(𝑡). (26) The deduction of relation (14) for 𝑆₂(𝑡) follows now from direct calculations. Finally, we prove the upper bound formula for the error. The bound (15) for the error 𝐸_1,𝑁(𝑡) at time 𝑡 follows easily from [18]. With respect to sum 𝑆₂, the error at 𝑡 is bounded by𝐸_2,𝑁(𝑡≤𝛼)⁽¹⁾(𝑡)(𝑡−𝑎)^{−𝛼(𝑡)}Γ((𝑡×()1−𝛼))−ln𝑡−𝑎^∞_𝑘=𝑁+1(𝑡)𝑘−𝛼(−1)^𝑘×𝑘+1_𝑡𝑎𝑥⁽¹⁾(𝜏)(𝜏−𝑎)^𝑘+1(𝑡−𝑎)^𝑘+𝑑𝜏_{∞𝑘=𝑁+1}(𝑡)𝑘()−𝛼−1^𝑘×^∞_𝑝=𝑁+11𝑝×𝑘+𝑝+1_𝑡𝑎𝑥⁽¹⁾(𝜏)(𝜏−𝑎)^𝑘+𝑝+1(𝑡−𝑎)^𝑘+𝑝𝑑𝜏. (27) Define the quantities 𝐼₁(𝑡)=^∞_𝑁1𝑘^{1−𝛼(𝑡)}()𝐼𝑘+1)(𝑘+2𝑑𝑘,₂(𝑡)=^∞_𝑁^∞_𝑁1𝑘^{1−𝛼(𝑡)}𝑝𝑘+𝑝+1𝑘+𝑝+2𝑑𝑝𝑑𝑘. (28) Inequality (16) follows from relation (𝑡)𝑘≤−𝛼exp𝛼²(𝑡)−𝛼(𝑡)𝑘^{1−𝛼(𝑡)} (29) and the upper bounds 𝐼₁(𝑡)<^∞_𝑁1𝑘^{2−𝛼(𝑡)}1𝑑𝑘=(1−𝛼(𝑡))𝑁^{1−𝛼(𝑡)},𝐼₂(𝑡)<^∞_𝑁^∞_𝑁1𝑘^{2−𝛼(𝑡)}𝑝²1𝑑𝑝𝑑𝑘=(1−𝛼(𝑡))𝑁^{2−𝛼(𝑡)} (30) for 𝐼₁ and 𝐼₂.

Similarly as done in Theorem 4 for the left Riemann-Liouville fractional derivative, an approximation formula can be deduced for the right Riemann-Liouville fractional derivative.

Theorem 5. Fix 𝑛∈ℕ and 𝑁≥𝑛+1, and let 𝑥(⋅)∈𝐶^𝑛+1([𝑎,𝑏],ℝ). Define the (right) moment of 𝑥 of order 𝑘 by 𝑊_𝑘(𝑡)=(𝑘−𝑛)^𝑏_𝑡(𝑏−𝜏)^{𝑘−𝑛−1}𝑥(𝜏)𝑑𝜏. (31) Then, _𝑡𝐷_𝑏^𝛼(𝑡)𝑥(𝑡)=𝑆₁(𝑡)+𝑆₂(𝑡)+𝐸_1,𝑁(𝑡)+𝐸_2,𝑁(𝑡) (32) with 𝑆₁(𝑡)=(𝑏−𝑡)^{−𝛼(𝑡)𝑛}_𝑘=0𝐴(𝛼(𝑡),𝑘)(𝑏−𝑡)^𝑘𝑥^(𝑘)(𝑡)+^𝑁_𝑘=𝑛+1𝐵(𝛼(𝑡),𝑘)(𝑏−𝑡)^𝑛−𝑘𝑊_𝑘(𝑡), (33) where 𝐴(𝛼(𝑡))=(),𝑘−1^𝑘Γ((𝑡×𝑘+1−𝛼))1+^𝑁_{𝑝=𝑛+1−𝑘}Γ𝑝−𝑛+𝛼(𝑡)Γ(𝛼(𝑡))!,𝐵(𝛼(𝑡))=()−𝑘𝑝−𝑛+𝑘𝑘=0,…,𝑛,,𝑘−1^𝑛+1Γ(𝑘−𝑛+𝛼(𝑡))Γ((𝑡Γ((𝑡)!,𝑆−𝛼))1+𝛼))(𝑘−𝑛₂(𝑡)=𝑥(𝑡)𝛼⁽¹⁾(𝑡)Γ((𝑡()1−𝛼))𝑏−𝑡^{1−𝛼(𝑡)}×ln(𝑏−𝑡)(𝑡)−11−𝛼(1−𝛼(𝑡))²−ln(𝑏−𝑡)^𝑁_𝑘=0(𝑡)𝑘−𝛼(−1)^𝑘+𝑘+1_𝑁𝑘=0(𝑡)𝑘()−𝛼−1^𝑘𝑁_𝑝=11𝑝+𝛼𝑘+𝑝+1₍₁₎(𝑡)Γ((𝑡()1−𝛼))𝑏−𝑡^{1−𝛼(𝑡)}×ln(𝑏−𝑡)^𝑁+𝑛+1_𝑘=𝑛+1(𝑡)×()−𝛼𝑘−𝑛−1−1^{𝑘−𝑛−1}()𝑘−𝑛𝑏−𝑡^𝑛−𝑘𝑊_𝑘(𝑡)−^𝑁+𝑛+1_𝑘=𝑛+1(𝑡)()−𝛼𝑘−𝑛−1−1^{𝑘−𝑛−1}×^𝑁_𝑝=11𝑝()𝑘+𝑝−𝑛𝑏−𝑡^{𝑛−𝑘−𝑝}𝑊_𝑘+𝑝(𝑡). (34) The error of the approximation _𝑡𝐷_𝑏^𝛼(𝑡)𝑥(𝑡)≈𝑆₁(𝑡)+𝑆₂(𝑡) is given by 𝐸_1,𝑁(𝑡)+𝐸_2,𝑁(𝑡), where 𝐸_1,𝑁(𝑡) and 𝐸_2,𝑁(𝑡) are bounded by 𝐸_1,𝑁(𝑡)≤𝐿_𝑛+1(𝑡)×exp(𝑛−𝛼(𝑡))²+𝑛−𝛼(𝑡)Γ(𝑛+1−𝛼(𝑡))(𝑛−𝛼(𝑡))𝑁^{𝑛−𝛼(𝑡)}×(𝑏−𝑡)^{𝑛+1−𝛼(𝑡)},𝐸_2,𝑁(𝑡≤𝐿)¹(𝑡)𝛼⁽¹⁾(𝑡)(𝑏−𝑡)^{2−𝛼(𝑡)}exp𝛼²(𝑡)−𝛼(𝑡)Γ(2−𝛼(𝑡))𝑁^{1−𝛼(𝑡)}×(+1|ln𝑏−𝑡)|𝑁 (35) with 𝐿_𝑗(𝑡)=max_{[]𝜏∈𝑎,𝑡}𝑥^(𝑗)(𝜏),𝑗∈{1,𝑛+1}. (36)

Using the techniques presented in [20], similar formulas as the ones given by Theorems 4 and 5 can be proved for the left and right Riemann-Liouville fractional integrals of order 𝛼(⋅). For example, for the left fractional integral one has the following result.

Theorem 6. Fix 𝑛∈ℕ and 𝑁≥𝑛+1, and let 𝑥(⋅)∈𝐶^𝑛+1([𝑎,𝑏],ℝ). Then, _𝑎𝐼_𝑡^𝛼(𝑡)𝑥(𝑡)=(𝑡−𝑎)^{𝛼(𝑡)𝑛}_𝑘=0𝐴(𝛼(𝑡),𝑘)(𝑡−𝑎)^𝑘𝑥^(𝑘)(𝑡)+^𝑁_𝑘=𝑛+1𝐵(𝛼(𝑡),𝑘)(𝑡−𝑎)^𝑛−𝑘𝑉_𝑘(𝑡)+𝐸_𝑁(𝑡), (37) where 𝐴(𝛼(𝑡))=1,𝑘Γ((𝑡×𝑘+1+𝛼))1+^𝑁_{𝑝=𝑛+1−𝑘}Γ𝑝−𝑛−𝛼(𝑡)Γ((𝑡))!,𝐵(𝛼(𝑡))=Γ((𝑡−𝛼−𝑘𝑝−𝑛+𝑘𝑘=0,…,𝑛,,𝑘𝑘−𝑛−𝛼))Γ(𝛼(𝑡Γ((𝑡)!,𝑉))1−𝛼))(𝑘−𝑛_𝑘(𝑡)=(𝑘−𝑛)^𝑡_𝑎(𝜏−𝑎)^{𝑘−𝑛−1}𝑥(𝜏)𝑑𝜏,𝑘=𝑛+1,…,𝑁. (38) A bound for the error 𝐸_𝑁(𝑡) is given by 𝐸_𝑁(𝑡)≤𝐿_𝑛+1(𝑡)×exp(𝑛+𝛼(𝑡))²+𝑛+𝛼(𝑡)Γ(𝑛+1+𝛼(𝑡))(𝑛+𝛼(𝑡))𝑁^{𝑛+𝛼(𝑡)}×(𝑡−𝑎)^{𝑛+1+𝛼(𝑡)}. (39)

We now focus our attention on the left Marchaud fractional derivative _𝑎𝔻_𝑡^𝛼(𝑡)𝑥(𝑡). Splitting the integral (5), we deduce that _𝑎𝔻_𝑡^𝛼(𝑡)𝑥(𝑡)𝛼(𝑡)=−Γ(1−𝛼(𝑡))^𝑡_𝑎𝑥(𝜏)(𝑡−𝜏)^1+𝛼(𝑡)𝑑𝜏. (40)Integrating by parts,_𝑎𝔻_𝑡^𝛼(𝑡)𝑥(𝑡)=1Γ((𝑡𝑥(𝑎)1−𝛼))(𝑡−𝑎)^𝛼(𝑡)+^𝑡_𝑎(𝑡−𝜏)^{−𝛼(𝑡)}𝑥⁽¹⁾(𝜏)𝑑𝜏, (41) which is a representation for the left Riemann-Liouville fractional derivative when the order is constant, that is, when 𝛼(𝑡)≡𝛼 [17, Lemma 2.12]. For this reason, the Marchaud fractional derivative is more suitable as the inverse operation for the Riemann-Liouville fractional integral. With (41) and Theorem 4 in mind, it is not difficult to obtain the corresponding formula for _𝑎𝔻_𝑡^𝛼(𝑡)𝑥(𝑡).

Theorem 7. Fix 𝑛∈ℕ and 𝑁≥𝑛+1, and let 𝑥(⋅)∈𝐶^𝑛+1([𝑎,𝑏],ℝ). Then, _𝑎𝔻_𝑡^𝛼(𝑡)𝑥(𝑡)=𝑆₁(𝑡)+𝐸_1,𝑁(𝑡), (42) where 𝑆₁(𝑡) and 𝐸_1,𝑁(𝑡) are as in Theorem 4.

Similarly, having into consideration that _𝑡𝔻_𝑏^𝛼(𝑡)𝑥(𝑡)=1Γ((𝑡𝑥(𝑏)1−𝛼))(𝑏−𝑡)^𝛼(𝑡)−^𝑏_𝑡(𝜏−𝑡)^{−𝛼(𝑡)}𝑥⁽¹⁾(𝜏)𝑑𝜏, (43) the following result holds.

Theorem 8. Fix 𝑛∈ℕ and 𝑁≥𝑛+1, and let 𝑥(⋅)∈𝐶^𝑛+1([𝑎,𝑏],ℝ). Then, _𝑡𝔻_𝑏^𝛼(𝑡)𝑥(𝑡)=𝑆₁(𝑡)+𝐸_1,𝑁(𝑡), (44) where 𝑆₁(𝑡) and 𝐸_1,𝑁(𝑡) are as in Theorem 5.

4. Examples

For illustrative purposes, we consider the left Riemann-Liouville fractional integral and the left Riemann-Liouville and Marchaud fractional derivatives of order 𝛼(𝑡)=(𝑡+1)/4. Similar results as the ones presented here are easily obtained for the other fractional operators and for other functions 𝛼(⋅). All computations were done using the Computer Algebra System Maple.

4.1. Test Function

We test the accuracy of our approximations with an example.

Example 1. Let 𝑥 be the function 𝑥(𝑡)=𝑡⁴ with 𝑡∈[0,1]. Then, for 𝛼(𝑡)=(𝑡+1)/4, it follows from Example 3 that ₀𝐼_𝑡^𝛼(𝑡)𝑥(𝑡)=24Γ))𝑡((𝑡+21/4^(𝑡+17)/4₀,𝐷_𝑡^𝛼(𝑡)𝑥(𝑡)=24Γ))𝑡((19−𝑡/4^{(15−𝑡)/4}−6Γ))𝑡((23−𝑡/4^{(19−𝑡)/4}×ln(𝑡)−𝜓23−𝑡4+𝜓3−𝑡4₀,𝔻_𝑡^𝛼(𝑡)𝑥(𝑡)=24Γ))𝑡((19−𝑡/4^{(15−𝑡)/4}. (45)(46)(47) In Figures 1, 2, and 3 one can compare the exact expressions of the fractional operators of variable order (45), (46), and (47), respectively, with the approximations obtained from our results of Section 3 with 𝑛=2 and 𝑁∈{3,5}. The error 𝐸 is measured using the norm 𝐸𝑓,𝑔(𝑡)=¹₀𝑓(𝑡)−𝑔(𝑡)²𝑑𝑡. (48)

Figure 1 

Exact (45) and numerical approximations of the left Riemann-Liouville integral 0𝐼𝑡𝛼(𝑡)𝑥(𝑡) with 𝑥(𝑡)=𝑡4 and 𝛼(𝑡)=(𝑡+1)/4 obtained from Theorem 6 with 𝑛=2 and 𝑁∈{3,5}. The error (48) is 𝐸≈0.02169 for 𝑁=3 and 𝐸≈0.00292 for 𝑁=5.

Figure 2 

Exact (46) and numerical approximations of the left Riemann-Liouville derivative 0𝐷𝑡𝛼(𝑡)𝑥(𝑡) with 𝑥(𝑡)=𝑡4 and 𝛼(𝑡)=(𝑡+1)/4 obtained from Theorem 4 with 𝑛=2 and 𝑁∈{3,5}. The error (48) is 𝐸≈0.03294 for 𝑁=3 and 𝐸≈0.003976 for 𝑁=5.

Figure 3 

Exact (47) and numerical approximations of the left Marchaud derivative 0𝔻𝑡𝛼(𝑡)𝑥(𝑡) with 𝑥(𝑡)=𝑡4 and 𝛼(𝑡)=(𝑡+1)/4 obtained from Theorem 7 with 𝑛=2 and 𝑁∈{3,5}. The error (48) is 𝐸≈0.04919 for 𝑁=3 and 𝐸≈0.01477 for 𝑁=5.

4.2. Fractional Differential Equations of Variable Order

Consider the following fractional differential equation of variable order: ₀𝔻_𝑡^𝛼(𝑡)𝑥(𝑡)(𝑡)=1+𝑥Γ))𝑡((7−𝑡/4^(3−𝑡)/4𝑥(0)+𝑡,=0, (49) with 𝛼(𝑡)=(𝑡+1)/4. It is easy to check that 𝑥(𝑡)=𝑡 is a solution to (49). We exemplify how Theorem 7 may be applied in order to approximate the solution of such type of problems. The main idea is to replace all the fractional operators that appear in the differential equation by a finite sum up to order 𝑁, involving integer derivatives only, and, by doing so, to obtain a new system of standard ordinary differential equations that is an approximation of the initial fractional variable order problem. As the size of 𝑁 increases, the solution of the new system converges to the solution of the initial fractional system. The procedure for (49) is the following. First, we replace ₀𝔻_𝑡^𝛼(𝑡)𝑥(𝑡) by ₀𝔻_𝑡^𝛼(𝑡)𝑥(𝑡)≈𝐴(𝛼(𝑡),𝑁)𝑡^{−𝛼(𝑡)}𝑥(𝑡)+𝐵(𝛼(𝑡),𝑁)𝑡^{1−𝛼(𝑡)}𝑥⁽¹⁾(𝑡)+^𝑁_𝑘=2𝐶(𝛼(𝑡),𝑘)𝑡^{1−𝑘−𝛼(𝑡)}𝑉_𝑘(𝑡), (50) where 𝐴(𝛼(𝑡))=1,𝑁Γ(1−𝛼(𝑡))1+^𝑁_𝑝=2Γ𝑝−1+𝛼(𝑡)Γ(𝛼(𝑡!,𝐵(𝛼(𝑡))=1))𝑝−1,𝑁Γ(2−𝛼(𝑡))1+^𝑁_𝑝=1Γ𝑝−1+𝛼(𝑡)Γ(𝛼(𝑡)),𝐶(𝛼(𝑡))=Γ((𝑡−1𝑝!,𝑘𝑘−1+𝛼))Γ((𝑡Γ((𝑡)!,−𝛼))1+𝛼))(𝑘−1 (51) and 𝑉_𝑘(𝑡) is the solution of the system 𝑉_𝑘⁽¹⁾(𝑡)=(𝑘−1)𝑡^𝑘−2𝑥(𝑡),𝑉_𝑘(0)=0,𝑘=2,3,…,𝑁. (52) Thus, we get the approximated system of ordinary differential equations 𝐴(𝛼(𝑡),𝑁)𝑡^{−𝛼(𝑡)}+1𝑥(𝑡)+𝐵(𝛼(𝑡),𝑁)𝑡^{1−𝛼(𝑡)}𝑥⁽¹⁾(𝑡)+^𝑁_𝑘=2𝐶(𝛼(𝑡),𝑘)𝑡^{1−𝑘−𝛼(𝑡)}𝑉_𝑘(𝑡)=1Γ))𝑡((7−𝑡/4^(3−𝑡)/4𝑉+𝑡,_𝑘(1)(𝑡)=(𝑘−1)𝑡^𝑘−2𝑥(𝑡)𝑥(0)𝑉,𝑘=2,3,…,𝑁,=0,_𝑘(0)=0,𝑘=2,3,…,𝑁. (53) Now we apply any standard technique to solve the system of ordinary differential equations (53). We used the command dsolve of Maple. In Figure 4 we find the graph of the approximation 𝑥₃(𝑡) to the solution of problem (49), obtained by solving (53) with 𝑁=3. Table 1 gives some numerical values of such approximation, illustrating numerically the fact that the approximation 𝑥₃(𝑡) is already very close to the exact solution 𝑥(𝑡)=𝑡 of (49). In fact the plot of 𝑥₃(𝑡) in Figure 4 is visually indistinguishable from the plot of 𝑥(𝑡)=𝑡.

Table 1 

Some numerical values of the solution 𝑥₃(𝑡) of (53) with 𝑁=3, very close to the values of the solution 𝑥(𝑡)=𝑡 of the fractional differential equation of variable order (49).

Figure 4 

Approximation 𝑥3(𝑡) to the exact solution 𝑥(𝑡)=𝑡 of the fractional differential equation (49), obtained from the application of Theorem 7, that is, obtained by solving (53) with 𝑁=3.

4.3. Fractional Variational Calculus of Variable Order

We now exemplify how the expansions obtained in Section 3 are useful to approximate solutions of fractional problems of the calculus of variations [21]. The fractional variational calculus of variable order is a recent subject under strong current development [15, 16, 22, 23]. So far, only analytical methods to solve fractional problems of the calculus of variations of variable order have been developed in the literature, which consist in the solution of fractional Euler-Lagrange differential equations of variable order [15, 16, 22, 23]. In most cases, however, to solve analytically such fractional differential equations is extremely hard or even impossible, so numerical/approximating methods are needed. Our results provide two approaches to this issue. The first was already illustrated in Section 4.2 and consists in approximating the necessary optimality conditions proved in [15, 16, 22, 23], which are nothing else than fractional differential equations of variable order. The second approach is now considered. Similar to Section 4.2, the main idea here is to replace the fractional operators of variable order that appear in the formulation of the variational problem by the corresponding expansion of Section 3, which involves only integer-order derivatives. By doing it, we reduce the original problem to a classical optimal control problem, whose extremals are found by applying the celebrated Pontryagin maximum principle [24]. We illustrate this method with a concrete example. Consider the functional 𝐽(𝑥)=¹₀₀𝔻_𝑡^𝛼(𝑡)1𝑥(𝑡)−𝑡Γ((7−𝑡)/4)^{(3−𝑡)/42}𝑑𝑡, (54) with fractional order 𝛼(𝑡)=(𝑡+1)/4, subject to the boundary conditions 𝑥(0)=0,𝑥(1)=1. (55) Since 𝐽(𝑥)≥0 for any admissible function 𝑥 and taking 𝑥(𝑡)=𝑡, which satisfies the given boundary conditions (55), gives 𝐽(𝑥)=0, we conclude that 𝑥 gives the global minimum to the fractional problem of the calculus of variations that consists in minimizing functional (54) subject to the boundary conditions (55). The numerical procedure is now explained. Since we have two boundary conditions, we replace ₀𝔻_𝑡^𝛼(𝑡)𝑥(𝑡) by the expansion given in Theorem 7 with 𝑛=1 and a variable size 𝑁≥2. The approximation becomes ₀𝔻_𝑡^𝛼(𝑡)𝑥(𝑡)≈𝐴(𝛼(𝑡),𝑁)𝑡^{−𝛼(𝑡)}𝑥(𝑡)+𝐵(𝛼(𝑡),𝑁)𝑡^{1−𝛼(𝑡)}𝑥⁽¹⁾(𝑡)+^𝑁_𝑘=2𝐶(𝛼(𝑡),𝑘)𝑡^{1−𝑘−𝛼(𝑡)}𝑉_𝑘(𝑡). (56) Using (56), we approximate the initial problem (54)-(55) by the following one: to minimize 𝐽̃(𝑥)=¹₀𝐴(𝛼(𝑡),𝑁)𝑡^{−𝛼(𝑡)}𝑥(𝑡)+𝐵(𝛼(𝑡),𝑁)𝑡^{1−𝛼(𝑡)}𝑥⁽¹⁾(𝑡)+^𝑁_𝑘=2𝐶(𝛼(𝑡),𝑘)𝑡^{1−𝑘−𝛼(𝑡)}𝑉_𝑘(𝑡)−1𝑡Γ((7−𝑡)/4)^{(3−𝑡)/42}𝑑𝑡 (57) subject to 𝑉_𝑘⁽¹⁾(𝑡)=(𝑘−1)𝑡^𝑘−2𝑥(𝑡),𝑉_𝑘(0)𝑥(0)(1)=0,𝑘=2,…,𝑁,=0,𝑥=1, (58) where 𝛼(𝑡)=(𝑡+1)/4. This dynamic optimization problem has a system of ordinary differential equations as a constraint, so it is natural to solve it as an optimal control problem. For that, define the control 𝑢 by 𝑢(𝑡)=𝐴(𝛼(𝑡),𝑁)𝑡^{−𝛼(𝑡)}𝑥(𝑡)+𝐵(𝛼(𝑡),𝑁)𝑡^{1−𝛼(𝑡)}𝑥⁽¹⁾(𝑡)+^𝑁_𝑘=2𝐶(𝛼(𝑡),𝑘)𝑡^{1−𝑘−𝛼(𝑡)}𝑉_𝑘(𝑡). (59) We then obtain the control system 𝑥⁽¹⁾(𝑡)=𝐵⁻¹𝑡^{𝛼(𝑡)−1}𝑢(𝑡)−𝐴𝐵⁻¹𝑡⁻¹𝑥(𝑡)−^𝑁_𝑘=2𝐵⁻¹𝐶_𝑘𝑡^−𝑘𝑉_𝑘(𝑡)∶=𝑓(𝑡,𝑥(𝑡),𝑢(𝑡),𝑉(𝑡)), (60) where, for simplification, (𝛼(𝑡)),(𝛼(𝑡)),𝐶𝐴=𝐴,𝑁𝐵=𝐵,𝑁_𝑘(𝛼(𝑡)),𝑉(𝑡)=𝑉=𝐶,𝑘₂(𝑡),…,𝑉_𝑁(𝑡). (61) In conclusion, we wish to minimize the functional 𝐽̃(𝑥,𝑢,𝑉)=¹₀𝑢(𝑡)−1𝑡Γ((7−𝑡)/4)^{(3−𝑡)/42}𝑑𝑡 (62) subject to the first-order dynamic constraints 𝑥⁽¹⁾(𝑡)(),𝑉=𝑓𝑡,𝑥,𝑢,𝑉_𝑘(1)(𝑡)=(𝑘−1)𝑡^𝑘−2𝑥(𝑡),𝑘=2,…,𝑁, (63) and the boundary conditions 𝑥(0)𝑥(1)𝑉=0,=1,_𝑘(0)=0,𝑘=2,…,𝑁. (64) In this case, the Hamiltonian is given by 𝐻()=1𝑡,𝑥,𝑢,𝑉,𝜆𝑢−𝑡Γ((7−𝑡)/4)^{(3−𝑡)/42}+𝜆₁𝑓(𝑡,𝑥,𝑢,𝑉)+^𝑁_𝑘=2𝜆_𝑘(𝑘−1)𝑡^𝑘−2𝑥 (65) with the adjoint vector 𝜆=(𝜆₁,𝜆₂,…,𝜆_𝑁) [24]. Following the classical optimal control approach of Pontryagin et al. [24], we have the following necessary optimality conditions: 𝜕𝐻𝜕𝑢=0,𝑥⁽¹⁾=𝜕𝐻𝜕𝜆₁,𝑉_𝑝⁽¹⁾=𝜕𝐻𝜕𝜆_𝑝,𝜆₁⁽¹⁾=−𝜕𝐻,𝜆𝜕𝑥_𝑝(1)=−𝜕𝐻𝜕𝑉_𝑝. (66) That is, we need to solve the system of differential equations 𝑥⁽¹⁾(𝑡)=𝐵⁻¹Γ))−1((7−𝑡/42𝐵⁻²𝑡^{2𝛼(𝑡)−2}𝜆₁(𝑡)−𝐴𝐵⁻¹𝑡⁻¹𝑥(𝑡)−^𝑁_𝑘=2𝐵⁻¹𝐶_𝑘𝑡^−𝑘𝑉_𝑘(𝑡),𝑉_𝑘⁽¹⁾(𝑡)=(𝑘−1)𝑡^𝑘−2𝑥(𝑡)𝜆,𝑘=2,…,𝑁,₁₍₁₎(𝑡)=𝐴𝐵⁻¹𝑡⁻¹𝜆₁−^𝑁_𝑘=2(𝑘−1)𝑡^𝑘−2𝜆_𝑘(𝑡),𝜆_𝑘⁽¹⁾(𝑡)=𝐵⁻¹𝐶_𝑘𝑡^−𝑘𝜆₁,𝑘=2,…,𝑁, (67) subject to the boundary conditions 𝑥(0)𝑉=0,_𝑘(0)𝑥(1)𝜆=0,𝑘=2,…,𝑁,=1,_𝑘(1)=0,𝑘=2,…,𝑁. (68)

Figure 5 plots the numerical approximation 𝑥₂(𝑡) to the global minimizer 𝑥(𝑡)=𝑡 of the variable order fractional problem of the calculus of variations (54)-(55), obtained by solving (67)-(68) with 𝑁=2. The approximation 𝑥₂(𝑡) is already visually indistinguishable from the exact solution 𝑥(𝑡)=𝑡, and we do not increase the value of 𝑁. The effectiveness of our approach is also illustrated in Table 2, where some numerical values of the approximation 𝑥₂(𝑡) are given.

Table 2 

Some numerical values of the solution 𝑥₂(𝑡) of (67)-(68) with 𝑁=2, close to the values of the global minimizer 𝑥(𝑡)=𝑡 of the fractional variational problem of variable order (54)-(55).

Figure 5 

Approximation 𝑥2(𝑡) to the exact solution 𝑥(𝑡)=𝑡 of the fractional problem of the calculus of variations (54)-(55), obtained from the application of Theorem 7 and the classical Pontryagin maximum principle, that is, obtained by solving (67)-(68) with 𝑁=2.

Acknowledgments

This work is supported by FEDER Funds through COMPETE—Operational Programme Factors of Competitiveness (“Programa Operacional Factores de Competitividade”) and by Portuguese Funds through the Center for Research and Development in Mathematics and Applications (University of Aveiro) and the Portuguese Foundation for Science and Technology (“FCT—Fundação para a Ciência e a Tecnologia”), within Project PEst-C/MAT/UI4106/2011 with COMPETE no. FCOMP-01-0124-FEDER-022690. Delfim F.M. Torres was also supported by EU Funding under the 7th Framework Programme FP7-PEOPLE-2010-ITN, Grant Agreement no. 264735-SADCO.