Examples XLV from G. H. Hardy, ``A Course on Pure Mathematics'' 10th ed, p. 229 (orig. ed. 1908) restart; read "C:/Data/talks/mocaa/SymbolicOrder.mpl"; macro( sdiff=SymbolicOrder:-Differentiate ): In this worksheet we look at symbolic order derivatives---that is, differentiate f(x) a number n times, but without specifying n in advance. This is a quite typical first year calculus problem, that the students are expected to solve by guessing the formula and then proving it by induction.

phi[m] := x^m; The following command used to return unevaluated in Maple: now it works! diff( phi[m], x$n ); All that used to be available was repeated (iterated) differentiation. diff( f(x), x$20 ); lprint(%); diff( phi[m], x ); diff( phi[m], x$2 ); diff( phi[m], x$3 ); simplify(%); factor(%); phi := 'phi'; From the original talk: Assume diff( phi[m], x$k ) = x^(m-k) * falling( m, k ). Then diff( phi[m], x$(k+1) ) = diff( x^(m-k)*falling(m,k), x ) = falling(m,k) * (m-k)* x^(m-k-1) = falling( m, k+1) * x^(m-(k+1)). QED. How does Maple help with that, in any way? Maybe gfun can help? with(gfun): From an anonymous ISSAC referee (probably Bruno, because it seems to work :-) f := x^m; temp := n!*rsolve( diffeqtorec( holexprtodiffeq( subs( x=a+x, f ), y(a) ), y(a), u(n)), u(n)); What does "holonomic" mean? (the routine holexprtodiffeq means "holonomic expression to differential equation") What does Google say? Of course, the word occurs 52 times in the Maple help pages (on 9 pages in help for various routines for gfun, one page for SumTools:-Hypergeometric, and one page in Groebner (hilbertdim)), but nowhere is it defined. My copy of Marsden and Ratiu (Introduction to Mechanics and Symmetry) says "In geometry, when an orthonormal frame returns after traversing a closed path to its original position but rotated, the rotation is referred to as holonomy". Or, "A holonomic constraint can be defined for our purposes as the specification of a submanifold N in Q of a given configuration manifold Q. (More generally a holonomic constraint is an integrable subbundle of Q)." Eh? Of course the word is just being used in a different sense. A "holonomic function" (or D-finite function) is a function that satisfies a linear differential equation with polynomial coefficients (presumably over some closure of the rationals). It would be nice if that was somewhere in the help pages (Maple glossary, anyone?) LUklaGVscEc2IjYjUSpob2xvbm9taWNGJA== maybe := simplify( temp ) assuming x > 0; seq( factor( simplify( limit( maybe, n=k ) ) ), k=0..3 ) ; HOORAY! Maple can now compute the nth derivative of x^m! QyQtSSZzZGlmZkc2IjYlKUkieEdGJUkibUdGJUYoSSJuR0YlIiIi diff( x^m, x$n ); Is that correct? seq( simplify( % ), n=0..4 ); _EnvFallingNotation := factorial; QyQtSSZzZGlmZkc2IjYlKUkieEdGJUkibUdGJUYoSSJuR0YlIiIi diff( x^m, x$n ); _EnvFallingNotation := GAMMA; QyQtSSZzZGlmZkc2IjYlKUkieEdGJUkibUdGJUYoSSJuR0YlIiIi diff( x^m, x$n ) ; QyQtSSlhc3N1bWluZ0dJKF9zeXNsaWJHNiI2JDcjLUkmc2RpZmZHRiY2JSlJInhHRiZJIm1HRiZGLUkibkdGJjcjMkYuIiIhIiIi diff( x^m, x$n ) assuming m < 0; _EnvFallingNotation := pochhammer;

