FRACTIONAL OSCILLATIONS AND MITTAG-LEFFLER FUNCTIONS

Rudolf GORENFLO^{1} and Francesco MAINARDI^{2}

1 Department of Mathematics and Computer Science, Berlin Free University, D-14195 Berlin, Germany e-mail: gorenflo@math.fu-berlin.de

2 Department of Physics, University of Bologna, I-40126 Bologna, Italy e-mail: mainardi@bologna.infn.it

Abstract

The fractional oscillation equation is obtained from the classical equation for linear oscillations by replacing the second-order time derivative by a fractional derivative of order α with 1 < α < 2 . Using the method of the Laplace transform, it is shown that the fundamental solutions can be expressed in terms of Mittag-Leffler functions, and exhibit a finite number of damped oscillations with an algebraic decay. For completeness we also discuss both the cases 0 < α < 1 (fractional relaxation) and 2 < α ≤ 3 (growing oscillations), showing the key role of the Mittag-Leffler functions.

1991 Mathematics Subject Classification: 26A33, 33E20, 33E30, 45E10, 45J05, 70J99. Key Words: Fractional differential equations, fractional calculus, Mittag-Leffler functions.

1. Introduction

For real α > 0 (later only for 0 < α ≤ 3) we consider the fractional differential equation

D^{α} (u(t) − ∑_{k=0}^{m−1} \frac{t^{k}}{k!} u^{(k)}(0^{+})) = −u(t) + q(t) , t > 0 ,  (1.1)

where q(t) is a given continuous function. Here m is an integer uniquely defined by m − 1 < α ≤ m , which provides the number of the prescribed initial values u^{(k)}(0^{+}) = c_{k} , k = 0, 1, 2, …, m − 1 . Implicit in the form of (1.1) is our desire to obtain solutions u(t) for which the u^{(k)}(t) are continuous for t ≥ 0, k = 0, 1, …, m−1.

The purpose of this paper is to evaluate the solutions of (1.1) in analytical form, providing a detailed discussion of the main properties of interest in mechanics. In our analysis we point out the key role of the Mittag-Leffler functions, improving the results recently obtained by Gorenflo and Rutman [1] and by Mainardi [2], who had considered some mathematical aspects of (1.1).

A numerical treatment of (1.1) has been recently given by Blank [3], who has provided reliable results by means of the collocation approximation with polynomial splines. The plots of our analytical solutions confirm the numerical results in [3].

Several applications of fractional differential equations in physics and engineering have also appeared, see e.g. [4-15].

Before starting our analysis of (1.1), we have to explain the meaning of the operator D^{α} of fractional differentiation. Hereafter, for the sake of completeness, we provide the essentials of Fractional Calculus according to the Riemann-Liouville approach. Let us first define the operator of fractional integration J^{β} , β > 0 ,

J^{β} u(t) := \frac{1}{Γ(β)} ∫_{0}^{t}(t − τ)^{β−1} u(τ) dτ , t > 0 .  (1.2)

Then, taking the positive integer m with m − 1 < α ≤ m , we define, with D = \frac{d}{dt} ,

D^{α} u(t) := D^{m} J^{m−α} u(t) , t > 0 .  (1.3)

We also define J^{0} = D^{0} = I , where I denotes the Identity operator. We refer to the recent books on Fractional Calculus [16-18] for general treatment and discussion of the operators of fractional integration and derivation. Here we limit ourselves to report the simple, but relevant results for the power function, with t > 0 ,

J^{p} t^{q} = \frac{Γ(q + 1)}{Γ(q + p + 1)} t^{q+p} , D^{p} t^{q} = \frac{Γ(q + 1)}{Γ(q − p + 1)} t^{q−p} , p ≥ 0 , q > −1 . (1.4)

We also observe that an alternative definition of fractional derivative, which is due to Caputo [4], is

D̃^{α} u(t) := J^{m−α} D^{m} u(t) .  (1.5)

We note in general that

D^{m} J^{m−α} u(t) ≠ J^{m−α} D^{m} u(t) .  (1.6)

In fact, assuming that the passage of the m-derivative under the integral is legitimate, one recognizes that, for m − 1 < α < m , and t > 0

D^{α} u(t) = D̃^{α} u(t) + ∑_{k=0}^{m−1} \frac{t^{k−α}}{Γ(k − α + 1)} u^{(k)}(0^{+}) ,  (1.7)

and therefore, recalling the fractional derivative of the power function (1.4),

D^{α} (u(t) − ∑_{k=0}^{m−1} \frac{t^{k}}{k!} u^{(k)}(0^{+})) = D̃^{α} u(t) .  (1.8)

We also note that while

α → m^{∓} , D^{α} u(t) → D^{m} u(t) ;  (1.9a)

we have

