What's NiNJI251RzYi about the Derivative? A Reprise Robert M. Corless, ORCCA & Applied Math, UWO Acknowledgements: Mhenni M. Benghorbal, David Jeffrey, Elena S. Smirnova, Ryan Morris, Matt Davison, Chris Essex, Bruno Salvy, Richard Fateman, Edgardo Cheb-Terrab, Jacques Carette, Stephen M. Watt system( "yap C:/Data/talks/mocaa/whatsnu.dvi" ); The ability to do fractional and symbolic differentiation was installed in 2004 on the MONET symbolic services page, which can be seen at http://ptibonum.scl.csd.uwo.ca:16661/MonetServiceClient/ . This includes code for derivatives that was written by Mhenni M. Benghorbal, polished by Ryan Morris, and wrapped by Elena S. Smirnova. This was the subject of the MOCAA talk in 2004. This code has now been incorporated (and improved somewhat) into Maple 10, courtesy Edgardo Cheb-Terrab. We will look at both the original and the incorporated versions. Symbolic Order Differentiation (some examples of G.H. Hardy)

An easy proof by hand system("C:/Program Files/Adobe/Acrobat 5.0/Reader/AcroRd32.exe C:/Data/talks/mocaa/iteratedintegrals.pdf" ); In Maple, the formula is (Aye@@n)(f)(x) = Int( (x-t)^(n-1)/(n-1)!*f(t),t=a..x ); Let's try some examples of that formula. Eh := unapply( rhs(%), n, f ); Eh(1,sin); Eh(2,sin); int( int( int( sin(teh), teh=a..tau), tau=a..t ), t=a..x ); value( Eh(3,sin) ); Eh( 0, sin ); We'll come back to that, later. But, the formula seems to work.