QyQtSSZzZGlmZkc2IjYlKSwmKiZJImFHRiUiIiJJInhHRiVGK0YrSSJiR0YlRitJIm1HRiVGLEkibkdGJUYr diff( (a*x+b)^m, x$n ); Whoops, it appears that the "improvements" are not all in the positive direction. The original formula is right. seq( eval( factor(%%-simplify(diff( (a*x+b)^m, x$n ))), m=13), n=1..4 ) assuming m > 0;

One problem that the referee didn't tell us about: gfun:-holexprtodiffeq( A/(x-alpha)^p, y(x) ); Darn. But, we can avoid the problem here. sdiff( A/(x-alpha)^(p), x, n ); seq(simplify(%-diff(A/(x-alpha)^p,x$n)),n=1..3); sdiff( 1/x^p, x, n );

phi := 1/(1-x^2); diff( phi, x ); diff( phi, x$5 ); diff( phi, x$10 ); phi2 := convert( phi, parfrac, x ); diff( phi2, x ); diff( phi2, x$10 ); dnphi := sdiff( phi, x, n ); QyQ+SSZkbnBoaUc2Ii1JJWRpZmZHJSpwcm90ZWN0ZWRHNiRJJHBoaUdGJS1JIiRHRig2JEkieEdGJUkibkdGJSIiIg==

phi1 := (x+1)/(x^2-4); phi2 := x^4/(x-1)/(x-2); phi3 := 4*x/((x-1)^2*(x+2)); sdiff( phi1, x, n ); QyQtSSVkaWZmRyUqcHJvdGVjdGVkRzYkSSVwaGkxRzYiLUkiJEdGJTYkSSJ4R0YoSSJuR0YoIiIi sdiff( phi2, x, n ); QyQtSSVkaWZmRyUqcHJvdGVjdGVkRzYkSSVwaGkyRzYiLUkiJEdGJTYkSSJ4R0YoSSJuR0YoIiIi trial := sdiff( phi3, x, n ); QyQtSSVkaWZmRyUqcHJvdGVjdGVkRzYkSSVwaGkzRzYiLUkiJEdGJTYkSSJ4R0YoSSJuR0YoIiIi seq( normal(trial-diff(phi3,x$n)), n=1..5 ); normal( eval(trial,n=0)-phi3 );

phi := x^3/(x^2-1); _EnvFallingNotation := GAMMA; dn := sdiff( phi, x, n ); QyQtSSVkaWZmRyUqcHJvdGVjdGVkRzYkSSRwaGlHNiItSSIkR0YlNiRJInhHRihJIm5HRigiIiI= eval( dn, n=1 ); eval( %, x=0 ); I said in 2004 that "The following is not true if n=1, because x^(1-n) = 1 if n=1 (usual double limit problem)", and this is still correct. Substituting x=0 then n=1 gives a number that is not zero, whereas substituting n=1 and then n=0 gives zero. simplify(eval(dn,x=0)); QyQtSSVldmFsRyUqcHJvdGVjdGVkRzYkLCQqKCMiIiIiIiNGKiwmKSEiIkkibkc2IkYqRipGLkYqLUkmR0FNTUFHNiRGJUkoX3N5c2xpYkdGMDYjLCZGL0YqRipGKkYqRiovRi9GKkYq simplify( ) assuming n::even; simplify( eval(dn,x=0) ) assuming n::odd, n>1; _EnvFallingNotation := pochhammer; QyQ+SSNkbkc2Ii1JJWRpZmZHJSpwcm90ZWN0ZWRHNiRJJHBoaUdGJS1JIiRHRig2JEkieEdGJUkibkdGJSIiIg== QyQtSSlhc3N1bWluZ0dJKF9zeXNsaWJHNiI2JDcjLUkpc2ltcGxpZnlHNiQlKnByb3RlY3RlZEdGJTYjLUklZXZhbEdGLDYkSSNkbkdGJi9JInhHRiYiIiE3JCdJIm5HRiZJJG9kZEdGLDIiIiJGN0Y6 LUklbXJvd0c2Iy9JK21vZHVsZW5hbWVHNiJJLFR5cGVzZXR0aW5nR0koX3N5c2xpYkdGJzYjLUkjbWlHRiQ2JVEhRicvJSdpdGFsaWNHUSV0cnVlRicvJSxtYXRodmFyaWFudEdRJ2l0YWxpY0Yn

