[The following webpage is an HTML export of fractional_derivatives.mw]
The function "fracD" computes the fractional derivative of a function using a combination of
the composition law and fractional integration. To compute the fractional integral-derivative
we first (fractional) integrate up to an integer level above our target (using the Riemann-Liouville
definition of fractional integration):
\(\displaystyle {}_aD_t^{-q}[f(t)] = \frac{1}{\Gamma(q)}\int_a^t (t-s)^{q-1}f(s)\,ds \)
which is valid for \( q >0 \). Then we use regular derivatives to move down to the desired result.
For example: \(\displaystyle {}_aD_t^{0.5}[f(t)] = \frac{d}{dt}\left[ {}_aD_t^{-0.5}[f(t)]\right] = \frac{d}{dt}\left[ \frac{1}{\Gamma(1/2)}\int_a^t (t-s)^{-(-0.5)-1}f(s)\,ds\right]\)
Or again, briefly, \(D^{8.4}=D^9\;D^{-0.6} \) where \( D^{-0.6} \) is computed using an integral and \( D^9 \) is the regular derivative taken 9 times.
| > | restart;
with(plots): fracD := proc(f,t,mu,b) # This computes the fractional derivative of f(t) # of order mu based at t=b. local m; # fractional integrals can be computed using an integral... if mu<0 then RETURN(1/GAMMA(-mu)*int((t-u)^(-mu-1)*subs(t=u,f),u=b..t)); end if; # finding "m" such that "m" is a positive integer and "mu-m < 0". if type(mu,integer) then m := mu+1; else m := ceil(mu); end if; # integrate up the fractional part "mu-m" then differentiate # down using regular derivatives "m" times. RETURN(diff(fracD(f,t,mu-m,b),t$m)); end proc: fracD3D := proc(f,x,x0,a,b,c,d,N,hl,twoD) # Assuming f is an expression depending on x. # Plot x0_D_x^t[f] for c < t < d and a < x < b. # N is the number of steps. # Of course, a < b and c < d also N is a positive integer. # hl is a list of highlighted derivatives. local i,j,pts,ft,t,xV,plotPts,plotCrvs,plothl,plothl2d,anim2D; plotPts := Array(1..N+1); plotCrvs := Array(1..N); pts := Array(1..N,1..N); plothl := Array(1..nops(hl)); plothl2d := Array(1..nops(hl)); anim2D := Array(1..N); t := c; for j from 1 to N do t := t + (d-c)/N; ft := fracD(f,x,t,x0); plotCrvs[j] := spacecurve([x,t,ft],x=a+(b-a)/N..b,color=blue); anim2D[j] := plot(ft,x=a+(b-a)/N..b,color=red,thickness=2); xV := a; for i from 1 to N do xV := xV + (b-a)/N; pts(i,j) := evalf(subs(x=xV,ft)); end do: end do: for i from 1 to N do plotPts[i] := pointplot3d([seq([a+(b-a)/N*i,c+(d-c)/N*j,pts[i][j]],j=1..N)],connect=true,color=black); end do: # plotPts := pointplot3d([seq(seq([a+(b-a)/N*i,c+(d-c)/N*j,pts[i][j]],i=1..N),j=1..N)],connect=true,color=black); for i from 1 to nops(hl) do ft := fracD(f,x,hl[i],x0); plothl[i] := spacecurve([x,hl[i],ft],x=a+(b-a)/N..b,color=red,thickness=3); plothl2d[i] := plot(ft,x=a+(b-a)/N..b,color=blue,thickness=2); end do: for i from 1 to N do anim2D[i] := display(anim2D[i],seq(plothl2d[j],j=1..nops(hl))); end do: if twoD then display(seq(anim2D[i],i=1..N),insequence=true); else display3d(seq(plotCrvs[i],i=1..N),seq(plotPts[i],i=1..N),seq(plothl[i],i=1..nops(hl))); end if; end proc: |
To compute the fractional integral/derivative \( {}_aD_t^q[f(t)] \) (the \( q^{th} \) integral/derivative of \( f(t) \) with respect to \( t \) based at \( t=a \) ), we use "fracD(f(t),t,q,a);"
| > | fracD(exp(3*t),t,3/2,-infinity); |
| \( \color{blue}{3\sqrt{3} e^{3t}} \) | (1) |
| > | fracD(sin(t),t,1/2,-infinity); |
| \( \displaystyle \color{blue}{\frac{1}{2}\left(\sin(t)+\cos(t)\right)\sqrt{2}} \) | (2) |
Note that \( \displaystyle \frac{\sqrt{2}}{2}\sin(t)+\frac{\sqrt{2}}{2}\cos(t)=\sin(t)\cos\left(\frac{\pi}{4}\right)+\cos(t)\sin\left(\frac{\pi}{4}\right)=\sin\left(t+\frac{\pi}{4}\right) \)
In general, \( \displaystyle {}_{-\infty}D_t^k[\sin(t)]=\sin\left(t+\frac{\pi}{2}k\right) \).
| > | fracD(t^2,t,1/2,0); |
| \( \displaystyle \color{blue}{\frac{8}{3} \frac{t^{3/2}}{\sqrt{\pi}} }\) | (3) |
These are the derivatives of orders 0 through 4 of \( f(x)=x^3 \) based at \( x=-2 \). The animation and 3D graphs show how the fractional derivatives interpolate between the regular derivatives.
| > | graph1 := fracD3D(x^3,x,-2,0,4,-1/5,6,60,[0,1,2,3],true):
display(graph1,insequence=true,view=[0..4,-2..15]); |
|
| > | graph2 := fracD3D(x^3,x,-2,0,4,-1/5,6,29,[0,1,2,3],false):
display3d(graph2,orientation=[-90,90,0],viewpoint=circleleft); |
 |
These are the derivatives of orders 0 through 4 of \( f(x)=\sin(x) \). Here derivatives amount to nothing more than translations.
| > | graph1 := fracD3D(sin(x),x,0,Pi,3*Pi,0,4,30,[0,1,2,3],true):
display(graph1,insequence=true,view=[Pi..3*Pi,-1.2..1.2]); |
 |
| > | graph2 := fracD3D(sin(x),x,0,Pi,3*Pi,-1/5,4,30,[0,1,2,3],false):
display3d(graph2,orientation=[-90,90,0],viewpoint=circleleft,scaling=constrained); |
 |