Take the 2-minute tour ×
Mathematica Stack Exchange is a question and answer site for users of Mathematica. It's 100% free, no registration required.

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?

share|improve this question
    
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
2  
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

2 Answers 2

up vote 3 down vote accepted

Using NonlinearModelFit

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]}]]
  }]

fit

share|improve this answer

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]];

Notes:

  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
   O[x]^23

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

share|improve this answer

Your Answer

 
discard

By posting your answer, you agree to the privacy policy and terms of service.

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