sdiff( sin(x), x, r ); QyQtSSVkaWZmRyUqcHJvdGVjdGVkRzYkLUkkc2luRzYkRiVJKF9zeXNsaWJHNiI2I0kieEdGKy1JIiRHRiU2JEYtSSJyR0YrIiIi combine(%%,trig); That is the answer returned by Macsyma 2.4 (Richard Fateman, personal communication) sdiff( x*sin(x), x, n ); QyQtSSVkaWZmRyUqcHJvdGVjdGVkRzYkKiZJInhHNiIiIiItSSRzaW5HNiRGJUkoX3N5c2xpYkdGKTYjRihGKi1JIiRHRiU2JEYoSSJuR0YpRio= sdiff( x^2/(x^4-x^3+x-1), x, n ); QyQtSSVkaWZmRyUqcHJvdGVjdGVkRzYkKiZJInhHNiIiIiMsKiokRigiIiUiIiIqJEYoIiIkISIiRihGLkYxRi5GMS1JIiRHRiU2JEYoSSJuR0YpRi4= sdiff( exp(x), x, r ); QyQtSSVkaWZmRyUqcHJvdGVjdGVkRzYkLUkkZXhwRzYkRiVJKF9zeXNsaWJHNiI2I0kieEdGKy1JIiRHRiU2JEYtSSJyR0YrIiIi QyQtSSVkaWZmRyUqcHJvdGVjdGVkRzYkKUkiZUc2IkkieEdGKS1JIiRHRiU2JEYqSSJyR0YpIiIi QyQtSSZzZGlmZkc2IjYlKUkiZUdGJUkieEdGJUYpSSJyR0YlIiIi sdiff( exp(3*x), x, r ); QyQtSSVkaWZmRyUqcHJvdGVjdGVkRzYkLUkkZXhwRzYkRiVJKF9zeXNsaWJHNiI2IywkSSJ4R0YrIiIkLUkiJEdGJTYkRi5JInJHRisiIiI= sdiff( exp( 5*x+7 )^3, x, q ); QyQtSSVkaWZmRyUqcHJvdGVjdGVkRzYkKiQtSSRleHBHNiRGJUkoX3N5c2xpYkc2IjYjLCZJInhHRiwiIiYiIigiIiIiIiQtSSIkR0YlNiRGL0kicUdGLEYy Aha! Another place where functionality was removed. sdiff( sin(x)^3, x, n ); QyQtSSVkaWZmRyUqcHJvdGVjdGVkRzYkKiQtSSRzaW5HNiRGJUkoX3N5c2xpYkc2IjYjSSJ4R0YsIiIkLUkiJEdGJTYkRi5JIm5HRiwiIiI= And another. Well, is the first answer right? LUklbXJvd0c2Iy9JK21vZHVsZW5hbWVHNiJJLFR5cGVzZXR0aW5nR0koX3N5c2xpYkdGJzYjLUkjbWlHRiQ2JVEhRicvJSdpdGFsaWNHUSV0cnVlRicvJSxtYXRodmFyaWFudEdRJ2l0YWxpY0Yn trm := simplify(combine(,trig)); trm0 := evalc(simplify( eval(trm,x=0) )) assuming n::odd; sum( eval(trm0,n=2*k+1)*x^(2*k+1)/(2*k+1)!, k=0..infinity ); expand(%); combine( % - sin(x)^3, trig ); This one is a bit disappointing: gfun:-holexprtodiffeq fails here sdiff( log(x), x, n ); But we may get around it by using an interesting trick: the nth derivative is the n-1st derivative of the 1st derivative... sdiff( diff(log(x),x), x, n-1 ); QyQtSSVkaWZmRyUqcHJvdGVjdGVkRzYkLUkkbG9nRzYkRiVJKF9zeXNsaWJHNiI2I0kieEdGKy1JIiRHRiU2JEYtSSJuR0YrIiIi And again here! sdiff( Ei(x), x, n ); ein := sdiff( diff(Ei(x),x), x, n-1 ); Something different happens in the incorporated version: QyQtSSVkaWZmRyUqcHJvdGVjdGVkRzYkLUkjRWlHNiRGJUkoX3N5c2xpYkc2IjYjSSJ4R0YrLUkiJEdGJTYkRi1JIm5HRisiIiI= QyQtSSlzaW1wbGlmeUc2JCUqcHJvdGVjdGVkR0koX3N5c2xpYkc2IjYjSSIlR0YoIiIi But the original trick still works :-) QyQtSSVkaWZmRyUqcHJvdGVjdGVkRzYkLUYkNiQtSSNFaUc2JEYlSShfc3lzbGliRzYiNiNJInhHRi1GLy1JIiRHRiU2JEYvLCZJIm5HRi0iIiIhIiJGNUY1 But the trick fails here: sdiff( sin(x^2), x, n ); sdiff( diff(sin(x^2),x),x,n-1); QyQtSSVkaWZmRyUqcHJvdGVjdGVkRzYkLUkkc2luRzYkRiVJKF9zeXNsaWJHNiI2IyokSSJ4R0YrIiIjLUkiJEdGJTYkRi5JIm5HRisiIiI= A genuine improvement in the incorporation: that complicated answer is genuinely useful. The trick fails in the following, though. sdiff( DESol(x*diff(y(x),x) + y(x), y(x) ), x, n ); sdiff( int( DESol(x*diff(y(x),x) + y(x), y(x) ), x ), x, n+1 ); LUklbXJvd0c2Iy9JK21vZHVsZW5hbWVHNiJJLFR5cGVzZXR0aW5nR0koX3N5c2xpYkdGJzYjLUkjbWlHRiQ2JVEhRicvJSdpdGFsaWNHUSV0cnVlRicvJSxtYXRodmFyaWFudEdRJ2l0YWxpY0Yn QyQtSSVkaWZmRyUqcHJvdGVjdGVkRzYkLUkmREVTb2xHNiRGJUkoX3N5c2xpYkc2IjYkLCYqJkkieEdGKyIiIi1GJDYkLUkieUdGKzYjRi9GL0YwRjBGM0YwRjMtSSIkR0YlNiRGL0kibkdGK0Yw QyQtSSdmYWN0b3JHNiQlKnByb3RlY3RlZEdJKF9zeXNsaWJHNiI2Iy1JJWRpZmZHRiY2JC1JKUxhbWJlcnRXR0YlNiNJInhHRigtSSIkR0YmNiRGMCIiJSIiIg== QyQtSSdmYWN0b3JHNiQlKnByb3RlY3RlZEdJKF9zeXNsaWJHNiI2Iy1JJWRpZmZHRiY2JC1JLFdyaWdodG9tZWdhR0YlNiNJInhHRigtSSIkR0YmNiRGMCIiJSIiIg== QyQtSSVkaWZmRyUqcHJvdGVjdGVkRzYkLUksV3JpZ2h0b21lZ2FHNiRGJUkoX3N5c2xpYkc2IjYjSSJ4R0YrLUkiJEdGJTYkRi1JIm5HRisiIiI= QyQtSSlzaW1wbGlmeUc2JCUqcHJvdGVjdGVkR0koX3N5c2xpYkc2IjYjLCYtSSxXcmlnaHRvbWVnYUdGJTYjSSJ4R0YoIiIiLUkjbG5HRiU2I0YrRi9GLw== QyQtJkkmZXZhbGZHJSpwcm90ZWN0ZWRHNiMiI102Iy1JLFdyaWdodG9tZWdhRzYkRiZJKF9zeXNsaWJHNiI2IyQiIz8hIiIiIiI=

