Fractional Relaxation Equations and a Cauchy Formula for Repeated Integration of the Resolvent

Cauchy’s formula for repeated integration is shown to be valid for the function R(t) = λΓ(q)tEq,q(−λΓ(q)t) where λ and q are given positive constants with q ∈ (0, 1), Γ is the Gamma function, and Eq,q is a MittagLeffler function. The function R is important in the study of Volterra integral equations because it is the unique continuous solution of the so-called resolvent equation R(t) = λtq−1 − λ ∫ t 0 (t− s)q−1R(s) ds on the interval (0,∞). This solution, commonly called the resolvent, is used to derive a formula for the unique continuous solution of the Riemann-Liouville fractional relaxation equation Dx(t) = −ax(t) + g(t) (a > 0) on the interval [0,∞) when g is a given polynomial. This formula is used to solve a generalization of the equation of motion of a falling body. The last example shows that the solution of a fractional relaxation equation may be quite elementary despite the complexity of the resolvent.


Introduction
In classical physics, the ordinary differential equation x (t) = −ax(t) + g(t) (1.1) is sometimes called the relaxation equation (cf.[6, p. 138], [11]) when the constant a is positive.A generalization of (1.1) in a fractional calculus setting is the fractional relaxation equation where D q denotes a fractional differential operator of order q with q ∈ (0, 1) (cf.[6, p. 138]; [12, p. 292]; [17, p. 224]).This paper is a study of (1.2) when g(t) is a polynomial.For given a > 0 and g, we will prove that this equation has a unique continuous solution on the half-closed interval [0, ∞) and that necessarily x(0) = 0. Furthermore, in Section 7, we will derive a formula that expresses this solution as a sum involving twoparameter Mittag-Leffler functions (cf.(7.13)).Moreover, we will show that each term of this sum can also be expressed as a convolution integral involving the solution of the integral equation where λ is a positive constant related to the value of the constant a.In fact, it is well-established that (R λ ) has a unique continuous solution on the interval (0, ∞) whenever λ and q are positive constants with q ∈ (0, 1).A proof of this for a more general version of equation (R λ ) can be found in the 1971 monograph by Miller [14,Ch. IV].
There is also the recent paper [3] that investigates (R λ ) directly.Not only is the existence and uniqueness of a continuous solution of (R λ ) on (0, ∞) proven there but also a formula for it is derived, namely R(t) = λΓ(q)t q−1 E q,q (−λΓ(q)t q ), ( where E α,β (α, β > 0) denotes the two-parameter Mittag-Leffler function: We use the following established terminology: resolvent equation refers to equation (R λ ) and resolvent refers to its solution (1.3).

Riemann-Liouville Operators
For a function f that is (Riemann) integrable, we employ the integral operator J defined by Furthermore, for n ∈ N (set of natural numbers), let the operator J n denote the n-fold iterate of J; that is, where J 0 := I, the identity operator.For example, taking n = 2 and applying J 2 to an integrable function f , we have In general, This particular iterated integral can be expressed in terms of a single integral with a weighted integrand as in (2.2) below.It is known as Cauchy's formula for repeated integration (cf.[16, p. 38]).This formula is found in Abramowitz and Stegun's handbook [1, (25.4.58)]; and in some textbooks, such as [8, p. 487], it appears as an exercise.We omit its proof here because it is basically the mathematical induction argument that is used in the proof of Theorem 3.2 in Section 3.
Theorem 2.1 (Cauchy's formula for repeated integration).Let n ∈ N. If f is integrable on [0, T ], then Now let us extend the values of n in (2.2) from N to R + (set of strictly positive real numbers) by replacing (n − 1)! with Γ(n), where Γ denotes the Gamma function.This leads to the well-known definition of the Riemann-Liouville integral operator of order n.Definition 2.2.For any n ∈ R + , J n denotes the integral operator where f denotes a function for which the integral exists.J n is called the Riemann-Liouville fractional integral operator of order n.Furthermore, J and J 0 denote the operators where I denotes the identity operator.
Just as the integral operator J n can be defined for all values of n ∈ R + , the same is true of D n , namely, the classical ordinary differential operator of order n ∈ N.That is, for an n-times differentiable function f , This can be expressed recursively as follows: where D 1 := D and D 0 := I, the identity operator.
In the following extension of the definition of D n , we employ the floor function • , where n denotes the largest integer less than or equal to n. Definition 2.3.For a given n ∈ R + , D n denotes the differential operator where m = n + 1 and f denotes a function for which the right-hand side exists.For n = 0, D n f := f .D n is called the Riemann-Liouville fractional differential operator of order n (cf.[6, p. 27]).
Remark 2.4.The symbol D n on the left-hand side of (2.5) denotes the fractional differential operator of order n whereas D m on the right-hand side denotes the ordinary differential operator Thus the definition of the operator D n is well-defined.
We now use this theorem to show that Cauchy's formula for repeated integration can be applied to the resolvent R(t), notwithstanding the singularity at t = 0. Theorem 3.2.For n ∈ N, for t ≥ 0.
Proof.For a given t ≥ 0, it is well-known (cf.[14, Ch.IV]) that the resolvent integral function exists.Furthermore, it is proven in [3, Thm.9.5] that where E α , for α ∈ R + , denotes the one-parameter Mittag-Leffler function, which is defined by Since E q (z) is an entire function of z (cf.[6, Thm.4.1]) in the complex plane, it follows that exists and is continuous on [0, ∞) for each n ∈ N.
Note that (3.1) simplifies to (3.2) when n = 1.We complete the proof using mathematical induction to establish that (3.1) is true for all n ∈ N. Suppose that (3.2) is also true when n = k for some k ∈ N. Then Interchanging the order of integration and appealing to Theorem 3.1, we have which is precisely (3.1) when n = k + 1.Thus, as (3.1) is true for n = 1, it must be true for all n ∈ N.
The following result will be needed in the next section to prove Lemma 4.2.Although the proof is straightforward, it can be found in a number of places (e.g., [6, p. 28]).Lemma 3.3.Let q ∈ (0, 1) and p > −1.If p = q − 1, then for t > 0. If p = q − 1, then D q t p = 0 for t > 0.

Solution of a fractional relaxation equation
The following proof is an adaptation of a proof in [6, Thm.2.14].
If a function f is continuous and absolutely integrable on an interval (0, T ], then for all t ∈ (0, T ].
Proof.This is trivially true for n = 0 since by definition J 0 := I and D 0 := I.It is also true for n = 1 because by the Fundamental Theorem of Calculus for 0 < t ≤ T .It follows from this and an induction argument that (4.1) is true for all n ∈ N 0 , where N 0 := N ∪ {0}.Now consider (4.1) for a given n > 0 when it is not a positive integer.Then, by Definition 2.3, where m = n + 1.Since f by hypothesis is continuous and absolutely integrable on (0, T ] and m + n ≥ 1, we have The following result relates solutions of (1.2) to those of a Volterra integral equation when g(t) ≡ b, a constant.It will be extended to all polynomials in Section 7.
Lemma 4.2.Let a, b, and q be constants with a > 0, b ∈ R, and q ∈ (0, 1).If there is a continuous solution of the fractional relaxation equation on the interval [0, ∞), then it is also a solution of on [0, ∞) when β and λ have the values and λ = a Γ(q) .(4.4) Conversely, let β ∈ R and λ > 0 and suppose there is a continuous solution of the integral equation Then it is also a continuous solution of (4.2) on [0, ∞) when a = λΓ(q) and b = βΓ(q + 1).(4.5) Proof.Let β ∈ R and λ > 0 be given constants.Then let a and b be defined by (4.5).Suppose there is a continuous function x(t) that satisfies the integral equation ( 4.3) on [0, ∞).Expressing this in terms of the Riemann-Liouville integral operator (2.3), we obtain Applying the Riemann-Liouville differential operator D q to this and using Theorem 4.1, we get since D q t q = Γ(1 + q) (cf.Lemma 3.3).In other words, the function x(t) must also be a solution of (4.2) on [0, ∞).Note from (4.3) that x(0) = 0. Now let a > 0 and b ∈ R be given constants.Then define constants β ∈ R and λ > 0 by (4.4) and suppose x(t) is a continuous function satisfying (4.2) on [0, ∞).And so Because of this and the continuity of x on [0, ∞), we obtain x(s) ds + bt = −aJx(t) + bt upon taking the limit of both sides of (4.6) as η → 0 + .Then the application of D 1−q yields D 1−q J 1−q x(t) = −aD 1−q Jx(t) + bD 1−q t, which because of Theorem 4.1 and (2.5) simplifies to It follows from [4, Lemma 4.8] that and from Lemma 3.3 that Therefore, we conclude from (4.8) that any continuous solution of (4.2) on [0, ∞) must also be a solution of Moreover, we see from this integral equation that x(0) = 0. Remark 4.3.Observe in the statement of Lemma 4.2 that no initial condition accompanies the fractional differential equation (4.2).At first this may appear to be an oversight until we realize from the proof that positing the existence of a continuous solution x(t) of (4.2) for t ≥ 0 implies x(0) = 0.
With the next theorem we complete what was initiated with Lemma 4.2 and that is to show that (4.2) and (4.3) do in fact share the same continuous solution on [0, ∞).But first let us dispose of the special case b = 0. Lemma 4.4.There is one and only one continuous solution of on [0, ∞); it is the trivial solution x(t) ≡ 0.
Proof.It follows from Lemma 4.2 that any continuous solution of (4.9) on [0, ∞) must also be a continuous solution of where λ = a/Γ(q).But the only solution of this integral equation is x(t) ≡ 0 (cf.[3, p. 15]).
Proof.First consider the integral equation (4.3)where β ∈ R and λ > 0 are given constants.From [3,Thm. 8.3] we know that if a function f is continuous on the interval [0, ∞), then has a unique continuous solution on [0, ∞).Moreover, this solution is given by the linear variation of parameters formula: Taking f (t) = βt q , this becomes In other words, this is the unique continuous solution of (4.3) on [0, ∞).But we can simplify (4.11) as follows: integrating the resolvent equation (R λ ) and interchanging the order of integration (cf.Thm.3.1), we obtain Substituting this into (4.11) and defining a and b by (4.5), we get In [3, Thm.9.5], we find the formula Therefore (4.10) is the unique continuous solution of (4.3).Moreover, Lemma 4.2 implies that it is also the unique continuous solution of (4.2) on [0, ∞).
Remark 4.6.We have shown that there is one and only one continuous solution x(t) of the fractional relaxation equation (4.2) on the half-closed interval [0, ∞).Moreover, from (4.10) we see that x(0) = 0. Thus, the initial value problem D q x(t) = −ax(t) + b, x(0) = x 0 has no continuous solution on [0, ∞) unless x 0 = 0. Also, observe that if we formally let q = 1 in (4.10), then it simplifies to Proof.Properties (i) and (ii) follow from (4.10) from which we see that x(0) = 0 and Since the derivative of (4.10) is it follows that x (t) > 0 if b > 0. And so x(t) is strictly increasing on [0, ∞).This together with x(0) = 0 implies that x(t) > 0 for t > 0. Likewise, if b < 0, then x(t) is strictly decreasing on [0, ∞) and x(t) < 0 for t > 0. This concludes the proof of (iii).
To prove (iv), we use the result stated in Section 1 that the resolvent R is a completely monotone function on (0, ∞).Thus, And so Finally, (v) follows from (4.13) and the complete monotonicity of R.
Example 4.8.We illustrate some of the properties of solutions of (4.2) that are enumerated in Corollary 4.7 by choosing two different values of q and graphing the corresponding solutions (4.10).For both values, let a = b = 1.
First let q = 1/2.Then (4.10) is According to [10, (1.8.6)], where erf(z) is the error function: Thus, is the unique continuous solution of on [0, ∞).The graph of (4.15) is the solid concave-downward curve in Figure 1.(All the graphs in this paper were created with Maple TM 17.) Now let q = 1/3.By (4.10) the unique continuous solution of The second term can be calculated with the help of formula (1.8.5) in [10]: for m = 2, 3, 4, . . ., Consequently, Changing the variable of integration to z = u 3 , we obtain t 0 e −u 3 du = 1 3 Likewise, the same change of variable yields t 0 u e −u 3 du = 1 3 Thus, Therefore, an alternative form of (4.18) is

.21)
Its graph is the concave-downward dashed curve in Figure 1.

Generalization of the equation of motion of a falling body
Consider the vertical downward motion of a body of mass m when it is released from rest above the ground.Besides the force due to gravity and the resisting force (drag force) exerted on the body as it moves through the air, assume that all other forces acting on the body are negligible.Also, let us suppose that the drag force is proportional to the velocity v of the body.Since the direction of the drag force is opposite that of the velocity of the body (downward, which we take as negative), the drag force F d points upward.Thus F d = −kmv, where the proportionality constant k > 0. Since the gravitational force acting on the body is F g = −mg, Newton's second law of motion yields Hence, we obtain the familiar classical equation of motion that is found in most undergraduate physics textbooks, such as [13, p. 68].Solving (5.1) by either separating variables or using the integrating factor e kt , we obtain the solution Now suppose we generalize the equation of motion (5.1) by replacing the classical differential operator d/dt with the Riemann-Liouville operator D q .But this does not make complete sense due to the dimensional inconsistency of the units, where we see from (5.1) that k has the dimension of inverse time.However, we can rectify this with the replacement d dt → k 1−q D q suggested by Rosales et al. in [18, p. 519].(Actually they use the Caputo fractional derivative; however, since q ∈ (0, 1) and the initial condition is v(0) = 0, the Caputo and Riemann-Liouville derivatives are equivalent.) Consequently, the fractional generalization of (5.1) is From Theorem 4.5 we see that the unique continuous solution of (5.3) is This agrees with the velocity formula in [18].Note that by formally letting q = 1, (5.4) simplifies to (5.2) because of (3.3).
Proof.Integrating (R λ ) and interchanging the order of integration (cf.Thm.3.1), we get Theorem 6.2.Let n ∈ N. The nth repeated integral of the resolvent R(t), namely is given by the formula for t ≥ 0.
Proof.It follows from Lemma 6.1 that formula (6.2) holds for n = 1.Let us show via a proof by induction that it holds for all n ∈ N. Suppose for some k ∈ N that (6.2) holds for n = k.Then Interchanging the order of integration as in Theorem 3.1, we obtain This shows that (6.2) holds for n = k + 1 if it holds for n = k.Therefore, by induction, (6.2) holds for all n ∈ N.

Solution of (1.2) for a given polynomial g
The first result of this section generalizes Lemma 4.2.
Lemma 7.1.Let a > 0, b ∈ R, m ∈ N 0 , and q ∈ (0, 1).If there is a continuous solution of on [0, ∞), then it is also a solution of and λ = a Γ(q) .(7.3) Conversely, let β ∈ R and λ > 0 and suppose there is a continuous solution of (7.2) on [0, ∞).Then it is also a continuous solution of (7.1) on [0, ∞) when a = λΓ(q) and b = β m! Γ(q + m + 1).(7.4) Proof.Equation (7.2), written in terms of the Riemann-Liouville integral operator, is Applying the differential operator D q , we get D q x(t) = βD q t q+m − λΓ(q)D q J q x(t) = β Γ(q + m + 1) Γ(m + 1) where we have used Theorem 4.1, Lemma 3.3, and (7.4).Thus, if a continuous solution of (7.2) exists for t ≥ 0, it must also be a solution of (7.1) when a and b have the values given by (7.4).Conversely, suppose there exists a continuous solution x(t) of (7.1) on [0, ∞); hence Integrating, as in the proof of Lemma 4.2, we have Taking the limit of both sides as η → 0 + , we obtain since J 1−q x(η) → 0 (cf.(4.7)).Applying D 1−q , we get t m+q because of Lemma 3.3 and Theorem 4.1.Since m ∈ N 0 , this simplifies to We conclude that if a continuous solution x(t) of (7.1) exists, then it must also be a solution of In the next theorem we prove that (7.1) does have a unique continuous solution on [0, ∞).Moreover, with the following formula, which is found in [17, p. 25], we show how to express it in terms of a Mittag-Leffler function.
Proof.Let us use (1.4) to write the integrand as the sum where Using an integration formula in [4, (4.4)] (or [6, p. 229 ]), we find for t > 0 that It then follows from a generalization of Levi's theorem for series ([2, p. 269]) that With the integration formulas involving the resolvent that we found in Section 6.1 and the variation of parameters formula that was used earlier in the proof of Theorem 4.5, we can now establish the existence of continuous solutions of equations (7.1) and (7.2) on [0, ∞) and their uniqueness.Theorem 7.3.Let a > 0, b ∈ R, m ∈ N 0 , and q ∈ (0, 1).The fractional relaxation equation (7.1) has the unique continuous solution on [0, ∞), where R denotes the resolvent corresponding to λ = a/Γ(q).It is also the unique continuous solution of the integral equation (7.2) on [0, ∞) when β and λ have the values given by (7.3).
Proof.First consider the integral equation (7.2) for given values of β ∈ R and λ > 0. By the variation of parameters formula, the function is the unique continuous solution of (7.2) on [0, ∞).Because of Corollary 6.3, this solution can be simplified as follows: first define a and b by (7.4).Then In short, we have shown that (7.8) is the unique continuous solution of (7.2) on [0, ∞) if a and b have the values given by (7.4).Furthermore, we can see from Lemma 7.1 that it is also the unique continuous solution of the fractional differential equation (7.1) on [0, ∞).
Setting n = 2 + m and x = −at, we obtain Thus, Writing out some of the terms of (7.9), we have Now note that when we evaluate this at t = 0, we get x(0) = (−1) m bm! a m+1 + (−1) m+1 bm! a m+1 = 0, which agrees with the value of (7.7) at t = 0. Also note that the fractional differential operator D q is defined to be the ordinary first-order differential operator D when q = 1 (cf.Def.2.3).So the formal substitution q = 1 in (7.7) suggests that (7.9) is the solution of the classical initial value problem Let us see if this is truly the case.By the classical variation of parameters formula, the solution of (7. With an appropriate change in the index of summation, we see that this is equivalent to (7.9).In sum, we have shown that (7.9) does in fact solve the classical initial value problem (7.10).
Because of Corollary 6.3 and Theorem 7.3, we can express the convolution of t p for p > −1 and the resolvent R in terms of two-parameter Mittag-Leffler functions.This is the content of the next result.
Our final result employs Theorem 7.3 to obtain the unique continuous solution of on the interval [0, ∞) when g(t) is a given polynomial.
Example 7.8.The equation has the unique continuous solution on the interval [0, ∞).
Each of these four terms can be expressed as a finite sum of powers of t and a constant multiple of e t erf( √ t).For instance, consider the second term.From (1.4) we have . (7.18) By the Cauchy-Hadamard formula for convergence and Stirling's formula for the Gamma function, the power series (1.4) defining E α,β (z) converges absolutely for all z in the complex plane (cf.[6, p. 68]), and a fortiori for all real values of z.Consequently, we can rearrange the terms of (7.18) as follows: .
It then follows that .
Changing the index of summation, we have As in Remark 7.5, let us compare the solution (7.16) of the fractional relaxation equation (7.15) with the solution of the initial value problem y (t) = −y(t) + 1 − 3t − 2t 2 + t 3 , y(0) = 0. (7.23) Applying the variation of parameters formula or simply multiplying the differential equation by the integrating factor e t and then integrating by parts and using the initial condition y(0) = 0, we obtain the solution y(t) = −6 + 7t − 5t 2 + t 3 + 6e −t .(7.24) The graph of the solution y(t) (dashed curve) is shown in Figure 2. The solid curve is the graph of the solution (7.16) of the fractional relaxation equation (7.15).The curve that begins at (0, 1) (dotted curve ) is the graph of the polynomial (7.17).

0 e
10) is x(t) = e −at x(0) + t −a(t−s) bs m ds, (7.11) which simplifies to x(t) = be −at t 0 s m e as ds since x(0) = 0. Integrating by parts or consulting a table of integrals, such as [9, 2.321], we find that s m e as ds = e as m k=0 (−1) k k! a k+1 m k s m−k .Hence, x(t) = be −at t 0 s m e as ds = b m k=0