- Email: [email protected]
Analysis and numerical approximation of tempered fractional calculus of variations problems
Accepted Manuscript Analysis and numerical approximation of tempered fractional calculus of variations problems Ricardo Almeida, M. Luísa Morgado
PII: DOI: Reference:
S0377-0427(19)30190-6 https://doi.org/10.1016/j.cam.2019.04.010 CAM 12230
To appear in:
Journal of Computational and Applied Mathematics
Received date : 27 June 2018 Revised date : 28 February 2019 Please cite this article as: R. Almeida and M.L. Morgado, Analysis and numerical approximation of tempered fractional calculus of variations problems, Journal of Computational and Applied Mathematics (2019), https://doi.org/10.1016/j.cam.2019.04.010 This is a PDF file of an unedited manuscript that has been accepted for publication. As a service to our customers we are providing this early version of the manuscript. The manuscript will undergo copyediting, typesetting, and review of the resulting proof before it is published in its final form. Please note that during the production process errors may be discovered which could affect the content, and all legal disclaimers that apply to the journal pertain.
Analysis and Numerical Approximation of Tempered Fractional Calculus of Variations problems Ricardo Almeida1 and M. Lu´ısa Morgado2 1
Center for Research and Development in Mathematics and Applications (CIDMA) Department of Mathematics, University of Aveiro, 3810–193 Aveiro, Portugal, [email protected] 2 Center for Computational and Stochastic Mathematics, Instituto Superior T´ecnico and Department of Mathematics, University of Tr´as-os-Montes e Alto Douro,UTAD 5000-801, Vila Real, Portugal [email protected] Abstract In this paper, we study variational problems involving the tempered Caputo fractional derivative. The cost functional involves the tempered Caputo fractional derivative. Several important optimization conditions are derived to find the optimal solution. Sufficient and necessary conditions are presented for different variational problems. For example, the cases of integral (isoperimetric problem) and holonomic constraints are considered, as well as problems with high order derivatives. A numerical scheme is proposed to determine approximations of the solution and it is illustrated through some examples.
Mathematics Subject Classification 2010: 26A33, 49K05, 49M05. Keywords: Tempered fractional derivative, Euler–Lagrange equation, Numerical methods.
1
Introduction
Since its very beginning that Calculus has been used to describe many physical phenomena. Although Fractional Calculus (FC), which consists of the study of integral and differential operators whose order is not restricted to integer values, is as much as old as the classical (integer order) one, only recently it has been intensively used in the modelling of many physical and engineering processes. FC may be regarded as a generalization of the integer order one. However, several (and nonequivalent) definitions of fractional order derivatives exist. We refer the reader to the books [1–4] to see the potential applications of FC and also to the variety of different definitions for fractional derivatives. We dare to say that the most popular fractional derivatives are the Riemann-Liouville, the Gr¨ unwald-Letnikov and the Caputo derivatives, being the last one the most commonly found in applicatons. There are two reasons for such choice. First, the Caputo fractional derivative of a constant is the null function. Secondly, the Laplace transform of the Caputo fractional derivative involves integer-order initial conditions. In this paper we are concerned with the theoretical analysis and numerical approximation of functionals involving tempered fractional derivatives in the Caputo sense. These new fractional operators are a generalization of the classical ones, in the sense that they depend on a new parameter λ, and when λ = 0, we obtain the well-know fractional integrals/derivatives. Left and right sided tempered fractional integrals and derivatives can be given trough the definition of the fractional integrals and derivatives (see [5, 6] for example).
1
Definition 1. For any α > 0 and fixed parameter λ ≥ 0, the left and right-sided tempered Riemann-Liouville fractional integrals of a function x on (a, b), are given by Z t 1 α,λ α λt −λt (t − s)α−1 e−λ(t−s) x(s) ds, a It e x(t) = a It x(t) = e Γ(α) a Z b 1 α,λ −λt λt α x(t) = (s − t)α−1 e−λ(s−t) x(s) ds, t Ib x(t) = e t Ib e Γ(α) t
respectively, where a Itα and t Ibα are the left and right Riemann-Liouville fractional integrals ( [7]). Note that when λ = 0, the tempered Riemann-Liouville fractional integrals reduce to the (usual) Riemann-Liouville integrals. Definition 2. For any α > 0, let n ∈ N be given by n = [α] + 1 if α ∈ / N, or n = α otherwise, and let λ ≥ 0 be a fixed parameter. The left and right-sided tempered fractional Caputo derivatives of a function x are given by: n Z t e−λt d C α,λ −λt C α λt n−α−1 (t − s) eλs x(s) ds, (1) a Dt x(t) = e a Dt e x(t) = Γ(n − α) a ds n Z b d eλt C α,λ n−α−1 λt C α λt D (s − t) − e−λs x(s) ds, x(t) = e D −e x(t) = t t b b Γ(n − α) t ds
α C α respectively, where C a Dt and t Db denote the left and right Caputo fractional derivatives, respectively ( [7]).
Once again if λ = 0 then the Caputo tempered fractional derivatives reduce to the Caputo fractional derivatives, and therefore, tempered fractional integrals and derivatives can be regarded as a generalization of the usual (non-tempered) fractional integrals and derivatives. Finally, we present the equivalent definitions for the Riemann–Liouville fractional derivatives. Definition 3. For any α > 0, let n = [α] + 1, and let λ ≥ 0 be a fixed parameter. The left and right-sided tempered fractional Riemann–Liouville derivatives of a function x are given by: n Z t d e−λt α,λ −λt α λt D x(t) = e D e x(t) = (t − s)n−α−1 eλs x(s)ds, a t a t Γ(n − α) dt a n Z b d eλt α,λ λt α λt D x(t) = e D −e x(t) = − (s − t)n−α−1 e−λs x(s)ds, t b t b Γ(n − α) dt t
respectively, where a Dtα and t Dbα denote the left and right Riemann–Liouville fractional derivatives, respectively ( [7]). Tempered fractional derivatives have appeared as one of the solutions to overcome the problem of the divergence of the second moment (or variance) of the probability density function of −1−α L´evy flights, which decays as |x| . This L´evy distribution originates the substancial (usual) fractional Riemann-Liouville derivatives of order α. By exponentially tempering the probability of −1−α −λ|x| large jumps of L´evy flights with a parameter λ > 0, the density function decays as |x| e and then finite second moments are obtained and the definitions of tempered fractional derivatives arise ( [8]). In this work we are concerned with Fractional Calculus of Variations problems. This subject has been extensively developed in the last decades, since the pioneering works of Riewe [9]. Since then, numerous works have appeared, dealing with different fractional operators. To mention a few of them, we refer to the case with the Riemann–Liouville derivative [10–18], the Caputo derivative [19–23], the Riesz derivative [24, 25], symmetric fractional derivative [26, 27], fuzzy fractional variational problems [28], optimal control problems [29], distributed-order derivatives [30], etc. We also mention the recent books [31, 32], where analytical and numerical methods were developed to deal with fractional variational problems. 2
Here we will be concerned with functionals of the type: → R Z b α,λ 7 → L(t, x(t), C a Dt x(t)) dt,
J : C 1 [a, b] x
(2)
a
where L : [a, b] × R2 → R is continuously differentiable with respect to the second and third variables. The fractional order α is a real between 0 and 1, and λ ≥ 0. The following notation will be used throughout the text. Given a function of several variables L ≡ L(a1 , . . . , am ), we use the symbol ∂L . ∂i L(a1 , . . . , am ) := ∂ai (a1 ,...,am )
α,λ Also, the vector (t, x(t), C a Dt x(t)) will be abbreviated by [x](t). A function L(a1 , a2 , . . . , am ) is convex with respect to the variables a2 , . . . , am if ∂i L exists for any i ∈ {2, . . . , m} and the following inequality holds:
L(a1 , a2 +a02 . . . , am +a0m )−L(a1 , a2 , . . . , am ) ≥ ∂2 L(a1 , a2 , . . . , am )a02 +. . .+∂m L(a1 , a2 , . . . , am )a0m , for all vectors (a1 , a2 + a02 . . . , am + a0m ) and (a1 , a2 , . . . , am ). The paper is organized as follows. In Section 2 we present our theoretical study. Several necessary and sufficient conditions are proven, for different variational problems. We consider the fundamental problem and deduce the respective Euler–Lagrange equation and the transversality conditions. Two types of constraints are after considered: the integral and holonomic problems. Then, a generalization of the fundamental problem is studied, where the Lagrange function depends on high order derivatives. Then, in Section 3, a numerical procedure is provided, following a usual direct approach as in [33], where variational problems with (non-tempered) Caputo derivatives have been considered. Here, we use different approximation formulas, namely, we use a higher order quadrature formula for the approximation of the integral in the functional. In this way, we not only extend the results in [33] (which here correspond to the case where the tempered parameter λ is zero), but we also improve the numerical scheme. Some examples for which the analytical solution is known are considered to illustrate the efficiency of the method.
2
Necessary and sufficient conditions of optimality
To start, we will present and prove a fractional integration by parts formula needed for the sequel. Theorem 1. Let x ∈ C[a, b] and y ∈ C n [a, b] be two functions. Then, for α > 0 and n ∈ N such that n = [α] + 1, the following formula holds: Z
a
b
α,λ x(t)C a Dt y(t) dt
=
Z
a
b
y(t)t Dα,λ b x(t) dt+
# "n−1 k b X dn−k−1 d λt −λt n−α,λ e y(t) − e t Ib x(t) . dtn−k−1 dt k=0
a
Proof. By definition of the tempered Caputo fractional derivative, Z
a
b
α,λ x(t)C a Dt y(t) dt =
Z
a
b
Z
a
t
x(t)e−λt dn (t − s)n−α−1 n eλs y(s) ds dt. Γ(n − α) ds
If we apply the Dirichlet’s formula, the last formula is equal to Z
a
b
dn λt 1 e y(t) dtn Γ(n − α)
Z
b
t
3
(s − t)n−α−1 e−λs x(s) ds dt.
Using now integration by parts, we obtain Z
b
a
"
#b Z b dn−1 λt 1 n−α−1 −λs = (s − t) e y(t) e x(s) ds dtn−1 Γ(n − α) t a ! Z b Z b n−1 d 1 d λt n−α−1 −λs e y(t) − (s − t) + e x(s) ds dt. n−1 dt Γ(n − α) t a dt
α,λ x(t)C a Dt y(t) dt
Integrating once again by parts, we get Z
b
a
# k 1 b X dn−k−1 λt d −λt n−α,λ x(t) = e t Ib e y(t) − dtn−k−1 dt k=0 a ! 2 Z b n−2 Z b d 1 d λt n−α−1 −λs + e y(t) − (s − t) e x(s) ds dt. n−2 dt Γ(n − α) t a dt
α,λ x(t)C a Dt y(t) dt
"
Repeating the procedure, we deduce the following: Z
a
b
# "n−1 k b X dn−k−1 d λt −λt n−α,λ = e y(t) − e t Ib x(t) dtn−k−1 dt k=0 a ! n Z b Z b 1 d + (s − t)n−α−1 e−λs x(s) ds dt eλt y(t) − dt Γ(n − α) t a "n−1 # k b Z b X dn−k−1 d λt −λt n−α,λ e y(t) − = e t Ib x(t) + y(t)t Dα,λ b x(t) dt, dtn−k−1 dt a
α,λ x(t)C a Dt y(t) dt
k=0
a
proving the desired formula. In particular, when the fractional order α is between 0 and 1, we obtain the simple formula: Z
a
b
α,λ x(t)C a Dt y(t) dt =
Z
a
b
h ib 1−α,λ y(t)t Dα,λ x(t) dt + y(t) I x(t) . t b b a
(3)
We are now in conditions to present a necessary condition that allows us to find the candidates of minimizers for functionals J as in (2). This type of conditions is known in the literature as Euler–Lagrange equations. Theorem 2. Let x? be a local minimizer of functional J as in (2), subject to the boundary conditions x(a) = β and x(b) = γ, where β, γ ∈ R. (4) If the map
? t 7→ t Dα,λ b (∂3 L[x ](t)),
t ∈ [a, b],
is continuous, then x? is a solution of the fractional differential equation ∂2 L[x](t) + t Dα,λ b (∂3 L[x](t)) = 0,
t ∈ [a, b].
(5)
Proof. Let η ∈ C 1 [a, b] be an arbitrary function satisfying the boundary conditions η(a) = 0 = η(b) and a real number. Extrema are characterized by setting the first variation of the functional equal to zero. Since x? + η is an admissible variation for the problem, J (x? + η) − J (x? ) = 0, →0
δJ (x? , η) := lim
4
we obtain
Z
b
a
α,λ ∂2 L[x? ](t)η(t) + ∂3 L[x? ](t)C a Dt η(t) dt = 0,
and using the integration by parts formula (3), Z
b
a
h
i h ib 1−α,λ ? ? ∂2 L[x? ](t) + t Dα,λ (∂ L[x ](t)) η(t) dt + η(t) I (∂ L[x ](t)) = 0. 3 t 3 b b a
Since η is an arbitrary function vanishing at t = a and t = b, we get that x? is a solution of Eq. (5). Definition 4. Function x? is said to be an extremal for functional J in (2) if ? ∂2 L[x? ](t) + t Dα,λ b (∂3 L[x ](t)) = 0,
t ∈ [a, b].
Theorem 3. Suppose that L is convex with respect to the variables x and solution x? of Eq. (5) is a solution of the problem J (x) =
Z
b
a
α,λ L(t, x(t), C a Dt x(t)) dt → min
s.t. x(a) = β
C α,λ a Dt x.
Then, every
and x(b) = γ.
Proof. Consider a variation of x? of form x? + η with η(a) = 0 = η(b). Using the convexity assumption and applying integration by parts formula (3), we get ?
?
J (x + η) − J (x ) ≥ =
Z
b
a
Z
a
b
α,λ ∂2 L[x? ](t)η(t) + ∂3 L[x? ](t)C a Dt η(t) dt
h
i h ib 1−α,λ ? ∂2 L[x? ](t) + t Dα,λ (∂3 L[x? ](t)) b (∂3 L[x ](t)) η(t) dt + η(t)t Ib
a
which is equal to zero since x? satisfies Eq. (5). This proves that J (x? ) is a local minimum of J. Theorem 4. Let x? be a local minimizer of functional J as in (2). If the map ? t 7→ t Dα,λ b (∂3 L[x ](t)),
t ∈ [a, b],
is continuous, then x? is a solution of the tempered fractional differential equation ∂2 L[x](t) + t Dα,λ b (∂3 L[x](t)) = 0,
t ∈ [a, b].
(6)
at t = a and t = b.
(7)
Moreover, the two next transversality conditions hold: 1−α,λ (∂3 L[x? ](t)) t Ib
=0
Proof. As it was seen in the proof of Theorem 2, Z
a
b
h
i h ib 1−α,λ ? ∂2 L[x? ](t) + t Dα,λ (∂3 L[x? ](t)) = 0. b (∂3 L[x ](t)) η(t) dt + η(t)t Ib a
If we assume that η(a) = 0 = η(b), then we obtain Eq. (6). Therefore, h ib η(t)t I1−α,λ (∂3 L[x? ](t)) = 0, b a
and by the arbitrariness of η at t = a and t = b Eqs. (7) are obtained.
5
Theorem 5. Suppose that L is convex with respect to the variables x and solution x? of Eqs. (6)-(7) is a solution of the problem J (x) =
Z
b
a
C α,λ a Dt x.
Then, every
α,λ L(t, x(t), C a Dt x(t)) dt → min .
Proof. Similar to the proof of Theorem 3. We proceed our work considering the isoperimetric problem. In this case, besides the boundary conditions (4), we have an integral constraint. In what follows, M : [a, b]×R2 → R is a continuously differentiable function with respect to the second and third variables and k ∈ R is a constant. Theorem 6. Suppose that x? is a solution of the problem: Z b α,λ J (x) = L(t, x(t), C a Dt x(t)) dt → min a
subject to the constrains x(a) = β
and x(b) = γ,
and I (x) := If the maps
Z
a
b
where β, γ ∈ R,
α,λ M (t, x(t), C a Dt x(t)) dt = k.
? t 7→ t Dα,λ b (∂3 L[x ](t)), ?
? t 7→ t Dα,λ b (∂3 M [x ](t)),
(8)
(9) t ∈ [a, b],
are continuous, and x is not an extremal for I, then there exists a real number µ such that x? is a solution of the equation ∂2 H[x](t) + t Dα,λ b (∂3 H[x](t)) = 0,
t ∈ [a, b],
(10)
where H is the function defined by H := L + µM . Proof. Consider a variation of x? with two parameters x? + 1 η1 + 2 η2 , where ηi ∈ C 1 [a, b] and ηi (a) = 0 = ηi (b), for i = 1, 2. Consider the two functions j, i (defined in a neighborhood of zero) as j(1 , 2 ) :=J (x? + 1 η1 + 2 η2 ) i(1 , 2 ) :=I(x? + 1 η1 + 2 η2 ) − k. Since ∂i (0, 0) = ∂2
Z
a
b
h
i ? ∂2 M [x? ](t) + t Dα,λ b (∂3 M [x ](t)) η2 (t) dt,
and x? is not an extremal for I, we conclude that there must exist a function η2 for which ∂i/∂2 (0, 0) 6= 0. Thus, by the Implicit Function Theorem, there exists some function 2 (·), defined in a neighborhood of zero, such that i(1 , 2 (1 )) = 0. From now on, we will consider such subfamily of variations that fulfills the isoperimetric constraint (9). Since ∇i(0, 0) 6= (0, 0), and (0, 0) is a solution of the problem: minimize j
subject to i(1 , 2 ) = 0,
we can ensure the existence of a real µ with ∇(j + µi)(0, 0) = (0, 0). Fom the relation ∂(j + µi) (0, 0) = ∂1
Z
a
b
h i ? ∂2 (L + µM )[x? ](t) + t Dα,λ (∂ (L + µM )[x ](t)) η1 (t) dt, 3 b
and from the arbitrariness of η1 , we obtain Eq. (10). 6
Theorem 7. Suppose that functions L and M are convex with respect to the variables x and and let µ ≥ 0 be a real. Define the function H as H = L + µM . Then, every solution x? of Eq. (10) is a solution of the problem C α,λ a Dt x,
J (x) =
Z
b
a
α,λ L(t, x(t), C a Dt x(t)) dt → min
under the constraints (8)–(9). Proof. Since function H is convex, by Theorem 3, we conclude that x? minimizes that is, Z
b
L[x? + η](t) dt + µ
a
Z
b
a
Once
Z
M [x? + η](t) dt ≥
b ?
M [x + η](t) dt =
a
Z
Z
b
L[x? ](t) dt + µ
a
Z
Rb a
H[x](t) dt,
b
M [x? ](t) dt.
a
b
M [x? ](t) dt = k,
a
we conclude that x? minimizes J .
For our next problem, instead of an integral constraint, we assume a restriction of form g(t, x1 (t), x2 (t)) = 0,
(11)
where g : [a, b] × R2 → R is a continuous differentiable function with respect to x1 and x2 . For simplicity, we will write g[x](t) := g(t, x1 (t), x2 (t))
and x(t) := (x1 (t), x2 (t)).
The functional that we deal now depends on two functions: J2 :
C 1 [a, b] × C 1 [a, b] → 7→
(x1 , x2 )
R Z b a
α,λ C α,λ L(t, x1 (t), x2 (t), C a Dt x1 (t), a Dt x2 (t)) dt,
α,λ where L : [a, b] × R4 → R is continuously differentiable with respect to x1 , x2 , C a Dt x1 and C α,λ a Dt x 2 .
Theorem 8. Suppose that x? = (x?1 , x?2 ) is a solution of the problem J2 (x1 , x2 ) =
Z
a
b
α,λ C α,λ L(t, x1 (t), x2 (t), C a Dt x1 (t), a Dt x2 (t)) dt → min
subject to the constraints x(a) = β
and x(b) = γ,
where β, γ ∈ R2 ,
(12)
and the holonomic constraint (11). If the maps ? t 7→ t Dα,λ b (∂4 L[x ](t)),
? t 7→ t Dα,λ b (∂5 L[x ](t)),
t ∈ [a, b],
are continuous, and if ∂3 g[x? ](t) 6= 0,
∀t ∈ [a, b],
then there exists a continuous function µ : [a, b] → R such that x? is a solution of the equations ∂i L[x](t) + t Dα,λ b (∂i+2 L[x](t)) + µ(t)∂i g[x](t) = 0, 7
t ∈ [a, b], i = 2, 3.
(13)
Proof. Let η := (η1 , η2 ) ∈ C 1 [a, b] × C 1 [a, b] be a function such that η(a) = (0, 0) = η(b). By (11) we can appeal the Implicit Function Theorem to ensure the existence of a subfamily of variations fulfilling g(t, x?1 (t) + η1 (t), x?2 (t) + η2 (t)) = 0, t ∈ [a, b]. (14) Differentiating (14) with respect to and then setting = 0, we get ∂2 g[x? ](t)η1 (t) + ∂3 g[x? ](t)η2 (t) = 0,
t ∈ [a, b].
(15)
Define the continuous function µ : [a, b] → R by the rule µ(t) := −
? ∂3 L[x? ](t) + t Dα,λ b (∂5 L[x ](t)) . ? ∂3 g[x ](t)
(16)
Then Eq. (13) is proven for the case i = 3. It remains to prove the case i = 2. By (15)–(16), we conclude that ? ? [∂3 L[x? ](t) + t Dα,λ (17) b (∂5 L[x ](t))]η2 (t) = µ(t)∂2 g[x ](t)η1 (t). On the other hand, since x? minimizes J2 , the first variation must vanish: Z
b
a
α,λ ? C α,λ ∂2 L[x? ](t)η1 (t) + ∂3 L[x? ](t)η2 (t) + ∂4 L[x? ](t)C a Dt η1 (t) + ∂5 L[x ](t)a Dt η2 (t) dt = 0,
and performing integration by parts, Z
a
b
α,λ ? ? ? [∂2 L[x? ](t) + t Dα,λ b (∂4 L[x ](t))]η1 (t) + [∂3 L[x ](t) + t Db (∂5 L[x ](t))]η2 (t) dt = 0.
(18)
Finally, using Eqs. (17)–(18), we obtain Z
a
b
? ? [∂2 L[x? ](t) + t Dα,λ b (∂4 L[x ](t)) + µ(t)∂2 g[x ](t)]η1 (t) dt = 0,
and the case i = 2 is proven. α,λ Theorem 9. Suppose that L is convex with respect to x1 , x2 , C a Dt x1 and be given by (16). If ∂3 g[x? ](t) 6= 0, ∀t ∈ [a, b]
C α,λ a Dt x 2 ,
and let µ
holds and x? satisfies Eq. (13) with i = 2, then x? minimizes J2 , subject to the restrictions (11)–(12). Proof. Given a variation x? + η with η(a) = (0, 0) = η(b), we have ?
?
J2 (x + η) − J2 (x ) ≥
=
Z
a
Z
b
∂2 L[x? ](t)η1 (t) + ∂3 L[x? ](t)η2 (t)
a
α,λ ? C α,λ + ∂4 L[x? ](t)C a Dt η1 (t) + ∂5 L[x ](t)a Dt η2 (t) dt b
α,λ ? ? ? [∂2 L[x? ](t) + t Dα,λ b (∂4 L[x ](t))]η1 (t) + [∂3 L[x ](t) + t Db (∂5 L[x ](t))]η2 (t) dt.
Since g(t, x?1 (t) + η1 (t), x?2 (t) + η2 (t)) = 0,
t ∈ [a, b],
differentiating with respect to and setting = 0, we obtain ∂3 g[x? ](t)η2 (t) = −∂2 g[x? ](t)η1 (t).
8
(19)
By (16) and (19), ? ? ? [∂3 L[x? ](t) + t Dα,λ b (∂5 L[x ](t))]η2 (t) = −µ(t)∂3 g[x ](t)η2 (t) = µ(t)∂2 g[x ](t)η1 (t).
In conclusion, J2 (x? + η) − J2 (x? ) ≥
Z
a
b
? ? [∂2 L[x? ](t) + t Dα,λ b (∂4 L[x ](t)) + µ(t)∂2 g[x ](t)]η1 (t) dt = 0.
All the previous theoretical results were concerned with the case where 0 < α < 1. However, they can be extended to functionals involving higher order derivatives. In such case, the functional depends on derivatives of fractional orders αk ∈ (k − 1, k), for k = 1, . . . , n and n ∈ N. The functional is given by Jn : C n [a, b] → R Z b αn ,λ α1 ,λ x(t)) dt, x 7→ x(t), . . . , C L(t, x(t), C a Dt a Dt a
where L : [a, b] × Rn+1 → R is continuously differentiable with respect to x and k = 1, . . . , n.
C αk ,λ x, a Dt
for
Theorem 10. Let x? be a local minimizer of functional Jn , subject to the boundary conditions x(k) (a) = βk If the maps
and x(k) (b) = γk ,
where βk , γk ∈ R, k = 0, 1, . . . , n − 1.
k ,λ t 7→ t Dα (∂k+2 L[x? ](t)), b
(20)
t ∈ [a, b], k = 1, . . . , n,
are continuous, then x? is a solution of the fractional differential equation ∂2 L[x](t) +
n X
αk ,λ (∂k+2 L[x](t)) t Db
k=1
= 0,
t ∈ [a, b].
(21)
Proof. Given a function η ∈ C n [a, b] satisfying the boundary conditions η (k) (a) = 0 = η (k) (b), k = 1, . . . , n,, we obtain Z
b
∂2 L[x? ](t)η(t) +
a
n X
αk ,λ ∂k+2 L[x? ](t)C ηk (t) dt = 0. a Dt
k=1
Integrating by parts cf. Theorem 1, we get # Z b" n X αk ,λ ? ? ∂2 L[x ](t) + (∂k+2 L[x ](t)) η(t) dt = 0. t Db a
k=1
By the arbitrariness of η, Eq. (21) follows. Theorem 11. Suppose that L is convex with respect to the variables x and 1, . . . , n. Then, every solution x? of Eq. (21) is a solution of the problem Jn (x) =
Z
a
b
α1 ,λ αn ,λ L(t, x(t), C x(t), . . . , C x(t)) dt → min a Dt a Dt
subject to (20). Proof. Similar to the proof of Theorem 3. 9
C αk ,λ x, a Dt
for k =
3
Numerical Method
In this section we provide an Euler-type numerical method for the solution of problems described in the previous section. We will follow the approach in [33], where the authors proposed a discrete direct method for variational problems involving the usual Caputo derivatives, but we will use different approximation formulas. In [33], a rectangular rule was used to approximate the integral in the functional, which, as it is known, is a first order quadrature formula. Here we will use a second order one. Then, for the discretization of the derivative, we will use a modified L1 approximation formula. In order to discretize the functional (2), we first consider a partition of the time interval [a, b] into N subintervals of equal size τ = b−a N , and we define the midpoints of each one of these subintervals through 1 ti+ 21 = a + i + τ, i = 0, 1, . . . , N − 1. 2 Using the midpoint quadrature rule to approximate the integral in (2), we obtain: J (x) =
Z
b
a
α,λ L(t, x(t), C a Dt x(t)) dt ≈ τ
N −1 X i=0
α,λ ¯ L ti+ 12 , x(ti+ 21 ), C a Dt x(ti+ 21 ) = J (x).
In [34] the following modified L1 approximation for the Caputo derivative of a function y was provided: C α a Dt y(ti+ 12 )
i X h (α) (α) (α) (α) (α) bi−k − bi−k+1 y tk− 12 − bi − Bi y t 12 ≈ τ −α b0 y ti+ 12 − k=1
(α)
where the coefficients bn
(α)
and Bn b(α) n
=
Bn(α)
=
i (α) ¯α − Bi y(a) =C a Dt yi+ 21 ,
are given by: 1 (n + 1)1−α − nα Γ(2 − α) " # 1−α 2 1 n+ − nα . Γ(2 − α) 2
The authors of [34] also proved that if y ∈ C 2 ([a, b]) then C α 2−α ¯α , i = 0, 1, . . . , N − 1, a Dt y(ti+ 21 ) −C a Dt yi+ 21 ≤ cτ
for some positive constant c which is independent of N (please see Lemma 3.1 in [34]). To approximate the left tempered Caputo derivative at t = ti+ 12 , we take into account the relationship (1), and we use the modified L1 approximation formula to obtain: C α,λ a Dt x(ti+ 21 ) i X
k=1
=e
−λti+ 1 C α 2 a Dt
e
λti+ 1
2
h −λt 1 (α) λt 1 x(ti+ 21 ) ≈ e i+ 2 τ −α b0 e i+ 2 x ti+ 12 −
λt λt i 1 1 (α) (α) (α) (α) (α) bi−k − bi−k+1 e k− 2 x tk− 12 − bi − Bi e 2 x t 12 − Bi eλa x(a) =∗
C α,λ 1. a Dt xi+ 2
If x(t) ∈ C 2 ([a, b]), then ω(t) = eλt x(t) ∈ C 2 ([a, b]), and then, there exists a positive constant c1 , independent of N , such that −λti+ 1 C α C α,λ α,λ C ¯α −λa 2−α 2 D ω(t 1 = e 1) − 1 ≤ e D x D ω cτ ≤ c1 τ 2−α . a Dt x(ti+ 21 ) −∗ C i+ 2 i+ 2 a t a t a t i+ 2 10
Combining all these approximations, we obtain J (x) ≈ τ
N −1 X i=0
˜ (x) = τ Ψ x 1 , x 3 , . . . , x 1 , L ti+ 12 , xi+ 12 , Λi = J N − 2 2 2
where xi+ 12 denotes an approximate value of x ti+ 12 , and Λi
= e −
−λti+ 1
2
(α) bi
−
Next we show that In fact, we have:
i h X λt 1 (α) λt 1 (α) (α) τ −α b0 e i+ 2 x ti+ 12 − bi−k − bi−k+1 e k− 2 x tk− 12 (α) Bi
k=1
i (α) e 2 x t 12 − Bi eλa x(a) , λt 1
i = 0, 1, . . . , N − 1.
˜ (x) → J (x) , when N → +∞. J ˜ (x) ≤ J (x) − J ¯ (x) + J ¯ (x) − J ˜ (x) . J (x) − J
(22) (23)
If L(·) ∈ C 2 ([a, b]) then, from basic quadrature theory, there must exist a positive constant C1 such that: ¯ (x) ≤ C1 τ 2 . J (x) − J (24) On the other hand,
N −1 X ¯ α,λ ∗ C α,λ ˜ (x) ≤ τ J (x) − J L ti+ 21 , x(ti+ 12 ), C . a Dt x(ti+ 21 ) − L ti+ 12 , x(ti+ 21 ), a Dt xi+ 1 2 i=0
Using the mean value theorem, there must exist a positive constant C2 such that α,λ ∗ C α,λ 1) 1 , x(t 1 ), 1 D x(t − L t D x L ti+ 12 , x(ti+ 21 ), C i+ 2 i+ 2 i+ 2 i+ 2 a t a t α,λ ∗ C α,λ ≤ C2 C ≤ C2 c1 τ 2−α = C3 τ 2−α . a Dt x(ti+ 12 ) − a Dt xi+ 1 2 Hence,
¯ ˜ (x) ≤ τ N C3 τ 2−α = (b − a)C3 τ 2−α . J (x) − J
From (23), taking (24) and (25) into account, we can conclude that: ˜ (x) ≤ Cτ 2−α , J (x) − J
(25)
(26)
for some positive constant C independent of N . From (26) we see that if N → +∞, then τ → 0 ˜ (x) → 0. and therefore J (x) − J
Remark 1. Taking (26) into account, we conclude that the method proposed for the approximation of the functional has convergence order equal to (2 − α). Consequently, smaller values of α lead to ˜ to J , while for higher values of α, the convergence becomes higher velocities of convergence of J slower. On the other hand, the tempered parameter λ does not affect the convergence order but it influences the convergence factor of the proposed method. As it is usual in the Euler-type methods for these problems, the discrete solution of the variational problem corresponds to the solution of the system of equations ∂Ψ = 0, ∂xi+ 12
i = 0, 1, . . . , N − 1.
11
(27)
The first equation is given by ∂Ψ (α) = 0 ⇔ ∂2 L t 12 , x 21 , Λ0 + τ −α ∂3 L t 12 , x 21 , Λ0 B0 ∂x 12
+τ −α
N −1 X j=1
λ (α) (α) ∂3 L tj+ 12 , xj+ 21 , Λj −bj−1 + Bj e
t 1 −tj+ 1 2
2
= 0,
and the remaining N − 1 equations by (α) ∂Ψ −α 1,x 1 , Λm 1,x 1 , Λm = 0 ⇔ ∂ L t + τ ∂ L t b0 2 3 m+ m+ m+ m+ ∂x 1 2 2 2 2 m+
+τ
2
−α
N −1 X
l=m+1
λ (α) (α) ∂3 L tl+ 12 , xl+ 12 , Λl bl−m − bl−m−1 e
tm+ 1 −tl+ 1 2
(28)
2
= 0,
(29)
for m = 1, 2, . . . , N −1. Hence, the discrete solution of the variational problem, x 12 , x 32 , . . . , xN − 21 ,
is the solution of the system of N equations (28)-(29). For the numerical approximation of the isoperimetric problem, we proceed similarly by using the same numerical approach, but where L is replaced with the function H = L + µM , where µ is a real number. Since in this case, the resulting system of the form (27) is a system of N equations on the (N + 1) unknowns x 12 , x 23 , . . . , xN − 21 , µ, an extra equation is needed. This extra equation is obtained through the discretization of (9). Using the same quadrature formulas and the same approximation for the Caputo derivatives as before, we obtain: τ
N −1 X
M (ti+ 12 , xi+ 21 , Λi ) = k.
i=0
4
Examples
Example 1. Consider the following problem: J =
Z
0
1
−t 2 2
(x(t) − e t ) +
C 0.5,1 x(t) 0 Dt
−e
−t
subject to the conditions
2 Γ
5 2
t
3 2
!2
dt → min
1 . e It is easy to see that the function x? (t) = e−t t2 satisfies the Euler-Lagrange equation (5). Also, since J ≥ 0 and J (x? ) = 0, we conclude that x? is in a fact a solution to the variational problem. We used the numerical approach that we described in the previous section, to approximate the solution of this problem for several values of N . The numerical results are presented in Table 1, where the maximum of the absolute error, E = max xi+ 12 − x? ti+ 21 x(0) = 0
and
x(1) =
i=0,...,N −1
is displayed. From the results in this table, we observe that the error decreases when N increases, as it is expected. In Figure 1, we present the plot of the analytical solution as well as the discrete one, where a good agreement between them can be observed. Next, we will consider an example, from [33], with a usual (non-tempered) Caputo derivative. Example 2. J =
Z
0
1
C 0.5 0 Dt x(t)
−
3 2 t2 Γ(2.5)
12
2
dt → min,
N E
5 1.03 ×10−2
10 4.00 ×10−3
20 1.51 ×10−3
40 5.58 ×10−4
Table 1: Maximum of the committed absolute error in the approximation of the solution for Example 1. 0.35 0.30 0.25 0.20 0.15 0.10 0.05
* * * * * * 0.00 * * * * 0.0 0.2
*
*
*
*
*
*
*
*
*
*
*
*
0.4
*
*
*
0.6
*
*
*
*
*
*
*
*
*
*
*
0.8
*
*
*
*
1.0
Figure 1: Analytical (continuous line) and discrete (* points) solution for Example 1 with N=40. subject to the conditions x(0) = 0
and x(1) = 1.
It is easy to see that the analytical solution is x(t) = t2 . In this case, we must use our approach with λ = 0. Some numerical results are listed in Table 2, where the available results from [33] are also presented. N 5 10 30
E 3.58 ×10−2 1.42 ×10−2 2.98 ×10−3
E ( [33]) 0.03 0.02 0.006
Table 2: Maximum of the committed absolute error in the approximation of the solution for Example 2. Example 3. Let us now consider the following isoperimetric problem, J =
Z
0
1
2
C 0.5,1 x(t) 0 Dt
+
2 √ 15 π 2 −t t e dt → min, 16
subject to the conditions x(0) = 0
and
x(1) =
1 e
and to the integral constraint Z 1 √ 15 π 2 −tC 0.5,1 675π 7 t e 0 Dt x(t) dt = 1− 2 . 16 1024 e 0 We can easily verify that the function x? (t) = t5/2 e−t satisfies the necessary conditions of Theorem 6, with µ = −2. We have applied the method described previously to approximate the solution of this problem. The numerical results are presented in Table 3, where the maximum of the absolute error at the time discretization points are listed as well as the obtained value for µ, for several values of N .
13
N E µ
5 2.39 ×10−2 -2.02324
10 8.76 ×10−3 -2.00574
20 3.15 ×10−3 -2.00143
40 1.12×10−3 -2.00036
Table 3: Numerical values of µ and maximum of the committed absolute error in the approximation of the solution for Example 3.
5
Conclusions
Here we have extended previous theoretical results for variational problems with Caputo fractional derivatives to the case of functionals containing tempered Caputo fractional derivatives. All the obtained results reduce to known results for the usual Caputo derivatives, by considering λ = 0. We have also presented a numerical approach which is able to deal with tempered and non-tempered fractional variational problems. Moreover, we have improved the method in [33] by considering alternative approximation formulas. Since the ones we have used here correspond to higher order approximations, it is expected that the method proposed here produces an error that goes faster to zero, as N increases, than the one presented in [33]. The results in Table 2 indicate this, although the available values of N in [33] are still relatively low to properly illustrate this. We intend to provide this comparison in a forthcoming paper where we not only illustrate numerically the conclusions in remark 1, but we will also compare this method with other methods, such as indirect methods, which solve the Euler-Lagrange equations to provide a numerical solution to the variational problems.
Acknowledgments R. Almeida is supported by Portuguese funds through the CIDMA - Center for Research and Development in Mathematics and Applications, and the Portuguese Foundation for Science and Technology (FCT-Funda¸c˜ao para a Ciˆencia e a Tecnologia), within project UID/MAT/04106/2019. M. L. Morgado acknowledges the Portuguese Foundation for Science and Technology through project UID/Multi/04621/2013. The authors are very grateful to two anonymous referees, for valuable remarks and comments that improved this paper.
References [1] A. A. Kilbas, H. M. Srivastava and J. J. Trujillo, Theory and Applications of Fractional Differential Equations. North-Holland Mathematics Studies, 204. Elsevier Science B.V., Amsterdam, 2006. [2] K. S. Miller and B. Ross, An Introduction to the Fractional Calculus and Fractional Differential Equations. A Wiley-Interscience Publication, Wiley, New York, 1993. [3] I. Podlubny, Fractional differential equations, Mathematics in Science and Engineering, 198. Academic Press, Inc., San Diego, CA, 1999. [4] S. G. Samko, A. A. Kilbas and O. I. Marichev, Fractional Integrals and Derivatives, translated from the 1987 Russian original. Gordon and Breach, Yverdon, 1993. [5] F. Sabzikar, M. M. Meerschaert and J. H. Chen, Tempered fractional calculus. J. Comput. Phys. 293 (2015) 14–28. [6] J. W. Deng, L. J. Zhao and Y. J. Wu, Fast predictor-corrector approach for the tempered fractional ordinary differential equations. Numer. Algor. 74 (3) (2017) 717–754.
14
[7] K. Diethelm, The Analysis of Fractional Differential Equations: An Application-Oriented Exposition Using Differential Operators of Caputo Type. Springer, 2010. [8] B. Baeumer and M. M. Meerschaert, Tempered stable L´evy motion and transient superdiffusion. J. Comput. Appl. Math. 233 (10) (2010) 2438–2448. [9] F. Riewe, Nonconservative Lagrangian and Hamiltonian mechanics. Phys. Rev. E 53 (1996) 1890. [10] O. P. Agrawal, Formulation of Euler–Lagrange equations for fractional variational problems. J. Math. Anal. Appl. 272 (2002) 368–379. [11] T. M. Atanackovi´c, S. Konjik and S. Pilipovi´c, Variational problems with fractional derivatives: Euler–Lagrange equations. J. Phys. A: Math. Theor. 41 (2008) 095201. [12] D. Baleanu, S. I. Muslih and E. M. Rabei, On fractional Euler–Lagrange and Hamilton equations and the fractional generalization of total time derivative. Nonlinear Dynam. 53 (1–2) (2008) 67–74. [13] M. J. Lazo and D. F. M. Torres, The Legendre condition of the fractional calculus of variations. Optimization 63 (2014) 1157–1165. [14] D. Baleanu, S. I. Muslih and K. Ta¸s, Fractional Hamiltonian analysis of higher order derivatives systems. J. Math. Phys. 47 (10) (2006) 103503. [15] M. A. E. Herzallah and D. Baleanu, Fractional Euler–Lagrange equations revisited. Nonlinear Dynam. 69 (3) (2012) 977–982. [16] D. Baleanu and S. I. Muslih, Lagrangian Formulation of Classical Fields within Riemann– Liouville Fractional Derivatives. Phys. Scr. 72 (2–3) (2005) 119–121. [17] M. Klimek, Lagrangean and Hamiltonian fractional sequential mechanics. Czechoslovak J. Phys. 52 (11) (2002) 1247–1253. [18] A. Lotfi and , S. A. Yousefi, A numerical technique for solving a class of fractional variational problem. J. Comput. Appl. Math. 237 (1) (2013) 633–643. [19] O. P. Agrawal, Fractional variational calculus and the transversality conditions. J. Phys. A: Math. Gen. 39 (2006) 10375–10384. [20] O. P. Agrawal, Generalized Euler–Lagrange equations and transversality conditions for FVPs in terms of the Caputo derivative. J. Vib. Control 13 (2007) 1217–1237. [21] R. Almeida and D. F. M. Torres, Necessary and sufficient conditions for the fractional calculus of variations with Caputo derivatives. Commun. Nonlinear Sci. Numer. Simul. 16 (2011) 1490– 1500. [22] D. Baleanu, Fractional constrained systems and Caputo derivatives. J. Comput. Nonlinear Dyn. 3 (2008) 021102. [23] M. A. E. Herzallah and D. Baleanu, Fractional-order Euler–Lagrange equations and formulation of Hamiltonian equations. Nonlinear Dynam. 58 (1) (2009) 385. [24] O. P. Agrawal, Fractional variational calculus in terms of Riesz fractional derivatives. J. Phys. A 40 (24) (2007) 6287–6303. [25] R. Almeida, Fractional variational problems with the Riesz–Caputo derivative. Appl. Math. Lett. 25 (2012) 142–148. [26] M. Klimek, Fractional sequential mechanics – models with symmetric fractional derivative. Czechoslovak J. Phys. 51 (12) (2001) 1348–1354. 15
[27] D. Tavares, R. Almeida and D. F.M. Torres, Combined fractional variational problems of variable order and some computational aspects. J. Comput. Appl. Math. 339 (2018) 374–388. [28] O. S. Fard and M. Salehi, A survey on fuzzy fractional variational problems. J. Comput. Appl. Math. 271 (2014) 71–82. [29] I. Matychyn and V. Onyshchenko, On time-optimal control of fractional-order systems. J. Comput. Appl. Math. 339 (2018) 245–257. [30] R. Almeida and M. L. Morgado, The Euler–Lagrange and Legendre equations for functionals involving distributed-order fractional derivatives. Appl. Math. Comput. 331 (2018) 394–403. [31] R. Almeida, S. Pooseh and D. F. M. Torres, Computational Methods in the Fractional Calculus of Variations. Imperial College Press, London, 2015. [32] A. B. Malinowska and D. F. M. Torres, Introduction to the Fractional Calculus of Variations. Imperial College Press, London, 2012. [33] S. Pooseh, R. Almeida and D. F. M. Torres, Discrete Direct Methods in the Fractional Calculus of Variations. Comput. Math. Appl. 66 (2013) 668–676. [34] F. Zeng and C. Li, A new Crank–Nicolson finite element method for the time-fractional subdiffusion equation. Appl. Num. Mathem. 121 (2017) 82–95.
16
PII: DOI: Reference:
S0377-0427(19)30190-6 https://doi.org/10.1016/j.cam.2019.04.010 CAM 12230
To appear in:
Journal of Computational and Applied Mathematics
Received date : 27 June 2018 Revised date : 28 February 2019 Please cite this article as: R. Almeida and M.L. Morgado, Analysis and numerical approximation of tempered fractional calculus of variations problems, Journal of Computational and Applied Mathematics (2019), https://doi.org/10.1016/j.cam.2019.04.010 This is a PDF file of an unedited manuscript that has been accepted for publication. As a service to our customers we are providing this early version of the manuscript. The manuscript will undergo copyediting, typesetting, and review of the resulting proof before it is published in its final form. Please note that during the production process errors may be discovered which could affect the content, and all legal disclaimers that apply to the journal pertain.
Analysis and Numerical Approximation of Tempered Fractional Calculus of Variations problems Ricardo Almeida1 and M. Lu´ısa Morgado2 1
Center for Research and Development in Mathematics and Applications (CIDMA) Department of Mathematics, University of Aveiro, 3810–193 Aveiro, Portugal, [email protected] 2 Center for Computational and Stochastic Mathematics, Instituto Superior T´ecnico and Department of Mathematics, University of Tr´as-os-Montes e Alto Douro,UTAD 5000-801, Vila Real, Portugal [email protected] Abstract In this paper, we study variational problems involving the tempered Caputo fractional derivative. The cost functional involves the tempered Caputo fractional derivative. Several important optimization conditions are derived to find the optimal solution. Sufficient and necessary conditions are presented for different variational problems. For example, the cases of integral (isoperimetric problem) and holonomic constraints are considered, as well as problems with high order derivatives. A numerical scheme is proposed to determine approximations of the solution and it is illustrated through some examples.
Mathematics Subject Classification 2010: 26A33, 49K05, 49M05. Keywords: Tempered fractional derivative, Euler–Lagrange equation, Numerical methods.
1
Introduction
Since its very beginning that Calculus has been used to describe many physical phenomena. Although Fractional Calculus (FC), which consists of the study of integral and differential operators whose order is not restricted to integer values, is as much as old as the classical (integer order) one, only recently it has been intensively used in the modelling of many physical and engineering processes. FC may be regarded as a generalization of the integer order one. However, several (and nonequivalent) definitions of fractional order derivatives exist. We refer the reader to the books [1–4] to see the potential applications of FC and also to the variety of different definitions for fractional derivatives. We dare to say that the most popular fractional derivatives are the Riemann-Liouville, the Gr¨ unwald-Letnikov and the Caputo derivatives, being the last one the most commonly found in applicatons. There are two reasons for such choice. First, the Caputo fractional derivative of a constant is the null function. Secondly, the Laplace transform of the Caputo fractional derivative involves integer-order initial conditions. In this paper we are concerned with the theoretical analysis and numerical approximation of functionals involving tempered fractional derivatives in the Caputo sense. These new fractional operators are a generalization of the classical ones, in the sense that they depend on a new parameter λ, and when λ = 0, we obtain the well-know fractional integrals/derivatives. Left and right sided tempered fractional integrals and derivatives can be given trough the definition of the fractional integrals and derivatives (see [5, 6] for example).
1
Definition 1. For any α > 0 and fixed parameter λ ≥ 0, the left and right-sided tempered Riemann-Liouville fractional integrals of a function x on (a, b), are given by Z t 1 α,λ α λt −λt (t − s)α−1 e−λ(t−s) x(s) ds, a It e x(t) = a It x(t) = e Γ(α) a Z b 1 α,λ −λt λt α x(t) = (s − t)α−1 e−λ(s−t) x(s) ds, t Ib x(t) = e t Ib e Γ(α) t
respectively, where a Itα and t Ibα are the left and right Riemann-Liouville fractional integrals ( [7]). Note that when λ = 0, the tempered Riemann-Liouville fractional integrals reduce to the (usual) Riemann-Liouville integrals. Definition 2. For any α > 0, let n ∈ N be given by n = [α] + 1 if α ∈ / N, or n = α otherwise, and let λ ≥ 0 be a fixed parameter. The left and right-sided tempered fractional Caputo derivatives of a function x are given by: n Z t e−λt d C α,λ −λt C α λt n−α−1 (t − s) eλs x(s) ds, (1) a Dt x(t) = e a Dt e x(t) = Γ(n − α) a ds n Z b d eλt C α,λ n−α−1 λt C α λt D (s − t) − e−λs x(s) ds, x(t) = e D −e x(t) = t t b b Γ(n − α) t ds
α C α respectively, where C a Dt and t Db denote the left and right Caputo fractional derivatives, respectively ( [7]).
Once again if λ = 0 then the Caputo tempered fractional derivatives reduce to the Caputo fractional derivatives, and therefore, tempered fractional integrals and derivatives can be regarded as a generalization of the usual (non-tempered) fractional integrals and derivatives. Finally, we present the equivalent definitions for the Riemann–Liouville fractional derivatives. Definition 3. For any α > 0, let n = [α] + 1, and let λ ≥ 0 be a fixed parameter. The left and right-sided tempered fractional Riemann–Liouville derivatives of a function x are given by: n Z t d e−λt α,λ −λt α λt D x(t) = e D e x(t) = (t − s)n−α−1 eλs x(s)ds, a t a t Γ(n − α) dt a n Z b d eλt α,λ λt α λt D x(t) = e D −e x(t) = − (s − t)n−α−1 e−λs x(s)ds, t b t b Γ(n − α) dt t
respectively, where a Dtα and t Dbα denote the left and right Riemann–Liouville fractional derivatives, respectively ( [7]). Tempered fractional derivatives have appeared as one of the solutions to overcome the problem of the divergence of the second moment (or variance) of the probability density function of −1−α L´evy flights, which decays as |x| . This L´evy distribution originates the substancial (usual) fractional Riemann-Liouville derivatives of order α. By exponentially tempering the probability of −1−α −λ|x| large jumps of L´evy flights with a parameter λ > 0, the density function decays as |x| e and then finite second moments are obtained and the definitions of tempered fractional derivatives arise ( [8]). In this work we are concerned with Fractional Calculus of Variations problems. This subject has been extensively developed in the last decades, since the pioneering works of Riewe [9]. Since then, numerous works have appeared, dealing with different fractional operators. To mention a few of them, we refer to the case with the Riemann–Liouville derivative [10–18], the Caputo derivative [19–23], the Riesz derivative [24, 25], symmetric fractional derivative [26, 27], fuzzy fractional variational problems [28], optimal control problems [29], distributed-order derivatives [30], etc. We also mention the recent books [31, 32], where analytical and numerical methods were developed to deal with fractional variational problems. 2
Here we will be concerned with functionals of the type: → R Z b α,λ 7 → L(t, x(t), C a Dt x(t)) dt,
J : C 1 [a, b] x
(2)
a
where L : [a, b] × R2 → R is continuously differentiable with respect to the second and third variables. The fractional order α is a real between 0 and 1, and λ ≥ 0. The following notation will be used throughout the text. Given a function of several variables L ≡ L(a1 , . . . , am ), we use the symbol ∂L . ∂i L(a1 , . . . , am ) := ∂ai (a1 ,...,am )
α,λ Also, the vector (t, x(t), C a Dt x(t)) will be abbreviated by [x](t). A function L(a1 , a2 , . . . , am ) is convex with respect to the variables a2 , . . . , am if ∂i L exists for any i ∈ {2, . . . , m} and the following inequality holds:
L(a1 , a2 +a02 . . . , am +a0m )−L(a1 , a2 , . . . , am ) ≥ ∂2 L(a1 , a2 , . . . , am )a02 +. . .+∂m L(a1 , a2 , . . . , am )a0m , for all vectors (a1 , a2 + a02 . . . , am + a0m ) and (a1 , a2 , . . . , am ). The paper is organized as follows. In Section 2 we present our theoretical study. Several necessary and sufficient conditions are proven, for different variational problems. We consider the fundamental problem and deduce the respective Euler–Lagrange equation and the transversality conditions. Two types of constraints are after considered: the integral and holonomic problems. Then, a generalization of the fundamental problem is studied, where the Lagrange function depends on high order derivatives. Then, in Section 3, a numerical procedure is provided, following a usual direct approach as in [33], where variational problems with (non-tempered) Caputo derivatives have been considered. Here, we use different approximation formulas, namely, we use a higher order quadrature formula for the approximation of the integral in the functional. In this way, we not only extend the results in [33] (which here correspond to the case where the tempered parameter λ is zero), but we also improve the numerical scheme. Some examples for which the analytical solution is known are considered to illustrate the efficiency of the method.
2
Necessary and sufficient conditions of optimality
To start, we will present and prove a fractional integration by parts formula needed for the sequel. Theorem 1. Let x ∈ C[a, b] and y ∈ C n [a, b] be two functions. Then, for α > 0 and n ∈ N such that n = [α] + 1, the following formula holds: Z
a
b
α,λ x(t)C a Dt y(t) dt
=
Z
a
b
y(t)t Dα,λ b x(t) dt+
# "n−1 k b X dn−k−1 d λt −λt n−α,λ e y(t) − e t Ib x(t) . dtn−k−1 dt k=0
a
Proof. By definition of the tempered Caputo fractional derivative, Z
a
b
α,λ x(t)C a Dt y(t) dt =
Z
a
b
Z
a
t
x(t)e−λt dn (t − s)n−α−1 n eλs y(s) ds dt. Γ(n − α) ds
If we apply the Dirichlet’s formula, the last formula is equal to Z
a
b
dn λt 1 e y(t) dtn Γ(n − α)
Z
b
t
3
(s − t)n−α−1 e−λs x(s) ds dt.
Using now integration by parts, we obtain Z
b
a
"
#b Z b dn−1 λt 1 n−α−1 −λs = (s − t) e y(t) e x(s) ds dtn−1 Γ(n − α) t a ! Z b Z b n−1 d 1 d λt n−α−1 −λs e y(t) − (s − t) + e x(s) ds dt. n−1 dt Γ(n − α) t a dt
α,λ x(t)C a Dt y(t) dt
Integrating once again by parts, we get Z
b
a
# k 1 b X dn−k−1 λt d −λt n−α,λ x(t) = e t Ib e y(t) − dtn−k−1 dt k=0 a ! 2 Z b n−2 Z b d 1 d λt n−α−1 −λs + e y(t) − (s − t) e x(s) ds dt. n−2 dt Γ(n − α) t a dt
α,λ x(t)C a Dt y(t) dt
"
Repeating the procedure, we deduce the following: Z
a
b
# "n−1 k b X dn−k−1 d λt −λt n−α,λ = e y(t) − e t Ib x(t) dtn−k−1 dt k=0 a ! n Z b Z b 1 d + (s − t)n−α−1 e−λs x(s) ds dt eλt y(t) − dt Γ(n − α) t a "n−1 # k b Z b X dn−k−1 d λt −λt n−α,λ e y(t) − = e t Ib x(t) + y(t)t Dα,λ b x(t) dt, dtn−k−1 dt a
α,λ x(t)C a Dt y(t) dt
k=0
a
proving the desired formula. In particular, when the fractional order α is between 0 and 1, we obtain the simple formula: Z
a
b
α,λ x(t)C a Dt y(t) dt =
Z
a
b
h ib 1−α,λ y(t)t Dα,λ x(t) dt + y(t) I x(t) . t b b a
(3)
We are now in conditions to present a necessary condition that allows us to find the candidates of minimizers for functionals J as in (2). This type of conditions is known in the literature as Euler–Lagrange equations. Theorem 2. Let x? be a local minimizer of functional J as in (2), subject to the boundary conditions x(a) = β and x(b) = γ, where β, γ ∈ R. (4) If the map
? t 7→ t Dα,λ b (∂3 L[x ](t)),
t ∈ [a, b],
is continuous, then x? is a solution of the fractional differential equation ∂2 L[x](t) + t Dα,λ b (∂3 L[x](t)) = 0,
t ∈ [a, b].
(5)
Proof. Let η ∈ C 1 [a, b] be an arbitrary function satisfying the boundary conditions η(a) = 0 = η(b) and a real number. Extrema are characterized by setting the first variation of the functional equal to zero. Since x? + η is an admissible variation for the problem, J (x? + η) − J (x? ) = 0, →0
δJ (x? , η) := lim
4
we obtain
Z
b
a
α,λ ∂2 L[x? ](t)η(t) + ∂3 L[x? ](t)C a Dt η(t) dt = 0,
and using the integration by parts formula (3), Z
b
a
h
i h ib 1−α,λ ? ? ∂2 L[x? ](t) + t Dα,λ (∂ L[x ](t)) η(t) dt + η(t) I (∂ L[x ](t)) = 0. 3 t 3 b b a
Since η is an arbitrary function vanishing at t = a and t = b, we get that x? is a solution of Eq. (5). Definition 4. Function x? is said to be an extremal for functional J in (2) if ? ∂2 L[x? ](t) + t Dα,λ b (∂3 L[x ](t)) = 0,
t ∈ [a, b].
Theorem 3. Suppose that L is convex with respect to the variables x and solution x? of Eq. (5) is a solution of the problem J (x) =
Z
b
a
α,λ L(t, x(t), C a Dt x(t)) dt → min
s.t. x(a) = β
C α,λ a Dt x.
Then, every
and x(b) = γ.
Proof. Consider a variation of x? of form x? + η with η(a) = 0 = η(b). Using the convexity assumption and applying integration by parts formula (3), we get ?
?
J (x + η) − J (x ) ≥ =
Z
b
a
Z
a
b
α,λ ∂2 L[x? ](t)η(t) + ∂3 L[x? ](t)C a Dt η(t) dt
h
i h ib 1−α,λ ? ∂2 L[x? ](t) + t Dα,λ (∂3 L[x? ](t)) b (∂3 L[x ](t)) η(t) dt + η(t)t Ib
a
which is equal to zero since x? satisfies Eq. (5). This proves that J (x? ) is a local minimum of J. Theorem 4. Let x? be a local minimizer of functional J as in (2). If the map ? t 7→ t Dα,λ b (∂3 L[x ](t)),
t ∈ [a, b],
is continuous, then x? is a solution of the tempered fractional differential equation ∂2 L[x](t) + t Dα,λ b (∂3 L[x](t)) = 0,
t ∈ [a, b].
(6)
at t = a and t = b.
(7)
Moreover, the two next transversality conditions hold: 1−α,λ (∂3 L[x? ](t)) t Ib
=0
Proof. As it was seen in the proof of Theorem 2, Z
a
b
h
i h ib 1−α,λ ? ∂2 L[x? ](t) + t Dα,λ (∂3 L[x? ](t)) = 0. b (∂3 L[x ](t)) η(t) dt + η(t)t Ib a
If we assume that η(a) = 0 = η(b), then we obtain Eq. (6). Therefore, h ib η(t)t I1−α,λ (∂3 L[x? ](t)) = 0, b a
and by the arbitrariness of η at t = a and t = b Eqs. (7) are obtained.
5
Theorem 5. Suppose that L is convex with respect to the variables x and solution x? of Eqs. (6)-(7) is a solution of the problem J (x) =
Z
b
a
C α,λ a Dt x.
Then, every
α,λ L(t, x(t), C a Dt x(t)) dt → min .
Proof. Similar to the proof of Theorem 3. We proceed our work considering the isoperimetric problem. In this case, besides the boundary conditions (4), we have an integral constraint. In what follows, M : [a, b]×R2 → R is a continuously differentiable function with respect to the second and third variables and k ∈ R is a constant. Theorem 6. Suppose that x? is a solution of the problem: Z b α,λ J (x) = L(t, x(t), C a Dt x(t)) dt → min a
subject to the constrains x(a) = β
and x(b) = γ,
and I (x) := If the maps
Z
a
b
where β, γ ∈ R,
α,λ M (t, x(t), C a Dt x(t)) dt = k.
? t 7→ t Dα,λ b (∂3 L[x ](t)), ?
? t 7→ t Dα,λ b (∂3 M [x ](t)),
(8)
(9) t ∈ [a, b],
are continuous, and x is not an extremal for I, then there exists a real number µ such that x? is a solution of the equation ∂2 H[x](t) + t Dα,λ b (∂3 H[x](t)) = 0,
t ∈ [a, b],
(10)
where H is the function defined by H := L + µM . Proof. Consider a variation of x? with two parameters x? + 1 η1 + 2 η2 , where ηi ∈ C 1 [a, b] and ηi (a) = 0 = ηi (b), for i = 1, 2. Consider the two functions j, i (defined in a neighborhood of zero) as j(1 , 2 ) :=J (x? + 1 η1 + 2 η2 ) i(1 , 2 ) :=I(x? + 1 η1 + 2 η2 ) − k. Since ∂i (0, 0) = ∂2
Z
a
b
h
i ? ∂2 M [x? ](t) + t Dα,λ b (∂3 M [x ](t)) η2 (t) dt,
and x? is not an extremal for I, we conclude that there must exist a function η2 for which ∂i/∂2 (0, 0) 6= 0. Thus, by the Implicit Function Theorem, there exists some function 2 (·), defined in a neighborhood of zero, such that i(1 , 2 (1 )) = 0. From now on, we will consider such subfamily of variations that fulfills the isoperimetric constraint (9). Since ∇i(0, 0) 6= (0, 0), and (0, 0) is a solution of the problem: minimize j
subject to i(1 , 2 ) = 0,
we can ensure the existence of a real µ with ∇(j + µi)(0, 0) = (0, 0). Fom the relation ∂(j + µi) (0, 0) = ∂1
Z
a
b
h i ? ∂2 (L + µM )[x? ](t) + t Dα,λ (∂ (L + µM )[x ](t)) η1 (t) dt, 3 b
and from the arbitrariness of η1 , we obtain Eq. (10). 6
Theorem 7. Suppose that functions L and M are convex with respect to the variables x and and let µ ≥ 0 be a real. Define the function H as H = L + µM . Then, every solution x? of Eq. (10) is a solution of the problem C α,λ a Dt x,
J (x) =
Z
b
a
α,λ L(t, x(t), C a Dt x(t)) dt → min
under the constraints (8)–(9). Proof. Since function H is convex, by Theorem 3, we conclude that x? minimizes that is, Z
b
L[x? + η](t) dt + µ
a
Z
b
a
Once
Z
M [x? + η](t) dt ≥
b ?
M [x + η](t) dt =
a
Z
Z
b
L[x? ](t) dt + µ
a
Z
Rb a
H[x](t) dt,
b
M [x? ](t) dt.
a
b
M [x? ](t) dt = k,
a
we conclude that x? minimizes J .
For our next problem, instead of an integral constraint, we assume a restriction of form g(t, x1 (t), x2 (t)) = 0,
(11)
where g : [a, b] × R2 → R is a continuous differentiable function with respect to x1 and x2 . For simplicity, we will write g[x](t) := g(t, x1 (t), x2 (t))
and x(t) := (x1 (t), x2 (t)).
The functional that we deal now depends on two functions: J2 :
C 1 [a, b] × C 1 [a, b] → 7→
(x1 , x2 )
R Z b a
α,λ C α,λ L(t, x1 (t), x2 (t), C a Dt x1 (t), a Dt x2 (t)) dt,
α,λ where L : [a, b] × R4 → R is continuously differentiable with respect to x1 , x2 , C a Dt x1 and C α,λ a Dt x 2 .
Theorem 8. Suppose that x? = (x?1 , x?2 ) is a solution of the problem J2 (x1 , x2 ) =
Z
a
b
α,λ C α,λ L(t, x1 (t), x2 (t), C a Dt x1 (t), a Dt x2 (t)) dt → min
subject to the constraints x(a) = β
and x(b) = γ,
where β, γ ∈ R2 ,
(12)
and the holonomic constraint (11). If the maps ? t 7→ t Dα,λ b (∂4 L[x ](t)),
? t 7→ t Dα,λ b (∂5 L[x ](t)),
t ∈ [a, b],
are continuous, and if ∂3 g[x? ](t) 6= 0,
∀t ∈ [a, b],
then there exists a continuous function µ : [a, b] → R such that x? is a solution of the equations ∂i L[x](t) + t Dα,λ b (∂i+2 L[x](t)) + µ(t)∂i g[x](t) = 0, 7
t ∈ [a, b], i = 2, 3.
(13)
Proof. Let η := (η1 , η2 ) ∈ C 1 [a, b] × C 1 [a, b] be a function such that η(a) = (0, 0) = η(b). By (11) we can appeal the Implicit Function Theorem to ensure the existence of a subfamily of variations fulfilling g(t, x?1 (t) + η1 (t), x?2 (t) + η2 (t)) = 0, t ∈ [a, b]. (14) Differentiating (14) with respect to and then setting = 0, we get ∂2 g[x? ](t)η1 (t) + ∂3 g[x? ](t)η2 (t) = 0,
t ∈ [a, b].
(15)
Define the continuous function µ : [a, b] → R by the rule µ(t) := −
? ∂3 L[x? ](t) + t Dα,λ b (∂5 L[x ](t)) . ? ∂3 g[x ](t)
(16)
Then Eq. (13) is proven for the case i = 3. It remains to prove the case i = 2. By (15)–(16), we conclude that ? ? [∂3 L[x? ](t) + t Dα,λ (17) b (∂5 L[x ](t))]η2 (t) = µ(t)∂2 g[x ](t)η1 (t). On the other hand, since x? minimizes J2 , the first variation must vanish: Z
b
a
α,λ ? C α,λ ∂2 L[x? ](t)η1 (t) + ∂3 L[x? ](t)η2 (t) + ∂4 L[x? ](t)C a Dt η1 (t) + ∂5 L[x ](t)a Dt η2 (t) dt = 0,
and performing integration by parts, Z
a
b
α,λ ? ? ? [∂2 L[x? ](t) + t Dα,λ b (∂4 L[x ](t))]η1 (t) + [∂3 L[x ](t) + t Db (∂5 L[x ](t))]η2 (t) dt = 0.
(18)
Finally, using Eqs. (17)–(18), we obtain Z
a
b
? ? [∂2 L[x? ](t) + t Dα,λ b (∂4 L[x ](t)) + µ(t)∂2 g[x ](t)]η1 (t) dt = 0,
and the case i = 2 is proven. α,λ Theorem 9. Suppose that L is convex with respect to x1 , x2 , C a Dt x1 and be given by (16). If ∂3 g[x? ](t) 6= 0, ∀t ∈ [a, b]
C α,λ a Dt x 2 ,
and let µ
holds and x? satisfies Eq. (13) with i = 2, then x? minimizes J2 , subject to the restrictions (11)–(12). Proof. Given a variation x? + η with η(a) = (0, 0) = η(b), we have ?
?
J2 (x + η) − J2 (x ) ≥
=
Z
a
Z
b
∂2 L[x? ](t)η1 (t) + ∂3 L[x? ](t)η2 (t)
a
α,λ ? C α,λ + ∂4 L[x? ](t)C a Dt η1 (t) + ∂5 L[x ](t)a Dt η2 (t) dt b
α,λ ? ? ? [∂2 L[x? ](t) + t Dα,λ b (∂4 L[x ](t))]η1 (t) + [∂3 L[x ](t) + t Db (∂5 L[x ](t))]η2 (t) dt.
Since g(t, x?1 (t) + η1 (t), x?2 (t) + η2 (t)) = 0,
t ∈ [a, b],
differentiating with respect to and setting = 0, we obtain ∂3 g[x? ](t)η2 (t) = −∂2 g[x? ](t)η1 (t).
8
(19)
By (16) and (19), ? ? ? [∂3 L[x? ](t) + t Dα,λ b (∂5 L[x ](t))]η2 (t) = −µ(t)∂3 g[x ](t)η2 (t) = µ(t)∂2 g[x ](t)η1 (t).
In conclusion, J2 (x? + η) − J2 (x? ) ≥
Z
a
b
? ? [∂2 L[x? ](t) + t Dα,λ b (∂4 L[x ](t)) + µ(t)∂2 g[x ](t)]η1 (t) dt = 0.
All the previous theoretical results were concerned with the case where 0 < α < 1. However, they can be extended to functionals involving higher order derivatives. In such case, the functional depends on derivatives of fractional orders αk ∈ (k − 1, k), for k = 1, . . . , n and n ∈ N. The functional is given by Jn : C n [a, b] → R Z b αn ,λ α1 ,λ x(t)) dt, x 7→ x(t), . . . , C L(t, x(t), C a Dt a Dt a
where L : [a, b] × Rn+1 → R is continuously differentiable with respect to x and k = 1, . . . , n.
C αk ,λ x, a Dt
for
Theorem 10. Let x? be a local minimizer of functional Jn , subject to the boundary conditions x(k) (a) = βk If the maps
and x(k) (b) = γk ,
where βk , γk ∈ R, k = 0, 1, . . . , n − 1.
k ,λ t 7→ t Dα (∂k+2 L[x? ](t)), b
(20)
t ∈ [a, b], k = 1, . . . , n,
are continuous, then x? is a solution of the fractional differential equation ∂2 L[x](t) +
n X
αk ,λ (∂k+2 L[x](t)) t Db
k=1
= 0,
t ∈ [a, b].
(21)
Proof. Given a function η ∈ C n [a, b] satisfying the boundary conditions η (k) (a) = 0 = η (k) (b), k = 1, . . . , n,, we obtain Z
b
∂2 L[x? ](t)η(t) +
a
n X
αk ,λ ∂k+2 L[x? ](t)C ηk (t) dt = 0. a Dt
k=1
Integrating by parts cf. Theorem 1, we get # Z b" n X αk ,λ ? ? ∂2 L[x ](t) + (∂k+2 L[x ](t)) η(t) dt = 0. t Db a
k=1
By the arbitrariness of η, Eq. (21) follows. Theorem 11. Suppose that L is convex with respect to the variables x and 1, . . . , n. Then, every solution x? of Eq. (21) is a solution of the problem Jn (x) =
Z
a
b
α1 ,λ αn ,λ L(t, x(t), C x(t), . . . , C x(t)) dt → min a Dt a Dt
subject to (20). Proof. Similar to the proof of Theorem 3. 9
C αk ,λ x, a Dt
for k =
3
Numerical Method
In this section we provide an Euler-type numerical method for the solution of problems described in the previous section. We will follow the approach in [33], where the authors proposed a discrete direct method for variational problems involving the usual Caputo derivatives, but we will use different approximation formulas. In [33], a rectangular rule was used to approximate the integral in the functional, which, as it is known, is a first order quadrature formula. Here we will use a second order one. Then, for the discretization of the derivative, we will use a modified L1 approximation formula. In order to discretize the functional (2), we first consider a partition of the time interval [a, b] into N subintervals of equal size τ = b−a N , and we define the midpoints of each one of these subintervals through 1 ti+ 21 = a + i + τ, i = 0, 1, . . . , N − 1. 2 Using the midpoint quadrature rule to approximate the integral in (2), we obtain: J (x) =
Z
b
a
α,λ L(t, x(t), C a Dt x(t)) dt ≈ τ
N −1 X i=0
α,λ ¯ L ti+ 12 , x(ti+ 21 ), C a Dt x(ti+ 21 ) = J (x).
In [34] the following modified L1 approximation for the Caputo derivative of a function y was provided: C α a Dt y(ti+ 12 )
i X h (α) (α) (α) (α) (α) bi−k − bi−k+1 y tk− 12 − bi − Bi y t 12 ≈ τ −α b0 y ti+ 12 − k=1
(α)
where the coefficients bn
(α)
and Bn b(α) n
=
Bn(α)
=
i (α) ¯α − Bi y(a) =C a Dt yi+ 21 ,
are given by: 1 (n + 1)1−α − nα Γ(2 − α) " # 1−α 2 1 n+ − nα . Γ(2 − α) 2
The authors of [34] also proved that if y ∈ C 2 ([a, b]) then C α 2−α ¯α , i = 0, 1, . . . , N − 1, a Dt y(ti+ 21 ) −C a Dt yi+ 21 ≤ cτ
for some positive constant c which is independent of N (please see Lemma 3.1 in [34]). To approximate the left tempered Caputo derivative at t = ti+ 12 , we take into account the relationship (1), and we use the modified L1 approximation formula to obtain: C α,λ a Dt x(ti+ 21 ) i X
k=1
=e
−λti+ 1 C α 2 a Dt
e
λti+ 1
2
h −λt 1 (α) λt 1 x(ti+ 21 ) ≈ e i+ 2 τ −α b0 e i+ 2 x ti+ 12 −
λt λt i 1 1 (α) (α) (α) (α) (α) bi−k − bi−k+1 e k− 2 x tk− 12 − bi − Bi e 2 x t 12 − Bi eλa x(a) =∗
C α,λ 1. a Dt xi+ 2
If x(t) ∈ C 2 ([a, b]), then ω(t) = eλt x(t) ∈ C 2 ([a, b]), and then, there exists a positive constant c1 , independent of N , such that −λti+ 1 C α C α,λ α,λ C ¯α −λa 2−α 2 D ω(t 1 = e 1) − 1 ≤ e D x D ω cτ ≤ c1 τ 2−α . a Dt x(ti+ 21 ) −∗ C i+ 2 i+ 2 a t a t a t i+ 2 10
Combining all these approximations, we obtain J (x) ≈ τ
N −1 X i=0
˜ (x) = τ Ψ x 1 , x 3 , . . . , x 1 , L ti+ 12 , xi+ 12 , Λi = J N − 2 2 2
where xi+ 12 denotes an approximate value of x ti+ 12 , and Λi
= e −
−λti+ 1
2
(α) bi
−
Next we show that In fact, we have:
i h X λt 1 (α) λt 1 (α) (α) τ −α b0 e i+ 2 x ti+ 12 − bi−k − bi−k+1 e k− 2 x tk− 12 (α) Bi
k=1
i (α) e 2 x t 12 − Bi eλa x(a) , λt 1
i = 0, 1, . . . , N − 1.
˜ (x) → J (x) , when N → +∞. J ˜ (x) ≤ J (x) − J ¯ (x) + J ¯ (x) − J ˜ (x) . J (x) − J
(22) (23)
If L(·) ∈ C 2 ([a, b]) then, from basic quadrature theory, there must exist a positive constant C1 such that: ¯ (x) ≤ C1 τ 2 . J (x) − J (24) On the other hand,
N −1 X ¯ α,λ ∗ C α,λ ˜ (x) ≤ τ J (x) − J L ti+ 21 , x(ti+ 12 ), C . a Dt x(ti+ 21 ) − L ti+ 12 , x(ti+ 21 ), a Dt xi+ 1 2 i=0
Using the mean value theorem, there must exist a positive constant C2 such that α,λ ∗ C α,λ 1) 1 , x(t 1 ), 1 D x(t − L t D x L ti+ 12 , x(ti+ 21 ), C i+ 2 i+ 2 i+ 2 i+ 2 a t a t α,λ ∗ C α,λ ≤ C2 C ≤ C2 c1 τ 2−α = C3 τ 2−α . a Dt x(ti+ 12 ) − a Dt xi+ 1 2 Hence,
¯ ˜ (x) ≤ τ N C3 τ 2−α = (b − a)C3 τ 2−α . J (x) − J
From (23), taking (24) and (25) into account, we can conclude that: ˜ (x) ≤ Cτ 2−α , J (x) − J
(25)
(26)
for some positive constant C independent of N . From (26) we see that if N → +∞, then τ → 0 ˜ (x) → 0. and therefore J (x) − J
Remark 1. Taking (26) into account, we conclude that the method proposed for the approximation of the functional has convergence order equal to (2 − α). Consequently, smaller values of α lead to ˜ to J , while for higher values of α, the convergence becomes higher velocities of convergence of J slower. On the other hand, the tempered parameter λ does not affect the convergence order but it influences the convergence factor of the proposed method. As it is usual in the Euler-type methods for these problems, the discrete solution of the variational problem corresponds to the solution of the system of equations ∂Ψ = 0, ∂xi+ 12
i = 0, 1, . . . , N − 1.
11
(27)
The first equation is given by ∂Ψ (α) = 0 ⇔ ∂2 L t 12 , x 21 , Λ0 + τ −α ∂3 L t 12 , x 21 , Λ0 B0 ∂x 12
+τ −α
N −1 X j=1
λ (α) (α) ∂3 L tj+ 12 , xj+ 21 , Λj −bj−1 + Bj e
t 1 −tj+ 1 2
2
= 0,
and the remaining N − 1 equations by (α) ∂Ψ −α 1,x 1 , Λm 1,x 1 , Λm = 0 ⇔ ∂ L t + τ ∂ L t b0 2 3 m+ m+ m+ m+ ∂x 1 2 2 2 2 m+
+τ
2
−α
N −1 X
l=m+1
λ (α) (α) ∂3 L tl+ 12 , xl+ 12 , Λl bl−m − bl−m−1 e
tm+ 1 −tl+ 1 2
(28)
2
= 0,
(29)
for m = 1, 2, . . . , N −1. Hence, the discrete solution of the variational problem, x 12 , x 32 , . . . , xN − 21 ,
is the solution of the system of N equations (28)-(29). For the numerical approximation of the isoperimetric problem, we proceed similarly by using the same numerical approach, but where L is replaced with the function H = L + µM , where µ is a real number. Since in this case, the resulting system of the form (27) is a system of N equations on the (N + 1) unknowns x 12 , x 23 , . . . , xN − 21 , µ, an extra equation is needed. This extra equation is obtained through the discretization of (9). Using the same quadrature formulas and the same approximation for the Caputo derivatives as before, we obtain: τ
N −1 X
M (ti+ 12 , xi+ 21 , Λi ) = k.
i=0
4
Examples
Example 1. Consider the following problem: J =
Z
0
1
−t 2 2
(x(t) − e t ) +
C 0.5,1 x(t) 0 Dt
−e
−t
subject to the conditions
2 Γ
5 2
t
3 2
!2
dt → min
1 . e It is easy to see that the function x? (t) = e−t t2 satisfies the Euler-Lagrange equation (5). Also, since J ≥ 0 and J (x? ) = 0, we conclude that x? is in a fact a solution to the variational problem. We used the numerical approach that we described in the previous section, to approximate the solution of this problem for several values of N . The numerical results are presented in Table 1, where the maximum of the absolute error, E = max xi+ 12 − x? ti+ 21 x(0) = 0
and
x(1) =
i=0,...,N −1
is displayed. From the results in this table, we observe that the error decreases when N increases, as it is expected. In Figure 1, we present the plot of the analytical solution as well as the discrete one, where a good agreement between them can be observed. Next, we will consider an example, from [33], with a usual (non-tempered) Caputo derivative. Example 2. J =
Z
0
1
C 0.5 0 Dt x(t)
−
3 2 t2 Γ(2.5)
12
2
dt → min,
N E
5 1.03 ×10−2
10 4.00 ×10−3
20 1.51 ×10−3
40 5.58 ×10−4
Table 1: Maximum of the committed absolute error in the approximation of the solution for Example 1. 0.35 0.30 0.25 0.20 0.15 0.10 0.05
* * * * * * 0.00 * * * * 0.0 0.2
*
*
*
*
*
*
*
*
*
*
*
*
0.4
*
*
*
0.6
*
*
*
*
*
*
*
*
*
*
*
0.8
*
*
*
*
1.0
Figure 1: Analytical (continuous line) and discrete (* points) solution for Example 1 with N=40. subject to the conditions x(0) = 0
and x(1) = 1.
It is easy to see that the analytical solution is x(t) = t2 . In this case, we must use our approach with λ = 0. Some numerical results are listed in Table 2, where the available results from [33] are also presented. N 5 10 30
E 3.58 ×10−2 1.42 ×10−2 2.98 ×10−3
E ( [33]) 0.03 0.02 0.006
Table 2: Maximum of the committed absolute error in the approximation of the solution for Example 2. Example 3. Let us now consider the following isoperimetric problem, J =
Z
0
1
2
C 0.5,1 x(t) 0 Dt
+
2 √ 15 π 2 −t t e dt → min, 16
subject to the conditions x(0) = 0
and
x(1) =
1 e
and to the integral constraint Z 1 √ 15 π 2 −tC 0.5,1 675π 7 t e 0 Dt x(t) dt = 1− 2 . 16 1024 e 0 We can easily verify that the function x? (t) = t5/2 e−t satisfies the necessary conditions of Theorem 6, with µ = −2. We have applied the method described previously to approximate the solution of this problem. The numerical results are presented in Table 3, where the maximum of the absolute error at the time discretization points are listed as well as the obtained value for µ, for several values of N .
13
N E µ
5 2.39 ×10−2 -2.02324
10 8.76 ×10−3 -2.00574
20 3.15 ×10−3 -2.00143
40 1.12×10−3 -2.00036
Table 3: Numerical values of µ and maximum of the committed absolute error in the approximation of the solution for Example 3.
5
Conclusions
Here we have extended previous theoretical results for variational problems with Caputo fractional derivatives to the case of functionals containing tempered Caputo fractional derivatives. All the obtained results reduce to known results for the usual Caputo derivatives, by considering λ = 0. We have also presented a numerical approach which is able to deal with tempered and non-tempered fractional variational problems. Moreover, we have improved the method in [33] by considering alternative approximation formulas. Since the ones we have used here correspond to higher order approximations, it is expected that the method proposed here produces an error that goes faster to zero, as N increases, than the one presented in [33]. The results in Table 2 indicate this, although the available values of N in [33] are still relatively low to properly illustrate this. We intend to provide this comparison in a forthcoming paper where we not only illustrate numerically the conclusions in remark 1, but we will also compare this method with other methods, such as indirect methods, which solve the Euler-Lagrange equations to provide a numerical solution to the variational problems.
Acknowledgments R. Almeida is supported by Portuguese funds through the CIDMA - Center for Research and Development in Mathematics and Applications, and the Portuguese Foundation for Science and Technology (FCT-Funda¸c˜ao para a Ciˆencia e a Tecnologia), within project UID/MAT/04106/2019. M. L. Morgado acknowledges the Portuguese Foundation for Science and Technology through project UID/Multi/04621/2013. The authors are very grateful to two anonymous referees, for valuable remarks and comments that improved this paper.
References [1] A. A. Kilbas, H. M. Srivastava and J. J. Trujillo, Theory and Applications of Fractional Differential Equations. North-Holland Mathematics Studies, 204. Elsevier Science B.V., Amsterdam, 2006. [2] K. S. Miller and B. Ross, An Introduction to the Fractional Calculus and Fractional Differential Equations. A Wiley-Interscience Publication, Wiley, New York, 1993. [3] I. Podlubny, Fractional differential equations, Mathematics in Science and Engineering, 198. Academic Press, Inc., San Diego, CA, 1999. [4] S. G. Samko, A. A. Kilbas and O. I. Marichev, Fractional Integrals and Derivatives, translated from the 1987 Russian original. Gordon and Breach, Yverdon, 1993. [5] F. Sabzikar, M. M. Meerschaert and J. H. Chen, Tempered fractional calculus. J. Comput. Phys. 293 (2015) 14–28. [6] J. W. Deng, L. J. Zhao and Y. J. Wu, Fast predictor-corrector approach for the tempered fractional ordinary differential equations. Numer. Algor. 74 (3) (2017) 717–754.
14
[7] K. Diethelm, The Analysis of Fractional Differential Equations: An Application-Oriented Exposition Using Differential Operators of Caputo Type. Springer, 2010. [8] B. Baeumer and M. M. Meerschaert, Tempered stable L´evy motion and transient superdiffusion. J. Comput. Appl. Math. 233 (10) (2010) 2438–2448. [9] F. Riewe, Nonconservative Lagrangian and Hamiltonian mechanics. Phys. Rev. E 53 (1996) 1890. [10] O. P. Agrawal, Formulation of Euler–Lagrange equations for fractional variational problems. J. Math. Anal. Appl. 272 (2002) 368–379. [11] T. M. Atanackovi´c, S. Konjik and S. Pilipovi´c, Variational problems with fractional derivatives: Euler–Lagrange equations. J. Phys. A: Math. Theor. 41 (2008) 095201. [12] D. Baleanu, S. I. Muslih and E. M. Rabei, On fractional Euler–Lagrange and Hamilton equations and the fractional generalization of total time derivative. Nonlinear Dynam. 53 (1–2) (2008) 67–74. [13] M. J. Lazo and D. F. M. Torres, The Legendre condition of the fractional calculus of variations. Optimization 63 (2014) 1157–1165. [14] D. Baleanu, S. I. Muslih and K. Ta¸s, Fractional Hamiltonian analysis of higher order derivatives systems. J. Math. Phys. 47 (10) (2006) 103503. [15] M. A. E. Herzallah and D. Baleanu, Fractional Euler–Lagrange equations revisited. Nonlinear Dynam. 69 (3) (2012) 977–982. [16] D. Baleanu and S. I. Muslih, Lagrangian Formulation of Classical Fields within Riemann– Liouville Fractional Derivatives. Phys. Scr. 72 (2–3) (2005) 119–121. [17] M. Klimek, Lagrangean and Hamiltonian fractional sequential mechanics. Czechoslovak J. Phys. 52 (11) (2002) 1247–1253. [18] A. Lotfi and , S. A. Yousefi, A numerical technique for solving a class of fractional variational problem. J. Comput. Appl. Math. 237 (1) (2013) 633–643. [19] O. P. Agrawal, Fractional variational calculus and the transversality conditions. J. Phys. A: Math. Gen. 39 (2006) 10375–10384. [20] O. P. Agrawal, Generalized Euler–Lagrange equations and transversality conditions for FVPs in terms of the Caputo derivative. J. Vib. Control 13 (2007) 1217–1237. [21] R. Almeida and D. F. M. Torres, Necessary and sufficient conditions for the fractional calculus of variations with Caputo derivatives. Commun. Nonlinear Sci. Numer. Simul. 16 (2011) 1490– 1500. [22] D. Baleanu, Fractional constrained systems and Caputo derivatives. J. Comput. Nonlinear Dyn. 3 (2008) 021102. [23] M. A. E. Herzallah and D. Baleanu, Fractional-order Euler–Lagrange equations and formulation of Hamiltonian equations. Nonlinear Dynam. 58 (1) (2009) 385. [24] O. P. Agrawal, Fractional variational calculus in terms of Riesz fractional derivatives. J. Phys. A 40 (24) (2007) 6287–6303. [25] R. Almeida, Fractional variational problems with the Riesz–Caputo derivative. Appl. Math. Lett. 25 (2012) 142–148. [26] M. Klimek, Fractional sequential mechanics – models with symmetric fractional derivative. Czechoslovak J. Phys. 51 (12) (2001) 1348–1354. 15
[27] D. Tavares, R. Almeida and D. F.M. Torres, Combined fractional variational problems of variable order and some computational aspects. J. Comput. Appl. Math. 339 (2018) 374–388. [28] O. S. Fard and M. Salehi, A survey on fuzzy fractional variational problems. J. Comput. Appl. Math. 271 (2014) 71–82. [29] I. Matychyn and V. Onyshchenko, On time-optimal control of fractional-order systems. J. Comput. Appl. Math. 339 (2018) 245–257. [30] R. Almeida and M. L. Morgado, The Euler–Lagrange and Legendre equations for functionals involving distributed-order fractional derivatives. Appl. Math. Comput. 331 (2018) 394–403. [31] R. Almeida, S. Pooseh and D. F. M. Torres, Computational Methods in the Fractional Calculus of Variations. Imperial College Press, London, 2015. [32] A. B. Malinowska and D. F. M. Torres, Introduction to the Fractional Calculus of Variations. Imperial College Press, London, 2012. [33] S. Pooseh, R. Almeida and D. F. M. Torres, Discrete Direct Methods in the Fractional Calculus of Variations. Comput. Math. Appl. 66 (2013) 668–676. [34] F. Zeng and C. Li, A new Crank–Nicolson finite element method for the time-fractional subdiffusion equation. Appl. Num. Mathem. 121 (2017) 82–95.
16