# Implementing Riemann-Liouville fractional derivative

I'm working with fractional derivatives, primarily Riemann-Liouville and Caputo derivatives. In these definitions n=Ceiling[\alpha], and a is chosen to be 0. They are defined as:  I implemented the Caputo fractional derivative on Sech[x] like so

derivC[q_, x_] := Module[{s=Ceiling[q]}, If[FractionalPart[q]==0, D[Sech[x], {x, s}], 1/Gamma[s-q] NIntegrate[Derivative[s][Sech][x-u] u^(s-q-1), {u,0,x}]]]


Inside the integrand we change variables x->t in accordance with the definition then x-t->u for ease of computation. Inside the integrand we take an integer number of derivatives on Sech and then integrate with a power law. If the order of the derivative q is integer then it should reduce to an integer-order derivative. Then using NIntegrate I integrate up to a bound x that's a real number because it is given by the function wrapper.

Whew! Long backstory.

The problem with implementing the Riemann-Liouville derivative though is that you take the integer number of derivatives outside of the integral. Since I'm using NIntegrate I only have numbers, not a function, which would make the derivative 0.

I do not know of a way to implement this. Any help or hints or suggestions are appreciated.

To obtain derivatives of a numerical integral, one can try inactivating the integral, then differentiate, then activate. For example:

Inactive[NIntegrate][(x-ξ)^n f[ξ], {ξ, 0, x}];
D[%, x]


0^n f[x] + Inactive[NIntegrate][n (x - ξ)^(-1 + n) f[ξ], {ξ, 0, x}]

The problem with this approach is that in your application, $n<0$, and so differentiating this form of the integral is not useful. However, I again suggest that you transform the integration variable $\xi \to x - u$. Then:

Inactive[NIntegrate][u^n f[x-u], {u, 0, x}];
D[%, x]


x^n f + Inactive[NIntegrate][u^n Derivative1[f][-u + x], {u, 0, x}]

This version doesn't suffer from fatal divide by 0 issues. We can package up the above ideas into a function:

RLDerivative[α_][f_] := With[{s = 1/Gamma[Ceiling[α]-α]},
RLFunction[
α,
f,
Function[{x,ϵ},
Evaluate @ D[
s Inactive[NIntegrate][u^(Ceiling[α]-α-1)f[x-u], {u, 0, ϵ I, x}],
{x, Ceiling[α]}
]
]
]
]

RLFunction[α_, f_, func_][x_?NumericQ] := Replace[
Quiet @ Check[Activate[func[x,0]], $Failed],$Failed :> (Message[RLFunction::ipath]; Activate[func[x, .1]])
]

RLFunction::ipath = "Unable to achieve convergence along the real axis. Integral path was deformed";


First note that I changed the path of integration to include $i \epsilon$. The default is to use $\epsilon =0$. However, there can be situations where the numerical integral doesn't converge. In that case, deforming the path of the integral can avoid this convergence issue. One down side is that using a complex path can introduce a spurious machine epsilon imaginary part. The other issue is the function under consideration may have a pole near $i \epsilon$ which will cause the deformed path integration to return an incorrect result.

Next, I split up the computation into two pieces. RLDerivative[α][f] will return an RLFunction object that can be applied to x values. This way the derivative of the inactive NIntegrate object only needs to be done once.

With those caveats mentioned, let's see how it works. First a test case that can be performed exactly:

int[n_, x_] = Integrate[(x-ξ)^n Exp[ξ], {ξ, 0, x}];
rlexp[α_, x_] := 1/Gamma[Ceiling[α]-α] D[int[Ceiling[α]-α-1, \[FormalX]], {\[FormalX], Ceiling[α]}]/. \[FormalX]->x


Now, compare RLDerivative with the exact rlexp:

rf = RLDerivative[.9][Exp];
rf /@ Range
rlexp[.9, Range//N]

rf = RLDerivative[12.2][Exp];
rf /@ Range
rlexp[12.2, Range//N]


{2.75781, 7.40346, 20.0932, 54.603, 148.416}

{2.75781, 7.40346, 20.0932, 54.603, 148.416}

{1.12015*10^7, 2211.05, 34.6723, 55.0066, 148.438}

{1.12015*10^7, 2211.05, 34.6723, 55.0066, 148.438}

Now, for an example where the default integral path has convergence issues:

rf = RLDerivative[.9][Sin];
rf[.9]


RLFunction::ipath: Unable to achieve convergence along the real axis. Integral path was deformed

0.698049 + 2.99218*10^-14 I

Notice the warning message and the appearance of the small imaginary part resulting from the deformed path. Finally, the RLFunction looks like:

RLDerivative[1.9][Sin]

(*
RLFunction[1.9, Sin,
Function[{x$, ϵ$},
0.105114 (1/x$^0.9 + Inactive[NIntegrate][Sin[u - x$]/u^0.9, {u, 0, I ϵ$, x$}])]]
*)


This looks a bit ugly, so here is a MakeBoxes rules for RLFunction, intended to mimic the standard M11 appearance of summary boxes for these types of functions:

RLFunction /: MakeBoxes[RLFunction[a_, f_, _], StandardForm] := RowBox[{
"FractionalDerivative","[",
PanelBox[
RowBox[{
ToBoxes @ Plot[
f[x],
{x,0,1},
Axes->False,Frame->True,FrameTicks->None,
AspectRatio->1,ImageSize->20
],
"\" Order:\"",
}],
Background->GrayLevel[.95]
],
"]"
}]


Now, the RLFunction format hides the gory details:

RLDerivative[1.9][Sin] • A couple comments: 1. I don't think negative values for α work properly. 2. Functions with an integrable singularity at the origin won't work with this method. I could fix both issues with a bit of time. Jun 8 '17 at 2:57
• What a nice answer! Jun 8 '17 at 6:18