Everyone can interpolate polynomials. restart; CurveFitting:-PolynomialInterpolation( [a,b,c], [A,B,C], x, form=Lagrange ); seq( eval(%,x=i), i=[a,b,c] ); The factorial function is interpolated uniquely by a log-convex analytic function, the GAMMA function. (Bohr-Mollerup, 1922). Apparently, Euler invented the GAMMA function in order to interpolate a formula for derivatives... (4/3)! = evalf( (4/3)! ); convert(%,GAMMA); plotfact := plots[logplot]( [seq([n,n!],n=0..10)], style=POINT, symbolsize=30, colour=BLUE ): plotgam := plots[logplot]( GAMMA(x+1), x=0..10 ): plots[display]( {plotfact,plotgam} ); Charles M. Patton's observation Now, we can interpolate some discrete dynamical systems in an interesting way---this is called computing the log of the homomorphism, or "the method of modified equations" in numerical analysis. Let's interpolate the numerical solution of y' = lambda*y by Euler's method: u[n+1] = u[n] + h*u[n]' or u[n+1] = u[n] + h*lambda*u[n] rsolve( u(n+1) = (1+h*lambda)*u(n), u(n) ); exp( nh*log(1+h*lambda)/h ) = (1+h*lambda)^n; Here is a differential equation whose solution interpolates the Euler solution to y' = lambda*y. QyQ+SSZtb2RkZUc2Ii8tSSVkaWZmRyUqcHJvdGVjdGVkRzYkLUkiWUdGJTYjSSJ0R0YlRi4qKC1JI2xuRzYkRilJKF9zeXNsaWJHRiU2IywmIiIiRjYqJkkiaEdGJUY2SSdsYW1iZGFHRiVGNkY2RjZGOCEiIkYrRjZGNg== QyQtSSdzZXJpZXNHJSpwcm90ZWN0ZWRHNiUqJi1JI2xuRzYkRiVJKF9zeXNsaWJHNiI2IywmIiIiRi8qJkkiaEdGLEYvSSdsYW1iZGFHRixGL0YvRi9GMSEiIkYxIiIlRi8= dsolve( {modde, Y(0)=u(0)}, Y(t) ); lambda := -10; h := 0.07968121301; u[0] := 1; for i to 10 do u[i] := u[i-1] + h*lambda*u[i-1]; od: numsol := [seq([i*h,u[i]],i=0..10)]: exactsol := exp(-10*x); modifiedequationsol := exp(-20*x); solplot := plots[logplot]({exactsol,modifiedequationsol}, x=0..10*h): datplot := plots[logplot]( numsol, style=POINT, symbolsize=30 ): plots[display]({solplot,datplot} ); That is, the interpolant explains the numerical solution (in a very nice way). lambda := 'lambda'; h := 'h'; series( ln(1+h*lambda)/h, h ); That is O(h) different from lambda, and so Euler's method is 1st order (which we knew). lambda*(1-z/2+z^2/3-z^3/4+z^4/5); series( log(1+z)/z, z ); solve( log(1+x)/x=2, x ); evalf(%); Therefore, if h*lambda = -0.7968, the effective decay rate is 2*lambda, not lambda. This is an interesting way to look at the error behaviour of Euler's method. Now a harder problem: y'(t) = y(t)^2 (Error Backward, and the Beyn function) restart; diff( Y(t), t ) = B(h*Y(t))*Y(t)^2; Now, B(v) is chosen so that v' = B(v)v^2 interpolates v(tau+1) = v(tau) + v(tau)^2, cf Euler's method for the above problem u[n+1] = u[n] + h*u[n]^2; scale t and u by h to get t = h*tau, u = h*v. One can show that B(v) has a (divergent) series and a convergent infinite product... read "Beyn.mpl"; system( "C:/Program Files/WinEdt Team/WinEdt/WinEdt.exe C:/Data/talks/mocaa/Beyn.mpl" ); B(v) = add( c(n)*v^(n-1), n=1..10 ) + O(v^10); B(v) = Product( (1+v[k])^2/(1+2*v[k]), k=0..infinity ); B(0.2); Beyn := proc( n, t, v, vp ) vp[1] := B(v[1])*v[1]^2 end proc; sol := dsolve( procedure=Beyn, v(t), numeric, range=0..11, start=0, initial=array(1..1,[0.3]), procvars=[v(t)] ); u[0] := 0.3; for i to 11 do u[i] := u[i-1] + u[i-1]^2; od: pts := plot( [seq( [k,u[k]], k=0..7 )], style=POINT, symbolsize=30, colour=BLUE ): interplot := plots[odeplot]( sol, view=[0..7,0..20] ): plots[display]({pts,interplot}); exactsol := dsolve( {diff(y(t),t) = y(t)^2,y(0)=3/10}, y(t) ); exact := plot( rhs(exactsol), t=0..7, colour=BLACK, linestyle=DASH ): plots[display]({pts,interplot,exact}); plot( B(v), v=-1..1, y=-5..5 ); solve( v+v^2=u, v ); N := 12; poles := array(0..2^(N)-2); poles[0] := -0.5; last := 0; first := 1 ; for i to 2^(N-2) do u := poles[last]; poles[first] := (sqrt(1.0+4.0*u)-1.0)/2.0; first := first + 1; poles[first] := (-sqrt(1.0+4.0*u)-1.0)/2.0; first := first + 1; last := last + 1; end do: pts := map(t->[Re(t),Im(t)],[seq( poles[i], i=0..2^(N)-2 )]): plot( pts, style=POINT, axes=BOXED, colour=BLUE, symbol=POINT );

We define, analogous to our integer formula derived by interchange of order of integration, a formula valid for "fractional integrals". restart; Formula := Aye@@alpha = Int( (x-t)^(alpha-1)/GAMMA(alpha)*f(t), t=0..x ); limit( rhs(Formula), alpha=1, left ); That's good (what we expect).

limit( rhs(Formula), alpha=0, right ); That's not what we expect, and actually incorrect, because while the GAMMA function is going to infinity, so also is the singularity of the integral near t=x. "limit" is assuming continuity of the integrand, I think. Let us look a bit more carefully: kernel := (x-t)^(alpha-1)*f(t)/GAMMA(alpha); Int( kernel, t=0..x ) = Int( kernel, t=0..x-K*alpha^(1/2) ) + Int( kernel, t=x-K*alpha^(1/2)..x ); Suppose f is continuous on an interval containing [0,x]; in fact, Lipschitz continuous so that |f(x)-f(t)| \le L|x-t|. Then there exists a bound F such that |f(t)| \le F on the interval. Since (x-t)^(alpha-1) \le (K*alpha^(1/2))^(alpha-1) on the interval [0,x-K*alpha^(1/2)], then the first integral on the rhs above is bounded by F*Int( (K*alpha^(1/2))^(alpha-1)/GAMMA(alpha), t=0..x-K*alpha ); value(%) assuming alpha > 0; limit( %, alpha=0, right ); plot( sqrt(alpha)^(alpha-1)/GAMMA(alpha), alpha=0..0.2 ); Therefore in the limit as alpha -> 0, only the second integral contributes anything: Int( kernel, t=x-K*alpha^(1/2)..x ); Now, here, suppose f(t) is replaced by f(x) + e*(x-t); we know |e| \le L. Int( (x-t)^(alpha-1)*f(x)/GAMMA(alpha), t=x-K*alpha^(1/2)..x ) ; value(%) assuming alpha > 0; limit(%, alpha=0, right ); One remaining error term: L*Int( (x-t)^(alpha-1)*(x-t)/GAMMA(alpha), t=x-K*alpha^(1/2)..x ); value( simplify(%) ) assuming alpha > 0; limit(%,alpha=0,right); This completes the proof: as alpha -> 0, our integral tends to f(x).