\{ \begin{matrix}α & → m^{−} , D̃^{α} u(t) → D^{m} u(t) ; \\ α & → m^{+} , D̃^{α} u(t) → D^{m} u(t) − u^{(m)}(0^{+}) . \end{matrix}  (1.9b)

Therefore in (1.1) we have substantially adopted Caputo’s definition for the fractional derivative, which incorporates the initial values of the function and of its integer derivatives of lower order. This approach appears suitable in the applications, where the initial conditions are explicitly expressed in terms of the conventional integer derivatives. Observe, however, that in mathematical treatments of fractional differential equations, the initial conditions are usually given as limiting values of fractional integrals; see e.g. [17-18].

The equation (1.1) with the given initial conditions can be solved by elementary functions if α is an integer. For completeness of presentation we write down the corresponding solutions. In the case α = 1 we have the equation of relaxation

u′(t) = −u(t) + q(t) ,  (1.10)

whose solution, under the initial condition u(0) = c_{0} , is

u(t) = c_{0} e−t + ∫_{0}^{t}q(t − τ) e−τ dτ .  (1.11)

In the case α = 2 we have the equation of oscillation

u′′(t) = −u(t) + q(t) ,  (1.12)

whose solution, under the initial conditions u(0) = c_{0} and u′(0) = c_{1} , is

u(t) = c_{0} cos t + c_{1} sin t + ∫_{0}^{t}q(t − τ) sin τ dτ .  (1.13)

In the case α = 3 we have the equation

u′′′(t) = −u(t) + q(t) ,  (1.14)

whose solution, under the initial conditions u(0) = c_{0} , u′(0) = c_{1} , u′′(0) = c_{2} , is

u(t) = c_{0} u_{0}(t) + c_{1} u_{1}(t) + c_{2} u_{2}(t) − ∫_{0}^{t}q(t − τ) u′_{0}(τ) dτ .  (1.15)

Here

\{ \begin{matrix}u_{0}(t) & = \frac{1}{3} e −t + \frac{2}{3} e t/2 cos \frac{\sqrt{3}}{2} t , \\ u_{1}(t) & = J u_{0}(t) = −\frac{1}{3} e −t + \frac{2}{3} e t/2 cos (\frac{\sqrt{3}}{2} t − \frac{π}{3} ) , \\ u_{2}(t) & = J u_{1}(t) = J^{2} u_{0}(t) = \frac{1}{3} e −t + \frac{2}{3} e t/2 cos (\frac{\sqrt{3}}{2} t − 2\frac{π}{3} ) , \end{matrix}  (1.16)

with the operator J = J^{1} defined according to J u(t) = ∫_{0}^{t} u(τ) dτ .

We note that in the cases α = 1 and α = 2 we can use a notation in analogy with the case α = 3 , namely we write u_{0}(t) = e^{−t} and u_{0}(t) = cos t , u_{1}(t) = J u_{0}(t) = sin t = cos (t − π/2) , respectively.

Furthermore, we note that as for α = 3 [see (1.15)], also for α = 1 and α = 2 [see (1.11) and (1.13), respectively], the free term q(t) appears convoluted with the function

u_{δ}(t) := −u′_{0}(t) .  (1.17)

This function represents the so called impulse response solution for the corresponding equation [(1.10), (1.12), (1.14)], namely the particular solution of the inhomogeneous equation with all c_{k} ≡ 0 , k = 0, 1, …, m − 1 , m = 1, 2, 3 and with q(t) = δ(t) . Here δ(t) denotes the Dirac delta distribution.

R e m a r k : The more general equation

D^{α} (u(t) − ∑_{k=0}^{m−1} \frac{t^{k}}{k!} u^{(k)}(0^{+})) = −ρ^{α} u(t) + q(t) , ρ > 0 , t > 0 ,  (1.18)

can be reduced to (1.1) by a change of scale t → t/ρ . We prefer, for ease of notation, to discuss the ”dimensionless” form (1.1).

2. Formal solutions of the fractional differential equation

We prefer to write (1.1) with the prescribed initial conditions as an equivalent (fractional) integral equation and then treat the integral equation by the Laplace transform method (see e.g. [19] for this method). For this purpose we note that

J^{α} D̃^{α} u(t) = J^{α} J^{m−α} D^{m} u(t) = J^{m} D^{m} u(t) = u(t) − ∑_{k=0}^{m−1} c_{k} \frac{t^{k}}{k!} .  (2.1)

Applying the operator J^{α} to both sides of (1.1) we obtain

u(t) = ∑_{k=0}^{m−1} c_{k} \frac{t^{k}}{k!} − J^{α} u(t) + J^{α} q(t) ,  (2.2)

Using the sign ÷ for the juxtaposition of a function depending on t with its Laplace transform depending on s:

u(t) ÷ U(s) , q(t) ÷ Q(s) ,

we have

U(s) = ∑_{k=0}^{m−1} \frac{c_{k}}{s^{k+1}} − \frac{1}{s^{α}} U(s) + \frac{1}{s^{α}} Q(s) ,

hence

U(s) = ∑_{k=0}^{m−1} c_{k} \frac{s^{α−k−1}}{s^{α} + 1} + \frac{1}{s^{α} + 1} Q(s) .  (2.3)

Introducing the functions

e_{α}(t) ÷ \frac{s^{α−1}}{s^{α} + 1} ,  (2.4)

and

u_{k}(t) := J^{k}e_{α}(t) ÷ \frac{s^{α−k−1}}{s^{α} + 1} ,  (2.5)

we find, from inversion of the Laplace transforms in (2.3),

u(t) = ∑_{k=0}^{m−1} c_{k} u_{k}(t) − ∫_{0}^{t}q(t − τ) u′_{0}(τ) dτ .  (2.6)

For finding the last term in the R.H.S of (2.6) we have used the well-known rule for the Laplace transform of the derivative noting that u_{0}(0) = e_{α}(0) = 1 and

\frac{1}{s^{α} + 1} = − (s \frac{s^{α−1}}{s^{α} + 1} − 1) ÷ −u′_{0}(t) .  (2.7)

The formula (2.6) encompasses the solutions (1.11), (1.13) and (1.15) found for α = 1 , 2 , 3 , respectively.

When α is not integer, namely for m − 1 < α < m , we note that m − 1 represents the integer part of α (usually denoted by [α]) and m the number of initial conditions necessary and sufficient to ensure the uniqueness of the solution u(t). The m functions u_{k}(t) = J^{k}e_{α}(t) with k = 0, 1, …, m − 1 represent those particular solutions which satisfy the following m initial conditions

u_{k}^{(h)}(0^{+}) = δ_{h k} , h = 0, 1, …, m − 1 ,  (2.8)

and therefore they are referred to as the fundamental solutions of the fractional equation (1.1). In the next section we are going to compute the fundamental solutions for 0 < α < 3 (α ≠ 1, 2 ). For completeness of presentation we also exhibit the impulse response solution, u_{δ}(t) .

R e m a r k : The reader is invited to verify that the solution (2.6) has continuous derivatives u^{(k)}(t) for k = 0, 1, 2, …, m − 1 , which fulfill the m initial conditions u^{(k)}(0^{+}) = c_{k} . In fact, looking back at (2.2), one must recognize the smoothing power of the operator J^{α} . However, the so called impulse response solution of our equation (1.1), u_{δ}(t) , is expected to be not so regular like the ordinary solution (2.6). In fact, from (2.3) and (2.6-7), one obtains

u_{δ}(t) = −u′_{0}(t) ÷ \frac{1}{s^{α} + 1} ,  (2.9)

and therefore, using the limiting theorem for Laplace transforms, one can recognize that, being m − 1 < α < m ,

u_{δ}^{(k)}(0^{+}) = 0 , k = 0, 1, …, m − 2 ; u_{δ}^{(m−1)}(0^{+}) = ∞ .  (2.10)

3. Explicit solutions of the fractional differential equation

In order to obtain the inversion of the Laplace transform in (2.4) we need to recall the so-called Mittag-Leffler function E_{α}(z) defined in the whole complex plane as

E_{α}(z) := ∑_{n=0}^{∞} \frac{z^{n}}{Γ(αn + 1)} , z ∈ C ,  (3.1)

This function, introduced by Mittag-Leffler at the beginning of this century, turns out to be an entire function of order ρ = 1/α , on which the reader may get information from Erdélyi [20] and Dzherbashyan [21].

The Mittag-Leffler function is connected to the Laplace integral through the equation

∫_{0}^{∞} e^{−x} E_{α} (x^{α} z) dx = \frac{1}{1 − z} , α > 0 .  (3.2)

The integral was evaluated by Mittag-Leffler who showed that the region of convergence of (3.2) contains the unit circle and is bounded by the line Re z^{1/α} = 1. Putting in (3.2) x = st and x^{α} z = λt^{α} with t ≥ 0 , and λ ∈ C we get the following Laplace transform pair

E_{α} (λ t^{α}) ÷ \frac{s^{α−1}}{s^{α} − λ} .  (3.3)

As a consequence we obtain from (2.4) and (3.3)

e_{α}(t) = E_{α}(−t^{α}) .  (3.4)

Formally we can obtain the function e_{α}(t) in form of a series in powers of t^{α} by term-wise Laplace inversion of the series

\frac{s^{α−1}}{s^{α} + 1} = ∑_{n=0}^{∞}(−1)^{n} \frac{1}{s^{αn+1}} , |s^{α}| > 1 ,

as

e_{α}(t) = ∑_{n=0}^{∞}(−1)^{n} \frac{t^{αn}}{Γ(αn + 1)} = E_{α}(−t^{α}) , t ≥ 0 .  (3.5)

We here prefer to derive the relevant properties of our functions e_{α}(t), directly from their representation as a Laplace inversion integral

e_{α}(t) = \frac{1}{2πi} ∫_{Br}e st \frac{s^{α−1}}{s^{α} + 1} ds ,  (3.6)

in detail for 0 < α ≤ 3 , without detour on the general theory of Mittag-Leffler functions in the complex plane. In (3.6) Br denotes the Bromwich path, i.e. a line Re s = σ with a value σ ≥ 1.

The essential step consists in decomposing e_{α}(t) into two parts according to e_{α}(t) = f_{α}(t) + g_{α}(t) , as indicated below. The part f_{α}, present for α not integer, turns out to be a completely monotonic function, which tends to zero as t tends to infinity; the other part turns out to be of oscillatory character if 1 < α < 3 , but identically vanishing if 0 < α < 1 . For transparency reasons, we separately discuss the cases:

(a) 0 < α < 1 ,  (b) 1 < α < 2 ,  (c) 2 < α < 3 ,  (3.7)

recalling that for α = 1 , 2 , 3 we know e_{α}(t) as elementary function, namely

e_{1}(t) = e −t , e_{2}(t) = cos t , e_{3}(t) = \frac{1}{3} e −t + \frac{2}{3} e t/2 cos \frac{\sqrt{3}}{2} t .  (3.8)

For α not integer the power function s^{α} is uniquely defined as s^{α} = |s|^{α} e^{i arg s} , with −π < arg s < π , that is in the complex s-plane cut along the negative real axis.

In order to obtain the desired decomposition of e_{α} we bend the Bromwich path of integration Br into the equivalent Hankel path Ha(1^{+}), a loop which starts from −∞ along the lower side of the negative real axis, encircles the circular disc |s| = 1 in the positive sense and ends at −∞ along the upper side of the negative real axis. One obtains in all cases (a), (b), (c),

e_{α}(t) = f_{α}(t) + g_{α}(t) , t ≥ 0 ,  (3.9)

with

f_{α}(t) := \frac{1}{2πi} ∫_{Ha(ϵ)}e st \frac{s^{α−1}}{s^{α} + 1} ds ,  (3.10)

where now the Hankel path Ha(ϵ) denotes a loop constituted by a small circle |s| = ϵ with ϵ → 0 and by the two borders of the cut negative real axis, and

g_{α}(t) := ∑_{h} e s′h t Res [\frac{s^{α−1}}{s^{α} + 1}]_{s′_{h}} ,  (3.11)

where s′_{h} are the relevant poles of s^{α−1}/(s^{α} + 1). In fact the poles turn out to be s_{h} = exp [i(2h + 1)π/α] with unitary modulus; they are all simple but relevant are only those situated in the main Riemann sheet, i.e. the poles s′_{h} for which −π < arg s′_{h} < π .

If 0 < α < 1 , there are no such poles, since for all integer h we have |arg s_{h}| = |2h + 1| π/α > π; as a consequence,

g_{α}(t) ≡ 0 , hence e_{α}(t) = f_{α}(t) , if 0 < α < 1 .  (3.12)

If 1 < α < 3 , then there exist precisely two relevant poles, namely s′_{0} = exp(iπ/α) and s′_{−1} = exp(−iπ/α) = \overline{s_{0}}′ , and one obtains

g_{α}(t) = \frac{2}{α} et cos (π/α) cos [t sin (\frac{π}{α})] , if 1 < α < 3 .  (3.13)

We note that this function exhibits oscillations with circular frequency ω(α) = sin (π/α) and with an amplitude which exponentially decays or grows correspondingly to the case (b) (1 < α < 2) or (c) (2 < α < 3) with the rate λ(α) = | cos (π/α)| . Of course, in the classical case α = 2 , we recover the sinusoidal behaviour with unitary frequency, ω(2) = 1 .

R e m a r k : If and only if α is an odd integer there is a pole, namely −1, on the negative real axis on which the cut is no longer present.

It is now an exercise in complex analysis to show that the contribution from the Hankel path Ha(ϵ) as ϵ → 0 is provided by

f_{α}(t) := ∫_{0}^{∞}e −rt K_{α}(r) dr ,  (3.14)

with

K_{α}(r) = \frac{1}{π} \frac{r^{α−1} sin (απ)}{r^{2α} + 2 r^{α} cos (απ) + 1} .  (3.15)

This function K_{α}(r) vanishes identically if α is an integer, it is positive for all r if 0 < α < 1 or 2 < α < 3 , negative for all r if 1 < α < 2 . In fact in (3.15) the denominator is, for α not integer, always positive being > (r^{α} − 1)^{2} ≥ 0 . Hence f_{α}(t) has the aforementioned property of complete monotonicity: it is decreasing towards zero in cases (a) and (c), increasing towards zero in case (b). We also note that, in order to satisfy the initial condition e_{α}(0^{+}) = 1 , we find ∫_{0}^{∞} K_{α}(r) dr = 1 if 0 < α < 1 , ∫_{0}^{∞} K_{α}(r) dr = 1 − 2/α if 1 < α < 3 .

In order to compile all the fundamental solutions for 0 < α < 3 with α ≠ 1, 2 [cases (a), (b), (c)], we also need to compute the associated functions J^{k}e_{α}(t) for k = 1, 2 , precisely J^{1} e_{α}(t) if 1 < α < 3 , and J^{2} e_{α}(t) if 2 < α < 3 . This can be done in direct analogy to the computation of the functions e_{α}(t), the Laplace transform of J^{k}e_{α}(t) being given by (2.5). We obtain

J^{k}e_{α}(t) = ∫_{0}^{∞}e −rt K_{α,k}(r) dr + \frac{2}{α} et cos (π/α) cos [t sin (\frac{π}{α}) − k \frac{π}{α}] . (3.16)

Here

K_{α,k}(r) := (−1)^{k} r^{−k} K_{α}(r) = \frac{(−1)^{k}}{π} \frac{r^{α−1−k} sin (απ)}{r^{2α} + 2 r^{α} cos (απ) + 1} ,  (3.17)

with K_{α}(r) = K_{α,0}(r) , can be referred to as the spectral function of the corresponding solution J^{k} e_{α}(t) . This function turns out to be absolutely integrable (locally summable) in IR^{+} for any permissible integer k ≥ 0 , with

∫_{0}^{∞}K_{α,k}(r) dr = \{ \begin{matrix}1 , if 0 < α < 1 , k = 0 , \\ 1 − \frac{2}{α} , if 1 < α < 3 , k = 0 , \\ − \frac{2}{α} cos (k \frac{π}{α}) , if 1 < α < 3 , 1 ≤ k ≤ m − 1 = [α] , \end{matrix}

(3.18)

in order to satisfy the initial conditions e_{α}(0^{+}) = 1 and J^{k}e_{α}(t) → 0 as t → 0^{+} , where 1 ≤ k ≤ m − 1 < α ≤ m .

For the impulse response solution u_{δ}(t) = −u′_{0}(t) = −De_{α}(t) we similarly obtain

D e_{α}(t) = ∫_{0}^{∞}e −rt K_{α,−1}(r) dr + \frac{2}{α} et cos (π/α) cos [t sin (\frac{π}{α}) + \frac{π}{α}] . (3.19)

In Fig 1 we show the plots of the basic spectral functions K_{α}(r) for some values of α in the intervals (a) 0 < α < 1 , (b) 1 < α < 2 , (c) 2 < α < 3 .

We resume the fundamental solutions u_{k}(t) [including u_{δ}(t)] as follows:

(a) 0 < α < 1 ,

u(t) = c_{0} u_{0}(t) + ∫_{0}^{t} q(t − τ) u_{δ}(τ) dτ ,  (3.20a)

where

\{ \begin{matrix}u_{0}(t) & = ∫_{0}^{∞}e −rt K_{α,0}(r) dr , \\ u_{δ}(t) & = − ∫_{0}^{∞}e −rt K_{α,−1}(r) dr , \end{matrix}  (3.21a)

with

\{ \begin{matrix}u_{0}(0^{+}) & = 1 , \\ u_{δ}(0^{+}) & = ∞; \end{matrix}  (3.22a)

(b) 1 < α < 2 ,

u(t) = c_{0} u_{0}(t) + c_{1} u_{1}(t) + ∫_{0}^{t} q(t − τ) u_{δ}(τ) dτ ,  (3.20b)

where

\{ \begin{matrix}u_{0}(t) & = ∫_{0}^{∞}e −rt K_{α,0}(r) dr + \frac{2}{α} et cos (π/α) cos [t sin (\frac{π}{α})] , \\ u_{1}(t) & = ∫_{0}^{∞}e −rt K_{α,1}(r) dr + \frac{2}{α} et cos (π/α) cos [t sin (\frac{π}{α}) − \frac{π}{α}] , \\ u_{δ}(t) & = − ∫_{0}^{∞}e −rt K_{α,−1}(r) dr − \frac{2}{α} et cos (π/α) cos [t sin (\frac{π}{α}) + \frac{π}{α}] , \end{matrix} (3.21b)

with

\{ \begin{matrix}u_{0}(0^{+}) & = 1 , u′_{0}(0^{+}) = 0 , \\ u_{1}(0^{+}) & = 0 , u′_{1}(0^{+}) = 1 , \\ u_{δ}(0^{+}) & = 0 , u′_{δ}(0^{+}) = +∞; \end{matrix}  (3.22b)

(c) 2 < α < 3 ,

u(t) = c_{0} u_{0}(t) + c_{1} u_{1}(t) + c_{2} u_{2}(t) + ∫_{0}^{t} q(t − τ) u_{δ}(τ) dτ ,  (3.20c)

where

\{ \begin{matrix}u_{0}(t) & = ∫_{0}^{∞}e −rt K_{α,0}(r) dr + \frac{2}{α} et cos (π/α) cos [t sin (\frac{π}{α})] , \\ u_{1}(t) & = ∫_{0}^{∞}e −rt K_{α,1}(r) dr + \frac{2}{α} et cos (π/α) cos [t sin (\frac{π}{α}) − \frac{π}{α}] , \\ u_{2}(t) & = ∫_{0}^{∞}e −rt K_{α,2}(r) dr + \frac{2}{α} et cos (π/α) cos [t sin (\frac{π}{α}) − 2\frac{π}{α}] , \\ u_{δ}(t) & = − ∫_{0}^{∞}e −rt K_{α,−1}(r) dr − \frac{2}{α} et cos (π/α) cos [t sin (\frac{π}{α}) + \frac{π}{α}] , \end{matrix} (3.21c)

with

\{ \begin{matrix}u_{0}(0^{+}) & = 1 , u′_{0}(0^{+}) = 0 , u′′_{0}(0^{+}) = 0 , \\ u_{1}(0^{+}) & = 0 , u′_{1}(0^{+}) = 1 , u′′_{1}(0^{+}) = 0 , \\ u_{2}(0^{+}) & = 0 , u′_{2}(0^{+}) = 0 , u′′_{2}(0^{+}) = 1 , \\ u_{δ}(0^{+}) & = 0 , u′_{δ}(0^{+}) = 0 , u′′_{δ}(0^{+}) = +∞ . \end{matrix}  (3.22c)

We recall that in [3] plots of the solution u_{0}(t), which we agree to refer as the basic fundamental solution, are reported for some values of α . In Fig. 2 we quote the plots of this solution for the following cases : (a) α = 0.25 , 0.50 , 0.75 , 1 , and (b) α = 1.25 , 1.50 , 1.75 , 2 , obtained from the first formula in (3.21a) and (3.21b), respectively. We have verified that our present results confirm those obtained in [3] by a numerical treatment and those obtained in [2] by an analytical treatment, valid when α is a rational number.

For plotting we have not considered the case (c) since the corresponding solutions exhibit an exponential growth, and therefore are unstable. For the cases (a) (b) with α not integer, i.e. for the cases of fractional relaxation and fractional oscillation, we note an algebraic decay as t → ∞ , that will be investigated in the next section.

In Fig. 3 we show some plots of the basic fundamental solution u_{0}(t) = e_{α}(t) for α = 1.25 , 1.50 , 1.75 . Here the algebraic decay of the fractional oscillation can be recognized and compared with the two contributions provided by f_{α} (monotonic behaviour ) and g_{α}(t) (exponentially damped oscillation).

4. The decay of the fractional relaxations and fractional oscillations

We are going to describe the asymptotic behaviour as t → ∞ of the fundamental solutions for the cases (a) and (b), namely of u_{0}(t) = e_{α}(t) and for 0 < α < 1 , and u_{0}(t) = e_{α}(t) , u_{1}(t) = J^{1}e_{α}(t) for 1 < α < 2 . For completeness we consider the asymptotic behaviour of the corresponding impulse response solution u_{δ}(t).

Furthermore, for the case (b), we try to get some information about the position and the number of the zeros of the fundamental solution u_{0}, with special regard to the limiting cases α → 1^{+} and α → 2^{−} .

Let us start with the asymptotic behaviour of u_{0}(t) . To this purpose we first derive an asymptotic series for the function f_{α}(t), valid for t → ∞ . Using the identity

\frac{1}{s^{α} + 1} = 1 − s^{α} + s^{2α} − s^{3α} + … + (−1)^{N−1} s^{(N−1)α} + (−1)^{N} \frac{s^{Nα}}{s^{α} + 1} , (4.1)

in formula (3.10) and the Hankel representation of the reciprocal Gamma function, we (formally) obtain the following asymptotic expansion (for α non integer)

f_{α}(t) = ∑_{n=1}^{N}(−1)^{n−1} \frac{t^{−nα}}{Γ(1 − nα)} + O (t^{−(N+1)α}) , as t → ∞ .  (4.2)

The validity of this asymptotic expansion can be established rigorously using the generalized Watson lemma (see Bleistein and Handelsman [22]). We note that this asymptotic expansion coincides with that for u_{0}(t) = e_{α}(t), having assumed 0 < α < 2 (α ≠ 1). In fact the contribution of g_{α}(t) is identically zero if 0 < α < 1 and exponentially small as t → ∞ if 1 < α < 2 . The asymptotic expansion of the solutions u_{1}(t) and u_{δ}(t) are obtained from (4.2) integrating or differentiating term by term with respect to t .

In particular, taking the leading term in (4.2), we obtain the following asymptotic representations as t → ∞ ,

\{ \begin{matrix}u_{0}(t) & ∼ \frac{t^{−α}}{Γ(1 − α)} , (a) (b) \\ u_{δ}(t) & ∼ −\frac{t^{−α−1}}{Γ(−α)} , (a) (b) \\ u_{1}(t) & ∼ \frac{t^{1−α}}{Γ(2 − α)} , (b) \end{matrix}  (4.3)

that point out the algebraic decay of the fundamental and impulse response solutions.

R e m a r k : The asymptotic series (4.2) can also be derived from the existing formula in literature for E_{α}(z) when z is real negative and 0 < α < 2 (α ≠ 1), see e.g. Erdélyi [20].

Now we find it interesting to carry out some considerations about the zeros of e_{α}(t) in the case (b), namely about the solutions of the following equation

e_{α}(t) = f_{α}(t) + \frac{2}{α} e t cos (π/α) cos [t sin (\frac{π}{α})] = 0 .  (4.4)

Of course, in the case (a) e_{α}(t) does not exhibit any zero, while in the case (c) it exhibits infinitely many zeros since the oscillating part g_{α}(t) dominates by exponential increase of amplitude the decaying part f_{α}(t) .

We first note that the function e_{α}(t) exhibits an odd number of zeros, in that e_{α}(0) = 1 , and, for sufficiently large t, e_{α}(t) turns out to be permanently negative, as shown in (4.3) by the sign of Γ(1 − α) . The smallest zero lies in the first positivity interval of cos [t sin (π/α)] , hence in the interval 0 < t < π/[2 sin (π/α)] ; all other zeros can only lie in the succeeding positivity intervals of cos [t sin (π/α)] , in each of these two zeros are present as long as

\frac{2}{α} e t cos (π/α) ≥ |f_{α}(t)| .  (4.5)

When t is sufficiently large the zeros are expected to be found approximately from the equation

\frac{2}{α} e t cos (π/α) ≈ \frac{t^{−α}}{|Γ(1 − α)|} ,  (4.6)

obtained from (4.4) neglecting the oscillation factor of g_{α}(t) [see (3.13)] and taking the first term in the asymptotic expansion of f_{α}(t) [see (4.2-3)]. As we shall see, such an approximate equation turns out to be useful when α → 1^{+} and α → 2^{−} .

For α → 1^{+} , only one zero is present, which is expected to be very far from the origin in view of the large period of the function cos [t sin (π/α)] . In fact, since there is no zero for α = 1, and by increasing α more and more zeros arise, we are sure that only one zero exists for α sufficiently close to 1. Putting now α = 1 + ϵ the asymptotic position T_{∗} of this zero can be found from the relation (4.6) in the limit ϵ → 0^{+} . Assuming in this limit the first-order approximation, we get

cos(π/α) = cos[π/(1 + ϵ)] ∼ cos[π(1 − ϵ)] ∼ −1 ,

and

Γ(1 − α) = Γ(−ϵ) ∼ −1/ϵ .

Therefore

e−T^{∗} ∼ \frac{ϵ(1 + ϵ)}{2} T_{∗}^{−1−ϵ} ,

from which

T_{∗} ∼ log (\frac{2}{ϵ}) + log T_{∗} .  (4.7)

Neglecting log T_{∗} with respect to T_{∗} we obtain

T_{∗} ∼ log (\frac{2}{ϵ}) ,  (4.8a)

which shows that T_{∗} tends to infinity slower than 1/ϵ , as ϵ → 0 . Reversing this asymptotic expression we obtain

ϵ ∼ 2 e^{−T_{∗}} ,  (4.8b)

which shows that ϵ tends to zero exponentially as T_{∗} → ∞ .

For α → 2^{−}, there is an increasing number of zeros up to infinity since e_{2}(t) = cos t has infinitely many zeros [in t_{n}^{∗} = (n + 1/2)π , n = 0, 1, …]. We now are going to give an asymptotic estimate for the largest zero T_{∗} (and for the total number N_{∗} of zeros) when α = 2 − δ in the limit δ → 0^{+} . Assuming in this limit the first-order approximation, we get

cos(π/α) = cos[π/(2 − δ)] ∼ cos[(π/2) (1 + δ/2)] = sin(πδ/4) ∼ πδ/4 ,

and

Γ(1 − α) = Γ(−1 + δ) ∼ −1/δ .

Therefore

e−πδT^{∗}/4 ∼ \frac{(2 − δ)δ}{2} T_{∗}^{−2+δ} ,

from which

\frac{πδT_{∗}}{4} ∼ log (\frac{1}{δ}) + 2 log T_{∗} .  (4.9)

Since both the two terms in the R.H.S. of (4.9) diverge, δ T_{∗} → ∞, so that T_{∗} tends to infinity faster than 1/δ as δ → 0 , or, equivalently, that δ tends to 0 slower than 1/T_{∗} as T_{∗} → ∞ . As a matter of fact both terms turn out to be of the same order as we recognize putting in the L.H.S of (4.9) either T_{∗} ∼ (A/δ) log (1/δ) or δ ∼ B log (T_{∗})/T_{∗} , where a and B are positive constants to be determined. We easily obtain A = B = 12/π, namely the following equivalent asymptotic expressions

T_{∗} ∼ \frac{12}{π δ} log (\frac{1}{δ}) ,  (4.10a)

δ ∼ \frac{12}{π} \frac{log T_{∗}}{T_{∗}} .  (4.10b)

Now, for δ → 0^{+} the length of the positivity intervals of g_{α}(t) tends to π and, as long as t ≤ T_{∗} , there are two zeros in each positivity interval. Hence, in the limit δ → 0^{+} , there is in average one zero per interval of length π , so we expect that N_{∗} ∼ T_{∗}/π .

R e m a r k : For the above considerations we got inspiration from an interesting paper by Wiman [23] who, after having treated the Mittag-Leffler function in the complex plane, considered the position of the zeros of the function on the negative real axis (without providing any detail). While our expression of T_{∗} for α = 1 + ϵ is in agreement with that by Wiman (possibly, unless of a numerical factor!), that of T_{∗} for α = 2 − δ is completely different. In fact, using our notation, Wiman’s expression reads

α = 2 − δ ,  T_{∗} ∼ \frac{4}{π δ} log (\frac{1}{δ}) .  (4.11)

The discrepancy between Wiman’s and our results is (presumably) due to the fact that Wiman (arbitrarily) neglected log T_{∗} with respect to log (1/δ) in (4.9).

The comparison with the numerical results has proven the validity of our analysis, as shown in the plots of Figures 4 and 5, where we show the position of the unique zero when 1 < α ≤ 1.4 and the position of the last zero as 1.8 ≤ α < 2 , respectively.

Furthermore, for completeness we find it interesting to analyse the the phenomenon of the transition of the (odd) number of zeros as 1.4 ≤ α ≤ 1.8 . For this purpose, in Table I we report the intervals of amplitute Δα = 0.01 where these transitions occur, and the location T_{∗} (evaluated within a relative error of 0.1% ) of the largest zeros found at the two extreme values of the above intervals.

We recognize that the transition from 1 to 3 zeros occurs as 1.40 ≤ α ≤ 1.41, that one from 3 to 5 zeros occurs as 1.56 ≤ α ≤ 1.57, and so on. The last transition in the considered range of α is from 15 to 17 zeros, and it just occurs as 1.79 ≤ α ≤ 1.80 .

N_{∗}  α  T_{∗}
1 ÷ 3  1.40 ÷ 1.41  1.730 ÷ 5.726
3 ÷ 5  1.56 ÷ 1.57  8.366 ÷ 13.48
5 ÷ 7  1.64 ÷ 1.65  14.61 ÷ 20.00
7 ÷ 9  1.69 ÷ 1.70  20.80 ÷ 26.33
9 ÷ 11  1.72 ÷ 1.73  27.03 ÷ 32.83
11 ÷ 13  1.75 ÷ 1.76  33.11 ÷ 38.81
13 ÷ 15  1.78 ÷ 1.79  39.49 ÷ 45.51
15 ÷ 17  1.79 ÷ 1.80  45.51 ÷ 51.46

Table I N_{∗} = number of zeros, α = fractional order, T_{∗} location of the largest zero.

Acknowledgements

The authors like to thank Gianni de Fabritiis for help in numerical computations. The authors appreciate the support given by CNR for an extended visit of the first named author to University of Bologna when this research was started. From his side this work is also embedded in the research project ”Convolutions”, supported by the Research Commission of Free University of Berlin, and in the project ”Fractional Order Systems”, supported by NATO Collaborative Grant 940 843 (R. Gorenflo and R. Rutman).

References

[1] R. Gorenflo and R. Rutman: On ultraslow and on intermediate processes, in P. Rusev, I.

Dimovski and V. Kiryakova (eds): Transform Methods and Special Functions, Sofia

1994, Science Culture Technology, Singapore 1995, pp. 61-81.

[2] F. Mainardi: Fractional relaxation-oscillation and fractional diffusion-wave phenomena,

Chaos, Solitons and Fractals, in press (1996).

[3] L. Blank: Numerical treatment of differential equations of fractional order, MCCM Numerical Analysis Report No 287, The University of Manchester, 1996.

[4] M. Caputo: Elasticità e Dissipazione, Zanichelli, Bologna 1969 [in Italian].

[5] M. Caputo and F. Mainardi: Linear models of dissipation in anelastic solids, Riv. Nuovo Cimento (Ser. II), 1, 161-198 (1971).

[6] M. Caputo: A model for the fatigue in elastic materials with frequency independent Q, J. Acoust. Soc. Am., 66 (1), 176-179 (1979).

[7] P.J. Torvik and R.L. Bagley: On the appearance of the fractional derivatives in the behavior of real materials, J. Appl. Mech. (Trans. ASME), 51, 294-298 (1984).

[8] R. C. Koeller: Applications of fractional calculus to the theory of viscoelasticity, J. Appl. Mech. (Trans. ASME), 51, 299-307 (1984).

[9] W.R. Schneider and W. Wyss: Fractional diffusion and wave equations, J. Math. Phys., 30, 134–144 (1989).

[10] T.F. Nonnenmacher: Fractional integral and differential equations for a class of Levy-type probability densities, J. Phys A: Math. Gen., 23, L697-L700 (1990).

[11] T.F. Nonnenmacher and W.G. Glöckle: A fractional model for mechanical stress relaxation, Phil. Mag. Letters, 64 (2), 89-93 (1991).

[12] F. Mainardi, Fractional relaxation in anelastic solids, J. Alloys and Compounds, 211/212, 534-538 (1994).

[13] H. Beyer, S. Kempfle: Dämpfungsbeschreibung mittels gebrochener Ableitungen, ZAMM, 74, T657-T660 (1994).

[14] H. Beyer, S. Kempfle: Definition of physically consistent damping laws with fractional derivatives, ZAMM, 75, 623-635 (1995).

[15] I. Podlubny: Solutions of linear fractional differential equations, in P. Rusev, I. Dimovski

and V. Kiryakova (eds): Transform Methods and Special Functions, Sofia 1994,

Science Culture Technology, Singapore (1995), pp. 227-237.

[16] R. Gorenflo and S. Vessella: Abel Integral Equations: Analysis and Applications,

Lecture Notes in Mathematics # 1461, Springer-Verlag, Berlin 1991.

[17] S.G. Samko, A.A. Kilbas and O.I. Marichev: Fractional Integrals and Derivatives, Theory and Applications, Gordon and Breach, Amsterdam 1993.

[18] K.S. Miller and B. Ross: An Introduction to the Fractional Calculus and Fractional Differential Equations, Wiley, New York 1993.

[19] G. Doetsch: Introduction to the Theory and Application of the Laplace Transfor-

mation, Springer Verlag, Berlin 1974.

[20] A. Erdélyi (ed.): Higher Transcendental Functions, Vol. 3. Ch. 18, McGraw-Hill, New York 1955.

[21] M.M. Dzherbashyan: Integral Transforms and Representations of Functions in the

Complex Plane, Nauka, Moscow 1966 (in Russian).

[22] N. Bleistein and R. A. Handelsman: Asymptotic Expansions of Integrals, Ch. 4, p. 162, Dover, New York 1986.

[23] A. Wiman: Über die Nullstellen der Funktionen E_{α}(x) , Acta Math., 29, 217-234 (1905).