Knocking on Leibniz's Door

Fractional Calculus

The fractional calculus is implemented in Mathematica. The author discusses the theoretical background and demonstrates the practical application of the fractional calculus to the tautochrone problem formulated by Abel and a fractional relaxation problem in Weyl calculus.

by Gerd Baumann

The term fractional calculus is by no means new. It is a generalization of the ordinary differentiation by non-integer derivatives. The subject is as old as calculus itself and goes back to the time when Leibniz, Gauss, and Newton invented differentiation as a new tool. In a letter to L`Hopital in 1695 Leibniz raised the following question:

Gottfried Wilhelm von Leibniz; Born: 1 July 1646 in Leipzig,  Died: 14 November 1716 in Hannover.  Leibniz developed the present day notation for the differential and integral calculus. He never thought of the derivative as a limit.

Can the meaning of derivatives with integer order (((d^n) y(x))/(dx^n)) be generalized to derivatives with non-integer order; with, in general, n being a complex number? The story goes that L`Hopital was somewhat curious about that question and replied with another question to Leibniz. What if n  =  (1/2)? In a letter dated September 30, 1695 Leibniz replied: Il y a de l'apparence qu'on tirera un jour des consequences bien utiles de ces paradoxes, car il n'y a guere de paradoxes sans utilité. The translation reads: "It seems that useful consequences shall be one day drawn from these paradoxes, as there are no paradoxes that do not prove useful."

The question raised by Leibniz for a fractional derivative has been an ongoing topic during the last 300 years. Many mathematicians (e.g. Liouville, Riemann, and Weyl) made major contributions to the theory of fractional calculus.

Let us consider Leibniz's question in connection with a special set of functions, say powers. In this case a fractional derivative is useful and can be expressed again by powers. For example let us consider the nth derivative of (x^m). We know that the general expression for the nth derivative is given by

(((d^n) (x^m))/(dx^n)) =  ((m!)/((m - 
                  n)!)) (x^(m - n)).

We also know that the factorial is connected with Euler's Γ-function by n != Γ(n + 1). Replacing the factorials in (1) with Γ-functions, we are able to write

