<<RISC`HolonomicFunctions`

HolonomicFunctions Package version 1.7.3 (21-Mar-2017)
written by Christoph Koutschan
Copyright Research Institute for Symbolic Computation (RISC),
Johannes Kepler University, Linz, Austria

Basic Examples

Transformation of 2F1 functions

HoldForm[Hypergeometric2F1[-2n,2n+2a+1,a+1,(1-x)/2]==Hypergeometric2F1[-n,n+a+1/2,a+1,1-x^2]]//TraditionalForm

Mathematica does not know this transformation:

FullSimplify[Hypergeometric2F1[-2n,2n+2a+1,a+1,(1-x)/2]-Hypergeometric2F1[-n,n+a+1/2,a+1,1-x^2],Element[{a,n},Integers]&&n≥0&&a≥0]

Computing annihilating ideals for both expressions, we find that they agree.

lhs=Annihilator[Hypergeometric2F1[-2n,2n+2a+1,a+1,(1-x)/2],{S[n],S[a]}]

rhs=Annihilator[Hypergeometric2F1[-n,n+a+1/2,a+1,1-x^2],{S[n],S[a]}]

GBEqual[lhs,rhs]

We compare initial values. Which ones we have to consider, tells us the following command:

sing=AnnihilatorSingularities[lhs,{0,0}]

(Hypergeometric2F1[-2n,2n+2a+1,a+1,(1-x)/2]==Hypergeometric2F1[-n,n+a+1/2,a+1,1-x^2])/.(First/@sing)

This concludes the proof.

Strang's Integral  (∫_(-1)^1)(((((LegendreP[2 k+1, x])/x)))^2)dx==2

This integral can be found in the SFB proposal 2004, Section 7.2.2, Example 5.

TraditionalForm[HoldForm[Integrate[(LegendreP[2k+1,x]/x)^2,{x,-1,1}]==2]]

We start by computing an annihilating ideal for the integrand.

ann=Annihilator[(LegendreP[2k+1,x]/x)^2,{Der[x],S[k]}]

with Chyzak's algorithm

Next, compute a creative telescoping relation with Chyzak's algorithm:

ct=CreativeTelescoping[ann,Der[x],S[k]]

Check the delta part. It does not contribute any inhomogeneous part from telescoping.

delta=ApplyOreOperator[ct[[2,1]],(LegendreP[2k+1,x]/x)^2]//Simplify

(delta/.x->1)-(delta/.x->-1)

FullSimplify[%,Assumptions->Element[k,Integers]&&k>0]

Hence, S[k]-1 is a recurrence for the integral. The same result can be more easily (i.e., completely automatically) obtained by

Annihilator[Integrate[(LegendreP[2*k+1,x]/x)^2,{x,-1,1}],S[k],Assumptions->Element[k,Integers]&&k≥0]

One initial value suffices to obtain the evaluation of the integral:

Integrate[(LegendreP[1,x]/x)^2,{x,-1,1}]

with Takayama's algorithm

Takayama's algorithm delivers a bigger recurrence (which is a multiple of the above). But one has to be careful since the integral does not have natural boundaries.

Takayama[ann,{x}]

With the option Extended → True, we can force the command  Takayama to carry the delta parts through the Gröbner basis computation.

{tak,{deltas}}=Takayama[ann,{x},Extended->True]

But again, no inhomogeneous part arises.

ApplyOreOperator[deltas,(LegendreP[2k+1,x]/x)^2]//Simplify

(%/.x->1)-(%/.x->-1)

The recurrence obtained by Takayama is of order 2 and is a multiple of the one computed by Chyzak's algorithm.

OreReduce[tak,Annihilator[1,S[k]]]

with Zeilberger's slow algorithm (elimination by ansatz)

The command FindRelation can be used to compute an operator in the left ideal ann that is free of x. The methods used in FindRelation are ansatz ans coefficient comparison. Since the resulting operator is of considerable size, the computation takes a while.

Timing[rel=First[FindRelation[ann,Eliminate->x]];]

ByteCount[rel]

By substituting Der[x] → 0, we obtain the principal part. it is a recurrence of order 7, that again is a multiple of the minimal recurrence computed above:

pp=OrePolynomialSubstitute[rel,{Der[x]->0}];

Support[pp]

OreReduce[pp,Annihilator[1,S[k]]]

The delta part is obtained by subtracting the principal part from the original relation and divide it by Der[x]. This division cannot be done in an Ore algebra. That's why we switch to standard Mathematica polynomials for this step.

delta=Together[1/Der[x]*Normal[ChangeOreAlgebra[rel-pp,OreAlgebra[Der[x],x,S[k]]]]];

PolynomialQ[delta,Der[x]]

For the subsequent computations, it is of advantage to bring the delta part to normal form with respect to the annihilator ann. Then it is again seen that the delta part vanishes at the integration bounds (and hence the principal part computed above is a recurrence for the integral).

delta=OreReduce[delta,Together[ann]];

Together[ApplyOreOperator[delta,(LegendreP[2k+1,x]/x)^2]]/.x->1

FullSimplify[Together[ApplyOreOperator[delta,(LegendreP[2k+1,x]/x)^2]]/.x->-1,Element[k,Integers]]

with Zeilberger's slow algorithm (elimination by Gröbner bases)

The elimination of the indeterminate x with Gröbner bases does not work very well in this instance. The computation takes extremely long, and the recurrence that we obtain is of order 8 even.

TimeConstrained[gb=OreGroebnerBasis[ann,OreAlgebra[x,Der[x],S[k]],MonomialOrder->EliminationOrder[1]];,60]

With different options we can get at least some result, but still need a considerable amount of time. The option Incomplete → True has the effect that the Gröbner basis computation is interrupted as soon as an x-free element is found, and this element is returned. The method "elimination" is a greedy strategy towards an x-free operator, in the sense that always a pair is chosen whose S-polynomial has the least degree in x.

Timing[xfree=OreGroebnerBasis[ann,OreAlgebra[x,Der[x],S[k]],MonomialOrder->EliminationOrder[1],Method->"elimination",Incomplete->True];]

ByteCount[xfree]

pp=OrePolynomialSubstitute[xfree,{Der[x]->0}];

Support[pp]

OreReduce[pp,Annihilator[1,S[k]]]

Contiguous Relations for 2F1

We want to find all contiguous relations for the hypergeometric 2F1 function. We consider all relation between the function itself and two of its immediate neighbours, i.e., where one of the parameters a,b,c is decreased or increased by one. The possible combinations are (where we eliminated negative shifts):

