New in Fractional Differentiation

Oleg Marichev, Wolfram Research
Paco Jain, Wolfram Research

ABSTRACT

In the scientific literature, one sees a lot of approaches to fractional order integro-differentiation. In this talk, we present numerous essential details of the most natural Riemann–Liouville–Hadamard construction for fractional differentiation, which has been published in the Wolfram Function Repository as ResourceFunction[“FractionalOrderD”]. Together they form a new picture of fractional differentiation for analytic functions of complex variables. Elements of this picture, to be described in this talk, include:
  1) separation of symbolic positive integer n-th order simplest cases of fractional differentiation from generic result;
  2) developing formulas for about 500 basic single functions (like Sin[z], BesselJ[z,x], BesselJ[v,z],...);
  3) deriving generic rules for product, composition, ratio, inverse function, series, etc.;
  4) operation in the full complex plane with all functions involving so-called differential constants (like (Sqrt[(z^2)])/z) and    piecewise construction including logarithmic situations previously ignored in the literature;
  5) widely using representations of functions through the more general MeijerG and the new FoxH functions;
  6) developing formulas for new Mathematica Heun and Lamé functions, for Carlson elliptic integrals and for the FoxH function.

The talk is devoted to new operation FractionalOrderD[f[z],{z,α}]. FractionalOrderD is analytical by order α operator, which for integer positive α coincides with usual α-th order derivative D[f[z],{z,α}], for integer negative α coincides with repeatable -α times indefinite integral (Underscript[∫∫…∫, (Underscript[⏟, -α])])f[z]dz and for α==0 coincides with initial function f[z]. The building of such operator was not easy problem, because differentiation of basic for Taylor series power function (z^γ) always gives power function with constant coefficient and different power, but integration of (z^γ) can produce log(z) (for negative integer γ). We have solved this problem by extension well known from literature generic result (((Gamma[γ+1]) (z^(γ-α)))/(Gamma[γ-α+1]))for integro-differentiation of (z^γ) to value  (((((-1))^(γ-1)) ((PolyGamma[-γ])-(PolyGamma[γ-α+1])+log(z)))/((-γ-1)! (Gamma[γ-α+1])))(z^(γ-α)) , which should be applied for negative integer γ. More than 100 basic functions in reality are the infinite sums by k of  (Subscript[c, k])(z^(γ+k)) with constant (Subscript[c, k]). It allowed us to build fractional derivatives of order α for more than 85800 test-cases, including basic functions and their compositions. Their values coincide with results based on the most popular Riemann-Liouville approach to fractional integro differentiation.

Previously...

See our Wolfram Tech Conference 2020 talk notebook and screencast.

About 30 years ago we published a large (about 1000 pages with 3000 references) monograph
Fractional Integrals and Derivatives. Theory and Applications, 1993 by S.G.Samko, A.A.Kilbas & O.I.Marichev  which included descriptions of different approaches to fractional order integro-differentiation since Leibniz’ time (1695). We started  with the realization of the main Riemann-Liouville-Hadamard approach for analytic functions of complex variables in a neighbourhood of zero, which is not deeply described in literature. It gave many new results, which allow us to say that we have “new fractional differentiation”.

Preliminary

In our previous talk "Fractional integro-differentiation in Wolfram language" we presented some information about history of fractional  integro-differentiation and some our results in this direction (in particular our program for FractionalOrderD[f[z], {z,α}] was successfully tested with 85000+ cases).  Now we describe the most important aspects in our research and program, which demonstrate how many can be done if one use modern computer technology for building “Fractional calculus”.

We start with the “simplest” mathematical functions, involving only one or two letters   {(1/z),(Sqrt[z]),(z^b),(a^z),(e^z),(z^z)}. Next, we see “named functions” with head-titles like  {log(z),sin(z),(BesselJ[v, z]),(BesselJ[z, b])} . At last we see “composed functions” {(Sqrt[(z^2)]),((((z^a)))^b),(a^(z^c)),(sin^(-1))((z^3))}.

{(1/z),(Sqrt[z]),(z^b),(a^z),(E^z),(z^z)}  (*simplest functions*)

{Log[z], Sin[z], BesselJ[v,z], BesselJ[z,b]} (*named functions*)

{(Sqrt[(z^2)]),((((z^a)))^b),(a^(z^c)),ArcSin[(z^3)]}  (* composed functions *)

If we apply differentiation or indefinite integration, we find that not each integral and even derivative can be evaluated and in many cases we see special functions as result of integration.

List14={(1/z),(Sqrt[z]),(z^b),(a^z),(E^z),(z^z),Log[z], Sin[z], BesselJ[v,z], BesselJ[z,b],(Sqrt[(z^2)]),((((z^a)))^b),(a^(z^c)),ArcSin[(z^3)]};

{D[List14,z],Integrate[List14,z]}

Our new operation FractionalOrderD evaluates repeatable α-integer order derivatives and indefinite integrals and extends results for arbitrary complex or real order α:

{#,[#,{z,α}]}&/@List14//Simplify//TableForm

It was possible to make because we apply FractionalOrderD operation to each term of Taylor series expansions of all functions near zero:

[f[z],{z,α}]==(∑_(n=0)^∞)(Subscript[c, n])[(z^(ν+n)),{z,α}]/;f[z]==(z^ν)(∑_(n=0)^∞)(Subscript[c, n])(z^n)

f[z]==(∑_(n=0)^∞)((Derivative[n][f][0])/(n!))(z^n)         (*Exp[z]*)

f[z]==(z^ν)(∑_(n=0)^∞)(Subscript[c, n])(z^n) /;(Subscript[c, n])==((Derivative[n][f][0])/(n!))     (* BesselJ[ν,z] *)

Above formula works nicely for functions like Exp[z], BesselJ[ν,z], but sometimes series expansions include Log function (as in BesselK[0,z])

f[z]==(z^ν)Log[z](∑_(n=0)^∞)(Subscript[c, n])(z^n) +ϕ[z]/;(Subscript[c, n])==((Derivative[n][f][0])/(n!))  (* BesselK[0,z]*)

It means, that we must evaluate for arbitrary γ, α and integer k==0,1,2,… the following values

{[(z^γ),{z,α}],[(z^γ)((Log[z])^k),{z,α}]} (*/;k==0,1,2,…*)

This problem was solved, for example:

{[(z^γ),{z,α}],[(z^γ)Log[z],{z,α}]}

1) Definition of FractionalOrderD through Riemann-Liouville-Hadamard integral

By definition of FractionalOrderD[f[z],{z,α}] we have

[f[z],{z,α}]==
({{\[Piecewise], ({{Derivative[α][f][z], α∈Integers&&α>0}, {Integrate[(((((z-t))^(n-α-1)) Derivative[n][f][t])/(Gamma[n-α])),{t,0,z},GenerateConditions->False], n==1+Floor[α]&&Re[α]>0}, {Integrate[(((((z-t))^(-α-1)) f[t])/(Gamma[-α])),{t,0,z},GenerateConditions->False], α∈Integers&&α<0}, {f[z], α==0}, {D[Integrate[((((z-t))^(-α))/(Gamma[1-α]))f[t], {t,0,z},GenerateConditions->False],z], Re[α]==0&&Im[α]≠0}, {(∫_0^z)(((((z-t))^(-α-1)) f[t])/(Gamma[-α]))\[DifferentialD]t, Re[α]<0}, {0, True}})}})

Above we separated cases of symbolic positive integer n-th order derivatives from generic result integro-differentiation of fractional order.

In particular for α == -1, -2, … we have

Integrate[((((z-t))^(-α-1))/(Gamma[-α]))f[t], {t, 0,z},GenerateConditions->False]==(Underscript[∫∫…∫, (Underscript[⏟, -α])])f[z]\[DifferentialD]z/;α∈Integers&&α<0

With[{α=-3},(∫_0^z)((((z-t))^(-α-1))/(Gamma[-α]))Exp[t]\[DifferentialD]t==(∫_0^z)((∫_0^(Subscript[t, 3]))((∫_0^(Subscript[t, 2]))Exp[(Subscript[t, 1])]\[DifferentialD](Subscript[t, 1]))\[DifferentialD](Subscript[t, 2]))\[DifferentialD](Subscript[t, 3])==Exp[z]-1-z-((z^2)/2)]

With[{α=-1/2},(∫_0^z)((((z-t))^(-α-1))/(Gamma[-α]))Exp[t]\[DifferentialD]t==[Exp[z],{z,α}]==(E^z) Erf[(Sqrt[z])]]//FunctionExpand

(E^z) (1+((((2 (E^(-z)))/(Sqrt[z]))-2 (Sqrt[π]) (1-Erf[(Sqrt[z])]))/(2 (Sqrt[π]))))//Simplify

With[{α=1/2},Integrate[((((z-t))^(-α-1))/(Gamma[-α]))Exp[t],{t,0,z},GenerateConditions->False]==
[Exp[z],{z,α}]==(1/((Sqrt[π]) (Sqrt[z])))+(E^z) Erf[(Sqrt[z])]]//FunctionExpand//FunctionExpand//Simplify

2) Formulas for about 500 basic single functions (like Sin[z], BesselJ[z,x], BesselJ[v,z],...)