Well, we have a formula that works for all positive integers, and if we look carefully it does the right thing as alpha -> 0+. Let's declare it to be our "fractional" integral, and see what we can do. restart; Aye := (f,alpha) -> unapply( int( (x-t)^(alpha-1)/GAMMA(alpha)*f(t), t=0..x ), x ); Now let us derive a very useful formula: the fractional integral of a monomial. Aye( t->t^k, alpha ) assuming k::posint; simplify( %(x) ) assuming x > 0; In factorial notation, we have (Aye@@alpha)(x^k) = k!/(k+alpha)!*x^(k+alpha); simplify( limit( %, alpha=0, right ) ); Note the Maple bug on the left hand side... simplify( limit( %%, alpha=1, left ) ); Now let's try integrating some non-monomials. The following comes out either as AngerJ functions or Fresnel integrals, depending on Maple's mood. If we're lucky, it will be AngerJ, because then the 2nd whack works: Aye(sin,1/2); twice := Aye(%,1/2); simplify( twice(x) ); twice(0.3434); const := % +cos(0.3434); twice(0.889393 ); const - cos(0.889393 ); Well, either way it's clear that half-a-step, half-a-step gets us to 1-cos(x) = int( sin(t), t=0..x ). What about other functions? Aye( x->1/sqrt(x), 1/2 ); Aye( %, 1/3 ); Aye( %, 1/6 ) ; simplify( %(x) ) assuming x > 0; int( 1/sqrt(x), x ); 1/2 + 1/3 + 1/6; We have just demonstrated examples indicating that these iterated integrals compose with each other, at least up to constants. So, now we have some basis to believe that this interpolation for integrals will give us something useful (and indeed it does). Watch out for Maple nonrepeatabilities in the above, though---sometimes the integrals come out looking differently, and I have to admit I am nonplussed at that.

The main trick to fractional differentiation is to write the order of differentiation as an integer minus a fraction: restart; n = ceil(nu); alpha = n - nu, nu = n - alpha; n = ceil( 3.85 ); alpha = ceil(3.85) - 3.85; 3.85 = `4 - 0.15`; Since n >= nu, we have alpha >= 0. Since alpha is a fractional part, we have alpha < 1. Remember this: We will take the "nu"th derivative, by taking "n" ordinary derivatives and integrating alpha times! (D@@nu)(f)(x) = (Aye@@alpha)( (D@@n)(f) ) (x ); Yes, the order matters: put the (D@@n) on the inside, as we have done, and we get the Davison-Essex definition; put the alpha on the inside, and we get the Riemann-Liouville definition. We will see that these are (slightly) different. They are useful for different purposes. read "C:/Data/talks/mocaa/fracdiff.mpl": read "C:/Data/talks/mocaa/SymbolicOrder.mpl"; macro( sdiff=SymbolicOrder:-Differentiate ):

(D@@nu)(t->t^m)(x) = Int( (x-t)^(alpha-1)/GAMMA(alpha)*falling(m,n)*t^(m-n), t=0..x ); By our fractional integral formula for monomials, this can be simplified to (D@@nu)(x^m) = falling(m,n)*(m-n)!/(m-nu)!*x^(m-nu); # Davison-Essex (D@@nu)(x^m) = m!/(m-nu)!*x^(m-nu); # Riemann-Liouville Those differ ONLY in that the Davison-Essex derivative is zero if 0 <= m < n is an integer (because the ordinary derivative, n times, of x^m makes it zero. (One trick---we have extended the meaning of the derivative by analytic continuation) plots[display]( [seq( plot({x^3,3*x^2,3!/(3-k/51)!*x^(3-k/51)},x=0..1,colour=black),k=1..50 )], insequence=true );