supps=NormalizeVector[Append[#,1]]&/@Subsets[{S[a],1/S[a],S[b],1/S[b],S[c],1/S[c]},{2}]

Next we compute an annihilating ideal of 2F1 with respect to the shifts operators in question.

ann=Annihilator[Hypergeometric2F1[a,b,c,x],{S[a],S[b],S[c]}]

The command FindRelation searches for elements in the annihilating ideal that have the specified support.

ApplyOreOperator[Factor[FindRelation[ann,Support->supps]],F[a,b,c,x]]//TableForm

Module backwards shifting (and some typos) this list corresponds exactly to the one given in Gauß's paper "Disquisitiones generales circa series infinitam ...", 1813.

Recurrence obtained from creative telescoping can overshoot arbitrarily

The following sum is a classical example for the fact creative telescoping does not necessarily deliver the minimal recurrence for the sum/integral. By looping over the order of the principal part, it is guaranteed that the creative telescoping relation with the smallest principal part is found. However, it might well be, that the sum itself allows a shorter recurrence.

HoldForm[Sum[(-1)^k *Binomial[n,k]* Binomial[d*k,n],{k,0,n}]==(-d)^n]//TraditionalForm

Obviously, for any integer d, this sum satisfies a first-order recurrence. But for d>1 creative telescoping will deliver a recurrence of order d-1, showing that we can overshoot arbitrarily much.

Table[LeadingPowerProduct[CreativeTelescoping[(-1)^k *Binomial[n,k]* Binomial[d*k,n],S[k]-1,S[n]][[1,1]]],{d,2,10}]

Note: Also Zeilberger's slow or Takayama's algorithm do not improve this situation.

Second note: Also for integrals, similar things can happen, as the following example demonstrates:

Table[LeadingPowerProduct[First[Annihilator[Integrate[1/(x^2+1)^(n+1)/(x^d+1),{x,0,Infinity}],S[n],Assumptions->n≥0]]],{d,0,10}]

Takayama runs for ever due to extension/contraction

We want to integrate a simple rational function with Takayama's algorithm.

ann1=Annihilator[1/(x^2+y^2),{Der[x],Der[y]}]

Takayama's algorithm runs into an infinite loop (due to extension/contraction).

TimeConstrained[Takayama[ann1,{x}],10]

Singular's output of the full annihilator of 1/(x^2+y^2) in the polynomial Weyl algebra (using the command annfs of the library dmod.lib):

ann2=ToOrePolynomial[{y*Der[x]-x*Der[y],x*Der[x]+y*Der[y]+2,x^2*Der[y]+y^2*Der[y]+2*y},OreAlgebra[Der[x],Der[y]]]

We convince ourselves that these three elements generate the same ∂-finite ideal:

OreGroebnerBasis[ann2]

With the full annihilator, Takayama's algorithm is guaranteed to terminate.

Takayama[ann2,{x}]

(∑_(n=0)^∞)(((Subscript[H, n])(x) (Subscript[H, n])(y) (t^n))/(n!))==(1/(Sqrt[1-4 (t^2)])) (e^((y^2)-((((y-2 x t))^2)/(1-4 (t^2)))))

This example is taken from A=B, p. 197.

HoldForm[Sum[HermiteH[n,x]*HermiteH[n,y]*t^n/n!,{n,0,Infinity}]==(1-4t^2)^(-1/2)*Exp[y^2-(y-2*x*t)^2/(1-4t^2)]]//TraditionalForm

Annihilator[Sum[HermiteH[n,x]*HermiteH[n,y]*t^n/n!,{n,0,Infinity}],{Der[t],Der[x],Der[y]},Inhomogeneous->True]

If we believe that the occurring limits evaluate to 0, or even stronger, if we assume natural boundaries, we can just ignore the delta part, resp. use Takayama's algorithm:

lhs1=First[CreativeTelescoping[HermiteH[n,x]*HermiteH[n,y]*t^n/n!,S[n]-1,{Der[t],Der[x],Der[y]}]]

lhs2=Takayama[Annihilator[HermiteH[n,x]*HermiteH[n,y]*t^n/n!,{S[n],Der[t],Der[x],Der[y]}],{n}]

Alternatively, we could use the command DFiniteRE2DE. But it gives a smaller annihilating ideal (which is contained in the previously computed).

DFiniteRE2DE[Annihilator[HermiteH[n,x]*HermiteH[n,y]/n!,{S[n],Der[x],Der[y]}],{n},{t}]

OreReduce[%,lhs1]

For finding the right-hand side, we could now take either of the three first-order differential equations and solve it.

With[{var=First[Normal[LeadingPowerProduct[#]]]},DSolve[ApplyOreOperator[#,f[var]]==0,f[var],var]]&/@lhs1

However, concerning initial values, it is the simplest to consider t=0, because then the sum simplifies to 1.

solution=f[t]/.First[Last[%]]

Solve[(%/.t->0)==1,C[1]]

solution=FullSimplify[solution/.First[%]]

A hypergeometric double sum

This sum can be found in "Peter Paule: Parallel Algorithms in Computer Algebra - Zeilberger's Holonomic Systems Approach".

TraditionalForm[HoldForm[Sum[(-1)^(k1+k2)*Binomial[n1,k1]Binomial[n2,k2]Binomial[k1+k2,k1],{k1,0,n1},{k2,0,n2}]==KroneckerDelta[n1,n2]]]

F[n1_,n2_,k1_,k2_]:=(-1)^(k1+k2)Binomial[n1,k1]Binomial[n2,k2]Binomial[k1+k2,k1]

Such examples are a little bit problematic to do automatically with Chyzak's algorithm; since all computations are done in a rational Ore algebra, some polynomial factor that matters may be cancelled. We get a result that is true "nearly everywhere".

Annihilator[Sum[F[n1,n2,k1,k2],{k1,0,n1},{k2,0,n2}],{S[n1],S[n2]}]

The bad factor can be seen when looking at the certificate:

CreativeTelescoping[First[CreativeTelescoping[F[n1,n2,k1,k2],S[k1]-1,{S[k2],S[n1],S[n2]}]],S[k2]-1,{S[n1],S[n2]}]

Alternatively, we can employ other algorithms that are more robust against such misbehaviour. First, Zeilberger's slow algorithm:

ann=Annihilator[F[n1,n2,k1,k2],{S[n1],S[n2],S[k1],S[k2]}]

gb=OreGroebnerBasis[ann,OreAlgebra[k1,k2,S[k1],S[k2],S[n1],S[n2]],MonomialOrder->EliminationOrder[2]]

LeadingPowerProduct/@gb

OrePolynomialSubstitute[Take[gb,4],{S[k1]->1,S[k2]->1}]

lhs=Factor[OreGroebnerBasis[%,OreAlgebra[S[n1],S[n2]],MonomialOrder->DegreeLexicographic]]

Secondly, with Takayama's algorithm, which yields the same solution.

lhs=Factor[Takayama[ann,{k1,k2}]]

Checking for singularities reveals that on the diagonal something goes wrong.

AnnihilatorSingularities[%,{0,0}]

A closer inspection of the recurrences and initial values tells us that the sum is 0 at all non-singular points. The diagonal can be treated separately, to find that there we get 1.

Annihilator[Sum[F[n1,n1,k1,k2],{k1,0,n1},{k2,0,n1}],{S[n1]}]

Table[Sum[F[n1,n1,k1,k2],{k1,0,n1},{k2,0,n1}],{n1,0,1}]

Apagodu: Series evaluation...

Moa Apagodu: Series evaluation of a quartic integral and generalizations
The computations in this paper are easily reproduced with HolonomicFunctions.

Theorem 1:

reca=DFiniteDE2RE[Annihilator[Integrate[1/(x^4+2*a*x^2+1)^(m+1),{x,0,Infinity}],Der[a],Assumptions->m≥0],{a},{l}]

ApplyOreOperator[{1,Der[a]},1/(x^4+2*a*x^2+1)^(m+1)]/.a->0

Integrate[%,{x,0,Infinity},Assumptions->m>0]

RSolve[Prepend[Thread[{a[0],a[1]}==%],ApplyOreOperator[First[reca],a[l]]==0],a[l],l]

FullSimplify[%[[1,1,2]]]

Test:

Annihilator[1/4*Sum[(-1)^l*2^l*(l/2-3/4)!*(m+l/2-1/4)!/l!/m!*a^l,{l,0,Infinity}],{Der[a],S[m]},Assumptions->Abs[a]<1]

Annihilator[Integrate[1/(x^4+2*a*x^2+1)^(m+1),{x,0,Infinity}],{Der[a],S[m]},Assumptions->m>0]

Theorem 2:

Factor[CreativeTelescoping[1/(x^6+2*a*x^3+1)^(m+1),Der[x],Der[a]]]

Theorem 3:

Factor[CreativeTelescoping[1/(x^8+2*a*x^4+1)^(m+1),Der[x],Der[a]]]

Theorem 4:

Factor[CreativeTelescoping[1/(z^2+2*a*z+1)^(m+1)/n/z^( 1-1/n),Der[z],Der[a]]]

Theorem 5:

Factor[CreativeTelescoping[1/(z^2+n*a*z+1)^(m+1)/n/z^( 1-1/n),Der[z],Der[a]]]

Theorem 6:

Factor[CreativeTelescoping[1/(z^2+α*a*z+1)^(m+1)/n/z^( 1-1/n),Der[z],Der[a]]]

7 proofs of the trinomial identity

HoldForm[Sum[Binomial[n,j]Binomial[j,i]x^i y^(j-i) z^(n-j),{j,0,n},{i,0,j}]==(x+y+z)^n]//TraditionalForm

We precompute the annihilator of the summand.

ann=Annihilator[Binomial[n,j]Binomial[j,i]x^i y^(j-i) z^(n-j),{S[i],S[j],S[n]}]

Below we list seven commands which all deliver the first-order recurrence that proves the trinomial identity.

OrePolynomialSubstitute[First[OreGroebnerBasis[ann,OreAlgebra[i,j,S[i],S[j],S[n]],MonomialOrder->EliminationOrder[2]]],{S[i]->1,S[j]->1}]

OrePolynomialSubstitute[FindRelation[ann,Eliminate->{i,j}],{S[i]->1,S[j]->1}]

Takayama[ann,{i,j}]

First[CreativeTelescoping[First[CreativeTelescoping[ann,S[i]-1,{S[j],S[n]}]],S[j]-1,S[n]]]

First[Annihilator[Sum[Binomial[n,j]Binomial[j,i]x^i y^(j-i) z^(n-j),{j,0,n},{i,0,j}],S[n],Inhomogeneous->True]]

First[CreativeTelescoping[Annihilator[Sum[Binomial[n,j]Binomial[j,i]x^i y^(j-i) z^(n-j),{i,0,j}],{S[j],S[n]}],S[j]-1,S[n]]]

First[FindCreativeTelescoping[ann,{S[i]-1,S[j]-1}]]

Some identities for Jacobi polynomials

In Veronika Pillwein's thesis "Computer algebra tools for special functions in high order finite element methods", we find the following identities for Jacobi polynomials. We start with computing (actually, looking up) an annihilating ideal for Jacobi polynomials.

ann=Annihilator[JacobiP[n,a,b,x],{S[a],S[b],S[n],Der[x]}]

For some identities, we give several ways how to obtain them.

Formula (4.5)

TraditionalForm[HoldForm[D[JacobiP[n,a,b,x],x]==(n+a+b+1)/2*JacobiP[n-1,a+1,b+1,x]]]

FindRelation[ann,Support->{Der[x]*S[n],S[a]*S[b]}]

Formula (4.6)

TraditionalForm[HoldForm[(1-x)*D[JacobiP[a,b,n,x],x]==a*JacobiP[n,a,b,x]-(n+a)*JacobiP[n,a-1,b+1,x]]]

Take[OreGroebnerBasis[ann,OreAlgebra[b,S[n],Der[x],S[b],S[a]],MonomialOrder->EliminationOrder[2]],{2}]

FindRelation[ann,Support->{S[a],S[b],S[a]Der[x]}]

Formula (4.7)

TraditionalForm[HoldForm[(2 n+a+b) JacobiP[n,a-1,b,x]==(n+a+b) JacobiP[n,a,b,x]-(n+b) JacobiP[n-1,a,b,x]]]

First[OreGroebnerBasis[ann,OreAlgebra[x,Der[x],S[b],S[a],S[n]],MonomialOrder->EliminationOrder[3]]]

FindRelation[ann,Support->{S[n],S[a],S[n]S[a]}]

FindRelation[ann,Eliminate->x,Pattern->{_,0,_,0}]

Formula (4.8)

TraditionalForm[HoldForm[(2 n+a+b) JacobiP[n,a,b-1,x]==(n+a+b) JacobiP[n,a,b,x]+(n+a) JacobiP[n-1,a,b,x]]]

First[OreGroebnerBasis[ann,OreAlgebra[x,Der[x],S[a],S[b],S[n]],MonomialOrder->EliminationOrder[3]]]

FindRelation[ann,Support->{S[n],S[b],S[n]S[b]}]

FindRelation[ann,Eliminate->x,Pattern->{0,_,_,0}]

Recurrence for an elliptic integral

This problem was posed by Stefan Böttner. Find a recurrence in j for the following integral

TraditionalForm[HoldForm[Integrate[Sqrt[1-m*Sin[x]^2]^(2j-1),{x,0,Pi/2}]]]

The integrand is not holonomic:

Annihilator[Sqrt[1-m*Sin[x]^2]^(2j-1),{Der[x],S[j]}]

But after a simple substitution u=sin(x), the integral becomes holonomic:

TraditionalForm[HoldForm[Integrate[Sqrt[1-m*u^2]^(2j-1)/Sqrt[1-u^2],{u,0,1}]]]

Annihilator[Integrate[Sqrt[1-m*u^2]^(2j-1)/Sqrt[1-u^2],{u,0,1}],S[j]]

Hence the desired recurrence is:

ApplyOreOperator[Factor[First[%]],Int[j]]==0

Monthly Problem 11457

This problem has been proposed by M.L. Glasser in October 2009.
For  real numbers  0 ≤ a ≤ b:

TraditionalForm[HoldForm[Integrate[ArcCos[x/Sqrt[(a+b)*x-a*b]],{x,a,b}]=="?"]]

An annihilating ideal for the integral can be computed completely automatically:

ann=Annihilator[Integrate[ArcCos[x/Sqrt[(a+b)*x-a*b]],{x,a,b}],{Der[a],Der[b]},Assumptions->a≥0&&b≥a]

Fortunately, we get two first-order differential equations, which can easily be solved:

DSolve[ApplyOreOperator[First[ann],f[b]]==0,f[b],b]

solution=f[b]/.First[%]

The above expression contains all factors that depend on b. In principle, C[1] still depends on a. So we compute a system of differential equations for "the original solution divided by this part of the solution". This system then will give us the only missing non-constant factors that depend on a (we see that there are none, so we have already the full solution).

DFiniteTimes[ann,Annihilator[1/solution,{Der[a],Der[b]}]]

By plugging in a=0 and b=1, the constant C[1] can be determined:

Solve[(solution==int[ArcCos[x/Sqrt[(a+b)*x-a*b]],{x,a,b}])/.{a->0,b->1}/.int->Integrate,C[1]]

solution=solution/.First[%]

Alternative (recursive) approach

We first derive only the differential equation in a. To determine its initial condition (for a=0), we recursively solve a similar problem, this time by computing a differential equation in b.

de1=First[Annihilator[Integrate[ArcCos[x/Sqrt[(a+b)*x-a*b]],{x,a,b}],{Der[a]},Assumptions->a≥0&&b≥a]]

Now recursively solve the simpler problem (for a=0):

de2=First[Annihilator[Integrate[ArcCos[x/Sqrt[b*x]],{x,0,b}],Der[b],Assumptions->b≥0]]

init2=f[1]==Integrate[ArcCos[Sqrt[x]],{x,0,1}]

DSolve[{ApplyOreOperator[de2,f[b]]==0,init2},f[b],b]

Now the first differential equation can be solved:

DSolve[{ApplyOreOperator[de1,f[a]]==0,f[0]==b*Pi/4},f[a],a]

More advanced examples

Wallis's Integral

TraditionalForm[HoldForm[Integrate[Cos[x]^(2m),{x,0,Pi/2}]==Pi/2^(2m+1)*Binomial[2m,m]]]

Mathematica knows about this integral:

Integrate[Cos[x]^(2m),{x,0,Pi/2},Assumptions->Element[m,Integers]&&m≥0]

The integrand is not ∂-finite:

Annihilator[Cos[x]^(2m),{S[m],Der[x]}]

Annihilator[Cos[x]^(2m),S[m]]

Annihilator[Cos[x]^(2m),Der[x]]

In this special example, creative telescoping works, although by theory this is not guaranteed. Therefore we have to use the option Support in order to force some action.

CreativeTelescoping[Cos[x]^(2m),Der[x],S[m],Support->{1,S[m],S[m]^2,S[m]^3}]

Takayama-style solution (sin(x) and cos(x) in the coefficients)

We can easily read off two relations for the integrand (although it is not ∂-finite --- that's why we allow sin(x) and cos(x) to appear in the coefficients):

ann=ToOrePolynomial[{f[m+1,x]==Cos[x]^2*f[m,x],Cos[x]*D[f[m,x],x]==-2*m*Sin[x]*f[m,x]},f[m,x]]

ApplyOreOperator[ann,(Cos[x])^(2m)]//Simplify

For Takayama's algorithm we have to include certain multiples (by sin(x) and cos(x)) to the set of relations (we fill them up to degree 2 in sin(x) and cos(x):

ann=Join[ann,{Sin[x]**ann[[2]],Cos[x]**ann[[2]]}]

In order to separate the delta part (that has the form Der[x]*(...)) we commute Der[x] to the left and set it to 0:

ann=ChangeOreAlgebra[ann,OreAlgebra[Der[x],Sin[x],Cos[x],S[m]]]

ann=OrePolynomialSubstitute[ann,{Der[x]->0}]

The above steps somehow let Mathematica forget about the relations between sine and cosine since they are now treated as independent elements of the algebra. Hence we shall not forget to include the Pythagorean identity:

AppendTo[ann,Sin[x]^2+Cos[x]^2-1]

Now it is a simple exercise to eliminate sin(x) and cos(x), but without multiplying by them:

ann[[3]]-(2m+1)**ann[[5]]

%-(2m+2)**ann[[1]]

Alternatively, the elimination of course can be achieved by Gröbner bases:

OreGroebnerBasis[ann,OreAlgebra[Sin[x],Cos[x],S[m]],MonomialOrder->EliminationOrder[2],ModuleBasis->{1,2}]

All the above steps are done when calling the command Takayama:

Takayama[ToOrePolynomial[{-Cos[x]^2+S[m],Cos[x] Der[x]+2 m Sin[x],Sin[x]^2+Cos[x]^2-1}],{Sin[x],Cos[x]}]

Takayama[ToOrePolynomial[{-Cos[x]^2+S[m],Cos[x] Der[x]+2 m Sin[x],Sin[x]^2+Cos[x]^2-1}],{Sin[x],Cos[x]},Extended->True]

substitution

Performing the substitution u=cos(x) we obtain

TraditionalForm[HoldForm[Integrate[Cos[x]^(2m),{x,0,Pi/2}]==Integrate[u^(2m)/Sqrt[1-u^2],{u,0,1}]]]

Integrate[u^(2m)/Sqrt[1-u^2],{u,0,1},Assumptions->Element[m,Integers]&&m≥0]

Annihilator[Integrate[u^(2m)/Sqrt[1-u^2],{u,0,1}],S[m],Assumptions->m≥0]

The irresistible quartic integral

Evaluate the following integral (it can be found as formula (7.2.1) in the book "Irresistible Integrals" by Boros/Moll:

TraditionalForm[HoldForm[Integrate[1/(x^4+2 a x^2+1)^(m+1),{x,0,Infinity}]]]

Mathematica cannot evaluate this integral:

Integrate[1/(x^4+2 a x^2+1)^(m+1),{x,0,Infinity},Assumptions->m≥0]

With Takayama's algorithm

ann=Takayama[Annihilator[1/(x^4+2 a x^2+1)^(m+1),{S[m],Der[a],Der[x]}],{x},Saturate->True]

Hence we get the following differential equation:

de=ApplyOreOperator[ann[[2]],Int[a]]

However, the integral does not have natural boundaries, and therefore we should also care about the delta part:

{ann,{delta}}=Takayama[Annihilator[1/(x^4+2 a x^2+1)^(m+1),{S[m],Der[a],Der[x]}],{x},Saturate->True,Extended->True]

ApplyOreOperator[delta,1/(x^4+2 a x^2+1)^(m+1)]

Limit[%,x->Infinity,Assumptions->m≥0]-(%/.x->0)

The inhomogeneous part simplifies to zero, and hence everything is fine, i.e., ann indeed annihilates the integral.

Find a recurrence in m:

OreGroebnerBasis[ann,OreAlgebra[Der[a],S[m]],MonomialOrder->EliminationOrder[1]]

rec=ApplyOreOperator[Factor[First[%]],Int[m]]

With Chyzak's algorithm

{ann,delta}=CreativeTelescoping[1/(x^4+2 a x^2+1)^(m+1),Der[x],{S[m],Der[a]}]

ApplyOreOperator[delta,1/(x^4+2 a x^2+1)^(m+1)]

Limit[%,x->Infinity,Assumptions->m≥0]-(%/.x->0)

These steps can be performed automatically by the following command:

Annihilator[Integrate[1/(x^4+2 a x^2+1)^(m+1),{x,0,Infinity}],{S[m],Der[a]},Assumptions->m≥0]

Exact solution using DSolve with initial values

Let Int[a] be the integral. We compute initial values for Int[0] and Int'[0].

i0=Integrate[1/(x^4+1)^(m+1),{x,0,Infinity},Assumptions->Re[m]≥0]

i1=Integrate[D[1/(x^4+2 a x^2+1)^(m+1),a]/.a->0,{x,0,Infinity},Assumptions->Re[m]≥0]

Now we use Mathematica's DSolve on the differential equation that we derived in order to get a closed form.

solution=Int[a]/.First[DSolve[{de==0,Int[0]==i0,Int'[0]==i1},Int[a],a]]

Some simplifications...

solution=FullSimplify[solution/.LegendreP[a_,b_,0]->Limit[LegendreP[a,b,zz],zz->0]/.LegendreQ[a_,b_,0]->Limit[LegendreQ[a,b,zz],zz->0],Assumptions->Element[m,Integers]]

Mathematica 7 returned the following formula as the output of the previous command:

-(((1+I) I^-m 2^(-1-m) (-1+a^2)^(-(1/4)-m/2) Sqrt[π] LegendreQ[m,1/2+m,a])/Gamma[1+m])

Boros/Moll's evaluation of the quartic

We reproduce the steps written in the book for finding the evaluation of the integral.

dlm=Annihilator[2^(-2m)*Sum[2^k*Binomial[2m-2k,m-k]*Binomial[m+k,m]*Binomial[k,l],{k,l,m}],{Der[a],S[l],S[m]}]

ann1=DFiniteTimes[Annihilator[a^l,{Der[a],S[l],S[m]}],dlm]

pma=First[CreativeTelescoping[ann1,S[l]-1,{S[m],Der[a]}]]

eval=DFiniteTimes[Annihilator[Pi/2/(2*(a+1))^(m+1/2),{S[m],Der[a]}],pma]

GBEqual[%,ann]

Final solution of Boros and Moll

In the above mentioned book, the following evaluation of the integral is given.

TraditionalForm[mollsol=Pi*JacobiP[m,m+1/2,-m-1/2,a]/2^(m+3/2)/(a+1)^(m+1/2)]

We can immediately prove that it is correct (and hence equivalent to the evaluation that we found)

Annihilator[mollsol,{S[m],Der[a]}]

GBEqual[%,ann]

FullSimplify[ApplyOreOperator[UnderTheStaircase[ann],mollsol]/.a->0/.m->0]

{Integrate[1/(x^4+1),{x,0,Infinity}],Integrate[D[1/(x^4+2 a x^2+1),a]/.a->0,{x,0,Infinity}]}

Symeonidis identities

We found these two identities in the paper of Flavia Stan: "Computer-assisted proofs of special function identities related to Poisson integrals".

Theorem 1

For real x and t such that Abs[x]<1 and Abs[t]<1, and integer n>=3:

TraditionalForm[HoldForm[Sum[Binomial[k+n-2,k]/Binomial[k+n/2-2,k]*t^k*Hypergeometric2F1[k,1-n/2,k+n/2,t^2]*GegenbauerC[k,(n-2)/2,x],{k,0,Infinity}]==((1-t^2)/(1-2*t*x+t^2))^(n-1)]]

lhs=Annihilator[Sum[Binomial[k+n-2,k]/Binomial[k+n/2-2,k]*t^k*Hypergeometric2F1[k,1-n/2,k+n/2,t^2]*GegenbauerC[k,(n-2)/2,x],{k,0,Infinity}],{S[n],Der[x],Der[t]},Assumptions->Abs[x]<1&&Abs[t]<1&&n≥3&&Element[n,Integers]&&Element[k,Integers]]

The limit of  GegenbauerC[k, (n - 2)/2, x] for k -> Infinity cannot be computed by Mathematica, that's why the above evaluation does not work. Assuming that the sum has a natural upper bound, we find some annihilating relations for the sum on the left-hand side.

Timing[{lhs,inhom}=Annihilator[Sum[Binomial[k+n-2,k]/Binomial[k+n/2-2,k]*t^k*Hypergeometric2F1[k,1-n/2,k+n/2,t^2]*GegenbauerC[k,(n-2)/2,x],{k,0,Infinity}],{S[n],Der[x],Der[t]},Inhomogeneous->True];]

ReleaseHold[inhom/.HoldPattern[Limit[__]]->0]//FullSimplify

lhs

The annihilating ideal for the right-hand side can be obtained immediately:

rhs=Annihilator[((1-t^2)/(1-2*t*x+t^2))^(n-1),{S[n],Der[x],Der[t]}]

We see that they do not agree, but that lhs is contained in rhs.

OreReduce[lhs,rhs]

We have to compare two initial values. For any n, in particular for n=3 and n=4, we have:

Binomial[k+n-2,k]/Binomial[k+n/2-2,k]*t^k*Hypergeometric2F1[k,1-n/2,k+n/2,t^2]*GegenbauerC[k,(n-2)/2,x]/.t->0/.x->0//FullSimplify

%/.(0^k->1)/.k->0

((1-t^2)/(1-2*t*x+t^2))^(n-1)/.t->0/.x->0

Theorem 2

For real x and t such that Abs[x]<1 and Abs[t]<1, and integer n>=3:

TraditionalForm[HoldForm[Sum[Gamma[(k+n)/2]*Gamma[(k+1)/2]/Gamma[n/2+k-1]* t^k *Hypergeometric2F1[k,k+n-1,k+n/2,(1-Sqrt[1-t^2])/2]*GegenbauerC[k,(n-2)/2,x],{k,0,Infinity}]==(n-2)*Gamma[(n+1)/2]/2/Gamma[n/2+1]* Sqrt[1-t^2] *Hypergeometric2F1[n,1,n/2+1,(x*t+1)/2]]]

The computation is similar to that of Theorem 1, but a little more complex.

Timing[{lhs,inhom}=Annihilator[Sum[Gamma[(k+n)/2]*Gamma[(k+1)/2]/Gamma[n/2+k-1]* t^k *Hypergeometric2F1[k,k+n-1,k+n/2,(1-Sqrt[1-t^2])/2] *GegenbauerC[k,(n-2)/2,x],{k,0,Infinity}],{S[n],Der[t],Der[x]},Inhomogeneous->True];]

This time we get some inhomogeneous contribution:

inhom=ReleaseHold[inhom/.HoldPattern[Limit[__]]->0]//FunctionExpand//FullSimplify

lhs=OreGroebnerBasis[Flatten[MapThread[Function[{a,b},(#**b)&/@a],{Annihilator[#,{S[n],Der[t],Der[x]}]&/@inhom,lhs}]]]

rhs=Annihilator[(n-2)*Gamma[(n+1)/2]/2/Gamma[n/2+1]* Sqrt[1-t^2]*Hypergeometric2F1[n,1,n/2+1,(x*t+1)/2],{S[n],Der[t],Der[x]}]

This time we have that the annihilating ideals that we computed for the left-hand and the right-hand side are identical.

GBEqual[lhs,rhs]

uts=UnderTheStaircase[rhs]

ApplyOreOperator[uts,Gamma[(k+n)/2]*Gamma[(k+1)/2]/Gamma[n/2+k-1]* t^k *Hypergeometric2F1[k,k+n-1,k+n/2,(1-Sqrt[1-t^2])/2] *GegenbauerC[k,(n-2)/2,x]]/.t->0/.x->0/.(0^k)->1/.k->0/.n->3

ApplyOreOperator[uts,(n-2)*Gamma[(n+1)/2]/2/Gamma[n/2+1]* Sqrt[1-t^2] *Hypergeometric2F1[n,1,n/2+1,(x*t+1)/2]]/.t->0/.x->0/.n->3

Conversion of triangular Jacobi indices

This problem has been sent by Joachim Schöberl (see the document converttrig.pdf).

TableForm[TraditionalForm[HoldForm[{ψ[i,j]==JacobiP[i,a,0,y/(1-x)]*(1-x)^i*JacobiP[j,c+2i,0,x],ϕ[i,j]==JacobiP[i,a+1,0,y/(1-x)]*(1-x)^i*JacobiP[j,c+1+2i,0,x]}]]]

Task: Find relations between ψ(i, j) and ϕ(i, j)! Indeed, we need c=a+1.

ann=Annihilator[JacobiP[i,a,0,y/(1-x)]*(1-x)^i*JacobiP[j,c+2 i,0,x],{S[a],S[c],S[i],S[j]}];

First, the general case:

FindRelation[ann,Eliminate->{x,y},Pattern->{a_,a_,_,_}]

Now, the special case (c=a+1):

anna=Annihilator[JacobiP[i,a,0,y/(1-x)]*(1-x)^i*JacobiP[j,a+1+2 i,0,x],{S[a],S[i],S[j]}];

FindRelation[anna,Eliminate->{x,y}]

Factor[%]

ApplyOreOperator[First[%],f[a,i,j]]==0/.f[a,ii_,jj_]->psi[ii,jj]/.f[a+1,ii_,jj_]->phi[ii,jj]

FEM application: Basis functions in 2D

Joachim Schöberl (RWTH Aachen) defines the basis functions as follows:

TraditionalForm[HoldForm[phi[i_,j_,x_,y_]==LegendreP[i,2y/(1-x)-1] (1-x)^i JacobiP[j,2i+1,0,2x-1]]]

Task: For the above function, find a representation of the derivatives w.r.t. x and y in terms of shifts of this function in i and j, i.e., find a relation of the form ((∂phi(i,j,x,y))/(∂x))==(∑_(m=0)^M)(∑_(n=0)^N)c(m,n) phi(i+m,j+n,x,y). The coefficients have to be free of x and y.

ann=Annihilator[LegendreP[i,2y/(1-x)-1] (1-x)^i JacobiP[j,2i+1,0,2x-1],{S[i],S[j],Der[x],Der[y]}];

Timing[Factor[FindRelation[ann,Eliminate->{x,y},Pattern->{_,_,0|1,0}]]]

Timing[Factor[FindRelation[ann,Eliminate->{x,y},Pattern->{_,_,0,0|1}]]]

We see that it is (in this case) more efficient to use only these operators that are really involved:

Timing[Factor[FindRelation[Annihilator[LegendreP[i,2y/(1-x)-1] (1-x)^i JacobiP[j,2i+1,0,2x-1],{S[i],S[j],Der[x]}],Eliminate->{x,y},Pattern->{_,_,0|1}]]]

Timing[Factor[FindRelation[Annihilator[LegendreP[i,2y/(1-x)-1] (1-x)^i JacobiP[j,2i+1,0,2x-1],{S[i],S[j],Der[y]}],Eliminate->{x,y},Pattern->{_,_,0|1}]]]

Integral over 3 Gegenbauer polynomials

The following identity can be found as formula (6.8.10) in Adrews/Askey/Roy, p. 321. 
For λ>-1/2 and λ ≠ 0, for  l+m+n  even and the sum of any two of l,m,n is not less than the third. The integral is zero in all other cases.

TraditionalForm[HoldForm[Integrate[GegenbauerC[l,λ,x]*GegenbauerC[m,λ,x]*GegenbauerC[n,λ,x]*(1-x^2)^(λ-1/2),{x,-1,1}]==Pochhammer[λ,(m+n-l)/2]*Pochhammer[λ,(n+l-m)/2]*Pochhammer[λ,(l+m-n)/2]*2^(1-2*λ)*Pi*Gamma[(l+m+n)/2+2*λ]/
(((m+n-l)/2)!*((l+n-m)/2)!*((l+m-n)/2)!*Pochhammer[λ,(l+m+n)/2]*Gamma[λ]^2*((l+m+n)/2+λ))]]

This example is quite involved. First, if we try to apply Annihilator to the integral, some problems concerning the automatic simplification of the delta parts arise. In order to assist in this part by hand, we use the option Inhomogeneous. Still, the computation takes a while:

integrand=GegenbauerC[l,λ,x]*GegenbauerC[m,λ,x]*GegenbauerC[n,λ,x]*(1-x^2)^(λ-1/2);
expr=Pochhammer[λ,(m+n-l)/2]*Pochhammer[λ,(n+l-m)/2]*Pochhammer[λ,(l+m-n)/2]*2^(1-2*λ)*Pi*Gamma[(l+m+n)/2+2*λ]/
(((m+n-l)/2)!*((l+n-m)/2)!*((l+m-n)/2)!*Pochhammer[λ,(l+m+n)/2]*Gamma[λ]^2*((l+m+n)/2+λ));

(*With $NormalizeMethod="gcd",Mathematica crashes*)nm=$NormalizeMethod;$NormalizeMethod="ftl";

Timing[{lhs,inh}=Annihilator[Integrate[integrand,{x,-1,1}],{S[l],S[m],S[n]},Assumptions->λ>1/2,Inhomogeneous->True];]

$NormalizeMethod=nm;

ReleaseHold[inh/.HoldPattern[Limit[a_,b_,_]]:>(Expand[a]/.b)]/.(0^a_):>FullSimplify[0^a,λ>1/2]

So the delta parts vanish and lhs is an annihilating ideal for the integral. The right-hand side gives rise to exactly the same annihilating ideal.

Factor[rhs=Annihilator[expr,{S[l],S[m],S[n]}]]

GBEqual[lhs,rhs]

The creative telescoping relations can much faster be found by our ansatz that guesses the denominators in advance:

Timing[Factor[First[FindCreativeTelescoping[integrand,Der[x],{S[l],S[m],S[n]}]]]]

We see that in the leading coefficients many singularities may be hidden. To find them, the command AnnihilatorSingularities can be used.

sing=AnnihilatorSingularities[rhs,{0,0,0},Assumptions->λ>1/2&&Element[(l+m+n)/2,Integers]&&l+m≥n&&l+n≥m&&m+n≥l]

For the 3 single singular points, the identity is readily verified:

FullSimplify[Integrate[integrand/.First[#],{x,-1,1},Assumptions->Last[#]]==expr/.First[#],Last[#]]&/@Drop[sing,3]

The first element in sing can be ignored since the condition is never fulfilled. The next two entries describe lines of points where none of the recurrences can be applied. They cannot be checked by Mathematica, so we recursively apply HolonomicFunctions. Both conditions lead to the same integrand:

integrand/.(First/@Take[sing,{2,3}])

For both sides, we get in this special instance the same recurrence, which this time does not have any singularities.

Annihilator[Integrate[GegenbauerC[m,λ,x]^2*(1-x^2)^(λ-1/2),{x,-1,1}],S[m],Assumptions->λ>1/2]

Factor[%]

expr/.sing[[2,1]]

Annihilator[%,S[m]]

Hence one initial value suffices (actually, the case l=m=n=0 has already been checked).

FullSimplify[Integrate[integrand/.sing[[2,1]]/.m->0,{x,-1,1},Assumptions->λ>1/2]==expr/.sing[[2,1]]/.m->0]

Krawtchouk and Meixner polynomials

We find that all three definitions given on MathWorld are equivalent:

TraditionalForm[HoldForm[Superscript[Subscript[k,n],"(p)"][x,nn]==Sum[(-1)^(n-v)*Binomial[nn-x,n-v]*Binomial[x,v]*p^(n-v)*(1-p)^v,{v,0,n}]==(-1)^n*Binomial[nn,n]*p^n*Hypergeometric2F1[-n,-x,-nn,1/p]==(-1)^n*p^n/n!*Gamma[nn-x+1]/Gamma[nn-x-n+1]*Hypergeometric2F1[-n,-x,nn-x-n+1,(p-1)/p]]]

kra1=Annihilator[Sum[(-1)^(n-v)*Binomial[nn-x,n-v]*Binomial[x,v]*p^(n-v)*(1-p)^v,{v,0,n}],{Der[p],S[nn],S[x],S[n]}]

kra2=Annihilator[(-1)^n*Binomial[nn,n]*p^n*Hypergeometric2F1[-n,-x,-nn,1/p],{Der[p],S[nn],S[x],S[n]}]

kra3=Annihilator[(-1)^n*p^n/n!*Gamma[nn-x+1]/Gamma[nn-x-n+1]*Hypergeometric2F1[-n,-x,nn-x-n+1,(p-1)/p],{Der[p],S[nn],S[x],S[n]}]

GBEqual[kra1,kra2]

GBEqual[kra1,kra3]

UnderTheStaircase[kra1]

Table[{Sum[(-1)^(n-v)*Binomial[nn-x,n-v]*Binomial[x,v]*p^(n-v)*(1-p)^v,{v,0,n}],(-1)^n*Binomial[nn,n]*p^n*Hypergeometric2F1[-n,-x,-nn,1/p],(-1)^n*p^n/n!*Gamma[nn-x+1]/Gamma[nn-x-n+1]*Hypergeometric2F1[-n,-x,nn-x-n+1,(p-1)/p]},{n,0,1}]//FullSimplify

Thomas Stoll was interested in a differential equation. We easily can derive it (in many different ways):

Timing[OreGroebnerBasis[kra1,OreAlgebra[S[nn],S[x],S[n],Der[p]],MonomialOrder->EliminationOrder[3]]]

Timing[FGLM[kra1,OreAlgebra[S[nn],S[x],S[n],Der[p]],MonomialOrder->EliminationOrder[3]]]

Annihilator[(-1)^n*Binomial[nn,n]*p^n*Hypergeometric2F1[-n,-x,-nn,1/p],Der[p]]

The same for the Meixner polynomials of the first kind  (on MathWorld there are two definitions, but they deliver the same differential equation.

Annihilator[Hypergeometric2F1[-n,-x,b,1-1/c],{S[n],S[x],S[b],Der[c]}]

Annihilator[Pochhammer[b,n]*Hypergeometric2F1[-n,-x,b,1-1/c],{S[n],S[x],S[b],Der[c]}]

(∫_0^π)(LegendreP[2 n, cos(t)])dt

In Abramowitz/Stegun  (22.13.6,  S. 785)  we find:

HoldForm[Integrate[LegendreP[2n,Cos[t]],{t,0,Pi}]==Pi/16^n*Binomial[2n,n]^2]//TraditionalForm

This problem can be solved by using non-standard Ore algebras which have Sin[t] and Cos[t] in the coefficients. The only algorithm that works  properly in this setting is Takayama's algorithm. The annihilating ideal has to be generated by hand.

ann=NormalizeCoefficients/@OrePolynomialSubstitute[Annihilator[LegendreP[2n,x],{Der[x],S[n]}],{x->cos[x],Der[x]->-Der[x]/sin[x]}]/.cos->Cos/.sin->Sin

AppendTo[ann,ToOrePolynomial[Sin[x]^2+Cos[x]^2-1,OreAlgebra[ann]]];

ApplyOreOperator[ann,LegendreP[2n,Cos[x]]]/.n->8//FullSimplify

tak=Takayama[ann,{Sin[x],Cos[x]}]

rhs=Annihilator[Pi/16^n*Binomial[2n,n]^2,S[n]]

OreReduce[tak,rhs]

But a simple substitution (u=cos(x)) brings the problem to a perfectly ∂-finite integrand:

TraditionalForm[HoldForm[Integrate[LegendreP[2n,Cos[t]],{t,0,Pi}]==Integrate[LegendreP[2n,u]/Sqrt[1-u^2],{u,-1,1}]]]

Annihilator[Integrate[LegendreP[2n,u]/Sqrt[1-u^2],{u,-1,1}],S[n]]

Integrate[1/Sqrt[1-u^2],{u,-1,1}]

The initial value certainly agrees with the one for the right-hand side.

(∫_0^(2 π))sin(m x) sin(n x)dx==π (Subscript[δ, m,n])

TraditionalForm[HoldForm[Integrate[Sin[m*x]*Sin[n*x],{x,0,2Pi}]==Pi*KroneckerDelta[m,n]]]

for integers m>0 and n>0.

First we look what Mathematica does with this integral:

FullSimplify[Integrate[Sin[m*x]*Sin[n*x],{x,0,2Pi}],Assumptions->Element[{m,n},Integers]]

FullSimplify[Integrate[Sin[n*x]*Sin[n*x],{x,0,2Pi}],Assumptions->Element[n,Integers]]

Table[Integrate[Sin[n x]^2,{x,0,2Pi}],{n,0,3}]

Compute an annihilating ideal for the integrand:

ann=Annihilator[Sin[m*x]*Sin[n*x],{S[m],S[n],Der[x]}]

AnnihilatorDimension[ann]

In order to perform the integration, we need one more relation!

Annihilator[Sin[m*x]*Sin[n*x],{S[m],S[n]}]

gb=OreGroebnerBasis[%,OreAlgebra[Cos[x],S[m],S[n]],MonomialOrder->EliminationOrder[1]]

The first operator is free of x and Der[x], hence it trivially annihilates the whole integral. But we can do better by including this operator to the annihilating ideal from above and perform Chyzak's algorithm.

ApplyOreOperator[First[gb],Sin[m*x]*Sin[n*x]]//FullSimplify

ann=OreGroebnerBasis[Append[ann,First[gb]],OreAlgebra[Der[x],S[m],S[n]],MonomialOrder->DegreeLexicographic]

ct=CreativeTelescoping[ann,Der[x],{S[m],S[n]}]

Check the delta part (it vanishes):

With[{delta=ApplyOreOperator[ct[[2,1]],Sin[m*x]*Sin[n*x]]},FullSimplify[(delta/.x->2 Pi)-(delta/.x->0),Assumptions->Element[{m,n},Integers]]]

Hence the first part of ct is an annihilating operator for the integral. In recurrence notation it looks as follows.

ApplyOreOperator[ct[[1,1]],f[m,n]]/.m->m-1

f[m,n]==Factor[f[m,n]/.First[Solve[%==0,f[m,n]]]]

We investigate the leading coefficient:

Factor[LeadingCoefficient[ct[[1,1]]]/.m->m-1]

Since the leading monomial of the recurrence is S[m], we have to give initial values for m=1 and arbitrary n. The leading coefficient contains (m-1) as a factor, which is an additional reason to include initial values for m=1.

Similarly for the diagonal m=n:

FullSimplify[Integrate[Sin[n x]^2,{x,0,2Pi}],Assumptions->Element[n,Integers]]

And for m=1, n>1 the integral is 0 (symmetries around x=Pi/2 and x=3Pi/2).

The initial conditions match with the expression on the right hand side, and also the recurrence is fine:

ApplyOreOperator[ct[[1,1]],KroneckerDelta[m,n]]//FullSimplify

Identity involving binomial sums

It was Gabor Hedegüs who asked for a proof of the following formula.

binom0[n_Integer,k_Integer]:=If[n≥k,Binomial[n,k],0]

Define

TraditionalForm[HoldForm[h[n,m,t]==binom0[n+1+t,n+1]+Sum[(-1)^j*binom0[n,j]*binom0[n+t-2j+1,n+1],{j,1,n}]+Sum[(-1)^(j+1)*binom0[n,j]*binom0[n+1+t-m-2j,n+1],{j,0,n}]]]

and prove that

TraditionalForm[HoldForm[h[n,m,t]-h[n,m,t-1]==Sum[Binomial[n,t-i],{i,0,m-1}]]]

First we convince ourselves that the function binom0[n,k] satisfies the same recurrences as Binomial[n,k] (this can be seen by looking at the leading coefficients of the recurrences):

Annihilator[Binomial[n,k],{S[k],S[n]}]

Union[Flatten[Table[Evaluate[ApplyOreOperator[%,binom0[n,k]]],{k,-10,10},{n,-10,10}]]]

Now we compute an annihilating ideal for the expression h[n,m,t]:

annh=Annihilator[Binomial[n+1+t,n+1]+Sum[(-1)^j*Binomial[n,j]*Binomial[n+t-2j+1,n+1],{j,1,n}]+Sum[(-1)^(j+1)*Binomial[n,j]*Binomial[n+1+t-m-2j,n+1],{j,0,n}],{S[m],S[n],S[t]}];

It's quite big, that's why we do not display it here.

ByteCount[annh]

Support[annh]

Next we use closure properties to compute an annihilating ideal for the left-hand side:

annlhs=DFiniteOreAction[DFiniteSubstitute[annh,{t->t-1}],S[t]-1];

ByteCount[annlhs]

uts=UnderTheStaircase[annlhs]

Now the right-hand side. This is easier:

annrhs=Annihilator[Sum[Binomial[n,t-i],{i,0,m-1}],{S[m],S[n],S[t]}]

UnderTheStaircase[annrhs]

The two annihilating ideals for left- and right-hand side are not identical, but one is contained in the other:

OreReduce[annlhs,annrhs]

This means that we have to compare 4 initial values, according to the elements under the staircase of the smaller ideal (that one for the left-hand side):

h[n_Integer,m_Integer,t_Integer]:=binom0[n+1+t,n+1]+Sum[(-1)^j*binom0[n,j]*binom0[n+t-2j+1,n+1],{j,1,n}]+Sum[(-1)^(j+1)*binom0[n,j]*binom0[n+1+t-m-2j,n+1],{j,0,n}]

ApplyOreOperator[uts,h[n,m,t]-h[n,m,t-1]]/.{m->1,n->1,t->0}

Clear[rhs];rhs[n_Integer,m_Integer,t_Integer]:=Sum[Binomial[n,t-i],{i,0,m-1}]

ApplyOreOperator[uts,rhs[n,m,t]]/.{m->1,n->1,t->0}

Clear[rhs];

Since the initial values match, the identity is proven for all integers m>0, n>0, and t>=0.

To be completely rigorous, we should have a look at the leading coefficients of the generators of annlhs and annrhs (appropriately shifted). If they have common zeros in the area of N^3 (where we claim the identity to hold), then we would have a problem.

GroebnerBasis[(LeadingCoefficient[#]/.Thread[{m,n,t}->({m,n,t}-LeadingExponent[#])])&/@annlhs,MonomialOrder->Lexicographic]

The only common zero is (m=1,n=0,t=0) for which we do not claim the identity to hold.

GroebnerBasis[(LeadingCoefficient[#]/.Thread[{m,n,t}->({m,n,t}-LeadingExponent[#])])&/@annrhs,MonomialOrder->Lexicographic]

Stirling sum from Gabor Hegedüs

Find an evaluation of the expression

TraditionalForm[HoldForm[(-1)^(d-i)*Sum[(-1)^j*Binomial[i+1,j]*Product[n-l,{l,0,j-1}]*((-1)^(n-j+m)*StirlingS1[n-j,m]),{j,0,n-m}]+(-1)^(n-i-1)*Product[n-j,{j,0,i}]*((-1)^(n-i-1+m)*StirlingS1[n-i-1,m])]]

We deal with the sum only.

(-1)^(d-i)*Sum[(-1)^j*Binomial[i+1,j]*Pochhammer[n-j+1,j]*((-1)^(n-j+m)*StirlingS1[n-j,m]),{j,0,n-m}]

annSum=Annihilator[Sum[(-1)^(d-i)*(-1)^j*Binomial[i+1,j]*Pochhammer[n-j+1,j]*((-1)^(n-j+m)*StirlingS1[n-j,m]),{j,0,n-m}],{S[i],S[m],S[n]},Assumptions->Element[m,Integers]&&m>0]

The recurrences seem to be correct:

mySum[i_Integer,m_Integer,n_Integer]:=(-1)^(d-i)*Sum[(-1)^j*Binomial[i+1,j]*Product[n-l,{l,0,j-1}]*((-1)^(n-j+m)*StirlingS1[n-j,m]),{j,0,n-m}]

test=ApplyOreOperator[annSum,mySum[i,m,n]];

Union[Flatten[Table[test,{i,4},{m,4},{n,4},{d,2}]]]

(∫_0^1)(x^a) log(log((1/x)))dx==-((ℽ+log(a+1))/(a+1))

This integral was communicated to us by Victor Moll.

TraditionalForm[HoldForm[Integrate[x^a Log[Log[1/x]],{x,0,1}]==-(EulerGamma+Log[a+1])/(a+1)]]

Perform the substitution x=Exp[-u]:

TraditionalForm[HoldForm[Integrate[x^a Log[Log[1/x]],{x,0,1}]==Integrate[Exp[-u(a+1)]*Log[u],{u,0,Infinity}]]]

lhs=Annihilator[Integrate[Exp[-u(a+1)]*Log[u],{u,0,Infinity}],Der[a],Assumptions->a>-1]

Initial values:

Integrate[Log[Log[1/x]],{x,0,1}]

Integrate[D[x^a Log[Log[1/x]],a]/.a->0,{x,0,1}]

Find the right hand side:

DSolve[{ApplyOreOperator[First[lhs],f[a]]==0,f[0]==-EulerGamma,f'[0]==EulerGamma-1},f[a],a]

If the right hand is given, we can prove the identity by computing its annihilator (and comparing initial values):

Annihilator[-(EulerGamma+Log[a+1])/(a+1),Der[a]]

Experiment: using Log[x] in the coefficients of the algebra

We do not perform substitution, but allow Log[x] in the coefficients of our Ore polynomials (but be careful with such strange algebras).

ann1=ToOrePolynomial[{Der[a],(Log[x]^2 x^2)**Der[x]^2-(-Log[x]x-Log[x]^2x)**Der[x]},OreAlgebra[Der[x],Der[a]]]

FullSimplify[ApplyOreOperator[ann1,Log[Log[1/x]]],Assumptions->x>0]

ann2=ChangeOreAlgebra[Flatten[Annihilator[x^a,#]&/@{Der[a],Der[x]}],OreAlgebra[Der[x],Der[a]]]

ann=DFiniteTimes[ann1,ann2]

ann=ChangeOreAlgebra[ChangeOreAlgebra[ann,OreAlgebra[Der[x],Log[x],Der[a]]],OreAlgebra[Der[x],x,Log[x],Der[a]]]

By an easy computation we get a telescoping operator:

ct=ann[[2]]+((a+1)^2)**ann[[1]]

The principal part is

pp=OrePolynomialSubstitute[ct,{Der[x]->0}]

The delta part is

dp=ToOrePolynomial[Together[Normal[ct-pp]/Der[x]],OreAlgebra[Der[x],x,Log[x],Der[a]]]

delta=ApplyOreOperator[dp,x^a Log[Log[1/x]]]

Limit[delta,x->1]-Limit[delta,x->0,Assumptions->a>-1]

The delta part contributes an inhomogeneous part 1, hence we homogenize the differential equation by multiplying it by Der[a]

de=Der[a]**pp

Bilateral hypergeometric sum from Ramanujan's Lost Notebook

This identity appeared as Formula no. 10 in "Andrews/Berndt: Your Hit Parade: The Top Ten Most Fascinating Formulas in Ramanujan's Lost Notebook"

Define, for any real numbers θ and s and for any complex numbers a,b,c, and d, such that Re(a+b+c+d)>3,

TraditionalForm[HoldForm[Subscript[ϕ,s][θ]==Sum[E^((n+s)*I*θ)/(Gamma[a+s+n]*Gamma[b-s-n]*Gamma[c+s+n]*Gamma[d-s-n]),{n,-Infinity,Infinity}]]]

Theorem 1 states that

TraditionalForm[HoldForm[D[Subscript[ϕ,s][θ]/Subscript[ϕ,t][θ],θ]==I*Sin[Pi*(s-t)]*(2Sin[1/2*θ])^(a+b+c+d-4)*E^(I*(Pi-θ)*(a-b+c-d+2s+2t)/2)/(Pi*Subscript[ϕ,t][θ]^2*Gamma[a+b-1]*Gamma[b+c-1]*Gamma[c+d-1]*Gamma[d+a-1])]]

For convenience, we replace  Exp[I*theta]  by z. Further we multiply both sides by (((Subscript[ϕ, t])(θ))^2). Starting from the definition of ϕ, we are now going to compute an annihilating ideal for ((∂(Subscript[ϕ, s])(θ))/(∂θ)) (Subscript[ϕ, t])(θ)-((∂(Subscript[ϕ, t])(θ))/(∂θ)) (Subscript[ϕ, s])(θ).

phiSmnd[s_]:=z^(n+s)/Gamma[a+s+n]/Gamma[b-s-n]/Gamma[c+s+n]/Gamma[d-s-n]
dphiSmnd[s_]:=I*(n+s)*z^(n+s)/Gamma[a+s+n]/Gamma[b-s-n]/Gamma[c+s+n]/Gamma[d-s-n]

ops={S[a],S[b],S[c],S[d]};

Timing[ann=DFinitePlus[
DFiniteTimes[
First[CreativeTelescoping[dphiSmnd[s],S[n]-1,ops]],
First[CreativeTelescoping[phiSmnd[t],S[n]-1,ops]]],
DFiniteTimes[
First[CreativeTelescoping[dphiSmnd[t],S[n]-1,ops]],
First[CreativeTelescoping[phiSmnd[s],S[n]-1,ops]]]];]

ByteCount[ann]

Length[ann]

UnderTheStaircase[ann]

We now apply these operators to the right-hand side. In order to speed up the evaluation, we consider each summand of an operator in ann separately: we first apply the operator monomial, cancel away the Gammas, and then multiply by the relatively big coefficients. Then simplification basically boils down to rational function arithmetic.

myApply[OrePolynomial[data_List,alg_OreAlgebraObject,mon_]]:=
Module[{i,result=0,expr},
expr=I*Sin[Pi*(s-t)]*(2Sin[th/2])^(a+b+c+d-4)*Exp[I*(Pi-th)*(a-b+c-d+2s+2t)/2]/Pi/Gamma[a+b-1]/Gamma[b+c-1]/Gamma[c+d-1]/Gamma[d+a-1];
Do[
result+=(data[[i,1]]/.z->Cos[th]+I*Sin[th])*(FullSimplify[ApplyOreOperator[OrePolynomial[{{1,data[[i,2]]}},alg,mon],expr]/expr]/.Exp[I*x_]->Cos[x]+I*Sin[x]),
{i,Length[data]}];
Return[result];
];

FullSimplify[Together[myApply/@ann]]

Comparing initial values is not so easy here. We omit it.

Stoll/Zeng: differential equations related to ChebyshevU

The following examples can be found in:
Stoll/Zeng: "On the total weight of weighted matchings of segment graphs" (Electronic Journal of Combinatorics), 
Proposition 2.2

The task is to derive differential equations (in x) for U_n(x,a) where (according to Prop. 2.1) for n=2k+1 we have

TraditionalForm[HoldForm[U[n,x,a]=Sqrt[a]^k*x*ChebyshevU[k,(x^2-a-1)/(2Sqrt[a])]]]

and for n=2k we have

TraditionalForm[HoldForm[U[n,x,a]=Sqrt[a]^k*(ChebyshevU[k,(x^2-a-1)/(2Sqrt[a])]+Sqrt[a]*ChebyshevU[k-1,(x^2-a-1)/(2Sqrt[a])])]]

annUodd=Annihilator[Sqrt[a]^((n-1)/2)*x*ChebyshevU[(n-1)/2,(x^2-a-1)/(2Sqrt[a])],Der[x]]

Collect[#,x,Factor]&/@Reverse[OrePolynomialListCoefficients[First[annUodd]]]//TableForm

For even n, there are two options. Using the expression as it is given (standard closure properties like Plus and Times will be applied), delivers a DE of order 3:

Annihilator[Sqrt[a]^(n/2)*(ChebyshevU[n/2,(x^2-a-1)/(2Sqrt[a])]+Sqrt[a]*ChebyshevU[n/2-1,(x^2-a-1)/(2Sqrt[a])]),Der[x]]

Better is to slightly rewrite the expression in order to use the closure property "application of Ore operator". We now obtain the second order differential equation that is given in Prop. 2.2  (note that we have to include S[n], since it is used in ApplyOreOperator):

annUeven=Annihilator[Sqrt[a]^(n/2)*ApplyOreOperator[S[n]^2+Sqrt[a],ChebyshevU[n/2-1,(x^2-a-1)/(2Sqrt[a])]],{S[n],Der[x]}]

Collect[#,x,Factor]&/@Reverse[OrePolynomialListCoefficients[First[annUeven]]]//TableForm

Some remarks on the recurrences:

OreGroebnerBasis[annUeven,OreAlgebra[Der[x],S[n]],MonomialOrder->EliminationOrder[1]]

Annihilator[Sqrt[a]^((n-1)/2)*x*ChebyshevU[(n-1)/2,(x^2-a-1)/(2Sqrt[a])],S[n]]

Hence we have this recurrence for both the even and the odd case. The same recurrence can be obtained by combining the two given in formula (2):

recEven=ToOrePolynomial[S[n]^2-x*S[n]+a]

recOdd=ToOrePolynomial[S[n]^2-x*S[n]+1]

recEven**recOdd+(x*S[n])**(recEven+recOdd)

Note that recEven and recOdd commute:

recEven**recOdd-recOdd**recEven

Victor Moll: A Master Formula

This is called a "Master Formula" in Chapter 13 of "Victor Moll: Irresistible Integrals".

TraditionalForm[HoldForm[Integrate[(x^2/(x^4+2a x^2+1))^r*(x^2+1)/x^2/(x^s+1),{x,0,Infinity}]==Beta[r-1/2,1/2]/2^(r+1/2)/(a+1)^(r-1/2)]]

We can try to solve it right away, but the factor (x^s+1) in the denominator causes the integrand to be non-holonomic. The right-hand side, however, is holonomic.

lhs=Annihilator[Integrate[(x^2/(x^4+2a x^2+1))^r*(x^2+1)/x^2/(x^s+1),{x,0,Infinity}],{Der[a],S[r],S[s],Der[s]}]

rhs=Annihilator[Beta[r-1/2,1/2]/2^(r+1/2)/(a+1)^(r-1/2),{Der[a],S[r]}]

The relation that was found for the integral is contained in rhs. To conclude the proof, infinitely many initial values have to be compared (namely the cases a=0 and r>=1m as well as all higher derivatives with respect to a and r=1 (and of course, everything for symbolic s). This seems to be a cumbersome job, and we omit it.

OreReduce[lhs,rhs]

For the special case s=0, we obtain the same annihilating ideal as for the right-hand side

lhs1=Annihilator[Integrate[(x^2/(x^4+2a x^2+1))^r*(x^2+1)/x^2/(x^0+1),{x,0,Infinity}],{Der[a],S[r]},Assumptions->r≥1]

So this case can be made rigorous by comparing one initial value, which can be done by Mathematica:

Integrate[(x^2/(x^4+2a x^2+1))^r*(x^2+1)/x^2/2/.r->1/.a->0,{x,0,Infinity}]==Beta[r-1/2,1/2]/2^(r+1/2)/(a+1)^(r-1/2)/.r->1/.a->0

For other specific values of s, we still get a zero-dimensional ideal, but for s>0 only a subideal of rhs. Interestingly, its vector space dimension seems to be Floor[s/4]+2 for even s>=2 and Floor[s/2]+3 for odd s.

Table[Length[UnderTheStaircase[Annihilator[Integrate[(x^2/(x^4+2a x^2+1))^r*(x^2+1)/x^2/(x^s+1),{x,0,Infinity}],{Der[a],S[r]},Assumptions->r≥1]]],{s,0,15}]

Table[Which[s===0,1,EvenQ[s],Floor[s/4]+2,OddQ[s],Floor[s/2]+3],{s,0,15}]

Larry Glasser: Integral representations for exceptional univariate Lommel functions

This problem has been sent to us by Larry Glasser (see his paper of the same title, submitted to J.Phys.A:Math.Gen.).

We are considering the Lommel function (Subscript[S, m,n])(z) which is a particular solution to the inhomogeneous Bessel equation

BesselEq=z^2*w''[z]+z*w'[z]+(z^2-n^2)*w[z]==z^(m+1)

In the following computations we do not include treatment of initial values, which would be necessary in order to have rigorous proofs.

Formula (14) gives the following integral representation of this function:

TraditionalForm[HoldForm[Integrate[ArcSinh[x]/Sqrt[x^2+1]*E^(-z*x),{x,0,Infinity}]==Subscript[S,"-1,0"][z]]]

We compute an inhomogeneous differential equation for this integral and compare it to the defining Bessel equation.

Annihilator[Integrate[ArcSinh[x]/Sqrt[x^2+1]*E^(-z*x),{x,0,Infinity}],Der[z],
Assumptions->Re[z]>0,Inhomogeneous->True,Head->S]//ReleaseHold

FullSimplify[First[%]==(BesselEq/.{w->S,m->-1,n->0}),Re[z]>0]

Formula (15) states

TraditionalForm[HoldForm[Integrate[ArcSinh[x]^2*E^(-z*x),{x,0,Infinity}]==2/z*Subscript[S,"-1,0"][z]]]

We compute an inhomogeneous differential equation for this integral (we multiply the equation by z/2) and compare it to the defining Bessel equation.

Annihilator[Integrate[z/2*ArcSinh[x]^2*E^(-z*x),{x,0,Infinity}],Der[z],
Assumptions->Re[z]>0,Inhomogeneous->True,Head->S]//ReleaseHold

FullSimplify[First[%]==(BesselEq/.{w->S,m->-1,n->0}),Re[z]>0]

Formula (17):

TraditionalForm[HoldForm[(1/2)*Integrate[x*ArcSinh[x]/Sqrt[x^2+1]*E^(-z*x),{x,0,Infinity}]==Subscript[S,"-2,1"][z]]]

We compute an inhomogeneous differential equation for this integral and compare it to the defining Bessel equation.

Annihilator[Integrate[x*ArcSinh[x]/Sqrt[x^2+1]*E^(-z*x)/2,{x,0,Infinity}],Der[z],
Assumptions->Re[z]>0,Inhomogeneous->True,Head->S]//ReleaseHold

FullSimplify[First[%]==(BesselEq/.{w->S,m->-2,n->1})]

Formula (19):

TraditionalForm[HoldForm[Integrate[x^n*BesselJ[n,x]*Subscript[S,"-1,0"][x],{x,0,Infinity}]==2^n*Gamma[n+1/2]/Sqrt[Pi]*Integrate[ArcSinh[x]/(1+x^2)^(n+1),{x,0,Infinity}]]]

A recurrence for the right-hand side is readily derived:

Annihilator[2^n*Gamma[n+1/2]/Sqrt[Pi]*Integrate[ArcSinh[x]/(1+x^2)^(n+1),{x,0,Infinity}],S[n],Assumptions->n≥0]

For the left-hand side we have to do the calculation in several steps. First we write down an annihilating ideal for (Subscript[S, -1,0])(x) by homogenizing the above Bessel equation. Then use closure properties to get an annihilating ideal for the whole integrand.

annLommel=ToOrePolynomial[{S[n]-1,Der[x]**(x^2*Der[x]^2+x*Der[x]+x^2)},OreAlgebra[Der[x],S[n]]]

integrand=DFiniteTimes[Annihilator[x^n*BesselJ[n,x],{Der[x],S[n]}],annLommel]

Finally, we peform creative telescoping. Since the Lommel function is not given explicitely here, we omit checking the delta part and concentrate only on the principal part. It agrees with the recurrence for the right-hand side (all this is a little bit sloppy; the "proof" works also only for integer values of n).

First[CreativeTelescoping[integrand,Der[x],S[n]]]

Integrated Jacobi polynomials for tetrahedral p-FEM

In the paper by Beuchler/Pillwein: "Sparse shape functions for tetrahedral p-FEM using integrated Jacobi polynomials" (Computing 80(4), 2007), we find some formulas that can be found with HolonomicFunctions automatically. In the following, the function (Subscript[p, n]^a)(x) corresponds to the Jacobi polynomial with b=0, and (Subscript[(Overscript[p, "⋀"]), n]^a)(x) denotes its integrated version. We compute an annihilating ideal for the latter function:

phat=Annihilator[(1+x)/n JacobiP[n-1,a-1,1,x],{Der[x],S[n],S[a]}]

Formula (12)

... expresses (Subscript[(Overscript[p, "⋀"]), n]^a)(x) in terms of n-shifts of  (Subscript[p, n]^a)(x). Here are two ways to get it.

Factor[FindRelation[phat,Support->{S[n],Der[x],Der[x]S[n],Der[x]S[n]^2},Eliminate->x]]

Factor[FindRelation[phat,Eliminate->x,Pattern->{_,_,0}]]

Collect[ApplyOreOperator[First[%],ph[n,x]]/.Derivative[0,1][ph][k_,z_]->p[k-1,x]/.n->n-1,(p|ph)[_,_],Factor]

Formula (13)

... gives a recurrence in n for (Subscript[(Overscript[p, "⋀"]), n]^a)(x).

Collect[#,ph[_],Factor]&/@Solve[(ApplyOreOperator[First[FindRelation[phat,Pattern->{0,_,0}]],ph[n]]/.n->n-1)==0,ph[n+1]][[1,1]]/.Rule->Equal

Formula (14)

... expresses (Subscript[(Overscript[p, "⋀"]), n]^a)(x) in terms of (Subscript[p, n]^a)(x).

Factor[FindRelation[phat,Support->{S[a] S[n],S[n] Der[x],Der[x]}]]

ApplyOreOperator[First[FindRelation[phat,Eliminate->x]],ph[n,a,x]]/.a->a-1/.
Derivative[0,0,1][ph][nn_,aa_,x]->p[nn,aa,x]

Solve[%==0,ph[n,a,x]][[1,1]]/.Rule->Equal

Comparison with the paper reveals that we found a typo in this formula.

Triple sum from Bernd Sturmfels

In an e-mail dated from 2010-11-23, Bernd wrote:
"I wonder whether you would be willing to do me a quick "summation favor" of the "A=B" kind.
It concerns the sum in Theorem 5.15. Can it be simplified to something nicer ??"

smnd=(-4)^j*Binomial[d-2-k,j]*Multinomial[2i,k-i,n-k-i]*CatalanNumber[i+d-2-k-j];

TraditionalForm[HoldForm[Sum[s,{k,0,d-2},{i,0,n},{j,0,d-2-k}]]/.s->smnd]

After a few tries, one finds that the inner double sum can be simplified to a hypergeometric term. But we were not able to get a closed form for the entire expression.

Here is the annihilating ideal for the inner double sum (natural boundaries):

ann2=First[CreativeTelescoping[First[CreativeTelescoping[smnd,S[j]-1,{S[i],S[k],S[d],S[n]}]],S[i]-1,{S[k],S[d],S[n]}]]

Test whether this is correct:

SetDelayed@@(Hold[MySum2[k_Integer,d_Integer,n_Integer],Sum[Sum[s,{j,0,d-2-k}],{i,0,n}]]/.s->smnd)

test=ApplyOreOperator[ann2,MySum2[k,d,n]];
Position[Table[Total[test],{k,0,5},{d,0,5},{n,0,5}],a_Integer/;a≠0]-1

It seems that ann2 is not the correct annihilating ideal for MySum2:

test=ApplyOreOperator[ann2,MySum2[k,d,n]];
Position[Table[Total[test],{d,0,5},{n,0,5},{k,0,d-2}],a_Integer/;a≠0]-1

Indeed for k=d-2 we get a different ideal:

smnd2=smnd/.k->d-2;
ann2a=First[CreativeTelescoping[First[CreativeTelescoping[smnd2,S[j]-1,{S[i],S[d],S[n]}]],S[i]-1,{S[d],S[n]}]]

OrePolynomialSubstitute[Most[ann2],{k->d-2}]

However, we compute the hypergeometric solution of ann2 (step by step using one equation after the other):

sol1=RSolve[ApplyOreOperator[ann2[[1]],f[n]]==0,f[n],n][[1,1,2]]

This is not a good one:

Table[sol1,{k,3},{n,3},{d,3}]

We figure it out by hand:

TableForm[Table[Factor[RSolve[(1+n) (d-k+n) f[n]+(-1+k-n) (d-2 k+n) f[1+n]==0/.{d->dd,k->kk},f[n],n][[1,1,2]]],{dd,0,5},{kk,0,4}]/.C[1]->1]

TableForm[Table[Factor[RSolve[(1+n) (d-k+n) f[n]+(-1+k-n) (d-2 k+n) f[1+n]==0/.{d->dd,k->kk},f[n],n][[1,1,2]]]/Pochhammer[n-kk+1,kk]/Pochhammer[n-2kk+dd,kk],{dd,0,5},{kk,0,4}]/.C[1]->1]

sol1=Pochhammer[n-k+1,k]*Pochhammer[n-2k+d,k];

GBEqual[Annihilator[sol1,S[n]],{ann2[[1]]}]

Using the second equation, we compute the factors of the solution that depend on d but not on n:

DFiniteTimes[{ChangeOreAlgebra[ann2[[2]],OreAlgebra[S[d]]]},Annihilator[1/sol1,S[d]]]

RSolve[ApplyOreOperator[First[%],f[d]]==0,f[d],d][[1,1,2]]

sol2=sol1*%/C[1]

Using the third equation, we compute the remaining factors of the solution that depend on k only:

DFiniteTimes[{ChangeOreAlgebra[ann2[[3]],OreAlgebra[S[k]]]},Annihilator[1/sol2,S[k]]]

RSolve[ApplyOreOperator[First[%],f[k]]==0,f[k],k][[1,1,2]]

sol3=sol2*%/C[1]

Check whether this solution satisfies the original set of recurrences:

GBEqual[ann2,Annihilator[sol3,First[OreAlgebra[ann2]]]]

Find the constant:

MySum2[2,5,4]

sol3/.Thread[{k,d,n}->{2,5,4}]

sol3=FullSimplify[-sol3/4]

It turns out that this solution is also valid for k=d-2:

Union[Flatten[Table[sol3-MySum2[k,d,n],{n,0,10},{d,0,10},{k,0,d-2}]]]

Here is our result:

TraditionalForm[HoldForm[Sum[(-1)^d*2^(2d-3)*Pochhammer[1/2-k,d-1]*Pochhammer[d+n-2k,k]*Pochhammer[n+1-k,k]/(d-1)!/(2k)!,{k,0,d-2}]]]

A final check indicates that it is indeed correct:

MySol[d_Integer,n_Integer]:=Sum[(-1)^d*2^(2d-3)*Pochhammer[1/2-k,d-1]*Pochhammer[d+n-2k,k]*Pochhammer[n+1-k,k]/(d-1)!/(2k)!,{k,0,d-2}]

SetDelayed@@(Hold[MySum[d_Integer,n_Integer],Sum[Sum[Sum[s,{j,0,d-2-k}],{i,0,n}],{k,0,d-2}]]/.s->smnd)

Union[Flatten[Table[MySum[d,n]-MySol[d,n],{d,0,10},{n,0,10}]]]

Examples taken from Frédéric Chzyak's work

The shortcuts refer to the following literature:
98: Chyzak's thesis: Fonctions holonomes en Calcul formel
98b: Gröbner Bases, Symbolic Summation and Symbolic Integration  (Proceedings of the Conference 33 Years of Gröbner Bases)
00: An extension of Zeilberger's fast algorithm to general holonomic functions  (Discrete Math)

[98, p.140] or [98b, Section 3.2]: Indefinite Integral

We get a better result (annihilator of lower vector space dimension) by a simpler approach:

ann=Annihilator[Integrate[(1+x t+t^2)^(-2),{t,y,Infinity}],{Der[x],Der[y]}]

We try Mathematica on that...

Integrate[(1+x t+t^2)^(-2),{t,y,Infinity}]

[98, p.184]: Integral with four Bessel functions

HoldForm[Integrate[x BesselJ[1,a x] BesselI[1,a x] BesselY[0,x] BesselK[0,x],{x,0,Infinity}]==-Log[1-a^4]/(2Pi a^2)]//TraditionalForm

Takayama's algorithm works perfect on this example (if we believe that the integral has natural boundaries).

ann=Annihilator[x BesselJ[1,a x] BesselI[1,a x] BesselY[0,x] BesselK[0,x],{Der[x],Der[a]}]

tak=Takayama[ann,{x}]

ApplyOreOperator[First[tak],-Log[1-a^4]/(2Pi a^2)]//FullSimplify

CreativeTelescoping also works, but takes a while (because the relatively big dimension and because Mathematica takes very long in simplifying the delta part).

Length[UnderTheStaircase[ann]]

Timing[lhs=Annihilator[Integrate[x BesselJ[1,a x] BesselI[1,a x] BesselY[0,x] BesselK[0,x],{x,0,Infinity}],Der[a],Assumptions->0<a<1,Inhomogeneous->True];]

First[lhs]

Timing[lhs=Annihilator[Integrate[x BesselJ[1,a x] BesselI[1,a x] BesselY[0,x] BesselK[0,x],{x,0,Infinity}],Der[a],Assumptions->0<a<1]]

Using the command FindCreativeTelescoping, the creative telescoping relation can be found in much shorter time.

Timing[FindCreativeTelescoping[x BesselJ[1,a x] BesselI[1,a x] BesselY[0,x] BesselK[0,x],Der[x],Der[a]]]

rhs=Annihilator[-Log[1-a^4]/(2Pi a^2),Der[a]]

Compare two initial values  (we choose the point a=1/2):

Integrate[x BesselJ[1,a x] BesselI[1,a x] BesselY[0,x] BesselK[0,x]/.a->1/2,{x,0,Infinity}]

-Log[1-a^4]/(2Pi a^2)/.a->1/2

Here we cheat a little bit (Mathematica cannot evaluate the integral symbolically):

NIntegrate[D[x BesselJ[1,a x] BesselI[1,a x] BesselY[0,x] BesselK[0,x],a]/.a->1/2,{x,0,Infinity}]

D[-Log[1-a^4]/(2Pi a^2),a]/.a->1/2

N[%]

[98b]: Contiguous relations for the Appell (Subscript[F, 4]) hypergeometric function

We consider the Appell (Subscript[F, 4]) function that is defined as

TraditionalForm[HoldForm[Subscript[F,4][a,b,c,d,x,y]==Sum[Pochhammer[a,m+n]*Pochhammer[b,m+n]/m!/n!/Pochhammer[c,m]/Pochhammer[d,n]*x^m*y^n,{m,0,Infinity},{n,0,Infinity}]]]

These annihilating operators for the Appell (Subscript[F, 4]) function are given in Chyzak98b, formula (0.3):

ops=ToOrePolynomial[{x**Der[x]**(x**Der[x]+c-1)-x**(x**Der[x]+y**Der[y]+a)**(x**Der[x]+y**Der[y]+b),y**Der[y]**(y**Der[y]+d-1)-y**(x**Der[x]+y**Der[y]+a)**(x**Der[x]+y**Der[y]+b)}]

They do not form a Gröbner basis as the following computation shows:

gb=OreGroebnerBasis[ops,MonomialOrder->DegreeLexicographic]

In the Saito-Sturmfels-Takayama book we find on page 39 the following annihilating operators for the Appell (Subscript[F, 4]) function.

ops2=ChangeOreAlgebra[ToOrePolynomial[{Der[x]**(Euler[x]+c-1)-(Euler[x]+Euler[y]+a)**(Euler[x]+Euler[y]+b),Der[y]**(Euler[y]+d-1)-(Euler[x]+Euler[y]+a)**(Euler[x]+Euler[y]+b)}],OreAlgebra[ops]]

However, they do not agree with Chyzak's definition. We therefore use Takayama's algorithm in order to derive the relations ourselves (and find that Chyzak's results are correct).

GBEqual[gb,OreGroebnerBasis[ops2,MonomialOrder->DegreeLexicographic]]

tak=Takayama[Annihilator[Pochhammer[a,m+n]*Pochhammer[b,m+n]/m!/n!/Pochhammer[c,m]/Pochhammer[d,n]*x^m*y^n,{S[m],S[n],Der[x],Der[y]}],{m,n}];

GBEqual[gb,tak]

A "step-up" operator, i.e., an operator that maps (Subscript[F, 4])(a,b,c,d,x,y) to (Subscript[F, 4])(a+1,b,c,d,x,y), using the Ore operators (Subscript[D, x]) and (Subscript[D, y]) only, can be read off directly from the definition of the function:

stepUp=ToOrePolynomial[(x**Der[x]+y**Der[y]+a)/a]

The task is now to compute its inverse, i.e., a "step-down" operator that maps (Subscript[F, 4])(a,b,c,d,x,y) to (Subscript[F, 4])(a-1,b,c,d,x,y). This can be done by ansatz; the shape of the ansatz is determined by the monomials under the stairs of the Gröbner basis.

uts=UnderTheStaircase[gb]

stepDown=Array[c,4].uts

normalForm=Together[OreReduce[stepDown**stepUp-1,gb]];

matrix=Append[Coefficient[#,Array[c,4]],#/.c[_]->0]&/@OrePolynomialListCoefficients[normalForm];

solution=Together[First[NullSpace[matrix]]]

Since we started by reducing the relation stepDown*stepUp=1, we got the "step-down" operator that maps (Subscript[F, 4])(a+1,b,c,d,x,y) to (Subscript[F, 4])(a,b,c,d,x,y). Therefore we replace a by a-1.

stepDown=OrePolynomialSubstitute[Together[stepDown/.Thread[Array[c,4]->Most[solution]]],{a->a-1}];

We want to have some empirical correctness of our result.

test=ApplyOreOperator[stepDown,Pochhammer[a,m+n]*Pochhammer[b,m+n]/m!/n!/Pochhammer[c,m]/Pochhammer[d,n]*x^m*y^n]-Pochhammer[a-1,m+n]*Pochhammer[b,m+n]/m!/n!/Pochhammer[c,m]/Pochhammer[d,n]*x^m*y^n;

With[{bound=3},Drop[#,-2]&/@Drop[CoefficientList[Together[Sum[test,{m,0,bound},{n,0,bound}]],{x,y}],-2]]

[00, p.128]: hypergeometric double sum

TraditionalForm[HoldForm[Sum[(-1)^(n+r+s)Binomial[n,r]Binomial[n,s]Binomial[n+r,r]Binomial[n+s,s]Binomial[2n-(r+s),n],{r,0,Infinity},{s,0,Infinity}]==Sum[Binomial[n,k]^4,{k,0,Infinity}]]]

First observe that the sums have finite support, so we can safely write:

TraditionalForm[HoldForm[Sum[(-1)^(n+r+s)Binomial[n,r]Binomial[n,s]Binomial[n+r,r]Binomial[n+s,s]Binomial[2n-(r+s),n],{r,0,n},{s,0,n}]==Sum[Binomial[n,k]^4,{k,0,n}]]]

Trying the brute-force strategy is not a good idea here (it takes very long and gives a huge recurrence of order 7). The reasons are the delta parts that are relatively big expressions and require time-consuming simplification, and secondly the delta parts of the outer sum which involve the inner sum and which are not simplified to 0. Giving some assumptions improves the situation, but is still not satisfactory.

Timing[UnderTheStaircase[Annihilator[Sum[(-1)^(n+r+s)Binomial[n,r]Binomial[n,s]Binomial[n+r,r]Binomial[n+s,s]Binomial[2n-(r+s),n],{r,0,n},{s,0,n}],S[n],Assumptions->Element[n,Integers]&&n≥0]]]

However, the sum has finite support and therefore natural boundaries, so we need not to care about the delta parts. Using CreativeTelescoping, a recurrence is quickly obtained. Takayama's algorithm also works and delivers the same result, but takes a little bit longer. The right-hand side gives rise to the same recurrence, too.

First[CreativeTelescoping[First[CreativeTelescoping[(-1)^(n+r+s)Binomial[n,r]Binomial[n,s]Binomial[n+r,r]Binomial[n+s,s]Binomial[2n-(r+s),n],S[r]-1,{S[s],S[n]}]],S[s]-1,S[n]]]

Takayama[Annihilator[(-1)^(n+r+s)Binomial[n,r]Binomial[n,s]Binomial[n+r,r]Binomial[n+s,s]Binomial[2n-(r+s),n],{S[r],S[s],S[n]}],{r,s}]

Annihilator[Sum[Binomial[n,k]^4,{k,0,n}],S[n]]

The leading coefficient has no positive integer zeros:

Factor[LeadingCoefficient[First[%]]]

Hence we have to compare two initial values.

Table[Sum[(-1)^(n+r+s)Binomial[n,r]Binomial[n,s]Binomial[n+r,r]Binomial[n+s,s]Binomial[2n-(r+s),n],{r,0,n},{s,0,n}]==Sum[Binomial[n,k]^4,{k,0,n}],{n,0,1}]

[98, p.153]: Strehl identity (using Takayama's algorithm)

TraditionalForm[HoldForm[Sum[Binomial[n,k]^2*Binomial[n+k,k]^2,{k,0,n}]==Sum[Binomial[n,k] Binomial[n+k,k] Sum[Binomial[k,j]^3,{j,0,k}],{k,0,n}]]]

lhs=Takayama[Annihilator[Binomial[n,k]^2*Binomial[n+k,k]^2,{S[k],S[n]}],{k}]

We first do the inner sum:

rhs1=Takayama[Annihilator[Binomial[k,j]^3,{S[j],S[k],S[n]}],{j}]

Unfortunately, these are not enough relations for continuing. Using the option Saturate->True gives the desired result:

rhs1=Takayama[Annihilator[Binomial[k,j]^3,{S[j],S[k],S[n]}],{j},Saturate->True]

With the closure property product, we get a description of the summand of the outer sum:

rhs2=DFiniteTimes[rhs1,Annihilator[Binomial[n,k]*Binomial[n+k,k],{S[k],S[n]}]]

Again Takayama delivers a recurrence for the right-hand side.

rhs=Takayama[rhs2,{k}]

GBEqual[lhs,rhs]

Doing both summations in one stroke does not deliver the minimal recurrence, but a multiple of it:

Takayama[Annihilator[Binomial[n,k] Binomial[n+k,k]Binomial[k,j]^3,{S[j],S[k],S[n]}],{j,k}]

OreReduce[%,lhs]

We have to compare two initial values:

Table[Sum[Binomial[n,k]^2*Binomial[n+k,k]^2,{k,0,n}]==Sum[Binomial[n,k] Binomial[n+k,k] Sum[Binomial[k,j]^3,{j,0,k}],{k,0,n}],{n,0,1}]

[98b, p.18]: Takayama für q-Summation

TraditionalForm[HoldForm[Sum[q^((i+j)^2+j^2)/QPochhammer[q,q,n-i-j]/QPochhammer[q,q,i]/QPochhammer[q,q,j],{j,0,n},{k,1,j}]]]

In order to find a recurrence for this double sum, we first compute an annihilating ideal for the summand (it is q-hypergeometric).

ann=Annihilator[q^((i+j)^2+j^2)/QPochhammer[q,q,n-i-j]/QPochhammer[q,q,i]/QPochhammer[q,q,j],{QS[qi,q^i],QS[qj,q^j],QS[qn,q^n]}]

Takayama's algorithm in the first try gives an order 4 recurrence.

Printlevel=1;ann1=Takayama[ann,{i,j}]

But in this example, it pays off to extend the module in which the Takayama elimination takes place. Using the option StartDegree we can enforce that elements are included up to degree 7 at least (it started with degree 6 before).

Printlevel=0;ann2=Takayama[ann,{i,j},StartDegree->7]

Gradshteyn / Ryzhik

2.148   ∫(1/(((1+(x^2)))^n))dx

(p. 77) Formula 3 is stated as follows:

TraditionalForm[HoldForm[Integrate[1/(1+x^2)^n,x]==1/(2n-2)*x/(1+x^2)^(n-1)+(2n-3)/(2n-2)*Integrate[1/(1+x^2)^(n-1),x]]]

A slight rewriting gives

TraditionalForm[HoldForm[2n*Integrate[1/(1+x^2)^(n+1),x]-(2n-1)*Integrate[1/(1+x^2)^n,x]==x/(1+x^2)^n]]

We homogenize this recurrence by multiplying an annihilating operator for the inhomogeneous part to the operator on the left-hand side.

First[Annihilator[x/(1+x^2)^n,S[n]]]**((2n)**S[n]-(2n-1))

The recurrence that we directly obtain for the indefinite integral matches exactly:

Annihilator[Integrate[1/(1+x^2)^n,x],S[n]]

Formula 4 just unwinds the above inhomogeneous recurrence. However, it can be proven separately, too.

TraditionalForm[HoldForm[Integrate[1/(1+x^2)^n,x]==x/(2n-1)*Sum[Product[2n-2i+1,{i,k}]/2^k/Product[n-i,{i,k}]/(1+x^2)^(n-k),{k,n-1}]+Factorial2[2n-3]/2^(n-1)/(n-1)!*ArcTan[x]]]

We rewrite it into more convenient notation:

TraditionalForm[HoldForm[Integrate[1/(1+x^2)^n,x]==x/(2n-1)*Sum[Pochhammer[n-k+1/2,k]/Pochhammer[n-k,k]/(1+x^2)^(n-k),{k,1,n-1}]+Factorial2[2n-3]/2^(n-1)/(n-1)!*ArcTan[x]]]

Doing the same for the expression on the right-hand side delivers a recurrence of order 3 (which is not too bad; but we can do better).

Annihilator[x/(2n-1)*Sum[Pochhammer[n-k+1/2,k]/Pochhammer[n-k,k]/(1+x^2)^(n-k),{k,1,n-1}]+Factorial2[2n-3]/2^(n-1)/(n-1)!*ArcTan[x],S[n]]

By including the derivative w.r.t. x into the set of Ore operators, we find (besides mixed relations) a recurrence of order 2, that matches the one for the integral!

Annihilator[x/(2n-1)*Sum[Pochhammer[n-k+1/2,k]/Pochhammer[n-k,k]/(1+x^2)^(n-k),{k,1,n-1}]+Factorial2[2n-3]/2^(n-1)/(n-1)!*ArcTan[x],{S[n],Der[x]}]

Since the identity does not hold for n=0, we compare the initial values n=1 and n=2.

mySum[expr_,{k_,a_Integer,b_Integer}]:=Plus@@Table[expr/.k->kk,{kk,a,b}]

Table[x/(2n-1)*mySum[Pochhammer[n-k+1/2,k]/Pochhammer[n-k,k]/(1+x^2)^(n-k),{k,1,n-1}]+Factorial2[2n-3]/2^(n-1)/(n-1)!*ArcTan[x],{n,2}]

Factor[D[%,x]]

3.661.3  (∫_0^π)(((cos(t)+i cos(x) sin(t)))^n)dx==π (Subscript[P, n])(cos(t))

Gradshteyn/Ryzhik  3.661(3)  (S. 390):

TraditionalForm[HoldForm[Integrate[(Cos[t]+I Sin[t]Cos[x])^n,{x,0,Pi}]==Pi LegendreP[n,Cos[t]]]]

The integrand is not ∂-finite and we cannot find any relation automatically:

Annihilator[(Cos[t]+I Sin[t]Cos[x])^n,{S[n],Der[t],Der[x]}]

Performing the obvious substitution u=cos(x) helps. Further we replace cos(t)=c (in order to prevent Mathematica to do any simplification), and obtain the identity

TraditionalForm[HoldForm[Integrate[(c+I*u*Sqrt[1-c^2])^n/Sqrt[1-u^2],{u,-1,1}]==Pi LegendreP[n,c]]]

Now it is quite routine to prove the identity:

Annihilator[Integrate[(c+I*u*Sqrt[1-c^2])^n/Sqrt[1-u^2],{u,-1,1}],{Der[c],S[n]}]

Annihilator[LegendreP[n,c],{Der[c],S[n]}]

uts=UnderTheStaircase[%]

ApplyOreOperator[uts,(c+I*u*Sqrt[1-c^2])^n/Sqrt[1-u^2]]/.n->0/.c->0

Integrate[%,{u,-1,1}]

ApplyOreOperator[uts,Pi LegendreP[n,c]]/.n->0/.c->0

If we do not replace cos(t) by c, we can still do something, but this is less safe:

Annihilator[Integrate[(Cos[t]+I*u*Sqrt[1-Cos[t]^2])^n/Sqrt[1-u^2],{u,-1,1}],{S[n]}]

Annihilator[LegendreP[n,Cos[t]],{S[n]}]

3.953.1  (∫_0^∞)(x^(μ-1)) exp(-γ x-β (x^2)) sin(a x)dx==-(i/(2 (((2 β))^(μ/2))))exp((((γ^2)-(a^2))/(8 β))) (Gamma[μ]) (exp(-((i a γ)/(4 β))) (ParabolicCylinderD[-μ, ((γ-i a)/(Sqrt[2 β]))])-exp(((i a γ)/(4 β))) (ParabolicCylinderD[-μ, ((γ+i a)/(Sqrt[2 β]))]))

Gradshteyn/Ryzhik 3.953.1  (p. 503)

TraditionalForm[HoldForm[Integrate[x^(m-1)*Exp[-c*x-b*x^2]*Sin[a*x],{x,0,Infinity}]==(-I/(2*(2b)^(m/2)))*(Exp[(c^2-a^2)/(8b)]*Gamma[m]*(Exp[-I*a*c/(4b)]*ParabolicCylinderD[-m,(c-I*a)/Sqrt[2b]]-Exp[I*a*c/(4b)]*ParabolicCylinderD[-m,(c+I*a)/Sqrt[2b]]))]/.{b->β,c->γ,m->μ}]

For convenience, we replace the greek letters by roman ones.

rexpr=-I/(2*(2b)^(m/2))*Exp[(c^2-a^2)/(8b)]*Gamma[m]*(Exp[-I*a*c/(4b)]*ParabolicCylinderD[-m,(c-I*a)/Sqrt[2b]]-Exp[I*a*c/(4b)]*ParabolicCylinderD[-m,(c+I*a)/Sqrt[2b]]);

Mathematica can not do this integral:

Integrate[x^(m-1)*Exp[-c*x-b*x^2]*Sin[a*x],{x,0,Infinity},Assumptions->Re[m]>-1&&Re[b]>0&&a>0]

First observe that the identity holds for a=0 as well. We now compute annihilating ideal for both sides.

lhs=Annihilator[Integrate[x^(m-1)*Exp[-c*x-b*x^2]*Sin[a*x],{x,0,Infinity}],{S[m],Der[a],Der[b],Der[c]},Assumptions->Re[m]>-1&&Re[b]>0&&a≥0]

rhs=Annihilator[rexpr,{S[m],Der[a],Der[b],Der[c]}]

GBEqual[lhs,rhs]

Length[UnderTheStaircase[#]]&/@{lhs,rhs}

OreReduce[rhs,lhs]

We don't get the same ideal, but at least, rhs is a subideal of lhs. Why do we get less relations for the right-hand side?
The reason is the (Sqrt[b]) that appears in the argument of the parabolic cylinder function. Annihilator computes operators with coefficients  that are polynomials in the variables which correspond to the Ore operators in the algebra. This can increase the order of the resulting operators:

Annihilator[ParabolicCylinderD[-m,Sqrt[b]],S[m]]

Annihilator[ParabolicCylinderD[-m,Sqrt[b]],{Der[b],S[m]},MonomialOrder->Lexicographic]

However, the structure of our identity in question is special in the sense, that if we computed with (Sqrt[b]) in the coefficients, then all occurences of it would disappear in the final result. This can be reproduced by the following trick: We introduce a new variable b1=(Sqrt[b]) in the original expression, and in the end perform the corresponding ∂-finite substitution (and obtain the same annihilating ideal as for the left side):

rexpr1=rexpr/.b->b1^2/.(b1^2)^a1_->b1^(2a1)

rhs1=DFiniteSubstitute[Annihilator[rexpr1,{S[m],Der[a],Der[b1],Der[c]}],{b1->Sqrt[b]},Algebra->OreAlgebra[S[m],Der[a],Der[b],Der[c]]]

GBEqual[lhs,rhs1]

Alternatively, we can set Der[b] aside. In the following we will also not consider the shift in m, since this restricts us to the integers. And the initial values that we have to check are still simple enough with symbolic m and b.

lhs=Annihilator[Integrate[x^(m-1)*Exp[-c*x-b*x^2]*Sin[a*x],{x,0,Infinity}],{Der[a],Der[c]},Assumptions->Re[m]>-1&&Re[b]>0&&a≥0]

rhs=Annihilator[rexpr,{Der[a],Der[c]}]

We compare those initial  values that correspond to the monomials under the stairs.

uts=UnderTheStaircase[lhs]

ApplyOreOperator[uts,x^(m-1)*Exp[-c*x-b*x^2]*Sin[a*x]]/.{a->0,c->0}

Integrate[%,{x,0,Infinity},Assumptions->Re[m]>-1&&Re[b]>0]

FullSimplify[ApplyOreOperator[uts,rexpr]/.{a->0,c->0}]

4.535.1  (∫_0^1)((arctan(p x))/(1+(p^2) x))dx==(1/(2 (p^2)))arctan(p) log(1+(p^2))

Gradshteyn/Ryzhik 4.535.1  (p. 603)

TraditionalForm[HoldForm[Integrate[ArcTan[p*x]/(1+p^2*x),{x,0,1}]==(1/(2p^2))*ArcTan[p]*Log[1+p^2]]]

Mathematica takes quite long and gives a very complicated answer:

Integrate[ArcTan[p*x]/(1+p^2*x),{x,0,1}]

Compute a differential equation for either side:

lhs=Annihilator[Integrate[ArcTan[p*x]/(1+p^2*x),{x,0,1}],Der[p]]

rhs=Annihilator[1/(2p^2)*ArcTan[p]*Log[1+p^2],Der[p]]

OreReduce[lhs,rhs]

Why does the differential equation for the left side overshoot? The reason lies in the inhomogeneous part that remains after creative telescoping:

ct=CreativeTelescoping[ArcTan[p*x]/(1+p^2*x),Der[x],Der[p]]

delta=ApplyOreOperator[ct[[2,1]],ArcTan[p*x]/(1+p^2*x)]

FullSimplify[(delta/.x->1)-Limit[delta,x->0]]

Annihilator[%,Der[p]]

We get a third-order differential equation for the inhomogeneous part. Annihilator in this case uses the closure property "plus", but we can do better by observing that this expression can be written as an operator application (recall the derivative of the arctangent function!):

ApplyOreOperator[-2p*Der[p]-1,ArcTan[p]]

Now we get a second-order differential equation (and finally the same 4th-order DE as we obtained for the right-hand side)!

Annihilator[ApplyOreOperator[-2p*Der[p]-1,ArcTan[p]],Der[p]]

First[%]**ct[[1,1]]

GBEqual[%,rhs]

Now we could either compare initial values, or assume that we did not know the evaluation of the integral and recover it from solving the differential equation.

sol=DSolve[ApplyOreOperator[First[rhs],f[p]]==0,f[p],p][[1,1,2]]

inits=Integrate[Evaluate[FullSimplify[Table[D[ArcTan[p*x]/(1+p^2*x),{p,i}],{i,0,3}]/.p->0]],{x,0,1}]

From the fact that the integral is 0 for p=0, we can exclude the first two solutions, since they tend to infinity for p->0.

sol=sol/.C[1]->0/.C[2]->0

Limit[Simplify[Table[D[sol,{p,i}],{i,0,3}]],p->0]

Solve[Thread[%==inits],{C[3],C[4]}]

sol/.First[%]

4.539  (∫_0^∞)(x^(s-1)) (tan^(-1))(a exp(-x))dx==(2^(-s-1)) (Gamma[s]) a (LerchPhi[-(a^2), s+1, (1/2)])

4.539 (p. 604) states:

TraditionalForm[HoldForm[Integrate[x^(s-1)*ArcTan[a*Exp[-x]],{x,0,Infinity}]==2^(-s-1)*Gamma[s]*a*LerchPhi[-a^2,s+1,1/2]]]

We perform the obvious substitution in order to bring the integrand to a form that is closer to holonomic.

TraditionalForm[HoldForm[Integrate[(-Log[t])^(s-1)/t*ArcTan[a*t],{t,0,1}]]]

A relation that holds for both left-hand and right-hand side, is computed automatically:

Annihilator[Integrate[(-Log[t])^(s-1)/t*ArcTan[a*t],{t,0,1}],{Der[a],S[s]},Assumptions->s>0]

Annihilator[2^(-s-1)*Gamma[s]*a*LerchPhi[-a^2,s+1,1/2]]

The leading monomial  Der[a]*S[s]  suggests that we have to compare initial values for 
1)  a=0 and arbitrary integer s>0,
2)  s=1 and arbitrary n-th derivative with respect to a, evaluated at some point.

For a=0, the integrand becomes 0, as well as the right-hand side.

For s=1 we get:

Integrate[ArcTan[a*t]/t,{t,0,1},Assumptions->Element[a,Reals]]

If we rewrite a*LerchPhi[-a^2, 2, 1/2]/4 = a/4*Sum[(-a^2)^k / (k+1/2)^2, {k,0,Infinity}] (i.e., the right-hand side for s=1) in hypergeometric PFQ notation, it is easily seen, that it matches the above expression (even for symbolic a). 
(Mathematica versions earlier than 7 give a different output, using PolyLog).

Now we compare the first a-derivative of both sides for s=1 (this is not necessarily required):

Integrate[D[ArcTan[a*t]/t,a],{t,0,1},Assumptions->Element[a,Reals]&&a≠0]

Together[D[2^(-s-1)*Gamma[s]*a*LerchPhi[-a^2,s+1,1/2],a]/.s->1]

We get the same result for symbolic a. Hence all higher derivatives with respect to a, evaluated at 0, will also agree.

6.512.1  (∫_0^∞)(BesselJ[m, a x]) (BesselJ[n, b x])dx==(((b^n) (a^(-n-1)) (Gamma[(1/2) (m+n+1)]) (Hypergeometric2F1[(1/2) (m+n+1), (1/2) (n-m+1), n+1, ((b^2)/(a^2))]))/((Gamma[n+1]) (Gamma[(1/2) (m-n+1)])))
with Zeilberger's slow algorithm. CT does not work here.

6.512.1  on page 660. For a>0, b>0, Re(m+n)>-1, b<a:

TraditionalForm[HoldForm[Integrate[BesselJ[m,a*x]*BesselJ[n,b*x],{x,0,Infinity}]==b^n*a^(-n-1)*Gamma[(m+n+1)/2]/Gamma[n+1]/Gamma[(m-n+1)/2]*Hypergeometric2F1[(m+n+1)/2,(n-m+1)/2,n+1,b^2/a^2]]]

When we try this identity directly, we don't get a result, even when including all kinds of assumptions (depending on which Mathematica version is used, the limit computation may consume a lot of memory):

MemoryConstrained[Annihilator[Integrate[BesselJ[m,a*x]*BesselJ[n,b*x],{x,0,Infinity}],{Der[a],Der[b],S[m],S[n]},Assumptions->a>0&&b>0&&b<a&&m>0&&n>0&&Element[{m,n},Integers]],10^8]

Therefore, we have a closer look at what happens after creative telescoping:

{ann,delta}=CreativeTelescoping[BesselJ[m,a*x]*BesselJ[n,b*x],Der[x],{Der[a],Der[b],S[m],S[n]}]

The problem that arises here is that BesselJ[n,x] is asymptotically equivalent with Sqrt[2/(πx)] for x→∞  (see, e.g., 9.2.1 in Abramowitz/Stegun). This means that already with the delta part of the first relation, -x, we get into troubles, because x*Bessel*Bessel does not tend to zero:

Plot[x*BesselJ[3,x]*BesselJ[4,x],{x,0,1000}]

With[{h=10^13},Plot[x*BesselJ[5,x]*BesselJ[6,x],{x,h,h+10}]]

However, if we do not care about the delta part, we find that the principal parts agree with the annihilating ideal that is computed for the right-hand side:

pp=First[CreativeTelescoping[BesselJ[m,a*x]*BesselJ[n,b*x],Der[x],{Der[a],Der[b],S[m],S[n]}]]

rhs=Annihilator[b^n*a^(-n-1)*Gamma[(m+n+1)/2]/Gamma[n+1]/Gamma[(m-n+1)/2]*Hypergeometric2F1[(m+n+1)/2,(n-m+1)/2,n+1,b^2/a^2],{Der[a],Der[b],S[m],S[n]}]

GBEqual[pp,rhs]

Aiming at a rigorous proof, Zeilberger's original (slow) algorithm seems to be the perfect tool for this specific problem. Recall that this algorithm first completely eliminates x, and then does the splitting into principal part and delta part. Hence no x occurs in the delta part and we can expect the corresponding limits to converge. Additionally we would like to have no derivatives w.r.t. a or b in the delta parts, since they cause the same problems. 
Operators with these properties can be found by means of our command FindRelation, where the first condition can be encoded by the
option Eliminate, and the second by the option Support: we give all monomials up to total degree 2, but omit Der[a]*Der[x] and Der[b]*Der[x].

ops={Der[x],Der[a],Der[b],S[m],S[n]};
supp=Complement[Join[{1},ops,Flatten[Outer[Times,ops,ops]]],Der[x]*{Der[a],Der[b]}]

rels=FindRelation[Annihilator[BesselJ[m,a*x]*BesselJ[n,b*x],ops],Eliminate->x,Support->supp]

We now have to separate principal parts and delta parts by hand:

pps=OrePolynomialSubstitute[rels,{Der[x]->0}]

deltas=Together[(Normal/@(rels-pps))/Der[x]]

It is now obvious that all limits will tend to zero.

ass=Assumptions->a>0&&b>0&&m≥0&&n≥0;
inhom=((Limit[#,x->Infinity,ass]-FullSimplify[#/.x->0,ass]))&/@ApplyOreOperator[deltas,BesselJ[m,a*x]*BesselJ[n,b*x]]

A final Gröbner basis computation is needed to convert these into our original setting:

ann=OreGroebnerBasis[pps,OreAlgebra[Der[a],Der[b],S[m],S[n]],MonomialOrder->DegreeLexicographic]

Again, we get the same annihilating ideal as for the right-hand side:

GBEqual[rhs,ann]

Initial values:

uts=UnderTheStaircase[rhs]

ApplyOreOperator[uts,BesselJ[m,a*x]*BesselJ[n,b*x]]/.m->0/.n->0/.a->1/.b->1/2

In Mathematica 8 and 9, we get the following error message, while previously the last integral evaluated to 4/(3*Pi) (4*EllipticE[1/4] - 3*EllipticK[1/4]). Also later versions of Mathematica give this result.

Integrate[%,{x,0,Infinity}]

ApplyOreOperator[uts,b^n*a^(-n-1)*Gamma[(m+n+1)/2]/Gamma[n+1]/Gamma[(m-n+1)/2]*Hypergeometric2F1[(m+n+1)/2,(n-m+1)/2,n+1,b^2/a^2]]/.m->0/.n->0/.a->1/.b->1/2

FullSimplify[%-%%]

6.774  (∫_0^∞)(1/(Sqrt[(x^2)+(a^2)]))log((((Sqrt[(x^2)+(a^2)])+x)/((Sqrt[(x^2)+(a^2)])-x))) (BesselJ[0, b x])dx==((BesselK[0, ((a b)/2)])^2)

On page 747 we have for Re(a)>0 and b>0:

TraditionalForm[HoldForm[Integrate[Log[(Sqrt[x^2+a^2]+x)/(Sqrt[x^2+a^2]-x)]*BesselJ[0,b*x]/Sqrt[x^2+a^2],{x,0,Infinity}]==BesselK[0,a*b/2]^2]]

Annihilator[Integrate[Log[(Sqrt[x^2+a^2]+x)/(Sqrt[x^2+a^2]-x)]*BesselJ[0,b*x]/Sqrt[x^2+a^2],{x,0,Infinity}],{Der[a],Der[b]},Assumptions->Re[a]>0&&b>0]

So we have to do some work by hand.

{ann,delta}=CreativeTelescoping[Log[(Sqrt[x^2+a^2]+x)/(Sqrt[x^2+a^2]-x)]*BesselJ[0,b*x]/Sqrt[x^2+a^2],Der[x],{Der[a],Der[b]}]

Limit[ApplyOreOperator[delta[[1]],Log[(Sqrt[x^2+a^2]+x)/(Sqrt[x^2+a^2]-x)]*BesselJ[0,b*x]/Sqrt[x^2+a^2]],x->Infinity,Assumptions->Re[a]>0&&b>0]

Simplify[ApplyOreOperator[delta[[2]],Log[(Sqrt[x^2+a^2]+x)/(Sqrt[x^2+a^2]-x)]*BesselJ[0,b*x]/Sqrt[x^2+a^2]]]

Limit[%,x->Infinity,Assumptions->Re[a]>0&&b>0]

%%/.a->1/.b->1

Plot[%,{x,0,1000}]

So, although the limit of the delta part does not converge, the principal part makes perfect sense:

Annihilator[BesselK[0,a*b/2]^2,{Der[a],Der[b]}]

GBEqual[%,ann]

Fortunately, we can avoid this problem by reversing the order of the operators Der[a] and Der[b]. This will influence the creative telescoping relations in a way that the delta parts vanish:

Annihilator[Integrate[Log[(Sqrt[x^2+a^2]+x)/(Sqrt[x^2+a^2]-x)]*BesselJ[0,b*x]/Sqrt[x^2+a^2],{x,0,Infinity}],{Der[b],Der[a]},Assumptions->Re[a]>0&&b>0]

Annihilator[BesselK[0,a*b/2]^2,{Der[b],Der[a]}]

But still, the initial values are difficult to check:

uts=UnderTheStaircase[%]

FullSimplify[ApplyOreOperator[uts,Log[(Sqrt[x^2+a^2]+x)/(Sqrt[x^2+a^2]-x)]*BesselJ[0,b*x]/Sqrt[x^2+a^2]]/.a->1/.b->1]

Integrate[%,{x,0,Infinity}]

Unfortunately, Mathematica cannot do these integrals. So we content ourselves with a numerical check:

%/.Integrate->NIntegrate

ApplyOreOperator[uts,BesselK[0,a*b/2]^2]/.a->1/.b->1

N[%]

7.313  (∫_(-1)^1)(((1-(x^2)))^(nu-(1/2))) (Subscript[C, m]^((nu)))(x) (Subscript[C, n]^((nu)))(x)dx

The following identity (p. 795) holds for m ≠ n and Re[nu]>-1/2.

TraditionalForm[HoldForm[Integrate[(1-x^2)^(nu-1/2)*GegenbauerC[m,nu,x]*GegenbauerC[n,nu,x],{x,-1,1}]==0]]

The exception for m=n is not recognized automatically by our program, so we get 1 as an annihilator, whatever m and n are:

Annihilator[Integrate[(1-x^2)^(nu-1/2)*GegenbauerC[m,nu,x]*GegenbauerC[n,nu,x],{x,-1,1}],{S[nu],S[m],S[n]},Assumptions->nu>1/2]

To see that for m=n, something goes wrong, we have to look at the certificate of creative telescoping. Indeed, the factor (m-n) appears in the denominator:

CreativeTelescoping[(1-x^2)^(nu-1/2)*GegenbauerC[m,nu,x]*GegenbauerC[n,nu,x],Der[x],{S[nu],S[m],S[n]}]

In the case m=n we get (for Re[nu]>-1/2,  formula 7.313.2):

TraditionalForm[HoldForm[Integrate[(1-x^2)^(nu-1/2)*GegenbauerC[n,nu,x]^2,{x,-1,1}]==Pi*2^(1-2nu)*Gamma[2nu+n]/n!/(n+nu)/Gamma[nu]^2]]

This evaluation can be FOUND automatically! We compute a recurrence in n for the integral and one initial value for n=0 (for the automatic treatment, we have to sharpen the conditions on nu):

ann=Annihilator[Integrate[(1-x^2)^(nu-1/2)*GegenbauerC[n,nu,x]^2,{x,-1,1}],S[n],Assumptions->Re[nu]>1/2]

(1-x^2)^(nu-1/2)*GegenbauerC[n,nu,x]*GegenbauerC[n,nu,x]/.n->0

Integrate[%,{x,-1,1},Assumptions->Re[nu]>-1/2]/.Sec[nu*Pi]->Gamma[1/2+nu]*Gamma[1/2-nu]/Pi

FullSimplify[RSolve[{ApplyOreOperator[First[ann],f[n]]==0,f[0]==%},f[n],n]]

This evaluation agrees exactly with the formula given in Gradshteyn/Ryzhik.

7.322  (∫_0^(2 a))(((x (2 a-x)))^(nu-(1/2))) (Subscript[C, n]^((nu)))((x/a)-1) (e^(-b x))dx==(((((-1))^n) π Γ(n+2 nu) ((((a/(2 b))))^nu) (e^(-a b)) (Subscript[I, n+nu])(a b))/(n! Γ(nu)))
(misprint in the book!)

The formula 7.322  (S. 797) is not correct as we will show in the following.

TraditionalForm[HoldForm[Integrate[(x(2a-x))^(nu-1/2)*GegenbauerC[n,nu,x/a-1]*E^(-b*x),{x,0,2a}]==(-1)^n*Pi*Gamma[2nu+n]/n!/Gamma[nu]*(a/(2b))^nu*E^(-a*b)*BesselI[nu+n,a*b]]]

Annihilating ∂-finite ideal for the integrand

ann=Annihilator[(x(2a-x))^(nu-1/2)*GegenbauerC[n,nu,x/a-1]*E^(-b*x),{Der[x],Der[a],Der[b],S[n],S[nu]}]

We can use Zeilberger's slow algorithm or Takayama (assuming natural boundaries)

Timing[gb=OreGroebnerBasis[ann,OreAlgebra[x,Der[x],Der[a],Der[b],S[n],S[nu]],MonomialOrder->EliminationOrder[1]];]

LeadingPowerProduct/@gb

gb=OreGroebnerBasis[OrePolynomialSubstitute[Drop[gb,-1],{Der[x]->0}],OreAlgebra[Der[a],Der[b],S[n],S[nu]],MonomialOrder->DegreeLexicographic];

UnderTheStaircase[gb]

Timing[tak=Takayama[ann,{x}];]

UnderTheStaircase[tak]

We get the same annihilating ideal for the integral:

GBEqual[gb,tak]

But Chyzak's algorithm is preferable. The only problem is that in this instance the inhomogeneous part cannot be evaluated automatically. Hence we use the option Inhomogeneous->True and simplify it ourselves (it turns out to be 0 anyway).

{lhs,inh}=Annihilator[Integrate[(x(2a-x))^(nu-1/2)*GegenbauerC[n,nu,x/a-1]*E^(-b*x),{x,0,2a}],{Der[a],Der[b],S[n],S[nu]},Assumptions->nu≥1,Inhomogeneous->True];

Simplify[Simplify[ReleaseHold[inh/.Limit->myLimit]]/.myLimit->Limit,Assumptions->nu≥1]

rhs=Annihilator[(-1)^n*Pi*Gamma[2nu+n]/n!/Gamma[nu]*(a/(2b))^nu*E^(-a*b)*BesselI[nu+n,a*b],{Der[a],Der[b],S[n],S[nu]}]

GBEqual[lhs,rhs]

We see that with Zeilberger slow and Takayama we got only a subideal

OreReduce[tak,rhs]

Using the result from Chyzak's algorithm we have to compare two initial values:

uts=UnderTheStaircase[lhs]

ApplyOreOperator[uts,(x(2a-x))^(nu-1/2)*GegenbauerC[n,nu,x/a-1]*E^(-b*x)]/.a->1/.b->1/.n->1/.nu->0

Integrate[%,{x,0,2}]

ApplyOreOperator[uts,(-1)^n*Pi*Gamma[2nu+n]/n!/Gamma[nu]*(a/(2b))^nu*E^(-a*b)*BesselI[nu+n,a*b]]/.a->1/.b->1/.n->1/.nu->0

Obviously, on the right hand side there is a factor 2 missing. Hence we found a misprint in the Gradshteyn/Ryzhik book.

7.349  (∫_(-1)^1)(((Subscript[T, n])(1-(x^2) y))/(Sqrt[1-(x^2)]))dx==((π ((LegendreP[n, 1-y])+(LegendreP[n-1, 1-y])))/2)

On page 801, we find:

TraditionalForm[HoldForm[Integrate[ChebyshevT[n,1-x^2 y]/Sqrt[1-x^2],{x,-1,1}]==Pi/2*(LegendreP[n,1-y]+LegendreP[n-1,1-y])]]

We get two creative telescoping  relations for the integrand.

CreativeTelescoping[ChebyshevT[n,1-x^2 y]/Sqrt[1-x^2],Der[x],{S[n],Der[y]}]

We can check by hand that the delta parts vanish...

(Limit[#,x->1]-Limit[#,x->-1])&[ApplyOreOperator[Last[%],ChebyshevT[n,1-x^2 y]/Sqrt[1-x^2]]]

... or we can just use the command Annihilator, that does all these steps automatically:

Annihilator[Integrate[ChebyshevT[n,1-x^2 y]/Sqrt[1-x^2],{x,-1,1}],{S[n],Der[y]}]

Now the right-hand side:

Annihilator[Pi/2*(LegendreP[n,1-y]+LegendreP[n-1,1-y]),{S[n],Der[y]}]

This is a classical case where a little bit of human intervention is desired. When the closure property DFinitePlus is used, the annihilating ideal often overshoots, especially when the two functions to be combined are very "similar" (as it is the case here). We can get a better result by applying the closure property DFiniteOreAction. To this purpose, the sum of the two Legendre polynomials has to be reformulated as applying the operator S[n]+1 to the function LegendreP[n-1,1-y].

rhs=Annihilator[ApplyOreOperator[S[n]+1,LegendreP[n-1,1-y]],{S[n],Der[y]}]

Now we get the same annihilating ideal, and have to compare two initial values.

uts=UnderTheStaircase[rhs]

For n=1:

Integrate[ApplyOreOperator[uts,ChebyshevT[n,1-x^2 y]/Sqrt[1-x^2]]/.n->1/.y->0,{x,-1,1}]

Together[ApplyOreOperator[uts,Pi/2*(LegendreP[n,1-y]+LegendreP[n-1,1-y])]/.n->1]/.y->0

7.512.5  (∫_0^1)(x^(r-1)) (((1-x))^(s-1)) (Hypergeometric2F1[a, b, c, x])dx==(((Gamma[r]) (Gamma[s]) (Subscript[, 3])(Subscript[F, 2])(a,b,r;c,r+s;1))/(Gamma[r+s]))

HoldForm[Integrate[x^(r-1)(1-x)^(s-1)Hypergeometric2F1[a,b,c,x],{x,0,1}]==Gamma[r] Gamma[s]/Gamma[r+s] HypergeometricPFQ[{a,b,r},{c,r+s},1]]//TraditionalForm
(* Re[r]>0, Re[s]>0, Re[c+s-a-b]>0 *)

We prove this identity for integers a,b,c,r,s such that c + s - a - b > 0  and  a ≥ 0, b ≥ 0, c ≥ 1, r ≥ 1, s ≥ 2. The latter condition is because for s ≥ 1, the delta part in the following computation cannot be simplified. However, the special case s=1 could be treated separately in the very same manner.

lhs=Annihilator[Integrate[x^(r-1)*(1-x)^(s-1)*Hypergeometric2F1[a,b,c,x],{x,0,1}],{S[a],S[b],S[c],S[r],S[s]},Assumptions->s>1]

Now the right-hand side. Substituting 1 into a hypergeometric PFQ often leads to the following problem.

rhs=Annihilator[Gamma[r]*Gamma[s]/Gamma[r+s]*HypergeometricPFQ[{a,b,r},{c,r+s},1],{S[a],S[b],S[c],S[r],S[s]}]

We can compute an annihilating ideal for the PFQ by using its definition.

{pfq,inh}=Annihilator[Sum[Pochhammer[a,k]*Pochhammer[b,k]*Pochhammer[r,k]/Pochhammer[c,k]/Pochhammer[r+s,k]/k!,{k,0,Infinity}],{S[a],S[b],S[c],S[r],S[s]},Inhomogeneous->True]

Since we choose the parameters a,b,c,r,s in a way that everything goes well, we believe that the above limits go to 0, although Mathematica does not evaluate them (also not when giving assumptions). Alternatively, having ensured that the bounds of the sum are natural, we can call Takayama and get the same result:

tak=Takayama[Annihilator[Pochhammer[a,k]*Pochhammer[b,k]*Pochhammer[r,k]/Pochhammer[c,k]/Pochhammer[r+s,k]/k!,{S[k],S[a],S[b],S[c],S[r],S[s]}],{k}]

GBEqual[pfq,tak]

With closure properties an annihilating ideal for the right-hand side is now obtained, and it turns out that the previously one computed indeed was correct:

GBEqual[rhs,DFiniteTimes[Annihilator[Gamma[r]*Gamma[s]/Gamma[r+s],{S[a],S[b],S[c],S[r],S[s]}],pfq]]

So we found the same annihilating ideal for both sides.

GBEqual[rhs,lhs]

Initial values. As we can see, the leading coefficients do have a lot of zeros, and we have to ask ourselves where we cannot apply any of the recurrences (all such values have to be given as initial conditions).

Factor[lhs]//TableForm

AnnihilatorSingularities[lhs,{0,0,1,1,1}]

As we see here, there are infinitely many singular points . It turns out that they all lie in the hyperplane generated by the condition c + s - a - b > 0:

{a,b,c,r,s}/.Cases[%,{repl_,Except[True]}:>repl]

Apply[(#3+#5-#1-#2)&,%,{1}]

Hence we should include this condition when looking for singular points. Then only a finite collection of points is concerned, and they can easily be checked.

sing=AnnihilatorSingularities[lhs,{0,0,1,1,1},Assumptions->c+s-a-b>0]

MyInt[x^(r-1)(1-x)^(s-1)Hypergeometric2F1[a,b,c,x],{x,0,1}]==Gamma[r] Gamma[s]/Gamma[r+s] HypergeometricPFQ[{a,b,r},{c,r+s},1]/.(First/@sing)/.MyInt->Integrate

9.181  Hypergeometric functions of two variables

The definition is 9.180, the differential equation is 9.181 (both on p. 1018). For |x|<1 and |y|<1 we define

TraditionalForm[HoldForm[Subscript[F,1][a,b,bp,c;x,y]==Sum[Pochhammer[a,m+n]*Pochhammer[b,m]*Pochhammer[bp,n]/(Pochhammer[c,m+n]*m!*n!)*x^m*y^n,{m,0,Infinity},{n,0,Infinity}]]]

Since this double sum has natural boundaries, we can apply Takayama's algorithm in order to compute a system of partial differential equations for it.

ann=Annihilator[Pochhammer[a,m+n]*Pochhammer[b,m]*Pochhammer[bp,n]/(Pochhammer[c,m+n]*m!*n!)*x^m*y^n,{S[m],S[n],Der[x],Der[y]}]

pde=Takayama[ann,{m,n}]

When searching for a partial differential equation that does not contain bp or b, respectively, we obtain exactly the two PDEs given as 9.181.1:

Factor[First[OreGroebnerBasis[pde,OreAlgebra[bp,Der[x],Der[y]],MonomialOrder->EliminationOrder[1]]]]

Factor[First[OreGroebnerBasis[pde,OreAlgebra[b,Der[y],Der[x]],MonomialOrder->EliminationOrder[1]]]]

Olver's problem set  (from Abramowitz/Stegun)

Introduction

In an e-mail to Peter Paule, Frank Olver wrote:

"The writing of DLMF Chapter BS by Leonard Maximon and myself is now largely complete;... However, a problem has arisen in connection with about a dozen formulas from Chapter 10 of Abramowitz and Stegun for which we have not yet tracked down proofs,and the author of this chapter, Henry Antosiewiecz, died about a year ago. Since it is the editorial policy for the DLMF not to state formulas without indications of proofs, I am hoping that you will be willing to step into the breach and supply verifications by computer algebra methods... I will fax you the formulas later today..."

10.1.39  generating function for SphericalBesselY

HoldForm[Sin[Sqrt[z^2+2z t]]/z==Sum[(-t)^n/n! SphericalBesselY[n-1,z],{n,0,Infinity}]]//TraditionalForm

lhs=Annihilator[Sin[Sqrt[z^2+2z t]]/z,{Der[t],Der[z]}]

The right-hand side we can either compute with Takayama or with CreativeTelescoping:

smnd=Annihilator[(-t)^n/n!SphericalBesselY[n-1,z],{S[n],Der[t],Der[z]}]

rhs1=Takayama[smnd,{n}]

GBEqual[lhs,rhs1]

ct=CreativeTelescoping[smnd,S[n]-1,{Der[t],Der[z]}]

delta=ApplyOreOperator[ct[[2]],(-t)^n/n!SphericalBesselY[n-1,z]]

delta/.n->0

Assuming that the delta part vanishes for n -> Infinity, we get

rhs2=ct[[1]]

GBEqual[lhs,rhs2]

The same we can achive by saying:

Annihilator[Sum[(-t)^n/n! SphericalBesselY[n-1,z],{n,0,Infinity}],{Der[t],Der[z]},Inhomogeneous->True]

Compare initial values:

uts=UnderTheStaircase[lhs]

ApplyOreOperator[uts,Sin[Sqrt[z^2+2z t]]/z]/.t->0/.z->1

For the right-hand side, only the first summand is nonzero:

ApplyOreOperator[uts,1/n! SphericalBesselY[n-1,z]]/.t->0/.z->1/.n->0//FunctionExpand//FullSimplify

10.1.40  generating function for SphericalBesselJ

HoldForm[Cos[Sqrt[z^2-2*z*t]]/z==Sum[t^n/n!*SphericalBesselJ[n-1,z],{n,0,Infinity}]]//TraditionalForm

lhs=Annihilator[Cos[Sqrt[z^2-2*z*t]]/z,{Der[z],Der[t]}]

rhs1=Takayama[Annihilator[t^n/n!*SphericalBesselJ[n-1,z],{S[n],Der[z],Der[t]}],{n}]

rhs2=Annihilator[Sum[t^n/n!*SphericalBesselJ[n-1,z],{n,0,Infinity}],{Der[z],Der[t]},Inhomogeneous->True]

Assuming that the delta part vanishes for n -> Infinity, we get

rhs2=rhs2[[1]]

GBEqual[lhs,rhs2]

Compare initial values:

uts=UnderTheStaircase[lhs]

ApplyOreOperator[uts,Cos[Sqrt[z^2-2*z*t]]/z]/.t->0/.z->1

ApplyOreOperator[uts,Sum[t^n/n!*SphericalBesselJ[n-1,z],{n,0,Infinity}]]

In the first sum, only the summand for n=0 will be nonzero, in the second one only the summand n=1:

{1/n!*SphericalBesselJ[n-1,z]/.z->1/.n->0,D[t^n/n!*SphericalBesselJ[n-1,z],t]/.n->1/.z->1/.t->0}//FullSimplify

10.1.41 - 10.1.44  derivatives of SphericalBessel w.r.t. the first parameter

We use the following transformation rules (10.1.1, 9.1.64, 9.1.65)

transRules={SphericalBesselJ[n_,z_]->Sqrt[Pi/(2z)]*BesselJ[n+1/2,z],SphericalBesselY[n_,z_]->Sqrt[Pi/(2z)]*BesselY[n+1/2,z],Derivative[1,0][BesselJ][nu_,z_]->BesselJ[nu,z]*Log[z/2]-(z/2)^nu*Hold[Sum[(-1)^k*PolyGamma[nu+k+1]*(z^2/4)^k/Gamma[nu+k+1]/k!,{k,0,Infinity}]],Derivative[1,0][BesselY][arg_,z_]:>Cot[arg*Pi]*(D[BesselJ[arg,z],nu]-Pi*BesselY[arg,z])-Csc[arg*Pi]*D[BesselJ[-arg,z],nu]-Pi*BesselJ[arg,z]};

TableForm[TraditionalForm[HoldForm[Evaluate[transRules]]/.(Rule|RuleDelayed)->Equal/.Pattern->myP/.myP[n_,Blank[]]->n/.arg->nu]]

TraditionalForm[HoldForm[D[SphericalBesselJ[nu,x],nu]/.nu->0==1/x*(CosIntegral[2x]*Sin[x]-SinIntegral[2x]*Cos[x])]]

lhs=D[SphericalBesselJ[nu,x]/.transRules,nu]/.transRules/.nu->0

annLhs=Annihilator[lhs,Der[x]]

annRhs=Annihilator[1/x*(CosIntegral[2x]*Sin[x]-SinIntegral[2x]*Cos[x]),Der[x]]

OreGroebnerBasis[Join[annLhs,annRhs]]

annCommon=DFinitePlus[annLhs,annRhs];

uts=UnderTheStaircase[annCommon]

ApplyOreOperator[uts,ReleaseHold[lhs]-1/x*(CosIntegral[2x]*Sin[x]-SinIntegral[2x]*Cos[x])]/.x->1//Simplify

TraditionalForm[HoldForm[D[SphericalBesselJ[nu,x],nu]/.nu->-1==1/x*(CosIntegral[2x]*Cos[x]+SinIntegral[2x]*Sin[x])]]

lhs=D[SphericalBesselJ[nu,x]/.transRules,nu]/.transRules/.nu->-1

annLhs=Annihilator[lhs,Der[x]]

annRhs=Annihilator[1/x*(CosIntegral[2x]*Cos[x]+SinIntegral[2x]*Sin[x]),Der[x]]

OreGroebnerBasis[Join[annLhs,annRhs]]

uts=UnderTheStaircase[DFinitePlus[annRhs,annLhs]]

ApplyOreOperator[uts,ReleaseHold[lhs]-1/x*(CosIntegral[2x]*Cos[x]+SinIntegral[2x]*Sin[x])]/.x->1//Simplify

TraditionalForm[HoldForm[D[SphericalBesselY[nu,x],nu]/.nu->0==1/x*(CosIntegral[2x]*Cos[x]+(SinIntegral[2x]-Pi)*Sin[x])]]

lhs=D[SphericalBesselY[nu,x]/.transRules,nu]//.transRules/.nu->0

annLhs=Annihilator[lhs,Der[x]]

annRhs=Annihilator[1/x*(CosIntegral[2x]*Cos[x]+(SinIntegral[2x]-Pi)*Sin[x]),Der[x]]

uts=UnderTheStaircase[DFinitePlus[annLhs,annRhs]]

ApplyOreOperator[uts,ReleaseHold[lhs]-1/x*(CosIntegral[2x]*Cos[x]+(SinIntegral[2x]-Pi)*Sin[x])]/.x->1//FullSimplify

TraditionalForm[HoldForm[D[SphericalBesselY[nu,x],nu]/.nu->-1==-1/x*(CosIntegral[2x]*Sin[x]-(SinIntegral[2x]-Pi)*Cos[x])]]

lhs=D[SphericalBesselY[nu,x]/.transRules,nu]//.transRules/.nu->-1

annLhs=Annihilator[lhs,Der[x]]

uts=UnderTheStaircase[DFinitePlus[annLhs,annRhs]]

ApplyOreOperator[uts,ReleaseHold[lhs]+1/x*(CosIntegral[2x]*Sin[x]-(SinIntegral[2x]-Pi)*Cos[x])]/.x->1//Simplify

10.1.41  (revisited)

With version 1.2 of HolonomicFunctions, the treatment of derivatives of Bessel functions with respect to the order became much simpler. Those expressions are now known to the Annihilator command.

TraditionalForm[HoldForm[D[SphericalBesselJ[nu,x],nu]/.nu->0==1/x*(CosIntegral[2x]*Sin[x]-SinIntegral[2x]*Cos[x])]]

lhs=Annihilator[Derivative[1,0][SphericalBesselJ][0,x],Der[x]]

rhs=Annihilator[1/x*(CosIntegral[2x]*Sin[x]-SinIntegral[2x]*Cos[x]),Der[x]]

OreReduce[rhs,lhs]

Hence we have to compare 6 initial values, which in this situation is not so easy. But we can also verify that the right-hand side satisfies the fourth-order differential equation, and compare 4 initial values.

ApplyOreOperator[lhs,1/x*(CosIntegral[2x]*Sin[x]-SinIntegral[2x]*Cos[x])]//FunctionExpand

Mathematica cannot evaluate the derivatives of the integrand. We use the representation introduced in the previous section:
                     Derivative[1,0][SphericalBesselJ](n,x)==(Sqrt[(π/(2 x))]) ((BesselJ[n+(1/2), x]) log((x/2))-((((x/2)))^(n+(1/2))) (∑_(k=0)^∞)(((((-1))^k) (PolyGamma[n+k+(3/2)]) (((((x^2)/4)))^k))/((Gamma[n+k+(3/2)]) k!)))
and define a function that evaluates the x-derivatives of Derivative[1,0][SphericalBesselJ](n,x) at certain values.

(* gives the value of Derivative[1,d][SphericalBesselJ][n,xx] *)
dxj[d_Integer/;d≥0,n_,xx_]:=Module[{x},(D[Sqrt[Pi/(2x)]*Log[x/2]*BesselJ[n+1/2,x],{x,d}]/.x->xx)-Sum[(-1)^k*PolyGamma[n+k+3/2]/Gamma[n+k+3/2]/k!*Together[D[Sqrt[Pi/(2x)]*(x/2)^(n+2k+1/2),{x,d}]]/.x->xx,{k,0,Infinity}]]

Table[dxj[d,0,1]-(D[1/x*(CosIntegral[2x]*Sin[x]-SinIntegral[2x]*Cos[x]),{x,d}]/.x->1),{d,0,3}]//FunctionExpand

10.1.48  BesselJ[0, z*Sin[t]], SphericalBesselJ, and LegendreP

TraditionalForm[HoldForm[BesselJ[0,z*Sin[t]]==Sum[(4n+1)*(2n)!/2^(2n)/n!^2*SphericalBesselJ[2n,z]*LegendreP[2n,Cos[t]],{n,0,Infinity}]]]

lhs=Annihilator[BesselJ[0,z*Sin[t]],Der[z]]

rhs=Annihilator[Sum[(4n+1)*(2n)!/2^(2n)/n!^2*SphericalBesselJ[2n,z]*LegendreP[2n,Cos[t]],{n,0,Infinity}],Der[z],Inhomogeneous->True]

tak=Takayama[Annihilator[(4n+1)*(2n)!/2^(2n)/n!^2*SphericalBesselJ[2n,z]*LegendreP[2n,Cos[t]],{S[n],Der[z]}],{n}]

Plug in the left-hand side into the differential equation for the right-hand side:

ApplyOreOperator[First[tak],BesselJ[0,z*Sin[t]]]//FullSimplify

Check initial values:

uts=UnderTheStaircase[tak]

ApplyOreOperator[uts,BesselJ[0,z*Sin[t]]]/.z->0

Cases[ApplyOreOperator[uts,(4n+1)*(2n)!/2^(2n)/n!^2*SphericalBesselJ[2n,z]*LegendreP[2n,Cos[t]]],SphericalBesselJ[n_,_]:>n,Infinity]//Union

Limit[ApplyOreOperator[uts,Sum[(4n+1)*(2n)!/2^(2n)/n!^2*SphericalBesselJ[2n,z]*LegendreP[2n,Cos[t]],{n,0,2}]],z->0]

10.1.49  SphericalBesselJ and SphericalBesselY

TraditionalForm[HoldForm[SphericalBesselJ[n,2z]==-n!*z^(n+1)*Sum[(2n-2k+1)/k!/(2n-k+1)!*SphericalBesselJ[n-k,z]*SphericalBesselY[n-k,z],{k,0,n}]]]

lhs=Annihilator[SphericalBesselJ[n,2z],{S[n],Der[z]}]

rhs=Annihilator[-n!*z^(n+1)*Sum[(2n-2k+1)/k!/(2n-k+1)!*SphericalBesselJ[n-k,z]*SphericalBesselY[n-k,z],{k,0,n}],{S[n],Der[z]}]

Compare initial values:

uts=UnderTheStaircase[rhs]

Limit[ApplyOreOperator[uts,SphericalBesselJ[n,2z]]/.n->0,z->0]

Limit[ApplyOreOperator[uts,-n!*z^(n+1)*(2n-2k+1)/k!/(2n-k+1)!*SphericalBesselJ[n-k,z]*SphericalBesselY[n-k,z]]/.{k->0,n->0},z->0]

10.1.52  sum over SphericalBesselJ^2

TraditionalForm[HoldForm[Sum[SphericalBesselJ[n,z]^2,{n,0,Infinity}]==SinIntegral[2z]/(2z)]]

lhs=Annihilator[Sum[SphericalBesselJ[n,z]^2,{n,0,Infinity}],Der[z],Inhomogeneous->True]

We assume that the above limit is 0:

FullSimplify[ReleaseHold[lhs[[2]]/.HoldPattern[Limit[__]]->0]]

Annihilator[First[%],Der[z]]

lhs=First[%]**lhs[[1,1]]

rhs=Annihilator[SinIntegral[2z]/(2z),Der[z]]

OreReduce[lhs,rhs]

Compare initial values:

uts=UnderTheStaircase[{lhs}]

Limit[ApplyOreOperator[uts,SinIntegral[2z]/(2z)],z->0]

Table[SphericalBesselJ[n,0],{n,0,10}]

We want to avoid that Mathematica evaluates the Sum[...]

mySum[ApplyOreOperator[#,SphericalBesselJ[n,z]^2],{n,0,Infinity}]&/@uts

Only the first terms of the sums contribute:

Limit[Simplify[%/.Infinity->5/.mySum->Sum],z->0]

10.2.30  kind of generating function for BesselI

TraditionalForm[HoldForm[Sinh[Sqrt[z^2-2I z t]]/z==Sum[(-I t)^n/n!*Sqrt[1/2*Pi/z]*BesselI[-n+1/2,z],{n,0,Infinity}]]]

lhs=Annihilator[Sinh[Sqrt[z^2-2I z t]]/z,{Der[t],Der[z]}]

rhs=Annihilator[Sum[(-I t)^n/n!*Sqrt[1/2*Pi/z]*BesselI[-n+1/2,z],{n,0,Infinity}],{Der[t],Der[z]},Inhomogeneous->True]

We assume that the limits are 0:

%/.HoldPattern[Limit[__]]->0//ReleaseHold

rhs=rhs[[1]]

GBEqual[lhs,rhs]

uts=UnderTheStaircase[lhs]

Limit[ApplyOreOperator[uts,Sinh[Sqrt[z^2-2I z t]]/z]/.t->0,z->0]

ApplyOreOperator[uts,(-I t)^n/n!*Sqrt[1/2*Pi/z]*BesselI[-n+1/2,z]]/.t->0//Simplify

We see that only the first summand will contribute:

Limit[%/. 0^n->1/.n->0,z->0]

10.2.31  kind of generating function for BesselI

HoldForm[Cosh[Sqrt[z^2+2I z t]]/z==Sum[(I t)^n/n!*Sqrt[1/2*Pi/z]*BesselI[n-1/2,z],{n,0,Infinity}]]//TraditionalForm

lhs=Annihilator[Cosh[Sqrt[z^2+2I z t]]/z,{Der[t],Der[z]}]

rhs1=Takayama[Annihilator[(I t)^n/n!*Sqrt[1/2*Pi/z]*BesselI[n-1/2,z],{S[n],Der[t],Der[z]}],{n}]

GBEqual[lhs,{I,1}*rhs1]

rhs2=Annihilator[Sum[(I t)^n/n!*Sqrt[1/2*Pi/z]*BesselI[n-1/2,z],{n,0,Infinity}],{Der[t],Der[z]},Inhomogeneous->True]

We may assume that the limits go to zero.

rhs2=First[rhs2]

GBEqual[lhs,rhs2]

Compare initial values:

uts=UnderTheStaircase[lhs]

ApplyOreOperator[uts,Cosh[Sqrt[z^2+2I z t]]/z]/.t->0/.z->1//FullSimplify

ApplyOreOperator[uts,(I t)^n/n!*Sqrt[1/2*Pi/z]*BesselI[n-1/2,z]]/.t->0/.z->1

We see that in both cases only the first summand survives:

%/.(0^n)->1/.n->0//FullSimplify

10.2.38  BesselK

TraditionalForm[HoldForm[BesselK[n+1/2,2z]==n!*Pi^(-1/2)*z^(n+1/2)*Sum[(-1)^k*(2n-2k+1)/k!/(2n-k+1)!*BesselK[n-k+1/2,z]^2,{k,0,n}]]]

lhs=Annihilator[BesselK[n+1/2,2z],{S[n],Der[z]}]

rhs=Annihilator[n!*Pi^(-1/2)*z^(n+1/2)*Sum[(-1)^k*(2n-2k+1)/k!/(2n-k+1)!*BesselK[n-k+1/2,z]^2,{k,0,n}],{S[n],Der[z]}]

uts=UnderTheStaircase[rhs]

ApplyOreOperator[uts,BesselK[n+1/2,2z]]/.n->0/.z->1

ApplyOreOperator[uts,n!*Pi^(-1/2)*z^(n+1/2)*(-1)^k*(2n-2k+1)/k!/(2n-k+1)!*BesselK[n-k+1/2,z]^2]/.{k->0,n->0}/.z->1//FullSimplify

Examples from Semrad's diploma thesis

4.2.1  integral representation of BesselJ

TraditionalForm[HoldForm[BesselJ[n,z]==(z/2)^n/Gamma[n+1/2]/Gamma[1/2]*Integrate[(1-t^2)^(n-1/2)*Cos[z t],{t,-1,1}]]]

lhs=Annihilator[BesselJ[n,z],{Der[z],S[n]}]

rhs=Annihilator[(z/2)^n/Gamma[n+1/2]/Gamma[1/2]*Integrate[(1-t^2)^(n-1/2)*Cos[z t],{t,-1,1}],{Der[z],S[n]},Assumptions->n≥0]

Compare initial values:

uts=UnderTheStaircase[rhs]

ApplyOreOperator[uts,(z/2)^n/Gamma[n+1/2]/Gamma[1/2]*(1-t^2)^(n-1/2)*Cos[z t]]/.z->0/. 0^n->1/.n->0

Integrate[%,{t,-1,1}]

ApplyOreOperator[uts,BesselJ[n,z]]/.z->0/.n->0

4.2.2  integral representation of BesselJ

TraditionalForm[HoldForm[BesselJ[n,z]==(z/2)^n/(2Pi I)*Integrate[t^(-n-1)*Exp[t-z^2/(4t)],{t,c-I*Infinity,c+I*Infinity}]]]

lhs=Annihilator[BesselJ[n,z]]

int=Annihilator[Integrate[t^(-n-1)*Exp[t-z^2/(4t)],{t,c-I*Infinity,c+I*Infinity}],{S[n],Der[z]},Inhomogeneous->True]

DFiniteTimes[Annihilator[(z/2)^n/(2Pi I),{S[n],Der[z]}],First[int]]

uts=UnderTheStaircase[lhs]

ApplyOreOperator[uts,BesselJ[n,z]]/.n->0/.z->0

ApplyOreOperator[uts,(z/2)^n/(2Pi I)*t^(-n-1)*Exp[t-z^2/(4t)]]/.n->0/.z->0

Integrate[%,{t,c-I*n,c+I*n}]

Limit[%[[1,1]],n->Infinity]

4.3.1  integral representation of BesselY  (2 integrals)
(wrong)

TraditionalForm[HoldForm[BesselY[n,z]==2(z/2)^n/Gamma[n+1/2]/Gamma[1/2]*(Integrate[Sin[z t](1-t^2)^(n+1/2),{t,0,1}]-Integrate[Exp[-z t](1+t^2)^(n-1/2),{t,0,Infinity}])]]

lhs=Annihilator[BesselY[n,z],{S[n],Der[z]}]

rhs=Annihilator[2(z/2)^n/Gamma[n+1/2]/Gamma[1/2]*(Integrate[Sin[z t](1-t^2)^(n+1/2),{t,0,1}]-Integrate[Exp[-z t](1+t^2)^(n-1/2),{t,0,Infinity}]),{S[n],Der[z]},Assumptions->z>0&&n≥0]

The annihilating ideal that we computed for rhs is contained in the annihilating ideal of the left hand side:

OreReduce[rhs,lhs]

Check identity for n=0 and z=1, something is wrong here:

2/Gamma[1/2]/Gamma[1/2]*(Integrate[Sin[t](1-t^2)^(1/2),{t,0,1}]-Integrate[Exp[- t](1+t^2)^(-1/2),{t,0,Infinity}])

N[%]

N[BesselY[0,1]]

4.3.2  integral representation of BesselY
(improper use of Takayama's algorithm)

TraditionalForm[HoldForm[BesselY[n,z]==-2/Gamma[1/2-n]/Gamma[1/2]/(z/2)^n*Integrate[Cos[z t]/(t^2-1)^(n+1/2),{t,1,Infinity}]]]

This integral can be found in Abramowitz/Stegun 9.1.24, where also the conditions Re[n]<1/2 and z>0 are given (which have been copied incorrectly in the thesis).

The first condition implies that it does not make sense to consider the operator S[n].

rhs=Annihilator[-2/Gamma[1/2-n]/Gamma[1/2]/(z/2)^n*Integrate[Cos[z t]/(t^2-1)^(n+1/2),{t,1,Infinity}],Der[z],Assumptions->z>0&&Abs[Re[n]]<1/2]

Let's see why we get this error. We do the creative telescoping by hand.

ct=CreativeTelescoping[Cos[z t]/(t^2-1)^(n+1/2),Der[t],Der[z]]

test=Simplify[ApplyOreOperator[ct[[2,1]],Cos[z t]/(t^2-1)^(n+1/2)]]

Unfortunately, the limit for t→∞ does not vanish, but gives ∞. Hence the method of creative telescoping does not enable us to do a clean reasoning; a fortiori the same is true for Takayama's algorithm (as described in the thesis).

Even if we ignore these problems, the next difficulties arise when comparing initial values. When we consider the first derivative w.r.t. z, we encounter an integral that does not converge for any n subject to the given conditions:

D[Cos[z t]/(t^2-1)^(n+1/2),z]/.n->0

Integrate[%,{t,1,Infinity},Assumptions->z>0]

4.4.1  Exp[Sinh[...]], not D-finite

Annihilator[Exp[z Sinh[w]-n*w],{S[n],Der[z],Der[w]}]

4.5.1  integral representation of BesselI

TraditionalForm[HoldForm[BesselI[n,z]==(z/2)^n/Gamma[n+1/2]/Gamma[1/2]*Integrate[Exp[z t]*(1-t^2)^(n-1/2),{t,-1,1}]]]

lhs=Annihilator[BesselI[n,z]]

rhs=Annihilator[(z/2)^n/Gamma[n+1/2]/Gamma[1/2]*Integrate[Exp[z t]*(1-t^2)^(n-1/2),{t,-1,1}],{S[n],Der[z]},Assumptions->n≥0]

uts=UnderTheStaircase[lhs]

ApplyOreOperator[uts,BesselI[n,z]]/.n->0/.z->0

ApplyOreOperator[uts,(z/2)^n/Gamma[n+1/2]/Gamma[1/2]*Exp[z t]*(1-t^2)^(n-1/2)]/.n->0/.z->0

Integrate[%,{t,-1,1}]

4.6.1  integral representation of BesselK
(corrected)

TraditionalForm[HoldForm[BesselK[n,z]==(z/2)^n/Gamma[n+1/2]*Gamma[1/2]*Integrate[Exp[-z t]*(t^2-1)^(n+1/2),{t,1,Infinity}]]]

Typo! Correct version is:

TraditionalForm[HoldForm[BesselK[n,z]==(z/2)^n/Gamma[n+1/2]*Gamma[1/2]*Integrate[Exp[-z t]*(t^2-1)^(n-1/2),{t,1,Infinity}]]]

lhs=Annihilator[BesselK[n,z]]

rhs=Annihilator[(z/2)^n/Gamma[n+1/2]*Gamma[1/2]*Integrate[Exp[-z t]*(t^2-1)^(n-1/2),{t,1,Infinity}],{S[n],Der[z]},Assumptions->z>0&&n≥0]

GBEqual[lhs,rhs]

uts=UnderTheStaircase[lhs]

ApplyOreOperator[uts,BesselK[n,z]]/.n->0/.z->1

ApplyOreOperator[uts,(z/2)^n/Gamma[n+1/2]*Gamma[1/2]*Exp[-z t]*(t^2-1)^(n-1/2)]/.n->0/.z->1

Integrate[%,{t,1,Infinity}]

4.6.2  integral representation of BesselK

TraditionalForm[HoldForm[BesselK[n,z]==1/2*(z/2)^n*Integrate[Exp[-t-z^2/(4t)]/t^(n+1),{t,0,Infinity}]]]

lhs=Annihilator[BesselK[n,z]]

rhs=Annihilator[1/2*(z/2)^n*Integrate[Exp[-t-z^2/(4t)]/t^(n+1),{t,0,Infinity}],{S[n],Der[z]},Assumptions->Element[z,Reals]]

GBEqual[lhs,rhs]

uts=UnderTheStaircase[lhs]

ApplyOreOperator[uts,BesselK[n,z]]/.n->0/.z->1

ApplyOreOperator[uts,1/2*(z/2)^n*Exp[-t-z^2/(4t)]/t^(n+1)]/.n->0/.z->1

Integrate[%,{t,0,Infinity}]

4.7.1  product of BesselJ, trick to eliminate cos(theta) with ChebyshevT

TraditionalForm[HoldForm[BesselJ[n,z]*BesselJ[m,z]==2/Pi*Integrate[BesselJ[m+n,2z*Cos[theta]]*Cos[(m-n)*theta],{theta,0,Pi/2}]]]

Trick: First substitute t=cos(theta), then rewrite cos((m-n) arccos(t)) to ChebyshevT[m-n,t]

TraditionalForm[HoldForm[BesselJ[n,z]*BesselJ[m,z]==Integrate[(2/Pi)*BesselJ[m+n,2z t]*ChebyshevT[m-n,t]/Sqrt[1-t^2],{t,0,1}]]]

lhs=Annihilator[BesselJ[n,z]*BesselJ[m,z],{S[m],S[n],Der[z]}]

CreativeTelescoping takes a bit. We can also use Takayama in this case (see below).

Timing[rhs=Annihilator[Integrate[2/Pi*BesselJ[m+n,2z t]*ChebyshevT[m-n,t]/Sqrt[1-t^2],{t,0,1}],{S[m],S[n],Der[z]},Assumptions->n≥0&&m≥0]]

GBEqual[lhs,rhs]

Timing[rhs2=Takayama[Annihilator[2/Pi*BesselJ[m+n,2z t]*ChebyshevT[m-n,t]/Sqrt[1-t^2],{Der[t],S[m],S[n],Der[z]}],{t}]]

GBEqual[lhs,rhs2]

4.8.1  integral representation for BesselJ
(wrong!)

TraditionalForm[HoldForm[Integrate[BesselJ[n,b t]*Exp[-p^2*t^2]*t^(n+1),{t,0,Infinity}]==b^n/(2p^2)^(n+1)*Exp[-b^2/(2p^2)]]]

lhs=Annihilator[Integrate[BesselJ[n,b t]*Exp[-p^2*t^2]*t^(n+1),{t,0,Infinity}],{Der[b],S[n],Der[p]},Assumptions->Element[p,Reals]&&n≥0&&b≥0]

In other Mathematica versions, this works. Here we cheat a bit:

lhs=Annihilator[Integrate[BesselJ[n,b t]*Exp[-p^2*t^2]*t^(n+1),{t,0,Infinity}],{Der[b],S[n],Der[p]},Inhomogeneous->True][[1]]

rhs=Annihilator[b^n/(2p^2)^(n+1)*Exp[-b^2/(2p^2)],{Der[b],S[n],Der[p]}]

UnderTheStaircase/@{lhs,rhs}

The annihilators are not the same, but both have the vector space dimension 1. Hence the identity seems to be wrong. Numerical counterexample:

subs={n->13,b->3/2,p->4/5}

N[b^n/(2p^2)^(n+1)*Exp[-b^2/(2p^2)]/.subs]

NIntegrate[BesselJ[n,b t]*Exp[-p^2*t^2]*t^(n+1)/.subs,{t,0,Infinity},WorkingPrecision->10]

4.8.2  BesselY and BesselK

TraditionalForm[HoldForm[Integrate[t^n*BesselY[n,a t]/(t^2+k^2),{t,0,Infinity}]==-k^(n-1)*BesselK[n,a k]]]

lhs1=Annihilator[Integrate[t^n*BesselY[n,a t]/(t^2+k^2),{t,0,Infinity}],{S[n],Der[a],Der[k]},Inhomogeneous->True]

delta=Last[lhs1]/. 0^_->0

Are these limits zero? If yes, then we get:

lhs=First[lhs1]

rhs=Annihilator[-k^(n-1)*BesselK[n,a k],{S[n],Der[a],Der[k]}]

UnderTheStaircase/@{lhs,rhs}

OreReduce[lhs,rhs]

Compare initial values:

uts=UnderTheStaircase[lhs]

ApplyOreOperator[uts,t^n*BesselY[n,a t]/(t^2+k^2)]/.a->1/.k->1/.n->0

Integrate[%,{t,0,Infinity}]

ApplyOreOperator[uts,-k^(n-1)*BesselK[n,a k]]/.a->1/.k->1/.n->0//FullSimplify

In Mathematica 9 we get the output in terms of MeijerG functions. Other Mathematica versions give the latter result in terms of BesselK functions, so that equality becomes apparent.

N[%%-%]

4.8.3  BesselJ, BesselY, BesselK

TraditionalForm[HoldForm[Integrate[Exp[-p^2*t^2]*BesselJ[0,a t]*BesselY[0,a t]*t,{t,0,Infinity}]==-1/(2Pi*p^2)*Exp[-a^2/(2p^2)]*BesselK[0,a^2/(2p^2)]]]

lhs=Annihilator[Integrate[Exp[-p^2*t^2]*BesselJ[0,a t]*BesselY[0,a t]*t,{t,0,Infinity}],{Der[a],Der[p]},Assumptions->Element[p,Reals]&&p≠0&&a>0]

rhs=Annihilator[-1/(2Pi*p^2)*Exp[-a^2/(2p^2)]*BesselK[0,a^2/(2p^2)],{Der[a],Der[p]}]

Compare initial values:

uts=UnderTheStaircase[rhs]

ApplyOreOperator[uts,Exp[-p^2*t^2]*BesselJ[0,a t]*BesselY[0,a t]*t]/.a->1/.p->1

Integrate[%,{t,0,Infinity}]

ApplyOreOperator[uts,-1/(2Pi*p^2)*Exp[-a^2/(2p^2)]*BesselK[0,a^2/(2p^2)]]/.a->1/.p->1

4.9.1  indefinite integral for BesselI

TraditionalForm[HoldForm[Integrate[z^(n+1)*BesselI[n,z],z]==z^(n+1)*BesselI[n+1,z]]]

ct=CreativeTelescoping[z^(n+1)*BesselI[n,z],Der[z],{}]

ApplyOreOperator[-ct[[2,1]],z^(n+1)*BesselI[n,z]]//FullSimplify

5.1.2 (1)  Beta and Hypergeometric2F1

TraditionalForm[HoldForm[Beta[x,a,b]==x^a/a*Hypergeometric2F1[a,1-b,a+1,x]]]

lhs=Annihilator[Beta[x,a,b],{S[a],S[b],Der[x]}]

rhs=Annihilator[x^a/a*Hypergeometric2F1[a,1-b,a+1,x],{S[a],S[b],Der[x]}]

Compare initial values:

uts=UnderTheStaircase[rhs]

ApplyOreOperator[uts,Beta[x,a,b]]/.a->1/.b->1/.x->0

ApplyOreOperator[uts,x^a/a*Hypergeometric2F1[a,1-b,a+1,x]]/.a->1/.b->1/.x->0

5.1.2 (2)  Beta and indefinite integral

TraditionalForm[HoldForm[Beta[x,a,b]==Integrate[x^(a-1)*(1-x)^(b-1),x]]]

lhs=Annihilator[Beta[x,a,b],{S[a],S[b]}]

For this indefinite integral we have to do a trick:

rhs=Annihilator[Integrate[t^(a-1)*(1-t)^(b-1),{t,0,x}],{S[a],S[b]},Assumptions->a≥1]

GBEqual[lhs,rhs]

uts=UnderTheStaircase[lhs]

ApplyOreOperator[uts,Beta[x,a,b]]/.a->1/.b->1//Expand

Integrate[ApplyOreOperator[uts,x^(a-1)*(1-x)^(b-1)]/.a->1/.b->1,x]//Expand

5.2.2  StruveH (incorrect argument of "natural boundaries")

TraditionalForm[HoldForm[Integrate[Exp[-a t]*StruveH[0,t],{t,0,Infinity}]==2/Pi/Sqrt[1+a^2]*Log[(1+Sqrt[1+a^2])/a]]]

lhs=Annihilator[Integrate[Exp[-a t]*StruveH[0,t],{t,0,Infinity}],Der[a],Assumptions->a>0]

rhs=Annihilator[2/Pi/Sqrt[1+a^2]*Log[(1+Sqrt[1+a^2])/a],Der[a]]

uts=UnderTheStaircase[rhs]

ApplyOreOperator[uts,Exp[-a t]*StruveH[0,t]]/.a->1

Integrate[%,{t,0,Infinity}]

ApplyOreOperator[uts,2/Pi/Sqrt[1+a^2]*Log[(1+Sqrt[1+a^2])/a]]/.a->1

%-%%//FullSimplify

5.3.2  Hypergeometric1F1  (result can be obtained much easier!)

Prove the recurrence relation

TraditionalForm[HoldForm[b*Hypergeometric1F1[a,b,x]-b*Hypergeometric1F1[a-1,b,x]-x*Hypergeometric1F1[a,b+1,x]==0]]

ann=Annihilator[Hypergeometric1F1[a,b,x],{S[a],S[b]}]

OreReduce[b*S[a]-b-x*S[a]*S[b],ann]

FindRelation[ann,Support->{S[a]*S[b],S[a],1}]

ApplyOreOperator[First[%],F[a,b,x]]/.a->a-1

5.4.2  Hypergeometric2F1  (Mgfun takes quite long)

TraditionalForm[HoldForm[Hypergeometric2F1[a,a+1/2,2a+1,z]==(1/2+Sqrt[1-z]/2)^(-2a)]]

lhs=Annihilator[Hypergeometric2F1[a,a+1/2,2a+1,z],{Der[z],S[a]}]

rhs=Annihilator[(1/2+Sqrt[1-z]/2)^(-2a),{Der[z],S[a]}]

uts=UnderTheStaircase[rhs]

ApplyOreOperator[uts,Hypergeometric2F1[a,a+1/2,2a+1,z]]/.z->0/.a->0

ApplyOreOperator[uts,(1/2+Sqrt[1-z]/2)^(-2a)]/.z->0/.a->0

Comparison with Mgfun:

Timing[Annihilator[Hypergeometric2F1[a,a+1/2,2a+1,z],{Der[z],S[a]}];]

Mgfun needs 35 seconds.

5.4.3  Hypergeometric2F1

TraditionalForm[HoldForm[Hypergeometric2F1[a,b,c,z]==(1-z)^(-a)*Hypergeometric2F1[a,c-b,c,z/(z-1)]]]

lhs=Annihilator[Hypergeometric2F1[a,b,c,z],{S[a],S[b],S[c],Der[z]}]

rhs=Annihilator[(1-z)^(-a)*Hypergeometric2F1[a,c-b,c,z/(z-1)],{S[a],S[b],S[c],Der[z]}]

GBEqual[lhs,rhs]

uts=UnderTheStaircase[lhs]

ApplyOreOperator[uts,Hypergeometric2F1[a,b,c,z]]/.a->0/.b->0/.c->0/.z->0

ApplyOreOperator[uts,(1-z)^(-a)*Hypergeometric2F1[a,c-b,c,z/(z-1)]]/.a->0/.b->0/.c->0/.z->0

5.4.4  quadratic transformation for Hypergeometric2F1

TraditionalForm[HoldForm[Hypergeometric2F1[a,b,2b,z]==(1-z/2)^(-a)*Hypergeometric2F1[a/2,a/2+1/2,b+1/2,(z/(2-z))^2]]]

lhs=Annihilator[Hypergeometric2F1[a,b,2b,z],{S[a],S[b],Der[z]}]

rhs=Annihilator[(1-z/2)^(-a)*Hypergeometric2F1[a/2,a/2+1/2,b+1/2,(z/(2-z))^2],{S[a],S[b],Der[z]}]

GBEqual[rhs,lhs]

uts=UnderTheStaircase[rhs]

ApplyOreOperator[uts,Hypergeometric2F1[a,b,2b,z]]/.a->0/.b->0/.z->0

ApplyOreOperator[uts,(1-z/2)^(-a)*Hypergeometric2F1[a/2,a/2+1/2,b+1/2,(z/(2-z))^2]]/.a->0/.b->0/.z->0

5.4.5  cubic transformation for Hypergeometric2F1
 (wrong!)

TraditionalForm[HoldForm[Hypergeometric2F1[3a,3a+1/2,4a+2/3,z]==(1-9/8)^(-2a)*Hypergeometric2F1[a,a+1/2,2a+5/6,27z^2(z-1)/(9z-8)^2]]]

lhs=Annihilator[Hypergeometric2F1[3a,3a+1/2,4a+2/3,z],{S[a],Der[z]}]

rhs=Annihilator[(1-9/8)^(-2a)*Hypergeometric2F1[a,a+1/2,2a+5/6,27z^2(z-1)/(9z-8)^2],{S[a],Der[z]}]

OreReduce[rhs,lhs]

Test incorrectness:

{Hypergeometric2F1[3a,3a+1/2,4a+2/3,z],(1-9/8)^(-2a)*Hypergeometric2F1[a,a+1/2,2a+5/6,27z^2(z-1)/(9z-8)^2]}/.a->1/.z->0

5.5.1  Hypergeometric2F3 and BesselJ
(Mgfun takes very long)

TraditionalForm[HoldForm[HypergeometricPFQ[{(m+n+1)/2,(m+n+2)/2},{m+1,n+1,m+n+1},-z^2]==Gamma[m+1]*Gamma[n+1]*(2/z)^(m+n)*BesselJ[m,z]*BesselJ[n,z]]]

lhs=Annihilator[HypergeometricPFQ[{(m+n+1)/2,(m+n+2)/2},{m+1,n+1,m+n+1},-z^2],{S[m],S[n],Der[z]}]

rhs=Annihilator[Gamma[m+1]*Gamma[n+1]*(2/z)^(m+n)*BesselJ[m,z]*BesselJ[n,z],{S[m],S[n],Der[z]}]

GBEqual[lhs,rhs]

uts=UnderTheStaircase[lhs]

ApplyOreOperator[uts,Sum[Pochhammer[(m+n+1)/2,k]*Pochhammer[(m+n+2)/2,k]/Pochhammer[m+1,k]/Pochhammer[n+1,k]/Pochhammer[m+n+1,k]/k!*(-z^2)^k,{k,0,3}]]/.m->0/.n->0/.z->0

ApplyOreOperator[uts,Gamma[m+1]*Gamma[n+1]*(2/z)^(m+n)*BesselJ[m,z]*BesselJ[n,z]]/.m->0/.n->0

Limit[%,z->0]

Compare with Mgfun:

Timing[Annihilator[HypergeometricPFQ[{(m+n+1)/2,(m+n+2)/2},{m+1,n+1,m+n+1},-z^2],{Der[z]}];]

Mgfun needs 85 seconds.

Timing[Annihilator[HypergeometricPFQ[{(m+n+1)/2,(m+n+2)/2},{m+1,n+1,m+n+1},-z^2],{S[m],S[n],Der[z]}];]

Mgfun can not compute that:
Error (in Mgfun:-MG_Internals:-expression_to_system) expression is not d-finite, or not recognized as such

5.5.2  recurrence for Hypergeometric0F1
(strange use of Mgfun, result can be obtained much easier!)

TraditionalForm[HoldForm[z*Hypergeometric0F1[b+1,z]+b*(b-1)*Hypergeometric0F1[b,z]-b*(b-1)*Hypergeometric0F1[b-1,z]==0]]

Annihilator[Hypergeometric0F1[b,z],S[b]]

ApplyOreOperator[Factor[First[%]],F[b,z]]/.b->b-1

5.5.3  Hypergeometric2F1^2 and Hypergeometric3F2
(takes quite long)

TraditionalForm[HoldForm[Hypergeometric2F1[a,b,a+b+1/2,z]^2==HypergeometricPFQ[{2a,2b,a+b},{a+b+1/2,2a+2b},z]]]

Timing[lhs=Annihilator[Hypergeometric2F1[a,b,a+b+1/2,z]^2,{S[a],S[b],Der[z]}];]

The following is an example where removing contents does not work well with FactorTermsList (the computation took around 4 hours). That's why we switch to "gcd".

Timing[nm=$NormalizeMethod;$NormalizeMethod="gcd";rhs=Annihilator[HypergeometricPFQ[{2a,2b,a+b},{a+b+1/2,2a+2b},z],{S[a],S[b],Der[z]}];$NormalizeMethod=nm;]

Still, it takes quite long. The annihilator of the 3F2 can also be computed with Takayama's algorithm using the definition of the HypergeometricPFQ function:

Timing[tak=Takayama[Annihilator[Pochhammer[2a,k]*Pochhammer[2b,k]*Pochhammer[a+b,k]/Pochhammer[a+b+1/2,k]/Pochhammer[2a+2b,k]/k!*z^k,{S[k],S[a],S[b],Der[z]}],{k}];]

GBEqual[rhs,tak]

GBEqual[lhs,rhs]

uts=UnderTheStaircase[lhs]

ApplyOreOperator[uts,Hypergeometric2F1[a,b,a+b+1/2,z]^2]/.z->0/.a->0/.b->0

ApplyOreOperator[uts,HypergeometricPFQ[{2a,2b,a+b},{a+b+1/2,2a+2b},z]]/.z->0/.a->0/.b->0

5.6.2  recurrence for GegenbauerC

TraditionalForm[HoldForm[(n+2+m)*GegenbauerC[n+2,m,x]==m*(GegenbauerC[n+2,m+1,x]-GegenbauerC[n,m+1,x])]]

ann=Annihilator[GegenbauerC[n,m,x],{S[m],S[n]}]

First solution: We simply take the recurrence and verify that it lies in the left ideal ann:

ToOrePolynomial[(n+2+m)*GegenbauerC[n+2,m,x]==m*(GegenbauerC[n+2,m+1,x]-GegenbauerC[n,m+1,x]),GegenbauerC[n,m,x]]

OreReduce[%,ann]

Second solution: We want to find a recurrence of the above form (i.e. with the same support):

FindRelation[ann,Support->{S[n]^2,S[m]*S[n]^2,S[m]}]

Third solution: We want to find a recurrence whose coefficients do not depend on x:

FindRelation[ann,Eliminate->x]

5.6.3  Laguerre and BesselJ
(OreSys bug)

TraditionalForm[HoldForm[LaguerreL[n,a,x]==Exp[x]*x^(-a/2)/n!*Integrate[Exp[-t]*t^(n+a/2)*BesselJ[a,2(t*x)^(1/2)],{t,0,Infinity}]]]

lhs=Annihilator[LaguerreL[n,a,x],{S[n],S[a],Der[x]}]

rhs=Annihilator[Exp[x]*x^(-a/2)/n!*Integrate[Exp[-t]*t^(n+a/2)*BesselJ[a,2(t*x)^(1/2)],{t,0,Infinity}],{S[n],S[a],Der[x]},Assumptions->n+a/2>0&&x≥0]

We don't get the same annihilator, but lhs is contained in rhs:

GBEqual[lhs,rhs]

OreReduce[rhs,lhs]

Compare initial values:

uts=UnderTheStaircase[rhs]

ApplyOreOperator[uts,LaguerreL[n,a,x]]/.n->0/.a->0/.x->0

Limit[ApplyOreOperator[uts,Exp[x]*x^(-a/2)/n!*Exp[-t]*t^(n+a/2)*BesselJ[a,2(t*x)^(1/2)]]/.n->0/.a->0,x->0]

Integrate[%,{t,0,Infinity}]

OreSys causes a wrong result on this example:

<<RISC`OreSys`

Package OreSys version 1.2 written by Stefan Gerhold
Copyright Research Institute for Symbolic Computation (RISC),
Johannes Kepler University, Linz, Austria

CreativeTelescoping[Exp[-t]*t^(n+a/2)*BesselJ[a,2(t*x)^(1/2)],Der[t],{S[n],S[a],Der[x]},Method->Zuercher]

5.7 Dellanoy Numbers
(wrong!)

TraditionalForm[HoldForm["D(m,n)"==Sum[Binomial[n,k]*Binomial[m+n-k,n],{k,0,n}]==Sum[Binomial[m,k]*Binomial[n,k],{k,0,n}]]]

lhs=Annihilator[Sum[Binomial[n,k]*Binomial[m+n-k,n],{k,0,n}],{S[m],S[n]}]

rhs=Annihilator[Sum[Binomial[m,k]*Binomial[n,k],{k,0,n}],{S[m],S[n]}]

OreReduce[lhs,rhs]

Table[Sum[Binomial[n,k]*Binomial[m+n-k,n],{k,0,n}]-Sum[Binomial[m,k]*Binomial[n,k],{k,0,n}],{n,0,5},{m,0,5}]

5.7 Dellanoy Numbers
(corrected version)

TraditionalForm[HoldForm["D(m,n)"==Sum[Binomial[n,k]*Binomial[m+n-k,n],{k,0,n}]==Sum[2^k*Binomial[m,k]*Binomial[n,k],{k,0,n}]]]

lhs=Annihilator[Sum[Binomial[n,k]*Binomial[m+n-k,n],{k,0,n}],{S[m],S[n]}]

rhs=Annihilator[Sum[2^k*Binomial[m,k]*Binomial[n,k],{k,0,n}],{S[m],S[n]}]

GBEqual[lhs,rhs]

uts=UnderTheStaircase[lhs]

ApplyOreOperator[uts,Binomial[n,k]*Binomial[m+n-k,n]]/.k->0/.m->0/.n->0

ApplyOreOperator[uts,2^k*Binomial[m,k]*Binomial[n,k]]/.k->0/.m->0/.n->0

5.8.2  AiryAi
(complicated initial values)

TraditionalForm[HoldForm[AiryAi[x]^2==2^(2/3)/Pi*Integrate[AiryAi[2^(2/3)*(u^2+x)],{u,0,Infinity}]]]

lhs=Annihilator[AiryAi[x]^2,Der[x]]

rhs=Annihilator[Integrate[AiryAi[2^(2/3)*(u^2+x)],{u,0,Infinity}],Der[x]]

Compare initial values:

uts=UnderTheStaircase[lhs]

initsLhs=ApplyOreOperator[uts,AiryAi[x]^2]/.x->0

ApplyOreOperator[uts,2^(2/3)/Pi*AiryAi[2^(2/3)*(u^2+x)]]/.x->0

We did not succeed in proving symbolically with Mathematica that the initial values match (it does not evaluate the integrals over the above expressions).
 At least, let's do a numerical check:

NIntegrate[%,{u,0,Infinity}]

N[initsLhs]

Non-holonomic examples

The first six examples in this section are taken from a talk given by Manuel Kauers on the "Non-holonomic systems approach to Special Function Identities"

Stirling summation

HoldForm[Sum[Binomial[k,m]*Pochhammer[x-(k-m)+1,k-m]*StirlingS2[n,k],{k,0,n}]==Sum[Binomial[n,k]*x^(n-k)*StirlingS2[k,m],{k,0,n}]]//TraditionalForm

ann1=Annihilator[Binomial[k,m]*Pochhammer[x-(k-m)+1,k-m]*StirlingS2[n,k],{S[k],S[m],S[n]}]

ann2=Annihilator[Binomial[n,k]*x^(n-k)*StirlingS2[k,m],{S[k],S[m],S[n]}]

gb=OreGroebnerBasis[ann1,OreAlgebra[k,S[k],S[m],S[n]],MonomialOrder->EliminationOrder[1]]

lhs1=OrePolynomialSubstitute[First[gb],{S[k]->1}]

lhs2=Takayama[ann1,{k}]

lhs3=First[CreativeTelescoping[ann1,S[k]-1,{S[m],S[n]}]]

lhs4=Annihilator[Sum[Binomial[k,m]*Pochhammer[x-(k-m)+1,k-m]*StirlingS2[n,k],{k,0,n}],{S[m],S[n]},Assumptions->Element[n,Integers]&&n≥0]

OreReduce[lhs1,lhs2]

gb=OreGroebnerBasis[ann2,OreAlgebra[k,S[k],S[m],S[n]],MonomialOrder->EliminationOrder[1]]

rhs1=OrePolynomialSubstitute[First[gb],{S[k]->1}]

rhs2=Takayama[ann2,{k}]

rhs3=First[CreativeTelescoping[ann2,S[k]-1,{S[m],S[n]}]]

Mathematica does not simplify the StirlingS2 recurrence, hence we give it in the Assumptions in order to obtain a nice result:

rhs4=Annihilator[Sum[Binomial[n,k]*x^(n-k)*StirlingS2[k,m],{k,0,n}],{S[m],S[n]},Assumptions->-StirlingS2[n,m]+(-1-m) StirlingS2[n,1+m]+StirlingS2[1+n,1+m]==0]

Abel-type sum (1)

HoldForm[Sum[Binomial[n,k]*(k+r)^(k-1)*(s-k)^(n-k),{k,0,n}]==(r+s)^n/r]//TraditionalForm

ann=Annihilator[Binomial[n,k]*(k+r)^(k-1)*(s-k)^(n-k),{S[k],S[n],S[r],S[s]}]

OreGroebnerBasis[ann,OreAlgebra[k,S[k],S[n],S[r],S[s]],MonomialOrder->EliminationOrder[1]]

lhs1=OrePolynomialSubstitute[First[%],{S[k]->1}]

lhs2=Takayama[ann,{k}]

Takayama[ann,{k},Extended->True]

lhs3=First[CreativeTelescoping[ann,S[k]-1,{S[n],S[r],S[s]}]]

rhs=Annihilator[(r+s)^n/r,{S[n],S[r],S[s]}]

OreReduce[{lhs1,lhs2,lhs3},rhs]

Abel-type sum (2)

HoldForm[Sum[(-1)^k*Binomial[n,k]*(x-k)^n,{k,0,n}]==n!]//TraditionalForm

ann=Annihilator[(-1)^k*Binomial[n,k]*(x-k)^n,{S[k],S[n],S[x]}]

OreGroebnerBasis[ann,OreAlgebra[k,S[k],S[n],S[x]],MonomialOrder->EliminationOrder[1]]

lhs1=OrePolynomialSubstitute[First[%],{S[k]->1}]

lhs2=Takayama[ann,{k}]

lhs3=Annihilator[Sum[(-1)^k*Binomial[n,k]*(x-k)^n,{k,0,n}],{S[n],S[x]}]

FullSimplify[ApplyOreOperator[lhs1,n!]]

OreReduce[lhs1,Annihilator[n!,{S[n],S[x]}]]

Stirling number identity

HoldForm[Sum[Binomial[n,k]*StirlingS2[k,l]*StirlingS2[n-k,m],{k,0,n}]==Binomial[l+m,l]*StirlingS2[n,l+m]]//TraditionalForm

ann=Annihilator[Binomial[n,k]*StirlingS2[k,l]*StirlingS2[n-k,m],{S[k],S[l],S[m],S[n]}]

OreGroebnerBasis[ann,OreAlgebra[k,S[k],S[l],S[m],S[n]],MonomialOrder->EliminationOrder[1]]

lhs1=OrePolynomialSubstitute[First[%],{S[k]->1}]

lhs2=Takayama[ann,{k}]

lhs3=First[CreativeTelescoping[ann,S[k]-1,{S[l],S[m],S[n]}]]

rhs=Annihilator[Binomial[l+m,l]*StirlingS2[n,l+m],{S[l],S[m],S[n]}]

OreReduce[{lhs1,lhs2,lhs3},rhs]

Integral representation of StirlingS2

TraditionalForm[HoldForm[StirlingS2[n,k]==n!/(2Pi I k!)ContourIntegral[(Exp[z]-1)^k/z^(n+1),Abs[z]==1]]]

ann=Annihilator[(Exp[z]-1)^k/z^(n+1)*n!/k!,{Der[z],S[k],S[n]}]

OrePolynomialSubstitute[First[OreGroebnerBasis[ann,OreAlgebra[z,Der[z],S[k],S[n]],MonomialOrder->EliminationOrder[1]]],{Der[z]->0}]

Takayama[ann,{z}]

CreativeTelescoping[ann,Der[z],{S[k],S[n]}]

Annihilator[StirlingS2[n,k],{S[k],S[n]}]

Mixed type

HoldForm[Sum[StirlingS2[k,m]*(m+1)^(n-k),{k,0,n}]==StirlingS2[n+1,m+1]]//TraditionalForm

annSmnd=Annihilator[StirlingS2[k,m]*(m+1)^(n-k),{S[k],S[m],S[n]}]

AnnihilatorDimension[annSmnd]

Since the dimension is 2, theory does not guarantee a telescoping relation. However, by inspection we see that the only relation the we have found for the summand is already k-free, hence we can do the sum by hand.

A correction term is introduced by the action of S[n] on the summation bound.

corr=Annihilator[Evaluate[StirlingS2[k,m]*(m+1)^(n-k)/.n->n+1/.k->n+1],{S[m],S[n]}]

lhs=First[corr]**First[annSmnd]

rhs=Annihilator[StirlingS2[n+1,m+1],{S[m],S[n]}]

OreReduce[lhs,rhs]

Abel's identity

HoldForm[Sum[Binomial[n,k]*r*(k+r)^(k-1)*(n-k+s)^(n-k),{k,0,n}]==(n+r+s)^n]//TraditionalForm

ann=Annihilator[Binomial[n,k]*r*(k+r)^(k-1)*(n-k+s)^(n-k),{S[k],S[n],S[r],S[s]}]

gb=OreGroebnerBasis[ann,OreAlgebra[k,S[k],S[n],S[r],S[s]],MonomialOrder->EliminationOrder[1]]

LeadingPowerProduct/@gb

lhs1=OrePolynomialSubstitute[gb[[1]],{S[k]->1}]

lhs2=Takayama[ann,{k}]

lhs3=First[CreativeTelescoping[Binomial[n,k]*r*(k+r)^(k-1)*(n-k+s)^(n-k),S[k]-1,{S[n],S[r],S[s]}]]

rhs=Annihilator[(n+r+s)^n,{S[n],S[r],S[s]}]

OreReduce[{lhs1,lhs2,lhs3},rhs]

Zeta and Beta

This example can be found in "Espinosa/Moll:  On some definite integrals involving the Hurwitz zeta function. Part 1.".
It also appeared in Bruno Salvy's FPSAC'09 talk.

TraditionalForm[HoldForm[Integrate[x^(k-1) Zeta[n,a+b x],{x,0,Infinity}]==b^(-k)Beta[k,n-k] Zeta[n-k,a]]]

Annihilator[Integrate[x^(k-1) Zeta[n,a+b x],{x,0,Infinity}],{S[a],Der[a],Der[b],S[k],S[n]}]

In Mathematica 7, the following limit evaluated to 0:

Limit[x^k*Zeta[n,a+b*x],x->Infinity,Assumptions->a>0&&b>0&&k>0&&n>k+1&&Element[{k,n},Integers]]

Annihilator[Integrate[x^(k-1) Zeta[n,a+b x],{x,0,Infinity}],{S[a],Der[a],Der[b],S[k],S[n]},Assumptions->a>0&&b>0&&k>0&&n>k+1&&Element[{k,n},Integers]]

So we cheat a little bit:

Annihilator[Integrate[x^(k-1) Zeta[n,a+b x],{x,0,Infinity}],{S[a],Der[a],Der[b],S[k],S[n]},Assumptions->k>0&&Zeta[_,_]==0]

Annihilator[b^(-k)Beta[k,n-k] Zeta[n-k,a],{S[a],Der[a],Der[b],S[k],S[n]}]

GBEqual[%,%%]

Bernoulli polynomials (∫_0^1)(BernoulliB2[m, x]) (BernoulliB2[n, x])dx==(((((-1))^(n-1)) m! n! (Subscript[B, m+n]))/((m+n)!))=(((((-1))^m) m! n! cos((1/2) π (m+n)) (Zeta[m+n]))/((2^(m+n-1)) (π^(m+n))))

TraditionalForm[HoldForm[Integrate[BernoulliB[n,x]*BernoulliB[m,x],{x,0,1}]==(-1)^(n-1)*m!*n!/(m+n)!*BernoulliB[n+m]]]

TraditionalForm[HoldForm[Integrate[BernoulliB[m,x]*BernoulliB[n,x],{x,0,1}]==(-1)^m*m!*n!/2^(m+n-1)*Cos[Pi*(m+n)/2]*Zeta[m+n]/Pi^(m+n)]]

Theory does not guarantee that a relation for the integral exists, since both factors of the integrand have dimension 1. Looking at its annihilator, it is easily seen that it delivers a creative telescoping relation:

Annihilator[BernoulliB[n,x]*BernoulliB[m,x],{S[m],S[n],Der[x]}]

delta=BernoulliB[m+1,x]*BernoulliB[n+1,x]

(delta/.x->1)-(delta/.x->0)

inh=FullSimplify[%]

Hence we get the following inhomogeneous recurrence:

rec=(n+1)*Int[m+1,n]+(m+1)*Int[m,n+1]-inh==0

Table[rec,{m,0,3},{n,0,3}]/.Int[m_,n_]:>Integrate[BernoulliB[m,x]*BernoulliB[n,x],{x,0,1}]

However, if m and n are not zero, the recurrence simplifies to

homrec=FullSimplify[rec,m>0&&n>0]

This recurrence, together with the observation that the integral is 0 for either m or n being 0 (but not both at the same time), immediately implies the closed form (-1)^(n + 1)*m!*n!/(m + n)!*BernoulliB[m + n] for m>0 and n>0.

Table[Which[m==n==0,1,m==0||n==0,0,True,(-1)^(n+1)*m!*n!/(m+n)!*BernoulliB[m+n]],{m,0,6},{n,0,6}]

Table[Int[m,n],{m,0,6},{n,0,6}]/.Int[m_,n_]:>Integrate[BernoulliB[m,x]*BernoulliB[n,x],{x,0,1}]

The second evaluation follows from rewriting the Bernoulli numbers in terms of Zeta values.

Table[Which[m==n==0,1,m==0||n==0,0,True,(-1)^m*m!*n!/2^(m+n-1)*Cos[Pi*(m+n)/2]*Zeta[m+n]/Pi^(m+n)],{m,0,6},{n,0,6}]

Bernoulli polynomials (∑_(k=1)^n)(({{n}, {k}})) (Subscript[B, n-k])(x)==n (x^(n-1))

TraditionalForm[HoldForm[Sum[Binomial[n,k]*BernoulliB[n-k,x],{k,1,n}]==n*x^(n-1)]]

Annihilator[Sum[Binomial[n,k]*BernoulliB[n-k,x],{k,1,n}],{S[n],Der[x]}]

This error can be avoided by splitting off the last summand.

BernoulliB[0,x]

Annihilator[1+Sum[Binomial[n,k]*BernoulliB[n-k,x],{k,1,n-1}],{S[n],Der[x]}]

ApplyOreOperator[First[%],n x^(n-1)]//Together

(∫_0^1)(x^n) (log^k)(x)dx==(((-1))^k) (((1+n))^(-k-1)) (Gamma[1+k])

This example was communicated to us by Clemens Raab.

TraditionalForm[HoldForm[Integrate[x^n Log[x]^k,{x,0,1}]==(-1)^k*(1+n)^(-k-1)*Gamma[1+k]]]

Integrate[x^n Log[x]^k,{x,0,1},Assumptions->n≥0&&k≥0]

Annihilator[Integrate[x^n Log[x]^k,{x,0,1}],{S[k],S[n]},Assumptions->n≥0&&k≥0]

Annihilator[(-1)^k*(1+n)^(-k-1)*Gamma[1+k],S[k]]

Compare initial values for k=0 and arbitrary n>=0:

Integrate[x^n,{x,0,1},Assumptions->n≥0]

(-1)^k*(1+n)^(-k-1)*Gamma[1+k]/.k->0