# How can I numerically solve for fractional functions and fractional derivatives?

I would like to plot fractional functions. Say $f=(\sin)^{1/2}(x)$. By that I mean that $f(f(x)) = \sin(x)$.

Similarly I can define a half-derivative to be an operator such that $H[H[f(x)]] = \frac{df(x)}{dx}$.

What I'd like to do, in an ideal situation, is to evaluate the function at different points and then plot it.

How do I numerically solve such equations in Mathematica? Is it possible to get a Taylor series?

Side question: Does Mathematica treat operators as separate entities?

-
Is there any transform on functions where you can do simple operation then invert back and get f[f[x]]? (Similar to differentiation under fourier transform) –  ssch Sep 2 '13 at 19:21
This article (library.wolfram.com/infocenter/Conferences/6139) outlines an implementation in Mathematica. In this (library.wolfram.com/infocenter/MathSource/7524) I have implemented fractional calculus for lists. –  Mike Honeychurch Sep 2 '13 at 22:58
Open article notebook directly with NotebookOpen["http://www.internationalmathematicasymposium.org/IMS99/paper46/Fr‌​actionalCalculus.nb"] and Mikes with NotebookOpen["http://library.wolfram.com/infocenter/MathSource/7524/SemiIntegra‌​tion.nb?file_id=7063"] –  ssch Sep 2 '13 at 23:43

Start with some good guess for what the solution operator might look like, apply that to itself and feed into NonLinearModelFit:

(* Degree of polynomial approximation *)
k = 5;
(* Points to fit against {{x1, Sin[x1]}, ... } *)
pts = {#, Sin[#]}\[Transpose] &@Range[0, Pi/2, 2/(Pi 2 k)] // N;
(* Our guess for a decent model, c[1]x + c[2]x^2 + ... + c[k] x^k *)
operator[x_] = Plus @@ Array[c@# x^# &, {5}];
(* Apply it to itself *)
doubleOperator = operator[operator[x]]
(* and hopefully find a good fit *)
model = NonlinearModelFit[pts, doubleOperator, Array[c, {5}], x];

GraphicsRow[{
Plot[
Abs[Sin[x] - model[x]],
{x, 0, Pi/2},
PlotLabel -> "Absolute fit error",
PlotRange -> All],
Plot[
{Sin[x], operator[x] /. model["BestFitParameters"]},
{x, 0, Pi/2},
PlotLabel -> Defer[{Sin, Sqrt[Sin]}]]
}]


-