<<HolonomicFunctions.m

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

Define the coefficient sequence of the series under consideration

ClearAll[F1,F2,F3,F4];
F1[n_Integer]:=If[n==0,1/x/y,F1[n]=Expand[((x+y+1/x+1/y)#-Coefficient[#,x,-1]-Coefficient[#,y,-1])&[F1[n-1]]]];(*checked*)F2[n_Integer]:=If[n==0,0,F2[n]=Expand[((x+y+1/x+1/y)#+Coefficient[F1[n-1],y,-1]-Coefficient[#,x,-1]-1/y Coefficient[#,y,0])&[F2[n-1]]]];(*checked*)F3[n_Integer]:=F2[n]/.{x->y,y->x};
F4[n_Integer]:=If[n==0,0,F4[n]=Expand[((x+y+1/x+1/y)#+Coefficient[F2[n-1],x,-1]+Coefficient[F3[n-1],y,-1]-1/x Coefficient[#,x,0]-1/y Coefficient[#,y,0])&[F4[n-1]]]];(*checked*)

Table[F1[n]+F2[n]+F3[n]+F4[n],{n,0,5}]

PP[expr_,x_] := Expand[x*(Expand[1/x expr]/.(x^i_/;i<0)->0)]; (* positive part extraction tool *)

s=x+y+1/x+1/y;rat=(x y-x/y-y/x+1/x/y)/(1 - t s);

Derive a differential equation for F1(1,1,t)

F1== 1/x/y HoldForm[PP][x,HoldForm[PP][y,rat]]==1/x/y RES[x,RES[y,1/x/y rat x/(1-x) y/(1-y)]]//TraditionalForm

deqF1=First[CreativeTelescoping[First[CreativeTelescoping[rat/(1-x)/(1-y), Der[x],{Der[y],Der[t]}]],Der[y],{Der[t]}]]

Derive a differential equation for F2(1,1,t)=F3(1,1,t)

F2 == t/y NP[x, RES[y, (x y - x/y-y/x+1/x/y)/(1 - t (x+y+1/x+1/y))] HoldForm[PP][y, (y-1/y)/(1 -t (x+y+1/x+1/y))]] // TraditionalForm

deqA1=First[CreativeTelescoping[rat,Der[y],{Der[t]}]]; (* differential equation for the residue *)
deqA2=First[CreativeTelescoping[Y/y^2/(1-Y/y)(y-1/y)/(1-t s),Der[y],{Der[t]}]] /. Y -> y; (* differential equation for the positive part *)
deqA = DFiniteTimes[deqA1,deqA2]; (* differential euation for the product of residue and positive part *)

aeqA=(-1+x)^2*x*(1+x)^2*y^2+(-1+x)*(1+x)*y*(t^2-2*t*x+x^2-2*t*x^3+t^2*x^4+t^2*x*y-t*x^2*y+t^2*x^3*y)*A+t^3*x^2*(t*x+t*y-x*y+t*x^2*y+t*x*y^2)*A^2;(* guessed algebraic equation for the same quantity *)
Together[OreReduce[First[deqA],Annihilator[Root[Function[A,Evaluate[aeqA]],1],{Der[t]}]]] (* algebraic function satisfies differential equation derived above if this returns zero *)

aeqAsol=Sum[a[n]t^n,{n,0,10}]/.First[Solve[CoefficientList[aeqA/.A->(Sum[a[n]t^n,{n,0,10}]+O[t]^11),t]==0,Table[a[n],{n,0,11}]]]; (* first terms of algebraic series solution of aeqA *)
truesol=Normal[Coefficient[Normal[Series[rat,{t,0,10}]],y,-1]PP[Normal[Series[(y-1/y)/(1-t s),{t,0,10}]],y]+O[t]^11]; (* first terms of F2(1,1,t) *)
aeqAsol-truesol//Together (* if this returns zero, F2(1,1,t) agrees with the algebraic function defined by aeqA *)

aux=Collect[Numerator[Together[Evaluate[aeqA/.A->A/(t/(1-x))/.y->1]]],A,Factor];
deqF2=First[CreativeTelescoping[Annihilator[Root[Function[A,Evaluate[aux]],1],{Der[x],Der[t]}],Der[x],{Der[t]}]] (* differential equation for F2(1,1,t) *)

Together[ApplyOreOperator[First[deqF2],(Sum[F2[n]t^n,{n,0,10}]+O[t]^11)/.{x->1,y->1}]] (* cross check, should give zero *)

Derive a Differential Equation for F4

F4==1/x/y t^2 RES[y,1/y^2/(1-1/y) RES[x, (y-1/y)/(1-t s) RES[y, rat]] HoldForm[PP][x,(x-1/x)/(1- t s)]]//TraditionalForm

deqf=First[CreativeTelescoping[rat,Der[y],{Der[x],Der[t]}]]; (* differential equations for innermost residue, as a function of x and t *)
deqf=OreGroebnerBasis[ChangeOreAlgebra[#,OreAlgebra[Der[x],Der[t],Der[y]]]&/@Append[deqf,ToOrePolynomial[Der[y]]]]; (* view the innermost residue also as a (constant) function with respect to y *)
deqf=First[CreativeTelescoping[DFiniteTimesHyper[deqf,(y-1/y)/(1-t s)],Der[x],{Der[y],Der[t]}]]; (* differential equations for the middle residue; this takes some time *)

The system of differential equations derived above is not minimal. Via guessing, we can find that the middle residue satisfies some additional differential equations that are not consequences of the derived equations (i.e.: not elements of the ideal generated by the corresponding operators). Using closure properties and comparison of initial values, we can certify that these guessed operators are correct, so we are entitled to add them to the ideal. Doing so will speed up some of the subsequent computations.

deqfguess=ToOrePolynomial[#,OreAlgebra[Der[y],Der[t]]]&/@(* additional guessed operators *){-2*t-3*y+96*t^2*y+18*t*y^2-18*y^3-192*t^2*y^3+18*t*y^4-3*y^5+96*t^2*y^5-2*t*y^6+t*(-1+y)^2*(1+y)^2*(-2*t-5*y+96*t^2*y-2*t*y^2)*Der[t]+t^2*(-1+4*t)*(1+4*t)*(-1+y)^2*y*(1+y)^2*Der[t]^2-(-1+y)*(1+y)*(2*t-3*y+2*t*y^2)*(4*t+3*y-8*t*y^2+y^3+4*t*y^4)*Der[y]-2*t^2*(-1+y)^3*(1+y)^3*(2*t-3*y+2*t*y^2)*Der[t]*Der[y]+(-1+y)^2*y^2*(1+y)^2*(2*t-3*y+2*t*y^2)*Der[y]^2,t^2-12*t*y+16*t^3*y+109*t^2*y^2-48*t*y^3-64*t^3*y^3+24*y^4-110*t^2*y^4-72*t*y^5+96*t^3*y^5+24*y^6-110*t^2*y^6-48*t*y^7-64*t^3*y^7+109*t^2*y^8-12*t*y^9+16*t^3*y^9+t^2*y^10+(t^3-4*t^2*y+8*t^4*y+37*t^3*y^2-4*t^2*y^3-32*t^4*y^3-38*t^3*y^4+16*t^2*y^5+48*t^4*y^5-38*t^3*y^6-4*t^2*y^7-32*t^4*y^7+37*t^3*y^8-4*t^2*y^9+8*t^4*y^9+t^3*y^10)*Der[t]+(-4*t^3+19*t^2*y-24*t*y^2+16*t^3*y^2+6*y^3+23*t^2*y^3-12*t*y^4-20*t^3*y^4+12*y^5-170*t^2*y^5+12*t*y^6-18*y^7+182*t^2*y^7+12*t*y^8+20*t^3*y^8-41*t^2*y^9+12*t*y^10-16*t^3*y^10-13*t^2*y^11+4*t^3*y^12)*Der[y]+(-2*t^4+7*t^3*y-5*t^2*y^2+8*t^4*y^2+5*t^3*y^3+10*t^2*y^4-10*t^4*y^4-50*t^3*y^5+50*t^3*y^7-10*t^2*y^8+10*t^4*y^8-5*t^3*y^9+5*t^2*y^10-8*t^4*y^10-7*t^3*y^11+2*t^4*y^12)*Der[t]*Der[y]+(6*t^2*y^2-12*t*y^3+6*y^4-24*t^2*y^4+12*t*y^5-12*y^6+36*t^2*y^6+12*t*y^7+6*y^8-24*t^2*y^8-12*t*y^9+6*t^2*y^10)*Der[y]^2+(t^2*y^3-2*t*y^4+y^5-5*t^2*y^5+4*t*y^6-3*y^7+10*t^2*y^7+3*y^9-10*t^2*y^9-4*t*y^10-y^11+5*t^2*y^11+2*t*y^12-t^2*y^13)*Der[y]^3};
UnderTheStaircase[DFinitePlus@@(DFiniteOreAction[deqf,#]&/@deqfguess)] (* initial values that have to be checked in order to confirm the correctness of the guessed operators. Since these terms were used for the guessing, the initial values match by construction so we do not perform a check here *)

deqf=OreGroebnerBasis[Join[deqf,deqfguess]];(* include the guessed equations into the system of differential equations for the middle residue *)

deqg=First[CreativeTelescoping[Annihilator[1/x^2/(1-1/x)(x-1/x)/(1-t s),{Der[y],Der[x],Der[t]}],Der[x],{Der[y],Der[t]}]]; (* differential equations for the positive part *)

deqF4=First[FindCreativeTelescoping[DFiniteTimesHyper[DFiniteTimes[deqf,deqg],t^2/y^2/(1-1/y)],{Der[y]},{Der[t]}]](* differential equation for the outer residue; this takes some time *)

ApplyOreOperator[First[deqF4],Sum[F4[n]t^n,{n,0,15}]+O[t]^15]/.{x->1,y->1} (* cross check,should give zero *)

Combine the derived differential operators to one for F.

deqF=DFinitePlus[deqF1,deqF2,deqF4]; (* certified operator; it has order 8 *)

deqFguessed=ToOrePolynomial[60+408 t+672 (t^2)+384 (t^3)+(-15+6 t+600 (t^2)+1584 (t^3)+1152 (t^4)) Der[t]+(-9 t-24 (t^2)+180 (t^3)+672 (t^4)+576 (t^5)) ((Der[t])^2)+(-(t^2)-4 (t^3)+12 (t^4)+64 (t^5)+64 (t^6)) ((Der[t])^3)]; (* guessed shorter operator *)

OreReduce[deqF,{deqFguessed}] (* if this gives zero, the guessed operator is correct *)

recF=DFiniteDE2RE[{deqFguessed},{t},{n}]; (* certified recurrence operator for the coefficients of F(1,1,t); it has order 4 *)

recFguessed = ToOrePolynomial[-192 - 736*n - 992*n^2 - 592*n^3 - 160*n^4 - 16*n^5 + (216 + 24*n - 124*n^2 - 60*n^3 - 8*n^4)*S[n] + (-48 + 52*n + 124*n^2 + 67*n^3 + 14*n^4 + n^5)*S[n]^2];

OreReduce[recF,{recFguessed}] (* if this is zero, the guessed operator is correct *)

<<Asymptotics.m

Asymptotics[ApplyOreOperator[recFguessed,f[n]],f[n]] (* negative integer exponents of n prove transcendence of generating functions *)