I would like to plot fractional functions. Say, $f(x)=\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 $H$ 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?

  • $\begingroup$ 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) $\endgroup$
    – ssch
    Commented Sep 2, 2013 at 19:21
  • 3
    $\begingroup$ 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. $\endgroup$ Commented Sep 2, 2013 at 22:58
  • $\begingroup$ Open article notebook directly with NotebookOpen["http://www.internationalmathematicasymposium.org/IMS99/paper46/FractionalCalculus.nb"] and Mikes with NotebookOpen["http://library.wolfram.com/infocenter/MathSource/7524/SemiIntegration.nb?file_id=7063"] $\endgroup$
    – ssch
    Commented Sep 2, 2013 at 23:43
  • $\begingroup$ For the equation f[f[x]]==Sin[x] it is interesting to consider solutions with period 2 Pi. These can be conveniently expressed as Fourier series. $\endgroup$
    – mikado
    Commented Jul 29, 2016 at 18:35

2 Answers 2


Use NonlinearModelFit.

Start with some good guess for what the solution operator might look like, apply that to itself and feed that 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];

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



Fractional Iterates

A way to obtain an approximate fractional iterate of a function is to use its Carleman matrix, which is formed from its Taylor coefficients, and then taking the appropriate $p$-th power of the matrix to obtain the series coefficients.

Note that I never said that $p$ had to be an integer; in the example given in the OP, then, we can take the square root of the Carleman matrix to get the coefficients of the semiiterate.

Using the routine CarlemanMatrix[] from this answer, we can do this:

x0 = 0; n = 21; (* expansion point and order *)
sinCM = N[CarlemanMatrix[Sin[x], {x, x0, n}], 30];
shalfCoeffs = MatrixPower[Transpose[sinCM], 1/2, UnitVector[n + 1, 2]];
shalf[x_] = Fold[(#1 x + #2) &, 0, Reverse[shalfCoeffs]];


  1. I used arbitrary precision to preserve some accuracy, as taking the square root of an unsymmetric matrix can be unstable.
  2. I used the "action" form of MatrixPower[], as I only need the coefficients of the semiiterate, which reside in the second row of the matrix square root (hence the Transpose[]).
  3. Even with the refinement in point 2, the computation takes a fair bit of time; more so if you're using a larger Carleman matrix.

To check that we now have the Taylor coefficients of the semiiterate, we can use ComposeSeries[]:

ComposeSeries[shalf[x] + O[x]^23, shalf[x] + O[x]^23] - Sin[x] // Chop

where we have exploited the fact that shalf[x] consists only of odd powers.

One can then also consider constructing a Padé approximant from the series, like so:

shalfpade[x_] = x PadeApproximant[shalf[x]/x, {x, x0, (n - 1)/2}];

(again exploiting the oddness of the function)

Unfortunately, both the Taylor and Padé approximants have a rather limited range of validity:

Plot[{shalf[x], shalfpade[x]}, {x, -3 π/4, 3 π/4}, 
     PlotStyle -> {RGBColor[7/19, 37/73, 22/31], RGBColor[59/67, 11/18, 1/7]}]

approximate semiiterates

There may be a way to construct better approximants from the Taylor coefficients of the semiiterate, but I haven't investigated them yet.

Fractional Derivatives

We can use the general formula of the $\alpha$-th derivative of $x^k$,

$$\frac{\mathrm d^\alpha}{\mathrm dx^\alpha}x^k=\frac{\Gamma(k+1)}{\Gamma(k+1-\alpha)}x^{k-\alpha}$$

and substitute this into the Maclaurin series for the sine. This yields

sinDFrac[x_, a_] = Gamma[a + 1] Sum[(-1)^k Binomial[2 k + 1, a]
                                    x^(2 k + 1 - a)/(2 k + 1)!, {k, 0, ∞}]
   x^(1 - a) Binomial[1, a] Gamma[a + 1]
   HypergeometricPFQ[{1}, {1 - a/2, 3/2 - a/2}, -x^2/4]

that is, a ${}_1 F_2$ hypergeometric function.

In the case $\alpha=\frac12$, we obtain a particularly nice form:

sinDHalf[x_] = FullSimplify[FunctionExpand[sinDFrac[x, 1/2]]]
   Sqrt[2] (Cos[x] FresnelC[Sqrt[2 x/π]] + FresnelS[Sqrt[2 x/π]] Sin[x])

Plot the sine and its semiderivative:

Plot[{Sin[x], sinDHalf[x]}, {x, 0, 2 π}, 
     PlotStyle -> {RGBColor[7/19, 37/73, 22/31], RGBColor[59/67, 11/18, 1/7]}]

sine and its semiderivative


Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.