\(TraditionalForm
      \`\(\(d\^n\) x\^m\)\/dx\^n = \ 
        \(\(Γ(m + 1)\)\/\(Γ(m - n + 1)\)\) 
          x\^\(m - n\)\).

This representation is equivalent to equation (1). However, equation (2) contains the potential of a generalization. We know that the Γ-function is defined for continuous arguments over the complex domain of numbers. If we change the integer value of n to a number q ∈  ℂ, we are able to generalize an integer differentiation to a noninteger one. We are even able to define a complex differentiation. Replacing n by q in (2) results into the general expression

\(TraditionalForm
      \`\(\(d\^q\) x\^m\)\/dx\^q = \ 
        \(\(Γ(m + 1)\)\/\(Γ(m - q + 1)\)\) 
          x\^\(m - q\)\).

From the mathematical point of view, relation (3) has a well defined meaning. However, this meaning is restricted to the special class of powers (x^m). If we try to fractionally differentiate such simple functions in Mathematica, we end up with the following result

\(∂\_{x, 1/2}\ x\^2\)

This means that Mathematica is not capable of dealing with fractional differentiation. However, the developers of Mathematica  designed the system in a way that the user can extend the definition of derivatives. The extension of the derivative D[] will be our next task. Telling Mathematica that fractional derivatives of powers are useful mathematical constructs is realized by the following lines

Unprotect[D];

\(D[x_\^m_. , {x_, q_}] := 
      \(Gamma[m + 1]\/Gamma[m - q + 1]\) x\^\(m - q\) /; 
        Head[q] == Real\  || \ Head[q] == Rational || Head[q] == Complex\)

Protect[D];

The definition of the fractional derivative of powers is based on equation (3) and restricts the order of differentiation to a rational, a real, or a complex number. A derivative for a negative rational number can now be calculated by

\(∂\_{x, \(-\(1\/2\)\)}\ x\)

If we set the order of differentiation q to a real number, we find

\(∂\_{x, 2.1}\ x\^2\)

If we use complex numbers as the differentiation order, we get the result

\(∂\_{x, 11.5 + I}\ x\^4\)

This kind of formulas were first discussed by Lacroix in 1819 [1]. In retrospect, formula (3) is the first analytical answer on Leibniz's question on fractional derivatives. The story on the fractional calculus was continued with contributions from Fourier, Abel, Liouville, Riemann and Weyl. For a historical survey, the reader may consult the books of  Oldham and Spanier [2] or Miller and Ross [3]. The historical developments culminated in two calculi which are based on the work of Riemann [4] and Liouville [5] on the one hand and on the work of Weyl [6] on the other. The formulations are connected and Weyl's calculus (W) is a subset of the Riemann-Liouville (RL) calculus. In the following section, we will introduce the RL calculus. We will also outline the connection between the W and RL calculi. The application of the two calculi to two physical problems will demonstrate the usefulness of the methods.

. The Riemann-Liouville and Weyl calculi

In the above discussion, we motivated the fractional derivative by heuristics using properties of Euler's Γ-function. Such a discussion is restricted to a special class of functions; i.e. powers. In this section we will motivate and introduce a calculus allowing the calculation of fractional derivatives for a more general class of functions. The results will be two calculi mainly conceived by Riemann and Liouville, and Weyl. Paradoxically the basis of the two calculi is not a derivative but an integral. However, we can understand integration as an inverse differentiation if we define the symbols

\(TraditionalForm
    \`\(d\^\(-1\)\/dx\^\(-1\)\) \(f(x)\)\  := \ 
      ∫\_0\%x\( f(t)\) dt\)

\(TraditionalForm
      \`\(d\^\(-2\)\/dx\^\(-2\)\) \(f(x)\)\  := \ 
        ∫\_0\%x
            \(∫\_0\%t\( f(s)\) ds 
              dt\)\) …

The negative order of differentiation means nothing more than an integration. Higher orders of negative differentiation are calculated by nesting the integrals on the right hand side. We will abbreviate this kind of recursion by the symbol (Subscript[𝒟, 0, x]^(-n)) where n is a positive integer, and 0, x denote lower and upper boundaries of the integral. Thus equation (4) is reduced to

\(TraditionalForm
      \`\(𝒟\_\(0, x\)\%\(-1\)\) \(f(x)\)\  = \ \ 
        ∫\_0\%x\( f(t)\) dt\).

The symbol (Subscript[𝒟, 0, x]^(-n)) contains the complete information to carry out the calculation of the negative differentiation. The lower two indices denote the lower and upper boundaries on the integral. The superscript indicates the order of "differentiation". A generalization of the notation is an arbitrary starting point  a at the lower boundary in the integral; meaning

\(TraditionalForm
      \`\(𝒟\_\(a, x\)\%\(-1\)\) \(f(x)\)\  = \ \ 
        ∫\_a\%x\( f(t)\) dt\).

If we consider the nth negative derivative (Subscript[𝒟, a, x]^(-n)) of an arbitrary function f(x), we write

\(TraditionalForm
      \`\(𝒟\_\(a, x\)\%\(-n\)\) \(f(x)\)\  = \ 
        ∫\_a\%x
            \(∫\_a\%x\_\(n - 1\)\(f(x\_0)\) 
              dx\_\(n - 1\) … 
              dx\_0\)\).

Now remember Cauchy's integral formula

((d^n)/(dx^n)) f(
            x) =  ((n!)/(2  πi)) (∫_C)(((ζ - 
                      z))^(-n - 
                    1)) f(ζ) dζ,

and applying (9) to (8), we can reduce equation (8) to

\(TraditionalForm
      \`\(𝒟\_\(a, x\)\%\(-n\)\) \(f(x)\)\  = \ 
        \(1\/\(\((n - 1)\)!\)\) 
          \(∫\_a\%x\(\((x - x\_0)\)\^\(n - 1\)\) \(f(x\_0)\) 
              dx\_0\)\).

Introducing again the Γ-functions for the factorials, we can generalize the result to an arbitrary order of fractional differentiation by n! -> Γ(n + 1). The general formula for arbitrary negative derivatives reads

(Subscript[𝒟, a, x]^(-q)) f(
            x)  =  (1/(Γ(
                  q))) (∫_a^x)(((x - (Subscript[x, 0])))^(q - 1)) f(
                  (Subscript[x, 0])) d(Subscript[x, 0])             q > 
          0,

where q is a fractional, real, or complex number. This kind of operator is known as the Riemann (R) version of the fractional integral [3]. The Liouville (L) version of this operator follows if we replace the lower boundary a in (11) by -∞; i.e. (Subscript[𝒟, -∞, 
            x]^(-q)) f(x)  is the Liouville fractional integral. The case with a = 0

\(TraditionalForm
    \`\(𝒟\_\(0, x\)\%\(-q\)\) \(f(x)\)\  = \ 
      \(1\/\(Γ(q)\)\) 
          \(∫\_0\%x\(\((x - x\_0)\)\^\(q - 1\)\) \(f(x\_0)\) 
              dx\_0\ \ \ \ \ \ \ \ \ \ \ \ \ q\) > 0\)

is called the Riemann-Liouville (RL) fractional integral. We observe that the various definitions of fractional integrals differ only in the lower boundary of the integral. You may think that this small difference is of minor importance. The following section will demonstrate that this assumption is not correct. The change of the lower boundary has far reaching consequences in the calculation of fractional derivatives.

So far we defined the notation of the fractional integral. A fractional derivative is connected with a fractional integral by introducing a positive order of differentiation in the operator  (Subscript[𝒟, a, x]^(-q)). This shift in the order of differentiation can be achieved by introducing an ordinary differentiation followed by a fractional integration. We thus define a fractional differentiation by

\(TraditionalForm
      \`\(𝒟\_\(a, x\)\%s\) \(f(x)\) := \ 
        \((d\^n\/dx\^n)\) \(𝒟\_\(a, x\)\%\(-\((n - s)\)\)\) 
          \(f(x)\)\)      with  n  ∈ ℕ,  s > 0,  
      n - s  > 0.

In this Riemann notation the fractional derivative depends on a lower boundary a of the integral. This dependence disappears if we consider only the RL operator with a = 0. Now we know in principle how general functions are treated in a fractional calculus. Before we apply the given notation to fractional differential equations, we state some general properties of the fractional derivative necessary to realize this operator in Mathematica.

The main properties needed in an implementation are linearity and the composition rule. These are two basic properties besides the Leibniz rule of differentiation and the chain rule. From a practical point of view only the property of linearity and the composition of derivatives are needed. The first property of linearity reads

\(TraditionalForm
    \`𝒟\_\(a, x\)\%s\ 
        \((α\ \(f(x)\) + \ β\ \(g(x)\))\)\  = 
      \  α\ \(𝒟\_\(a, x\)\%s\) 
          \(f(x)\) + \ 
        β\ \(𝒟\_\(a, x\)\%s\) \(g(x)\)\)

with α and β real constants. The composition rule of fractional derivatives combines two derivatives with different orders s and p in the form

\(TraditionalForm
    \`𝒟\_\(a, x\)\%s\ \(𝒟\_\(a, x\)\%p\) 
        \(f(x)\)\  = \ 
      \(𝒟\_\(a, x\)\%\(s + p\)\) \(f(x)\)\)

with p < 0 and f(x) finite at x = a. These properties are sufficient to define a function RiemannLiouville[] in Mathematica. The following lines contain the main steps of the definitions necessary to implement the RL operator.

\(

 (*\ \(\(--\(-\ 
              treatment\)\)\ of\ \(constants\ --\)\)\(-\)\ *) 
\(RiemannLiouville[f_, {x_, order_, a_:  0}] := 
        f\ RiemannLiouville[1, {x, order, a}] /; 
          FreeQ[f, x] && \ Not[SameQ[f, 1]];\)
	
 (*\ \(--\(-\ 
            remove\)\)\ constants\ from\ the\ \(argument\ --\)\ *) \)
\(\(RiemannLiouville[c_\ f_, {x_, order_, a_:  0}] := 
        c\ RiemannLiouville[f, {x, order, a}]\  /; \ 
          FreeQ[c, 
            x];\)

 (*\ \(--\(-\ \(linerity\ --\)\)\)\(-\)\ *) \)
\(RiemannLiouville[f_\  + \ g_, {x_, order_, a_:  0}] := 
      RiemannLiouville[f, {x, order, a}]\  + \ 
        RiemannLiouville[
          g, {x, order, 
            a}]
	
 (*\ \(--\(-\ \(composition\ --\)\)\)\(-\)\ *) \)
\(RiemannLiouville[\ 
        RiemannLiouville[f_, {x_, order1_, a_:  0}], {x_, order2_, a_:  0}] := 
      RiemannLiouville[f, {x, order1 + order2, a}]\  /; 
        order1 < 0

 (*\ \(--\(-\ \(Identity\ --\)\)\)\(-\)\ *) \)
\(RiemannLiouville[f_, {x_, 0, a_:  0}] := 
      f


 (*\(\(--\(-\ main\)\)\ \(function\ --\)\)\(-\)\ *) \)
\(\(\(RiemannLiouville[
        f_, {x_, order_, a_:  0}]\)\(:=\)\(Block[{n, int, $y1, q}, 
		If[
          NumericQ[order] && Simplify[order > 0], n = Floor[order]; \ 
          q = \ order - n, q = order; n = 0]; 
		int\  = \ 
          Integrate[\(\((x - $y1)\)\^\(\(-q\) - 1\)\) \((f /. 
                  x :> $y1)\), {$y1, a, x}, 
            GenerateConditions -> False]; 
		D[
            int\ /\ Gamma[\(-q\)], {x, n}]\  /; \ 
          FreeQ[int, $y1]
		]\)\(
\)\(			\)\)\)

The above program lines show that it is sufficient to use the properties and the definition of the RL operator given in equation s(11)-(15). The symbolic notation introduced in the discussion above can be generated in Mathematica by the utilities package Notation.

Needs["Utilities`Notation`"]

We define the symbolic pattern on the left hand side  and the function on the right hand side.

Notation[(Subscript[𝒟, a_, x_]^order)_ f_ ⟺ RiemannLiouville[f_, {x_, order_, a_}]]

A small palette allows us to generate the symbolic notion:

(Subscript[𝒟, \[Placeholder], 
            \[Placeholder]]^\[Placeholder]) \[Placeholder]

where the placeholders have the meaning given in the definition of the notation. Application of the symbol (Subscript[𝒟, 0, x]^(-1/2)) to the power function x gives the same result as calculated with the Lacroix definition (3)

(Subscript[𝒟, 0, x]^(-1/2)) x

An example frequently discussed in the literature is the differentiation of a constant [4]. From ordinary calculus we know that the differentiation of a constant vanishes. Surprisingly this is not the case in fractional calculus. Using the palette notation of the RL operator, we can write for a semi-differentiation

(Subscript[𝒟, 0, x]^(1/2)) 1

This result contradicts the general knowledge that the differentiation of a constant vanishes. This behavior can be understood if we go back to the definition of a fractional derivative based on an integral. Another interesting example of a fractional derivative occurs if we replace the lower limit of integration by a constant, say a

(Subscript[𝒟, a, x]^(1/2)) 1

The result shows that the lower boundary in the integral has some influence on the result. We also see that the results of the RL operator and the Riemann representation are different. The result for the Liouville definition is only accessible if we consider the limit of a -> -∞.

Limit[(Subscript[𝒟, a, x]^(1/2)) 1, 
      a -> -∞]

We observe that the result of the fractional differentiation mainly depends on the definition of the fractional operator. Another example demonstrating the exceptional behavior of the fractional derivative is given by

(Subscript[𝒟, 0, x]^q) (((a + x))^p)

The fractional differentiation of a general polynomial of order p results in a hypergeometric function (Subscript[F, 2, 1]). So far, we have demonstrated some properties of the Riemann-Liouville calculus. Another fractional operator with a broad physical application is the fractional Weyl operator. Weyl defined the fractional derivative in a way similar to Riemann and Liouville [6]. In 1917 Weyl stated his operator as

\(\(\(𝒲\&+\)\_t\%\(-q\)\) \(f(t)\)\) 
             =  (1/(Γ(q))) 
              (∫_t^∞)(((ζ - t))^(q - 1)) 
                  f(ζ) dζ      with    Re(q) > 0,  t > 0.

This definition is related to the Liouville definition by the relation

\(TraditionalForm
      \`\(\(𝒲\&+\)\_t\%\(-q\)\) \(f(t)\)\  = \ \ 
        \(\((\(-1\))\)\^q\) 
          \(𝒟\_\(∞, \ t\)\%\(-q\)\) \(f(t)\)\).

The Weyl and Liouville definitions differ in the change of limits in the integral and the factor (((-1))^q). Knowing this behavior, we can define the Weyl operator by

WeylPlus[f_, {t_, q_}] := 
      RiemannLiouville[(((-1))^q) f, {t, q, ∞}]

The related template is defined by

Notation[\(\(𝒲\&+\)\_t_\%q_[f_]\) ⟺ WeylPlus[f_, {t_, q_}]]

\(\(𝒲\&+\)\_\[Placeholder]\%\[Placeholder][
        \[Placeholder]]\)

Examples for some Weyl differentials and integrals follow. The first example is concerned with the qth differentiation of a negative power:

\(\(𝒲\&+\)\_x\%q[x\^\(-2\)]\)

The second example deals with a general polynomial of pth order:

\(\(𝒲\&+\)\_x\%q[\((x + a)\)\^p]\)

The Weyl fractional derivative is again a polynomial of order p - q. A semi-Weyl integration of the power -1/2 of 1/x results in

\(\(𝒲\&+\)\_x\%\(\(-1\)/2\)[1\/x]\)

The general Weyl derivative of νth order of (x^(-μ)) gives us again a power relation of order -μ-ν.

\(\(𝒲\&+\)\_x\%ν[x\^\(-μ\)]\)

The fractional integral of Weyl for the same function results in a power behavior with -μ+ν

\(\(𝒲\&+\)\_x\%\(-ν\)[x\^\(-μ\)]\)

Other examples that may be obtained by fractional Weyl integration are

\(\(𝒲\&+\)\_x\%\(-ν\)[
      \(x\^\(\(-1\)/2\)\) Exp[\(-α\)\ x\^\(1/2\)]]\)

containing the modified Bessel function of the second kind. More complicated integrals may be obtained by other functions in the argument, such as

\(\(𝒲\&+\)\_x\%\(-ν\)[
      \(x\^\(ν - 1\)\) Sin[α\ x]]\)

using the Bessel J function for representation of the result. Another non trivial example is

\(\(𝒲\&+\)\_x\%\(-ν\)[\(x\^\(-μ\)\) Log[x]]\)

where PolyGamma[n, z] gives n\(TraditionalForm\`th\^\*""""\) derivative of the digamma function (ψ^((n)))(z). PolyGamma[n,z] is given by (ψ^((n)))(
          z) = (d^n) ψ(z)/d  (z^n). Note that the digamma function corresponds to (ψ^((0)))(z). The general form (ψ^((n)))(z) is the (n + 1) th\(TraditionalForm\`\^\*""""\), not the n\(TraditionalForm\`th\^\*""""\) logarithmic derivative of the gamma function. The polygamma functions satisfy the relation (ψ^((n)))(
          z) = (((-1))^(n + 
                1)) n! (∑_(k = 0)^∞) 
              1/(((z + k))^(n + 1)).

So far, we implemented the basic properties of the Riemann–Liouville and Weyl fractional calculi. We already observed that application of the RL and W operators to well-known functions delivers extraordinary results. In the following section, we will discuss how these results are useful in connection with physical applications.

Applications of the RL and W calculi

In this section we discuss some examples demonstrating the physical relevance of the fractional calculus. We consider approaches and arguments used by early researchers in their attempts to apply the fractional calculus to physical problems. We discuss Abel's problem of the tautochrone as an example of the Riemann-Liouville calculus and fractional relaxation in connection with the Weyl calculus. Both examples demonstrate the usefulness of the fractional calculus.

The Tautochrone problem

Let us examine a particle of mass M moving in a vertical gravitational field along a path ℭ. The initial velocity of the particle is zero. Now let us find the shape of the path ℭ for which the time of descent τ from the starting point P to the origin is independent of the starting point. A curve ℭ describing this kind of movement is called a tautochrone.

The tautochrone problem should not be confused with the brachistochrone problem. That problem was to find the shape of the path ℭ such that the time of descent from P to the origin would be a minimum. However, as we will see, both problems  are solved by cycloids.

We now proceed to formulate the physical problem. Let s be the arc length measured along ℭ from the origin to an arbitrary point Q with coordinates (ξ,η) on ℭ. At this point (ξ,η), the particle of mass M increases the force -M g cos(β) with β an angle of inclination and g the gravitational constant. Newton's second law tells us

\(TraditionalForm
      \`M \(\( d\^2\) \(s(t)\)\)\/dt\^2 = \ \(-M\)\ g\ \(cos(β)\)\).

The change of the vertical coordinate η along the arc is given by cos(β)  =  dη/ds. Thus the equation of motion can be expressed by

(((d^2) s)/(dt^2))  =  -g d(η/ds)

Multiplying  equation (19) by ds/dt and integrating with respect to t gives

\(TraditionalForm
    \`\((ds\/dt)\)\^2\  = \ \(-2\)\ g\ η\  + \ c\_1\)

where (Subscript[c, 1]) is a first constant of integration. Since the particle started from rest, ds/dt is zero when η = y, and thus (Subscript[c, 1]) = 2  g y. The velocity of the particle becomes ds/dt = -(Sqrt[2  g(y - η)]). Thus the time of descent τ from P to the origin is

τ = -(1/(Sqrt[2  g])) 
                (∫_P^O) (1/(Sqrt[y - η])) 
                    ds.

It is obvious that the arc length s is a function of η, say s = h(η). h defines the shape of path ℭ. The line element ds in (21) is thus replaced by ds =  h' (η) dη which results in

\(TraditionalForm
    \`\@\(2\ g\)\ τ\  = \ 
      ∫\_0\%y\(\((y - η)\)\^\(\(-1\)/2\)\) h' \((η)\)\ 
          dη\)

Introducing the abbreviation f(η)  = h' (η) allows us to introduce the RL operator on the right hand side of (22)

(((Sqrt[2 g]) τ)/Gamma)[(1/2)] ==  
      (Subscript[𝒟, 0, y]^(-1/2)) f[y]

This relation was first used by Abel to solve the tautochrone problem. He attacked the equation by applying the RL operator  (Subscript[𝒟, 0,  y]^(1/2)) on both sides of the above equation

(Subscript[𝒟, 0, y]^(1/2)) 
        ((((Sqrt[2 g]) τ)/Gamma)[(1/2)]) == 
      (Subscript[𝒟, 0, y]^(1/2)) 
        ((Subscript[𝒟, 0, y]^(-1/2)) (f[y]))

The use of the composition relation allows to express the unknown function f(y) in an explicit way. The knowledge of f delivers the relation for h'. Since we know that  
      f = h'  =  ds/dy  =  (Sqrt[1 + (((dx/dy))^2)]),  we can write

equation  = 
      (Sqrt[1 + ((((∂_y) x[y]))^2)]) ==   
        (Subscript[𝒟, 0, y]^(1/2)) 
          (((Sqrt[2 g]) τ)/Gamma)[(1/2)]

The resulting first-order ordinary differential equation can be solved by DSolve[] to be

sol  = DSolve[equation, x, y]

The solutions represented above can be plotted for specific values of the parameters τ, g, and C[1]. The constant of integration C[1] can be chosen to be C[1]=0 since path  ℭ must intersect the origin. The gravitational constant takes the value of  9.81  m/(s^2) on earth. The time of descent τ is independent of the starting point and thus can be set to an arbitrary value (τ = 1).

sol1  =  
      x[y]  /.  sol  /. 
         {g -> 9.81, τ -> 1, C[1] -> 0}

The graphical representation of the two solution branches follows by

Plot[{sol1[[1]], 
        sol1[[2]]}, {y, 0, 1.8}, 
      PlotStyle 
        -> {RGBColor[0.996109,  0,  0], RGBColor[0,  0,  0.996109]}, 
      AxesLabel -> {"y", "x"}]

The two curves  show the solution of  the tautochrone problem for a specific choice of  parameters τ, g and C[1]; (τ = 1,  g = 9.81,  C[1] = 0). The upper branch is the solution realized in nature. The lower branch is a another mathematical possibility.

To find a more readable representation of the solution, let us replace coordinate y by the expression ((2 g (τ^(2 )) Sin[α]^2)/(π^2)), corresponding to a parametric representation of the path by means of angle α.

asol  =  
      (x[y]  /. sol )  /. 
               {y -> 
                  ((2 g (τ^(2 )) Sin[α]^2)/(π^2)), 
                C[1] -> 0} // Simplify // PowerExpand // Expand

The result is a parametric representation of a cycloid with coordinates

cycloid  =  
      Transpose[{
          asol, {((2 g (τ^(2 )) Sin[α]^2)/(π^2)), 
            ((2 g (τ^(2 )) Sin[α]^2)/(π^2))}}]

The two solutions are plotted for different times of descent τ by

ParametricPlot[
      Evaluate[Flatten[
          Table[cycloid /. {g -> 9.81, τ -> i}, {i,  .5, 
              1,  .1}], 1]], {α, 0, 2  π}, 
      AxesLabel -> {"x", "y"}, 
      PlotStyle 
        -> {RGBColor[0.996109,  0,  0], RGBColor[0,  0,  0.996109]}]

The figure shows the two solutions of  the tautochrone problem for different values of τ∈[0.5,1] in steps of 0.1. The right hand side of  the plot  represents the first solution,  the left side  of  the plot shows the second one.

At the end, we recognize that the tautochrone problem is solved by a cycloid.

A fractional relaxation equation

The second example discussed is related to the fractional calculus by Weyl. This example is connected with the problem of a fractional relaxation process. The standard Maxwell–Debye relaxation process, for example, is modeled by the initial value problem

τ (d/(d t)) ϕ(t)  =  -ϕ(t), with t > 0, and ϕ(0) = (Subscript[ϕ, 0])

where (Subscript[ϕ, 0]) is a constant and τ denotes the relaxation time. The solution is given by an exponential

ϕ(t) ={({{(Subscript[ϕ, 0]) (e^(-t/τ))     if
                     t ≥  0}, {0                   if t < 0. }})

The initial value problem (23) reads in an integral representation

\(TraditionalForm
      \`ϕ(t)\  - ϕ\_0\  = \ 
        \(\(-τ\^\(-1\)\) 
            \(∫\_0\%t\( ϕ(t')\)\ dt'\)\  = 
          \( : \ \(-τ\^\(-1\)\)\ 
            \(\(d\^\(-1\)\) \(ϕ(t)\)\)\/dt\^\(-1\)\)\)\).

The integral equation (25) incorporates the initial value (Subscript[ϕ, 0]) right from the beginning. This representation of a relaxation process contains the total information on the process in a nutshell. Relaxation processes deviating from the classical Maxwell–Debye behavior are referred to as non-Debye, nonexponential or anomalous relaxation processes. Anomalous relaxation processes have been observed in dielectrics, in diffusion-controlled relaxations, in liquid crystals, polymer melts, amorphous polymers, rubber, biopolymers, and other disordered systems. The deviation from the exponential decay is part of current problems in physics which have not been completely resolved. Typical non-Debye relaxation processes are described empirically either by a Kohlrausch–Williams–Watts (KWW) decay

\(TraditionalForm
    \`ϕ(t)\  = \ 
      \(ϕ\_0\) e\^\(-\((t/τ)\)\^α\)\)

with 0 < α < 1, or by an asymptotic power-law (Nutting law)

\(TraditionalForm
      \`ϕ(t)\  = \ 
        ϕ\_0\ 1\/\((1 + t/τ)\)\^γ\  ∼ \ 
          t\^\(-γ\)\)  if  t/τ -> ∞

with 0 < γ < 1. Both behaviors are connected to each other. In specific experiments, the transition between the KWW and the Nutting law can be observed, provided the range of the time domain extends over several decades. In fact, it is possible to find an analytic function interpolating between the KWW and the Nutting law. The interpolating function is a Mittag-Leffler function which is contained in the generalized function class of Fox's H-functions [7].

The Mittag-Leffler function, as a special member of the class of Fox's H-functions, describes anomalous relaxation processes. These occur as solutions of a generalized relaxation equation which we now discuss. The basis of the generalization is the integral representation (25) of a standard relaxation process. The generalization consists in introducing a fractional integral order in (25) by

\(TraditionalForm
    \`ϕ(t)\  + ϕ\_0\  = \ 
      \(\((\(-τ\))\)\^\(-q\)\ 
          \(\(d\^\(-q\)\) \(ϕ(t)\)\)\/dt\^\(-q\) = 
        \( : \ \(\((\(-τ\))\)\^\(-q\)\) 
            \(\(𝒲\&+\)\_t\%q\) \(ϕ(t)\)\)\)\)

where q, the fractional differentiation order, is a real constant  > 0. The symbol \(TraditionalForm\`\(𝒲\&+\)\_t\%\(-q\)\) denotes the Weyl  integral operator defined by \(TraditionalForm
      \`\(\(𝒲\&+\)\_t\%\(-q\)\) \(ϕ(t)\)\  = 
        ∫\_t\%∞\( ϕ(t')\)/\((t - t')\)\^\(1 - q\) 
            dt'/\(Γ(q)\)\). This generalization of a standard relaxation process assumes that anomalous relaxation processes are deeply connected with memory integrals.

Solution of such an equation is gained either by a power series ansatz or by Laplace–Mellin techniques. We prefer the second procedure because it directly connects the solution to Fox's H-functions. The first step to solve equation (28) is to apply the  Laplace transform (28) and gain an algebraic solution in Laplace space. The Laplace transform of equation (28) in Mathematica follows from

Assume[q > 0];

\(eq1 = 
      \(\((\(-τ\))\)\^q\) \(𝒲\&+\)\_t\%q[ϕ[t]] == 
        ϕ0 + ϕ[t]\)

under the assumption q > 0.

eql = (Subscript[ℒ, t]^s)[eq1] // PowerExpand

where C1[] is a constant introduced in the Laplace transform of the W operator. In the following we assume that this constant vanishes. The algebraic solution of (28) in the Laplace space is given by

sl = 
      Solve[eql, (Subscript[ℒ, t]^s)[ϕ[t]]] /. C1[] -> 0 // 
        Simplify

Since this expression contains terms with arbitrary powers q, we face the problem of calculating the inverse Laplace transform. At this point we note that Miller and Ross [3] developed a solution strategy for fractional differential equations that works well for the case if the fractional order q is a rational number. A way out of this restriction, applicable to arbitrary q > 0, is the additional transform to the Mellin space. The inclusion of the Mellin transform allows us to resolve the pole structure of the algebraic expression in the complex plane. The Mellin transform of the Laplace solution becomes

sm  = (Subscript[ℳ, s]^p)[sl]

where (Subscript[ℳ, s]^p)denotes the Mellin transform. A shift of the Mellin variable p allows us to rewrite the result

ssm = 
      Solve[sm /. Rule -> Equal /. p -> 1 - p, 
        (Subscript[ℳ, s]^p)[ϕ[s]]]

Since the inverse Mellin transform is closely connected with the definition of a Fox's H-function, we can use this link to derive the solution as

shh = 
      (Subscript[(((ℳ^(-1)))), p]^t)[ssm] // PowerExpand // Flatten

The result is a Mittag-Leffler function \(TraditionalForm\`\(E\_\(q, 1\)\)(
          z)\  = \ ∑\_\(k = 0\)\%∞ z\^k\/\(Γ(q\ k + 1)\)\). For q -> 1 this series approaches the exponential function. Derivation of the solution was possible since the Mellin representation of the solution is represented by Γ-functions. Whenever we can represent the algebraic solution in terms of Γ-functions, we are able to transform the solution from the Mellin space to the time domain. The mathematical link is a Barnes-type integral [8] defining either the inverse Mellin transform of the  Fox's H-functions. Knowing this connection it is straightforward to write down the Fox's H-function from the solution in Mellin space. The graphical representation of the result shows that a nonexponential decay occurs.

Plot[
      Evaluate[{Exp[-t], ϕ[t]} /. shh  /. {q -> 1/2, 
            ϕ0 -> 1, τ -> 1}], {t, 0, 13}, 
      AxesLabel -> {"t/τ", "ϕ/ϕ0"}, 
      PlotStyle 
        -> {RGBColor[0.996109,  0,  0], RGBColor[0,  0,  0.996109]}]

Comparing the exponential behavior with a power law behavior. The upper curve represents the solution of the fractional relaxation equation of Weyl type.

The figure above shows that an exponential (red) decays much faster than the solution of the Weyl equation (blue). The asymptotic behavior of the derived solution is given by

Series[ϕ[t] /. shh, {t, ∞, 1}]

confirming the power law behavior of anomalous relaxation process.

Conclusions

We demonstrated that Mathematica is able to handle fractional derivatives if the operations are defined in an appropriate way. We implemented with a minimal number of functions the Riemann–Liouville and the Weyl calculi. The achievement of both calculi allows us to calculate fractional integrals and fractional derivatives. We also demonstrated by several examples that the calculation of either fractional derivatives or integrals will result in special functions. This is true just for simple functions. We also discussed two examples applying the fractional calculi to solve physical problems. We solved the tautochrone problem originated by Abel and a fractional relaxation equation. Both physical problems enable us to solve the related integral equations.

REFERENCES

[1] S.F. Lacroix, Traité du Calcul Différentiel et du Calcul Intégral, 2nd ed., Vol.3,  409-410. Courcier, Paris (1819).

[2] J. Liouville, Mémoiresur le calcul des différentielles à indices quelconques, J. École Polytech., 13, 71-162, (1832).

[3] K.S. Miller and B. Ross, An Introduction to the Fractional Calculus and Fractional Differential Equations, John Wiley & Sons, Inc., New York, (1993).

[4] K.B. Oldham and J. Spanier, The Fractional Calculus, Academic Press, New York, (1974).

[5] G.F.B. Riemann, Gesammelte Werke, 353-366, Teubner, Leipzig, (1892).

[6] H. Weyl, Bemerkungen zum Begriff des Differentialquotienten gebrochener Ordnung, Vierteljahresschr. Naturforsch. Ges. Zürich, 62, 296-302, (1917).

[7] C. Fox, The G and H Functions as Symmetrical Fourier Kernels, Trans. Am. Math. Soc., 98, 395-429, (1961).

[8] B.L.J. Braaksma, Asymptotic Expansions and Analytic Continuations for a Class of Barnes-Integrals, Compos. Math. 15, 239-341, (1964).

ABOUT THE AUTHOR

Gerd Baumann is a faculty member in the Department of Mathematical Physics at the University of Ulm, Germany. He is the author of Mathematica in Theoretical Physics (1995) and Symmetry Analysis of Differential Equations with Mathematica (in press), both published by Springer-Verlag. 

Gerd Baumann
Department of Mathematical Physics
University of Ulm
Albert-Einstein-Allee 11
D-89069 Ulm
Germany
e-mail:Gerd.Baumann@physik.uni-ulm.de

ELECTRONIC SUBSCRIPTIONS

Included in the distribution for each electronic subscription is the file FractionalCalculus.nb,containing Mathematica code for the material described in this article.