This Mathematica notebook accompanies the paper 
"The integrals in Gradshteyn and Ryzhik. Part 18: Some automatic proofs"
by Christoph Koutschan and Victor Moll.

<<HolonomicFunctions.m

HolonomicFunctions package by Christoph Koutschan, RISC-Linz, Version 1.3 (25.01.2010)
⟶ Type  ?HolonomicFunctions  for help

The integrals in Gradshteyn and Ryzhik.
Part 18: Some automatic proofs

Section 3  (2.148.3)   ∫(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]]

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

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

inh=First[Annihilator[x/(1+x^2)^n,S[n]]]

Factor[inh**((2n)**S[n]-(2n-1))]

?Inhomogeneous

Annihilator[Integrate[1/(1+x^2)^n,x],S[n],Inhomogeneous->True]

Section 4  (9.181.1)  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[α,m+n]*Pochhammer[β,m]*Pochhammer[b,n]/(Pochhammer[γ,m+n]*m!*n!)*x^m*y^n,{S[m],S[n],Der[x],Der[y]}]

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

OreReduce[b*α+(b*x)**Der[x]+(-y+y^2)**Der[y]^2+(y+b*y+y*α-γ)**Der[y]+(-x+x*y)**Der[x]**Der[y],pde]

OreReduce[α*β+(-x+x^2)**Der[x]^2+(y*β)**Der[y]+(x+x*α+x*β-γ)**Der[x]+(-y+x*y)**Der[y]**Der[x],pde]

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

FindRelation[pde,Eliminate->β]

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

Section 5  (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])]]

Integrate[ChebyshevT[n,1-x^2 y]/Sqrt[1-x^2],{x,-1,1},Assumptions->Element[n,Integers]&&n≥0]

Annihilator[ChebyshevT[n,1-x^2 y]/Sqrt[1-x^2],{S[n],Der[x],Der[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]}]

Simplify[ApplyOreOperator[Last[%],ChebyshevT[n,1-x^2 y]/Sqrt[1-x^2]]]

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

(Limit[#,x->1]-Limit[#,x->-1])&/@%

... 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

Section 6  (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 *)

Integrate[x^(r-1)(1-x)^(s-1)Hypergeometric2F1[a,b,c,x],{x,0,1},Assumptions->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]}]

smnd=Pochhammer[a,k]*Pochhammer[b,k]*Pochhammer[r,k]/Pochhammer[c,k]/Pochhammer[r+s,k]/k!;

tak=Takayama[Annihilator[smnd,{S[k],S[a],S[b],S[c],S[r],S[s]}],{k}]

rhs2=DFiniteTimes[Annihilator[Gamma[r]*Gamma[s]/Gamma[r+s],{S[a],S[b],S[c],S[r],S[s]}],tak]

GBEqual[rhs,rhs2]

GBEqual[lhs,rhs]

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

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

UnderTheStaircase[lhs]

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

Section 7  (7.322)  (∫_0^(2 a))(((x (2 a-x)))^(ν-(1/2))) (Subscript[C, n]^((ν)))((x/a)-1) (e^(-b x))dx==(((((-1))^n) π Γ(n+2 ν) ((((a/(2 b))))^ν) (e^(-a b)) (Subscript[I, n+ν])(a b))/(n! Γ(ν)))

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

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

Integrate[(x(2a-x))^(ν-1/2)*GegenbauerC[n,ν,x/a-1]*E^(-b*x),{x,0,2a},Assumptions->Re[ν]>0&&Element[n,Integers]&&n≥0&&a≥0]

The 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))^(ν-1/2)*GegenbauerC[n,ν,x/a-1]*E^(-b*x),{x,0,2a}],{Der[a],Der[b],S[n],S[ν]},Assumptions->ν≥1,Inhomogeneous->True];

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

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

GBEqual[lhs,rhs]

UnderTheStaircase[lhs]

(x(2a-x))^(ν-1/2)*GegenbauerC[n,ν,x/a-1]*E^(-b*x)/.{a->1,b->1,n->1}/.{{ν->0},{ν->1}}

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

(-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},{nu->1}}

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

Section 8  (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)])))

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

Integrate[BesselJ[m,a*x]*BesselJ[n,b*x],{x,0,Infinity},Assumptions->a>0&&b>0&&b<a&&m>0&&n>0&&Element[{m,n},Integers]]

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

delta

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.

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:

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[ann,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.

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]

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

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

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:

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

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[lhs,rhs]

Initial values:

uts=UnderTheStaircase[rhs]

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

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[%-%%]

Section 9  (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=Simplify[rexpr/.b->s^2,s>0];

rhs1=Annihilator[rexpr1,{S[m],Der[a],Der[s],Der[c]}]

rhs1=DFiniteSubstitute[rhs1,{s->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}]

Section 10  (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]]

LeadingPowerProduct/@Flatten[{lhs,rhs}]

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

{{op},{inh}}=Annihilator[Integrate[ArcTan[p*x]/(1+p^2*x),{x,0,1}],Der[p],Inhomogeneous->True]

FullSimplify[inh]

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[%]**op

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[%]