Below you see 500 functions, analytical in variable z, for which we can evaluate FractionalOrderD operator. We should mention, that differentiation by parameters of these functions as rule can not produce well known functions. But we can say the following:
“ The first derivative with respect to an upper “parameter” (Subscript[a, k]), and all derivatives of symbolic integer order m with respect to a “lower” parameter (Subscript[b, k]) of the generalized hypergeometric function (Subscript[, p])(Subscript[F, q])((Subscript[a, 1]),…,(Subscript[a, p]);(Subscript[b, 1]),…,(Subscript[b, q]);z), can be expressed in terms of the Kampé de Fériet hypergeometric function of two variables” (see blog of Oleg Marichev & Yury Brychkov "New Derivatives of the Bessel Functions Have Been Discovered with the Help of the Wolfram Language!"

List500={(a^z),(E^z),(Sqrt[z]),(z^b),AiryAi[z],AiryAiPrime[z],AiryBi[z],AiryBiPrime[z],AlternatingFactorial[z],AngerJ[a,z],AngerJ[z,a],AngerJ[a,b,z],AngerJ[a,z,b],AngerJ[z,a,b],AppellF1[a,b1,b2,z,x,y],AppellF1[a,b1,z,d,x,y],AppellF1[a,b1,b2,c,x,z],AppellF1[a,b1,b2,c,z,y],AppellF1[a,z,c,d,x,y],AppellF1[z,b1,b2,d,x,y],ArcCos[z],ArcCosh[z],ArcCot[z],ArcCoth[z],ArcCsc[z],ArcCsch[z],ArcSec[z],ArcSech[z],ArcSin[z],ArcSinh[z],ArcTan[z],ArcTan[a,z],ArcTan[z,a],ArcTanh[z],ArithmeticGeometricMean[a,z],ArithmeticGeometricMean[b,z],BarnesG[z],BellB[ν,z],BernoulliB[a,z],BesselI[a,z],BesselI[z,a],BesselJ[a,z],BesselJ[z,a],BesselK[a,z],BesselK[z,a],BesselY[a,z],BesselY[z,a],Beta[a,z],Beta[z,b],Beta[a,b,z],Beta[a,z,c],Beta[z,a,b],Beta[a,b,c,z],Beta[a,b,z,d],Beta[c,z,a,b],Beta[z,c,a,b],BetaRegularized[a,b,z],BetaRegularized[a,z,c],BetaRegularized[z,a,b],BetaRegularized[a,b,c,z],BetaRegularized[a,b,z,d],BetaRegularized[c,z,a,b],BetaRegularized[z,c,a,b],Binomial[a,z],Binomial[z,b],CarlsonRC[x,z],CarlsonRC[z,y],CarlsonRD[x,z,a],CarlsonRD[z,y,a],CarlsonRD[x,y,z],CarlsonRE[x,z],CarlsonRE[z,y],CarlsonRF[x,z,a],CarlsonRF[z,y,a],CarlsonRF[x,y,z],CarlsonRG[z,y,a],CarlsonRG[x,z,a],CarlsonRG[x,y,z],CarlsonRJ[z,y,a,b],CarlsonRJ[x,z,a,b],CarlsonRJ[x,y,z,b],CarlsonRJ[x,y,a,z],CarlsonRK[x,z],CarlsonRK[z,y],CarlsonRM[x,z,a],CarlsonRM[z,y,a],CarlsonRM[x,y,z],CatalanNumber[z],ChebyshevT[a,z],ChebyshevT[z,b],ChebyshevU[a,z],ChebyshevU[z,b],Cos[z],Cosh[z],CoshIntegral[z],CosIntegral[z],Cot[z],Coth[z],Csc[z],Csch[z],DawsonF[z],EllipticE[z],EllipticE[a,z],EllipticE[z,b],EllipticF[a,z],EllipticF[z,b],EllipticK[z],EllipticNomeQ[z],EllipticPi[a,z],EllipticPi[z,b],EllipticPi[a,b,z],EllipticPi[a,z,b],EllipticPi[z,b,c],EllipticTheta[2,z],EllipticTheta[3,z],EllipticTheta[4,z],EllipticTheta[1,z,a],EllipticTheta[1,ν,z],EllipticTheta[2,z,a],EllipticTheta[2,ν,z],EllipticTheta[3,z,a],EllipticTheta[3,ν,z],EllipticTheta[4,z,a],EllipticTheta[4,ν,z],EllipticThetaPrime[1,z],EllipticThetaPrime[1,z,a],EllipticThetaPrime[1,ν,z],EllipticThetaPrime[2,z,a],EllipticThetaPrime[2,ν,z],EllipticThetaPrime[3,z,a],EllipticThetaPrime[3,ν,z],EllipticThetaPrime[4,z,a],EllipticThetaPrime[4,ν,z],Erf[z],Erf[a,z],Erf[z,b],Erfc[z],Erfi[z],EulerE[n,z],ExpIntegralE[a,z],ExpIntegralE[z,b],ExpIntegralEi[z],z!,z!!,FactorialPower[a,z],FactorialPower[z,b],FactorialPower[a,b,z],FactorialPower[a,z,c],FactorialPower[z,b,c],Fibonacci[z],Fibonacci[a,z],Fibonacci[z,b],FoxH[{{{(Subscript[a, 1]),(Subscript[α, 1])},{(Subscript[a, 2]),(Subscript[α, 1])}},{{(Subscript[a, 3]),(Subscript[α, 3])}}},{{{(Subscript[b, 1]),(Subscript[β, 1])}},{{(Subscript[b, 2]),(Subscript[β, 2])}}},z],FoxH[{Table[{(Subscript[a, i]),(Subscript[α, i])},{i,1,n}],Table[{(Subscript[a, i]),(Subscript[α, i])},{i,n+1,p}]},{Table[{(Subscript[b, i]),(Subscript[β, i])},{i,1,m}],Table[{(Subscript[b, i]),(Subscript[β, i])},{i,m+1,q}]},z],FresnelC[z],FresnelF[z],FresnelG[z],FresnelS[z],Gamma[z],Gamma[a,z],Gamma[z,a],Gamma[a,b,z],Gamma[a,z,b],Gamma[z,b,c],GammaRegularized[a,z],GammaRegularized[z,b],GammaRegularized[a,b,z],GammaRegularized[a,z,b],GammaRegularized[z,b,c],GegenbauerC[a,b,z],GegenbauerC[a,z,b],GegenbauerC[z,a,b],Gudermannian[z],HankelH1[a,z],HankelH1[z,a],HankelH2[a,z],HankelH2[z,a],HarmonicNumber[z],HarmonicNumber[a,z],HarmonicNumber[z,a],Haversine[z],HermiteH[a,z],HermiteH[z,a],HeunB[q,α,γ,δ,ϵ,z],HeunBPrime[q,α,γ,δ,ϵ,z],HeunC[q,α,γ,δ,ϵ,z],HeunCPrime[q,α,γ,δ,ϵ,z],HeunD[q,α,γ,δ,ϵ,z],HeunDPrime[q,α,γ,δ,ϵ,z],HeunG[a,q,α,β,γ,δ,z],HeunGPrime[a,q,α,β,γ,δ,z],HeunT[q,α,γ,δ,ϵ,z],HeunTPrime[q,α,γ,δ,ϵ,z],HurwitzLerchPhi[a,b,z],HurwitzLerchPhi[a,z,b],HurwitzLerchPhi[z,a,b],HurwitzLerchPhi[z,b,c],HurwitzZeta[a,z],HurwitzZeta[z,a],Hyperfactorial[z],Hypergeometric0F1[a,z],Hypergeometric0F1[z,a],Hypergeometric0F1Regularized[a,z],Hypergeometric0F1Regularized[z,a],Hypergeometric1F1[a,b,z],Hypergeometric1F1[a,z,b],Hypergeometric1F1[z,b,c],Hypergeometric1F1Regularized[a,b,z],Hypergeometric1F1Regularized[a,z,b],Hypergeometric1F1Regularized[z,b,c],Hypergeometric2F1[a,b,c,z],Hypergeometric2F1[a,b,z,c],Hypergeometric2F1[a,z,b,c],Hypergeometric2F1Regularized[a,b,c,z],Hypergeometric2F1Regularized[a,b,z,c],Hypergeometric2F1Regularized[a,z,b,c],HypergeometricPFQ[{},{(Subscript[b, 1]),(Subscript[b, 2]),(Subscript[b, 3]),(Subscript[b, 4])},z],HypergeometricPFQ[{(Subscript[a, 1]),(Subscript[a, 2])},{(Subscript[b, 1]),(Subscript[b, 2]),(Subscript[b, 3])},z],HypergeometricPFQ[Table[(Subscript[a, k]),{k,1,p}],Table[(Subscript[b, k]),{k,1,q}],z],HypergeometricPFQRegularized[{(Subscript[a, 1])},{(Subscript[b, 1]),(Subscript[b, 2]),(Subscript[b, 3])},z],HypergeometricPFQRegularized[Table[(Subscript[a, k]),{k,1,p}],Table[(Subscript[b, k]),{k,1,q}],z],HypergeometricU[a,b,z],HypergeometricU[a,z,b],HypergeometricU[z,a,b],InverseBetaRegularized[a,b,z],InverseBetaRegularized[a,z,c],InverseBetaRegularized[z,a,b],InverseBetaRegularized[a,b,c,z],InverseBetaRegularized[a,b,z,d],InverseBetaRegularized[x,z,a,b],InverseBetaRegularized[z,b,c,d],InverseEllipticNomeQ[z],InverseErf[z],InverseErf[x,z],InverseErf[z,b],InverseErfc[z],InverseGammaRegularized[a,z],InverseGammaRegularized[z,b],InverseGammaRegularized[a,x,z],InverseGammaRegularized[a,z,c],InverseGammaRegularized[z,b,c],InverseGudermannian[z],InverseHaversine[z],InverseJacobiCD[a,z],InverseJacobiCD[z,m],InverseJacobiCN[a,z],InverseJacobiCN[z,m],InverseJacobiCS[a,z],InverseJacobiCS[z,m],InverseJacobiDC[a,z],InverseJacobiDC[z,m],InverseJacobiDN[a,z],InverseJacobiDN[z,m],InverseJacobiDS[a,z],InverseJacobiDS[z,m],InverseJacobiNC[a,z],InverseJacobiNC[z,m],InverseJacobiND[a,z],InverseJacobiND[z,m],InverseJacobiNS[a,z],InverseJacobiNS[z,m],InverseJacobiSC[a,z],InverseJacobiSC[z,m],InverseJacobiSD[a,z],InverseJacobiSD[z,m],InverseJacobiSN[a,z],InverseJacobiSN[z,m],InverseWeierstrassP[z,{g2,g3}],JacobiAmplitude[a,z],JacobiAmplitude[z,m],JacobiCD[a,z],JacobiCD[z,m],JacobiCN[a,z],JacobiCN[z,m],JacobiCS[a,z],JacobiCS[z,m],JacobiDC[a,z],JacobiDC[z,m],JacobiDN[a,z],JacobiDN[z,m],JacobiDS[a,z],JacobiDS[z,m],JacobiEpsilon[a,z],JacobiEpsilon[z,b],JacobiNC[a,z],JacobiNC[z,m],JacobiND[a,z],JacobiND[z,m],JacobiNS[a,z],JacobiNS[z,m],JacobiP[a,b,c,z],JacobiP[a,b,z,c],JacobiP[a,z,b,c],JacobiP[z,a,b,c],JacobiSC[a,z],JacobiSC[z,m],JacobiSD[a,z],JacobiSD[z,m],JacobiSN[a,z],JacobiSN[z,m],JacobiZeta[a,z],JacobiZeta[z,b],JacobiZN[a,z],JacobiZN[z,b],KelvinBei[0,z],KelvinBei[a,z],KelvinBei[z,a],KelvinBer[0,z],KelvinBer[a,z],KelvinBer[z,a],KelvinKei[0,z],KelvinKei[a,z],KelvinKei[z,a],KelvinKer[0,z],KelvinKer[a,z],KelvinKer[z,a],KleinInvariantJ[z],LaguerreL[a,z],LaguerreL[z,a],LaguerreL[a,b,z],LaguerreL[a,z,b],LaguerreL[z,a,b],LameC[ν,j,z,m],LameCPrime[ν,j,z,m],LameS[ν,j,z,m],LameSPrime[ν,j,z,m],LegendreP[z,a],LegendreP[ν,z],LegendreP[a,b,z],LegendreP[a,z,b],LegendreP[z,a,b],LegendreP[a,b,2,z],LegendreP[a,b,3,z],LegendreP[a,z,2,b],LegendreP[a,z,3,b],LegendreP[z,a,2,b],LegendreP[z,a,3,b],LegendreQ[z,a],LegendreQ[ν,z],LegendreQ[a,b,z],LegendreQ[a,z,b],LegendreQ[z,a,b],LegendreQ[a,b,2,z],LegendreQ[a,b,3,z],LegendreQ[a,z,2,b],LegendreQ[a,z,3,b],LegendreQ[z,a,2,b],LegendreQ[z,a,3,b],LerchPhi[a,b,z],LerchPhi[a,z,b],LerchPhi[z,a,b],LerchPhi[z,b,c],((Log[a])/(Log[z])),Log[z],LogBarnesG[z],LogGamma[z],LogIntegral[z],LogisticSigmoid[z],LucasL[z],LucasL[a,z],LucasL[z,b],MarcumQ[a,b,z],MarcumQ[a,z,c],MarcumQ[a,b,c,z],MarcumQ[a,b,z,d],MarcumQ[a,z,c,d],MathieuC[a,q,z],MathieuCPrime[a,q,z],MathieuS[a,q,z],MathieuSPrime[a,q,z],MeijerG[{{(Subscript[a, 1])},{(Subscript[a, 2]),(Subscript[a, 3])}},{{(Subscript[b, 1]),(Subscript[b, 2])},{(Subscript[b, 3])}},z],MeijerG[{{(Subscript[a, 1])},{(Subscript[a, 2])}},{{(Subscript[b, 1]),(Subscript[b, 2])},{(Subscript[b, 3])}},z],MeijerG[{{(Subscript[a, 1])},{(Subscript[a, 2]),(Subscript[a, 3]),(Subscript[a, 4]),(Subscript[a, 5])}},{{(Subscript[b, 1]),(Subscript[b, 2])},{(Subscript[b, 3]),(Subscript[b, 4]),(Subscript[b, 5])}},z],MeijerG[{Table[(Subscript[a, i]),{i,1,n}],Table[(Subscript[a, n+j]),{j,1,-n+p}]},{Table[(Subscript[b, k]),{k,1,m}],Table[(Subscript[b, m+l]),{l,1,-m+q}]},z,r],MittagLefflerE[a,z],MittagLefflerE[z,b],MittagLefflerE[a,b,z],MittagLefflerE[a,z,b],MittagLefflerE[z,a,b],ModularLambda[z],Multinomial[a,b,z],Multinomial[a,b,c,z],NevilleThetaC[a,z],NevilleThetaC[z,b],NevilleThetaD[a,z],NevilleThetaD[z,b],NevilleThetaN[a,z],NevilleThetaN[z,b],NevilleThetaS[a,z],NevilleThetaS[z,b],NorlundB[ν,z],NorlundB[ν,z,b],NorlundB[ν,γ,z],OwenT[a,z],OwenT[z,a],ParabolicCylinderD[a,z],ParabolicCylinderD[z,a],Pochhammer[a,z],Pochhammer[z,b],PolyGamma[0,z],PolyGamma[a,z],PolyLog[2,z],PolyLog[a,z],PolyLog[z,a],PolyLog[2,b,z],PolyLog[a,b,z],PolyLog[z,a,b],ProductLog[z],ProductLog[a,z],ProductLog[k,z],QBinomial[a,z,c],QBinomial[z,b,c],QFactorial[z,b],QGamma[z,b],QHypergeometricPFQ[{(Subscript[a, 1]),…,(Subscript[a, r])},{(Subscript[b, 1]),…,(Subscript[b, s])},q,z],QPochhammer[z,b],QPochhammer[z,z],QPolyGamma[0,z,a],QPolyGamma[0,z,b],QPolyGamma[a,z,a],QPolyGamma[a,z,c],RamanujanTauL[z],RamanujanTauTheta[z],RamanujanTauZ[z],RiemannSiegelTheta[z],RiemannSiegelZ[z],ScorerGi[z],ScorerGiPrime[z],ScorerHi[z],ScorerHiPrime[z],Sec[z],Sech[z],Sin[z],Sinc[z],Sinh[z],SinhIntegral[z],SinIntegral[z],SphericalBesselJ[a,z],SphericalBesselJ[z,a],SphericalBesselY[a,z],SphericalHankelH1[a,z],SphericalHankelH1[z,a],SphericalHankelH2[a,z],SphericalHankelH2[z,a],SphericalHarmonicY[a,b,c,z],SphericalHarmonicY[a,b,z,c],SphericalHarmonicY[a,b,z,d],SphericalHarmonicY[a,z,c,d],SphericalHarmonicY[z,a,b,c],SpheroidalPS[ν,μ,γ,z],SpheroidalPS[ν,μ,2,γ,z],SpheroidalPSPrime[ν,μ,γ,z],SpheroidalPSPrime[ν,μ,2,γ,z],SpheroidalQS[ν,μ,γ,z],SpheroidalQS[ν,μ,2,γ,z],SpheroidalQSPrime[ν,μ,γ,z],SpheroidalQSPrime[ν,μ,2,γ,z],SpheroidalS1[ν,μ,γ,z],SpheroidalS1Prime[ν,μ,γ,z],SpheroidalS2[ν,μ,γ,z],SpheroidalS2Prime[ν,μ,γ,z],StruveH[z,a],StruveH[ν,z],StruveL[z,a],StruveL[ν,z],Subfactorial[z],(z^(1/3)),(z^(1/5)),(z^(1/7)),Tan[z],Tanh[z],WeberE[z,a],WeberE[ν,z],WeberE[a,z,b],WeberE[z,a,b],WeberE[ν,a,z],WhittakerM[a,b,z],WhittakerM[a,z,b],WhittakerM[z,a,b],WhittakerW[a,b,z],WhittakerW[a,z,b],WhittakerW[z,a,b],ZernikeR[a,b,z],Zeta[z],Zeta[a,z],Zeta[z,b],(InactiveHead[DirichletBeta])[z],(InactiveHead[DirichletEta])[z],(InactiveHead[DirichletLambda])[z],(InactiveHead[ErlangB])[a,z],(InactiveHead[ErlangB])[z,a],(InactiveHead[ErlangC])[a,z],(InactiveHead[ErlangC])[z,a],(InactiveHead[GegenbauerC])[a,z],(InactiveHead[GegenbauerC])[z,a],(InactiveHead[Log])[a,z],(InactiveHead[Log])[z,a],(InactiveHead[Log10])[z],(InactiveHead[Log2])[z],(InactiveHead[RiemannXi])[z]};

Length[List500]

data={#1,[#1,{z,α}]}&/@List500;

3) Generic rules for product, composition, ratio, inverse function, series, etc.

We should mention, that after making FractionalOrderD of above 500 functions user is able to use general abstract rules which are described by the following formulas for product, composition, ratio, inverse function & etc.

{#,[#,{z,α}]}&/@{(f^((-1)))[z],((f[z])^(g[z])),((f[z])^((g[z])^(h[z]))),z f[z],(z^2) f[z],(((a z))^b) ((Log[d f[z]])^c),(c^(b ((f[z])^a)+d)),(c^(b f[(z^a)]+d)),(c^(b f[(a^z)]+d)),(c^(b (a^(f[z]))+d)),(1/(f[z])),((f[z])^a),((f[z])/(g[z])),f[z] g[z],f[z] g[z] h[z],f[g[z]],f[g[h[z]]]}//Simplify//TableForm//Quiet

4) Operations with differential constants (like (Sqrt[(z^2)])/z, Log[(z^2)]-2Log[z])

As we see, in literature practically nobody try to investigate so called (by us) “differential constants”, which are not constants everywhere, but they are different constants on some domains of complex plane with discontinuity on some lines, for example,

{((Sqrt[(z^2)])/z)==({{\[Piecewise], ({{1, -(π/2)<Arg[z]≤(π/2)}, {-1, True}})}}),Log[(z^2)]-2 Log[z]==({{\[Piecewise], ({{0, -(π/2)<Arg[z]<=(π/2)}, {2π I, -π<Arg[z]≤(π/2)}, { -2π I, True}})}})};

Table[With[{z=Random[]Exp[Pi I k/4]},{((Sqrt[(z^2)])/z)-({{\[Piecewise], ({{1, -(π/2)<Arg[z]≤(π/2)}, {-1, True}})}}),Log[(z^2)]-2 Log[z]-({{\[Piecewise], ({{0, -(π/2)<Arg[z]<=(π/2)}, {2π I, -π<Arg[z]≤(π/2)}, { -2π I, True}})}})}],{k,0,7}]//Chop

Usual derivatives of these functions give zero:

D[{((Sqrt[(z^2)])/z),-2 Log[z]+Log[(z^2)]},z]=={0,0}//Simplify

FractionalOrderD provides the following results, which saves components (Sqrt[(z^2)]) &  (z^2) log==(-i z) log+(i z) log:

{#,[#,{z,α}]}&/@{((Sqrt[(z^2)])/z),-2 Log[z]+Log[(z^2)]}//FunctionExpand//Simplify//TableForm

FullSimplify[(Log[-I z]+Log[I z])-Log[(z^2)]]

5) Using representations through MeijerG and FoxH functions

For evaluation FractionalOrderD we widely use opportunity to present functions through “superfunction” MeijerG (now we have more general FoxH function). Such representations can be done for below 137 functions from mentioned above List500. Below you see corresponding tables, where we use existing version of MeijerGReduce, which already essential improved but still its improvement MeijerGForm2 is not implemented.

MeijerG137={(a^z),(E^z),(Sqrt[z]),(z^b),AiryAi[z],AiryAiPrime[z],AiryBi[z],AiryBiPrime[z],AngerJ[a,z],AngerJ[a,b,z],ArcCos[z],ArcCosh[z],ArcCot[z],ArcCoth[z],ArcCsc[z],ArcCsch[z],ArcSec[z],ArcSech[z],ArcSin[z],ArcSinh[z],ArcTan[z],ArcTan[a,z],ArcTan[z,a],ArcTanh[z],BesselI[a,z],BesselJ[a,z],BesselK[a,z],BesselY[a,z],Beta[z,a,b],Beta[c,z,a,b],Beta[z,c,a,b],BetaRegularized[z,a,b],BetaRegularized[c,z,a,b],BetaRegularized[z,c,a,b],CarlsonRC[x,z],CarlsonRC[z,y],CarlsonRE[x,z],CarlsonRE[y,z],CarlsonRK[x,z],CarlsonRK[y,z],ChebyshevT[a,z],ChebyshevU[a,z],Cos[z],Cosh[z],CoshIntegral[z],CosIntegral[z],DawsonF[z],EllipticE[z],EllipticK[z],Erf[z],Erf[a,z],Erf[z,b],Erfc[z],Erfi[z],ExpIntegralE[a,z],ExpIntegralEi[z],Fibonacci[z],Fibonacci[a,z],FresnelC[z],FresnelF[z],FresnelG[z],FresnelS[z],Gamma[a,z],Gamma[a,b,z],Gamma[a,z,b],GammaRegularized[a,z],GammaRegularized[a,b,z],GammaRegularized[a,z,b],GegenbauerC[a,b,z],HankelH1[a,z],HankelH2[a,z],Haversine[z],HermiteH[a,z],Hypergeometric0F1[a,z],Hypergeometric0F1Regularized[a,z],Hypergeometric1F1[a,b,z],Hypergeometric1F1Regularized[a,b,z],Hypergeometric2F1[a,b,c,z],Hypergeometric2F1Regularized[a,b,c,z],HypergeometricPFQ[{},{(Subscript[b, 1]),(Subscript[b, 2])},z],HypergeometricPFQ[Table[(Subscript[a, k]),{k,1,3}],Table[(Subscript[b, k]),{k,4}],z],HypergeometricPFQ[Table[(Subscript[a, k]),{k,1,p}],{(Subscript[b, 1]),(Subscript[b, 2]),(Subscript[b, 3]),(Subscript[b, 4])},z],HypergeometricPFQ[{(Subscript[a, 1]),(Subscript[a, 2])},Table[(Subscript[b, k]),{k,q}],z],HypergeometricPFQRegularized[{(Subscript[a, 1])},{(Subscript[b, 1]),(Subscript[b, 2]),(Subscript[b, 3])},z],HypergeometricPFQRegularized[Table[(Subscript[a, k]),{k,1,p}],Table[(Subscript[b, k]),{k,1,q}],z],HypergeometricPFQRegularized[Table[(Subscript[a, k]),{k,1,p}],{(Subscript[b, 1]),(Subscript[b, 2])},z],HypergeometricPFQRegularized[{(Subscript[a, 1])},Table[(Subscript[b, k]),{k,q}],z],HypergeometricU[a,b,z],InverseHaversine[z],KelvinBei[0,z],KelvinBei[a,z],KelvinBer[0,z],KelvinBer[a,z],KelvinKei[0,z],KelvinKei[a,z],KelvinKer[0,z],KelvinKer[a,z],LaguerreL[a,z],LaguerreL[a,b,z],LegendreP[ν,z],LegendreP[a,b,z],LegendreP[a,b,2,z],LegendreP[a,b,3,z],LegendreQ[ν,z],LegendreQ[a,b,z],LegendreQ[a,b,2,z],LegendreQ[a,b,3,z],Log[z],LucasL[z],LucasL[a,z],MeijerG[{{(Subscript[a, 1])},{(Subscript[a, 2]),(Subscript[a, 3])}},{{(Subscript[b, 1]),(Subscript[b, 2])},{(Subscript[b, 3])}},z],MeijerG[{{(Subscript[a, 1])},{(Subscript[a, 2])}},{{(Subscript[b, 1]),(Subscript[b, 2])},{(Subscript[b, 3])}},z],MeijerG[{{(Subscript[a, 1])},{(Subscript[a, 2]),(Subscript[a, 3]),(Subscript[a, 4]),(Subscript[a, 5])}},{{(Subscript[b, 1]),(Subscript[b, 2])},{(Subscript[b, 3]),(Subscript[b, 4]),(Subscript[b, 5])}},z],MeijerG[{Table[(Subscript[a, i]),{i,1,n}],Table[(Subscript[a, n+j]),{j,1,-n+p}]},{Table[(Subscript[b, k]),{k,1,m}],Table[(Subscript[b, m+l]),{l,1,-m+q}]},z,r],MeijerG[{{(Subscript[a, 1])},{(Subscript[a, 2])}},{{(Subscript[b, 1]),(Subscript[b, 2])},{(Subscript[b, 3])}},z,(1/3)],ParabolicCylinderD[a,z],PolyLog[2,z],PolyLog[a,z],ScorerGi[z],ScorerGiPrime[z],ScorerHi[z],ScorerHiPrime[z],Sin[z],Sinc[z],Sinh[z],SinhIntegral[z],SinIntegral[z],SphericalBesselJ[a,z],SphericalBesselY[a,z],SphericalHankelH1[a,z],SphericalHankelH2[a,z],StruveH[ν,z],StruveL[ν,z],WeberE[ν,z],WeberE[ν,a,z],WhittakerM[a,b,z],WhittakerW[a,b,z]};

{#,MeijerGReduce[#,z]}&/@ MeijerG137//TableForm

(({#,MeijerGForm2[#,z,op]}&/@MeijerG137//TableForm)/.Piecewise1->Piecewise/.ConditionalExpression1->ConditionalExpression/.Sum1->Inactive[Sum]/.MeijerM->Inactive[MeijerG])//TraditionalForm

6) Formulas for 10 Heun and 4 Lamé functions, for 8 Carlson elliptic integrals and for the FoxH function

In recent years Wolfram Language has added some new generic functions (10 Heun and 4 Lamé functions, 8 Carlson elliptic integrals, the FoxH function and 4 Coulomb functions). For these functions we have corresponding formulas of fractional order integro-differentiation, for example,

{#,[#,{z,α}]}&/@{HeunG[a,q,α1,β,γ,δ,z],LameC[ν,j,z,m],CarlsonRF[x,y,z],CarlsonRJ[,y,a,z],FoxH[{Table[{(Subscript[a, i]),(Subscript[α, i])},{i,1,n}],Table[{(Subscript[a, i]),(Subscript[α, i])},{i,n+1,p}]},{Table[{(Subscript[b, i]),(Subscript[β, i])},{i,1,m}],Table[{(Subscript[b, i]),(Subscript[β, i])},{i,m+1,q}]},z]}//Simplify//TableForm

Applications

Dimensions.Pi— publication frequency data

fractional conservation of mass equation ⧦ compressible flow

fractional advection dispersion equation ⧦ dispersion in porous media

Darcy’s law (groundwater flow)

anomalous dispersion ⧦ flow in crowded systems

protein movements in a cell

diffusive systems

diffusion of phosphorous/boron in silicon

fractional heat equation

viscoelastic damping

acoustic wave equations for complex media

fractional quantum mechanics

integral transforms

Solving Bessel’s equation

A common differential equation appearing in classical physics is Bessel’s equation, which can be written as follows:

(x^2) w''+ x w'+(x-(ν^2)/4)w==0

where ν is required to be an integer due to boundary conditions. Making the substitution u = (x^(±ν/2))w, this becomes

x u''+(1±ν) u'+u==0

Next, we make the ansatz that u can be written in terms of a fractional derivative: u = (D^(1/2±ν)) f for some function f(x). Some manipulation, including use of the Leibniz rule, then gives

(x f)''-3/2f'+f==0

where, remarkably, ν has dropped out entirely. Now, substituting the  independent variable by y = 2(Sqrt[x]), gives us (with primes being understood as x-derivatives everywhere):

f'==(((∂_y) f)/(Sqrt[x]))==(2/y)(∂_y) f

f''==-(4/(y^3))(∂_y) f+(4/(y^2))(∂_y^2) f

(x f)''==(f+x f')'=2f'+x f''=(3/y)(∂_y) f +(∂_y^2) f

Substituting these back above then gives

(3/y)(∂_y) f +(∂_y^2) f-3/2((2/y)(∂_y) f) +f==0

which simplifies down to (∂_y^2) f+f == 0, with the obvious solutions

(Subscript[f, 1])==Sin[y]==Sin[2(Sqrt[x])]

(Subscript[f, 2])==Cos[y]==Cos[2(Sqrt[x])]

Boundary conditions typically require that we take the first solution. Thus, we expect a family of solutions of the form (D^(1/2±ν)) Sin([2(Sqrt[x])) for any integer ν. Performing the requisite derivative, yields the Bessel functions, as expected:

[Sin[2Sqrt[x]],{x, 3/2}]

[Sin[2Sqrt[x]],{x, 5/2}]

Inverse scattering

For a collision between two particles with a spherically symmetric potential V(r), the scattering angle θ as a function of collision energy E is given by:

θ[E]= b (∫_(Subscript[r, c])^∞) ((\[DifferentialD]r)/((r^2)(Sqrt[1-((V[r])/E)-((b^2)/(r^2))])))

The radius of nearest approach, (Subscript[r, c]), and the impact parameter, b, are related through the potential via E=V((Subscript[r, c]))  + E (b^2)/(Subscript[r, c]^2) .  We are concerned with  finding the function V(r), given an experimentally measured θ[E].

With the changes of variables x==1/r and U(x)=V(x) + (b^2)E (x^2) and using above formulae, we can obtain

θ[E]==b(Sqrt[E])(∫_(Subscript[0, c])^E) (((E-U[x]))^(-3/2))((\[DifferentialD]x[U])/(\[DifferentialD]U))\[DifferentialD]U

From the definition of the Riemann-Liouville half-derivative, we recognize the right hand side, and can rewrite this equation as:

Gamma[1/2]b (Sqrt[E])((\[DifferentialD]x[U])/(\[DifferentialD]U))=(\[DifferentialD]/(\[DifferentialD]U))((D^(1/2)) θ[U])

or, integrating with respect to U up to E:

(1/r)==(1/(Gamma[1/2]b (Sqrt[E])))(D^(1/2)) θ[E]==((Sqrt[2μ])/(Gamma[1/2]L))(D^(1/2)) θ[E]

where L and μ are, respectively, the reduced mass and relative angular momentum of the two particles. In Rutherford scattering, the observed deflection angle involves a cotangent.

References

About 30 years ago we published large (about 1000 pages) monograph
Fractional Integrals and Derivatives. Theory and Applications, 1993 by S.G.Samko, A.A.Kilbas & O.I.Marichev
 with about 2000 references, which described different approaches to fractional differentiation and many related calculus directions.

There exists several magazines, devoted to fractional calculus, for example,
Fractional Calculus and Applied Analysis. Editor Virginia Kiryakova

Also we can mentioned several recently published books, including fractional calculus and its applications 
Theory and applications of fractional differential equations, A.A.Kilbas, H.M.Srivastava, J.J.Trujillo, 2006,
Theory and Numerical Approximations of Fractional Integrals and Derivatives, 2019 by Changpin Li & Min Cai