Now let us look at f := t -> sin(t); n = ceil(nu); nu = n - alpha; The nu-th = 1-alpha -th derivative of sin(t) is, by the Davison-Essex definition, obtainable: tadah := int( (x-t)^(alpha-1)*D(f)(t), t=0..x )/GAMMA(alpha); plots[display]( [seq( plot( {sin(x),cos(x),eval(tadah,alpha=1-k/51)}, x=0..4*Pi, colour=black), k=1..50 )],insequence=true ); fracdiff( sin(x), x, 1/3 ); fracdiff( sin(x), x, 1/3, method=direct ); series(%,x); fracdiff( sin(x), x, 1/3, method=laplace ); series(%,x); Now how about a hard one? fracdiff( tan(x), x, 1/3 ); fracdiff( tan(x), x, 1/3, method=direct ); fracdiff( tan(x), x, 1/3, method=laplace ); The following examples are particularly important in practice (Lavoie, Osler, Tremblay 1976) QyRJKHJlc3RhcnRHJSpwcm90ZWN0ZWRHIiIi fracdiff( z*log(z), z, 1/3 ); fracdiff( z*log(z), z, 1/3, method=direct ); fracdiff( z*log(z), z, 1/3, method=laplace ); fracdiff( z*log(z), z, -1/3, method=direct ); fracdiff( %, z, 1/3, method=direct ); fracdiff( z*log(z), z, 7/11, method=direct ); fracdiff( z^(2/3)*log(z), z, 7/13, method=direct ); fracdiff( z^(2/3)*log(z)^2, z, 1/3, method=direct ); fracdiff( %, z, -1/3, method=direct ); simplify(%); evalf(%);

restart; #read "C:/Data/talks/mocaa/fracdiff.mpl": #read "C:/Data/talks/mocaa/SymbolicOrder.mpl": #macro(sdiff=SymbolicOrder:-Differentiate); #sdiff( sin(x), x, r ); QyQtSSVkaWZmRyUqcHJvdGVjdGVkRzYkLUkkc2luRzYiNiNJInhHRiktSSIkR0YlNiRGK0kickdGKSIiIg== rth := combine(%,trig); plots[display]( [seq( plot({sin(x),cos(x),eval(rth,r=k/51)},x=0..2*Pi,colour=black), k=1..50)], insequence=true ); Notice that this function is always periodic in x, no matter what r is...we are tempted to say that the (real valued) rth derivative of sin(x) is this, but (beep) it isn't, except when r is an integer. Remember the graph in the previous section showing the alpha-th derivative of sin(x), and how the animation showed it was NOT periodic! In particular, that first local maxima decreased, but not for this formula. Therefore this is NOT the (same) rth derivative of sin(x) for real r. Here is an even more outrageous example of the failure of our intuition:

exp(x) = Sum( x^m/m!, m=0..infinity ); (D@@nu)(exp)(x) = Sum( m!/(m-nu)!*x^(m-nu)/m!, m=ceil(nu)..infinity ); value( rhs(%) ); convert( normal(% - exp(x)), GAMMA ); simplify(expand(%)); (D@@nu)(exp)(x) = exp(x) + %; Dnu := unapply( exp(x)*(1-GAMMA(ceil(nu)-nu,x)/GAMMA(ceil(nu)-nu)), nu, x ); limit( Dnu(nu,x), nu=3, left ); evalf(%); NumericEventHandler(division_by_zero=(()->args[3])); Dnu( -4, x ); expand(simplify(%)); Dnu( 4, x ); hah := Dnu( 4.9, x ); eval( hah, x=0.3 ); plot({hah,exp(x)}, x=-1..1,y=0..3 ); plot( GAMMA(alpha,0.3)/GAMMA(alpha), alpha=0..1, y=0..1, numpoints=101 ); limit( GAMMA(alpha,3/10)/GAMMA(alpha), alpha=1, left ); GAMMA(0,x); GAMMA(7,x); That derivative contains extra terms that turn into polynomials as alpha goes to 1, 2, 3, ..., which is weird until you realize that (D@@(-1))(exp)(x) = int( exp(t), t=0..x ) = exp(x) - 1, and similarly iteration of the integration operator knocks off 1+x, 1+x+x^2/2, and so on. But in between integer values, the derivative of exp(x) is NOT exp(x), because the incomplete Gamma function is not zero! But wait! If we integrated from -infinity instead of 0 in our fractional derivative definition, then all the incomplete Gamma functions go away...this extra term is an artifact of the lower limit. Putting the lower limit to -infinity gets us the Weyl definition. Here is the formula for fractional differentiation of exp(x) based at an arbitrary finite point instead of at zero: (exercise): (D[a]@@nu)(exp)(x) = exp(x)*(1-GAMMA(alpha,x-a)/GAMMA(alpha)); plot( GAMMA(0.4,0.22-a), a=-10..-1 ); For fixed x and nu (alpha), limit as a -> -infinity is zero.