We may differentiate other sums, such as Mittag-Leffler expansions (infinite partial fraction expansions for functions with an infinite number of poles) sdiff( Sum( 1/(x-j)^2, j=-infinity..infinity ), x, n ); QyQtSSVkaWZmRyUqcHJvdGVjdGVkRzYkLUkkU3VtRzYkRiVJKF9zeXNsaWJHNiI2JCokKSwmSSJ4R0YrIiIiSSJqR0YrISIiIiIjRjMvRjI7LCRJKWluZmluaXR5R0YlRjNGOC1JIiRHRiU2JEYwSSJuR0YrRjE= value(%); convert(%,GAMMA); convert(%,binomial); value(%); Not much happening. But "not much" is more than "nothing at all"! eval(%,n=3); eval(%%,n=6); eval(%%%,n=0); sum( 1/(x-j)^2, j=-infinity..infinity );

Here's another trick---based on extensive hard coding by Edgardo Cheb-Terrab, it seems! Thank you, Edgardo! sdiff( BesselJ(0,x), x, n ); QyQtSSVkaWZmRyUqcHJvdGVjdGVkRzYkLUkoQmVzc2VsSkc2JEYlSShfc3lzbGliRzYiNiQiIiFJInhHRistSSIkR0YlNiRGLkkibkdGKyIiIg== LUklbXJvd0c2Iy9JK21vZHVsZW5hbWVHNiJJLFR5cGVzZXR0aW5nR0koX3N5c2xpYkdGJzYlLUkjbWlHRiQ2JVEkc2VxRicvJSdpdGFsaWNHUSV0cnVlRicvJSxtYXRodmFyaWFudEdRJ2l0YWxpY0YnLUkobWZlbmNlZEdGJDYkLUYjNi0tSSNtb0dGJDYwUSJ+RicvRjNRJ25vcm1hbEYnLyUmZmVuY2VHUSZmYWxzZUYnLyUqc2VwYXJhdG9yR0ZCLyUpc3RyZXRjaHlHRkIvJSpzeW1tZXRyaWNHRkIvJShsYXJnZW9wR0ZCLyUubW92YWJsZWxpbWl0c0dGQi8lJ2FjY2VudEdGQi8lJWZvcm1HUSFGJy8lJ2xzcGFjZUdRJDBlbUYnLyUncnNwYWNlR0ZULyUobWluc2l6ZUdRIjFGJy8lKG1heHNpemVHUSlpbmZpbml0eUYnLUYsNiVRKXNpbXBsaWZ5RidGL0YyLUY2NiQtRiM2KC1GLDYlUSZ2YWx1ZUYnRi9GMi1GNjYkLUYjNiQtRiw2JVElZXZhbEYnRi9GMi1GNjYkLUYjNipGOi1JKG1hY3Rpb25HRiQ2JC1JKm12ZXJiYXRpbUdGJDYjRlEvJSthY3Rpb250eXBlR1E6bWFwbGVzb2Z0LmNvbTpsYWJlbChMMTk1KUYnLUY7NjBRIixGJ0Y+RkAvRkRGMUZFRkdGSUZLRk0vRlBRJmluZml4RidGUi9GVlEzdmVyeXRoaWNrbWF0aHNwYWNlRidGV0ZaRjotRiw2JVEibkYnRi9GMi1GOzYwUSI9RidGPkZARkNGRUZHRklGS0ZNRmlwL0ZTUS90aGlja21hdGhzcGFjZUYnL0ZWRmRxRldGWi1GLDYlUSJpRidGL0YyRjpGPkY+LUY7NjBRKiZ1bWludXMwO0YnRj5GQEZDRkVGR0ZJRktGTUZpcC9GU1EwbWVkaXVtbWF0aHNwYWNlRicvRlZGXXJGV0ZaLUYsNiVRJWRpZmZGJ0YvRjItRjY2JC1GIzYoLUYsNiVRKEJlc3NlbEpGJy9GMEZCRj4tRjY2JC1GIzYlLUkjbW5HRiQ2JFEiMEYnRj5GZXAtRiw2JVEieEYnRi9GMkY+RmVwRmJzLUY7NjBRIiRGJ0Y+RkBGQ0ZFRkdGSUZLRk1GT0ZSRlVGV0ZaRmZxRj5GOkY+RmVwRjpGZnFGYHEtRl9zNiRGWUY+LUY7NjBRIy4uRidGPkZARkNGRUZHRklGS0ZNL0ZQUShwb3N0Zml4RidGXHJGVUZXRlotRl9zNiRRIjNGJ0Y+RjpGPi1GOzYwUSI7RidGPkZARmhwRkVGR0ZJRktGTUZpcEZSRmVxRldGWg== The old way (using Edgardo's trick) was more complicated and less informative, but it also worked. S := convert( BesselJ(0,x), Sum ); _EnvFallingNotation:= GAMMA; sdiff( S, x, n ); The Pochhammer notation would have taken care of this automatically, but the GAMMA notation has a problem because we divide by infinity to get zero...but the following hack fixes it for this problem :-) subsop(2=(_k1=ceil(n/2)..infinity),%); value(%); wow := %: eval( wow, n=3 ); convert(%,StandardFunctions); %-diff(BesselJ(0,x),x$3); normal(%); _EnvFallingNotation := pochhammer; _EnvTrySum := true; notquite := sdiff( BesselJ(0,x), x, n ); convert(eval(%,n=3),StandardFunctions); NumericEventHandler(division_by_zero=(()->args[3])); convert( eval(notquite,n=3), StandardFunctions);