sdiff does not work for log(1+x) because gfun:-holexprtodiffeq doesn't recognize it; so we differentiate it once, and then work on it. restart; read "C:/Data/talks/mocaa/SymbolicOrder.mpl"; macro( sdiff = SymbolicOrder:-Differentiate ); read "C:/Data/talks/mocaa/fracdiff.mpl"; sdiff( diff(log(1+x),x), x, m-1 ); mlogder := %; seq( mlogder, m=1..5 ); It turns out to be better to write things using factorials. _EnvFallingNotation := factorial; sdiff( 1/(1+t), t, n-1 )/n!; simplify( eval(%,t=0) ); series( ln(1+x), x ); Well, we knew what ln(1+x) is as a sum, anyway. lnp1(x) = Sum( (-1)^(m-1)*x^m/m, m=1..infinity ); # we know 1st coeff is 0 Incidentally, here's a function not in the convert/Sum database: convert( log(1+x), Sum ); Now, remember our fractional derivative of a monomial: (D@@nu)(x^m) = m!*x^(m-nu)/(m-nu)!; # zero if ceil(nu) in {0,1,...m-1} forget(ceil); x := 'x'; lowlim := ceil(nu); # ceil(nu) or 1 missing := Sum( (-1)^(m-1)/m*m!/(m-nu)!*x^(m-nu), m=1..lowlim-1 ); Skip over this step to look at the Riemann-Liouville derivative; execute it to see the Davison-Essex derivative (because then falling( m, n ) = 0 for m < n, not the factorial formula above ). missing := 0; x := 'x'; nuthln := Sum( (-1)^(m-1)/m*m!/(m-nu)!*x^(m-nu), m=lowlim..infinity ); It is quite pleasant that the following gives an answer. value(%); convert(%,StandardFunctions); trial := unapply(%+missing,nu); trial(1); simplify(value(%)); trial(2); trial(3); value(simplify(%)); diff( ln(1+x), x$3 ); Wow! So the formula does indeed interpolate our integer-order derivatives...but it's supposed to do the fractional derivatives too. Let's try, at a random value of x. x := 0.3; interplot := plot( trial, 0..5): dataplot := plot( [seq([k,(D@@k)(t->log(1+t))(x)],k=0..5)], style=POINT, symbolsize=30, colour=BLUE ): plots[display]({interplot,dataplot},view=[0..5,-120..20]); If we didn't include the missing terms, then we would get discontinuities...the Davison-Essex definition is discontinuous (from the right) at integer values of the order (bigger than 1); but this turns out to be precisely what we need in order to solve initial-value problems. For integral equations, we need continuity in the order, and therefore we must use the Riemann-Liouville definition.

QyRJKHJlc3RhcnRHJSpwcm90ZWN0ZWRHIiIi fracdiff( hypergeom( [1,2], [3], x ), x, 1/3 ); fracdiff( hypergeom( [1,2], [3], x ), x, 1/3, method=direct ); sym := diff( hypergeom([1,2],[3],x), x$nu ); Is that true for noninteger nu? I don't think so! series( sym, x, 3 ); Sum( pochhammer(1,n)*pochhammer(2,n)/pochhammer(3,n)*x^n/n!, n=0..infinity ); diff( %, x$nu ); subsop(2=(n=ceil(nu)..infinity),%); synth := value(%); eval( synth-sym, [nu=rand()/1.0e11,x=rand()/1.0e11] ); x^(ceil(nu)-nu)*series( simplify(synth*x^(nu-ceil(nu))), x, 3 ); eval(%,nu=1/3); Th-th-th-that's all, folx.

References for fractional derivatives are a bit tricky: everyone thinks they "own" the area. But, references go at least back to Liouville, and even perhaps to Euler---whom Osler says invented the Gamma function in order to do fractional derivatives. But, here are a representative few. system("C:/Program Files/Adobe/Acrobat 5.0/Reader/AcroRd32.exe C:/Data/talks/mocaa/iteratedintegrals.pdf" );