sdiff( hypergeom([1,2],[3],x), x, n ); QyQtSSVkaWZmRyUqcHJvdGVjdGVkRzYkLUkqaHlwZXJnZW9tR0koX3N5c2xpYkc2IjYlNyQiIiIiIiM3IyIiJEkieEdGKi1JIiRHRiU2JEYxSSJuR0YqRi0= sdiff( hypergeom( [a,b], [d,e], x ), x, n ); convert( hypergeom( [a,b], [d,e], x ), Sum ); sdiff( %, x, n ); value(%); Sum( mul( pochhammer(a[i],k), i=1..3)/mul( pochhammer(b[i],k),i=1..3)*x^k/k!, k=0..infinity ); sdiff(%,x,n); value(%); _EnvFallingNotation := pochhammer; Sum( mul( pochhammer(a[i],k), i=1..3)/mul( pochhammer(b[i],k),i=1..3)*x^k/k!, k=0..infinity ); sdiff(%,x,n); subsop(2=(k=n..infinity),%); value(%); convert( pochhammer(1,n), GAMMA ); aha := %%; seq( eval( aha, [a[1]=1,a[2]=2,a[3]=3,b[1]=4,b[2]=1,b[3]=1,n=k] ), k=0..4 ); LUklbXJvd0c2Iy9JK21vZHVsZW5hbWVHNiJJLFR5cGVzZXR0aW5nR0koX3N5c2xpYkdGJzYjLUkjbWlHRiQ2JVEhRicvJSdpdGFsaWNHUSV0cnVlRicvJSxtYXRodmFyaWFudEdRJ2l0YWxpY0Yn LUklbXJvd0c2Iy9JK21vZHVsZW5hbWVHNiJJLFR5cGVzZXR0aW5nR0koX3N5c2xpYkdGJzYjLUkjbWlHRiQ2JVEhRicvJSdpdGFsaWNHUSV0cnVlRicvJSxtYXRodmFyaWFudEdRJ2l0YWxpY0Yn