Simplification with recurrence relations in Mathematica

Derivations and implementation

Oleksandr Pavlyk

Wolfram Research, Inc.

Univariate recurrence relations

We shall only consider fixed dispersion recurrence relations (fixed value of p), i.e.

c[n] == (Subscript[α, 1])c[n-1]+(Subscript[α, 2])c[n-2]+…+(Subscript[α, p])c[n-p]

An example

Many of these come from differential equations directly:

ToDifferentialEquation[HermiteH[n+1,x],y]==0

DSolve[%,y,x]

(%%/.%/.C[2]->0)//Simplify

Proving the recurrence equation thus implies proving the differential equation:

FullSimplify[%]

Restriction due to the assumption

This restriction rules out recurrence equations satisfied by, say, Bernoulli numbers, where p==n:

BernoulliB[n]==-(1/(((2^n)-1))) (∑_(k=0)^(n-1))((2^(k-1))-1) Binomial[n,k] BernoulliB[k]/; n>1

Further examples

ChebyshevU[n-1,z]+ChebyshevU[n+1,z]==2z ChebyshevU[n,z]//FullSimplify

ChebyshevT[n-1,z]+ChebyshevT[n+1,z]==2z ChebyshevT[n,z]//FullSimplify

n LaguerreL[n-1,z]+(n+1)LaguerreL[n+1,z]==(2n+1-z)LaguerreL[n,z]//FullSimplify

n LegendreP[n-1,z]+(n+1) LegendreP[n+1,z]==(2 n+1) z LegendreP[n,z]//FullSimplify

Simplification in univariate case is a matter of applying recurrence shifts to reduce all congruent instances to the same base.

n LegendreP[n-1,z]+(n+1) LegendreP[n+1,z]==(2 n+1) z LegendreP[n,z]/.{
LegendreP[m_,z]/;(m-n>0):>((2m-1)z LegendreP[m-1,z]-(m-1) LegendreP[m-2,z])/m
}//Simplify

Multivariate recurrence relations

Consider bi-parametric orthogonal polynomial (Subscript[L, n]^λ)(z):

ToDifferentialEquation[LaguerreL[n+2,λ,x],y]

%/.y->Function[x, LaguerreL[n+2,λ,x] ]

How does one prove this ?

Using a method similar to the the univariate case could work to prove the differential equation using simultaneous shifts of both n and λ, however ...

Multivariate recurrence relations

...  one should be also able to prove other classic recurrence equations, known for LaguerreL:

(n+λ) LaguerreL[n-1,λ,z]+(n+1) LaguerreL[n+1,λ,z]==(2 n+λ-z+1) LaguerreL[n,λ,z]

(n+λ) LaguerreL[n,λ-1,z]+z LaguerreL[n,λ+1,z]==(z+λ) LaguerreL[n,λ,z]

Since these seem like more basic, let us see if we can use them to prove the differential equation:

x LaguerreL[n,λ+2,x]+(x-1-λ) LaguerreL[n+1,λ+1,x]+(n+2) LaguerreL[n+2,λ,x]//.{
LaguerreL[ν_,μ_,z_]/;ν-n-1>0 :> ((2ν+μ-z-1)LaguerreL[ν-1,μ,z]-(ν+μ-1)LaguerreL[ν-2,μ,z])/ν,
LaguerreL[ν_,μ_,z_]/;μ-λ-1>0:> ((z+μ-1)LaguerreL[ν,μ-1,z]-(ν+μ-1)LaguerreL[ν,μ-2,z])/z
}//Simplify

Multivariate recurrence relations

Let us try to analyze what happens:

Manipulate[GraphicsRow[{
DifferenceGraphics[{{2,0},{1,1},{0,2}},{-1,3}],
DifferenceGraphics[{{1,0},{0,0},{1,1},{0,2}},{-1,3}],
DifferenceGraphics[{{1,0},{0,0},{1,1},{0,1}},{-1,3}]
}[[1;;k]],ImageSize->{200 k,200}],{{k,1,""},{1->"Diff.Eq.",2->"Shift by n", 3->"Shift by λ"},Appearance->"Row"},SaveDefinitions->True]

Multivariate recurrence relations

There must be an additional recurrence relations between these four LaguerreL functions.

((Subscript[α, 1])+(Subscript[α, 2])x+…) LaguerreL[n,λ,x]+((Subscript[β, 1])+(Subscript[β, 2])x+…) LaguerreL[n+1,λ,x]+((Subscript[γ, 1])+(Subscript[γ, 2])x+…) LaguerreL[n,λ+1,x]+((Subscript[δ, 1])+(Subscript[δ, 2])x+…)LaguerreL[n+1,λ+1,x]==0

Use series representaiton for LaguerreL polynomial:

c[n_,λ_,k_]:=((Gamma[λ+n+1] )/(Gamma[n+1]))((Pochhammer[-n,k])/(Gamma[λ+k+1] k!))

LaguerreL[n,λ,x]==(∑_(k=0)^n)c[n,λ,k](x^k)

Determining coefficients (Subscript[α, i]), (Subscript[β, i]),  (Subscript[γ, i]),  (Subscript[δ, i])

p=1;vars=Join[Table[(Subscript[α, k]),{k,p}],Table[(Subscript[β, k]),{k,p}],Table[(Subscript[γ, k]),{k,p}],Table[(Subscript[δ, k]),{k,p}]];

num=(1/(c[n,λ,k]))((∑_j^p)(Subscript[α, j]) c[n,λ,k-j+1]+(∑_j^p)(Subscript[β, j]) c[n+1,λ,k-j+1]+(∑_j^p)(Subscript[γ, j]) c[n,λ+1,k-j+1]+(∑_j^p)(Subscript[δ, j]) c[n+1,λ+1,k-j+1])//FunctionExpand//Together//Numerator;

sol=Solve[CoefficientList[num,k]==0,vars]//Simplify

(((∑_j^p)(Subscript[α, j]) (x^(j-1)))LaguerreL[n,λ,x]+((∑_j^p)(Subscript[β, j]) (x^(j-1)))LaguerreL[n+1,λ,x]+((∑_j^p)(Subscript[γ, j]) (x^(j-1)))LaguerreL[n,λ+1,x]+((∑_j^p)(Subscript[δ, j]) (x^(j-1)))LaguerreL[n+1,λ+1,x])/.sol

Outer[D,%,vars]//Flatten//Union//Simplify

The answer

sol1={0,LaguerreL[n,1+λ,x]+LaguerreL[1+n,λ,x]-LaguerreL[1+n,1+λ,x]};

sol2={0,LaguerreL[n,1+λ,x]+LaguerreL[1+n,λ,x]-LaguerreL[1+n,1+λ,x],x (LaguerreL[n,1+λ,x]+LaguerreL[1+n,λ,x]-LaguerreL[1+n,1+λ,x]),(((1+n+λ) LaguerreL[n,λ,x]+(-1-n+x) LaguerreL[1+n,λ,x]-x LaguerreL[1+n,1+λ,x])/(1+n+λ))};

Join[sol1,sol2]//KeepFactorsDependingOn[#,LaguerreL]&//Union//Simplify

We can prove this

FullSimplify[%]

Notice, that despite LaguerreL being bi-parametric, there just two independent parameters:

Manipulate[
Graphics[{PointSize[0.17],{col[[1]],Point[{{0,0}}]},{col[[2]],Point[{{1,1}}]},{col[[3]],Point[{{0,1}}]},{col[[4]],Point[{{1,0}}]},
{GrayLevel[0.75],Point[{{2,0},{2,1},{2,2},{1,2},{0,2},{-1,0},{-1,-1},{-1,1},{-1,2},{-1,3},{0,3},{0,-1},{1,3},{1,-1},{2,3},{2,-1},{3,0},{3,-1},{3,1},{3,2},{3,3}}]}
},GridLines->{Range[-1,3],Range[-1,3]},PlotRange->{{-1,3},{-1,3}},ImageSize->100],
{{col,{Red,Red,Gray,Gray},""},{
{Red,Red,Gray,Gray}->ArrayPlot[{{0,1},{1,0}},Mesh->True,ImageSize->20],
{Red,Gray,Red,Gray}->ArrayPlot[{{1,0},{1,0}},Mesh->True,ImageSize->20],
{Gray,Red,Red,Gray}->ArrayPlot[{{1,1},{0,0}},Mesh->True,ImageSize->20],
{Gray,Red,Gray,Red}->ArrayPlot[{{0,1},{0,1}},Mesh->True,ImageSize->20],
{Gray,Gray,Red,Red}->ArrayPlot[{{1,0},{0,1}},Mesh->True,ImageSize->20],
{Red,Gray,Gray,Red}->ArrayPlot[{{0,0},{1,1}},Mesh->True,ImageSize->20]
},ControlType->SetterBar,Appearance->"Row",ControlPlacement->Bottom},Alignment->Center
]

Machinery to find recurrence equations

Method of indeterminate coefficients is cumbersome.

For a large class of functions, namely hypergeometric functions, we can do better.

Shift operators

There two types of differential operators which shift a requested parameter:

(1/a)z D[#,z]+#&[ HypergeometricPFQ[{a,a2,a3},{b1,b2},z] ] == HypergeometricPFQ[{a+1,a2,a3},{b1,b2},z];

Series[%,{z,0,7}]//Simplify

(1/(b-1))z D[#,z]+#&[ HypergeometricPFQ[{a1,a2,a3},{b,b2},z] ] == HypergeometricPFQ[{a1,a2,a3},{b-1,b2},z];

Series[%,{z,0,7}]//Simplify

Differential equation

Hypergeometric function satisfies a simple differential equation:

(Underscript[(z (∂_z) +(Subscript[a, 1]))∘…∘(z (∂_z) +(Subscript[a, p])), (Underscript[⏟, p times])]) ∘f(z)    ==  (Underscript[(∂_z)∘(z (∂_z) +(Subscript[b, 1])-1)∘…∘(z (∂_z) +(Subscript[b, q])-1), (Underscript[⏟, q+1 times])])∘f(z)

Fold[z D[#1,z]+#2 #1&, f[z], {a1,a2,a3}]-D[Fold[z D[#1,z]+(#2-1) #1&,f[z],{b1,b2}],z]/.f->Function[z,HypergeometricPFQ[{a1,a2,a3},{b1,b2},z]]

FullSimplify[%]

Fold[z D[#1,z]+#2 #1&, f[z], {a1,a2,a3}]-D[Fold[z D[#1,z]+(#2-1) #1&,f[z],{b1,b2,b3}],z]/.f->Function[z,HypergeometricPFQ[{a1,a2,a3},{b1,b2,b3},z]]//FullSimplify

Recurrences

Now, if instead of default differentiation rule, we use identies noted above, we get different recurrence relations.

(Underscript[(z (∂_z) +(Subscript[a, 1]))∘…∘(z (∂_z) +(Subscript[a, p])), (Underscript[⏟, p times])]) ∘f(z)    ==  (Underscript[(∂_z)∘(z (∂_z) +(Subscript[b, 1])-1)∘…∘(z (∂_z) +(Subscript[b, q])-1), (Underscript[⏟, q+1 times])])∘f(z)

For  (Subscript[, 1])(Subscript[F, 2])(a;b,c;z), for example we rewrite the above differential equation as:

((Subscript[Op, 1]) +(Subscript[α, 1]))∘f(z)    == γ  ((Subscript[Op, 2])+(Subscript[α, 2]))∘((Subscript[Op, 3]) +(Subscript[α, 3]))∘((Subscript[Op, 3]) +(Subscript[α, 4]))∘f(z)

Simplify`GivePFQRecurrence[{a},{b,c},z,{1},{1,1,1}]

Manipulate[
pfqExpr=Collect[Simplify`GivePFQRecurrence[{a},{b,c},z,{p1},{p2,p3,p4}]//ReleaseHold//Factor//KeepFactorsDependingOn[#,HypergeometricPFQ]&,_HypergeometricPFQ,FullSimplify],{p1,1,3,None},{p2,1,3,None},{p3,1,3,None},{p4,1,3,None},
ExpressionCell[Row[{Grid[{{"(Subscript[Op, 1])",SetterBar[Dynamic[p1],{1->"(Subscript[A, 1])",2->"(Subscript[B, 1])",3->"(Subscript[B, 2])"},Appearance->"Row"],"            ","(Subscript[Op, 2])",SetterBar[Dynamic[p2],{1->"(Subscript[A, 1])",2->"(Subscript[B, 1])",3->"(Subscript[B, 2])"},Appearance->"Row"]},{"(Subscript[Op, 3])",SetterBar[Dynamic[p3],{1->"(Subscript[A, 1])",2->"(Subscript[B, 1])",3->"(Subscript[B, 2])"},Appearance->"Row"],"            ","(Subscript[Op, 4])",SetterBar[Dynamic[p4],{1->"(Subscript[A, 1])",2->"(Subscript[B, 1])",3->"(Subscript[B, 2])"},Appearance->"Row"]}}],"           ",Sequence@@If[Head[pfqExpr]===Plus,{
Dynamic@DifferenceGraphics3D[Cases[pfqExpr,HoldPattern[HypergeometricPFQ][w__]:>(Flatten[Most[{w}]]-{a,b,c}),Infinity]]}],"           ",Sequence@@If[Head[pfqExpr]===Plus,{
Style[Dynamic@Grid[Cases[pfqExpr,HoldPattern[HypergeometricPFQ][w__]:>(Flatten[Most[{w}]]-{a,b,c}),Infinity]],Larger]
},{}]}] 
],SaveDefinitions->True ]

Using recurrence equation in each of 3 independent directions, we can always reduce down to (3^3)==27 points. How many of these are independent ?

The answer seems to equal one less of the order of the differential equation, i.e. 3 in this case.

Finding the basis

Length[eqs=Union[Function[{p1,p2,p3,p4},Sequence@@DeleteCases[
Flatten[(ReleaseHold@Simplify`GivePFQRecurrence[{a},{b,c},z,{p1},{p2,p3,p4}])/.Flatten[Outer[Join,{{},{a->a+1},{a->a+2},{a->a+3}},{{},{b->b-1},{b->b-2},{b->b-3}},{{},{c->c-1},{c->c-2},{c->c-3}},1],2]],x_/;!FreeQ[x,HypergeometricPFQ[up_,lo_,z]/;Max[Abs[Join[up,lo]-{a,b,c}]]>2]]
] @@@Flatten[Outer[List,#,#,#,#]&[{1,2,3}],3]]]

Length[allInstances=Union[Cases[eqs,_HypergeometricPFQ,Infinity]]]

Graphics3D[{Red,PointSize[0.05],Point[allInstances/.HypergeometricPFQ[up_,lo_,z_]:>Join[up,lo]-{a,b,c}],Directive[Black,Opacity[0.3]],Polygon[{{2,0,0},{0,-2,0},{0,0,-2},{2,0,0}} ]  },FaceGrids->All,ImageSize->250]

baseFunctions={HypergeometricPFQ[{a},{b,c},z],HypergeometricPFQ[{a+1},{b,c},z],HypergeometricPFQ[{a+2},{b,c},z]};

rest=Complement[allInstances,baseFunctions];

sol1=Solve[Thread[eqs==0],rest]

Applications

Recurrence simplifiers are used in FullSimplify as we saw earlier.

Higher oder derivatives

Diferentiating Bessel function

D[BesselJ[ν,x],{x,4}]

Simplify`RecurrenceShift[%,TargetFunctions->BesselJ[ν,x]]

Differentiating orthogonal polynomial

D[Exp[-(x^2)]HermiteH[n,x],{x,5}]

Simplify`RecurrenceShift[%,TargetFunctions->HermiteH[n,x]]

Simplify`RecurrenceShift[%,TargetFunctions->HermiteH[n+5,x]]

Finding out a recurrence equation

d[n_]:=(((n!)^2)/((2n)!))LaguerreL[2n,n+(1/2),2]

ListLogPlot[d[Range[100]]]

What is the recurence equation satisfied by (Subscript[d, n]) ?

Simplify`InducedRecurrence[ Function[n,(((n!)^2)/((2n)!))], Function[n, {2n,n+(1/2),2}],LaguerreL,{{n},{n+1},{n+2}},y]//KeepFactorsDependingOn[#,y]&

Table[%,{n,10}]/.y->d

RecurrenceAsymptote[%%==0,y[n],{n,Infinity,2}]

ListLogPlot[d[#]/( (27/16)^#)&[Range[300]] ]

Applications

MathGroup post

Expansion of hypergeometric pFq near the pole

The question was taken from MathGroup post by Dimitris to find the limiting value at e == 0:

s[e_]:=HypergeometricPFQ[{(1/2)-e,(1/2)-e,(1/2)-e,1-e},{(1/2),1,1},1] Sin[e π]

Direct substitution does not work, because hypergeometric function diverges as e -> 0, because convergence condition 
		Total[{1/2-e,1/2-e,1/2-e,1-e}]-Total[{1,1,1/2}] ==-4e<0 
implies e>0.

s[0]

Try numerical extrapolation:

Table[s[(21-k)/2000`15],{k,20}]

Interpolation[%][21]//Rationalize

Symbolic analysis using recurrence relations

h[e_,k_,z_]:=e HypergeometricPFQRegularized[{(1/2)-e,(1/2)-e,(1/2)-e,1-e},{(1/2),1+k,1},z]

s[e]==(Sqrt[π])((Sin[π e])/e)h[e,0,1]//FullSimplify

eq[e_,k_,z_]=(Simplify`GivePFQRecurrence[{(1/2)-e,(1/2)-e,(1/2)-e,1-e},{(1/2),1+k+4,1},z,{6,6,6,6},{6,6,6,6}]/.{ Hold[HypergeometricPFQ][{(1/2)-e,(1/2)-e,(1/2)-e,1-e},{(1/2),m_,1},z]:>Pochhammer[k,m-k]hh[e,m-1,z] }
)//Factor//KeepFactorsDependingOn[#,hh]&//Collect[#,_hh,Simplify]&

eq[e,k,z]/.{hh->h}//Series[#,{z,0,5}]&//FullSimplify

Limit[(1-z)h[e,-1,z],z->1,Assumptions->e>0]

(Sqrt[π])((Sin[Pi e])/e)((hh[e,0,1]/.Solve[eq[e,-1,1]==0,hh[e,0,1]]/.hh->h//Together)/.e->0)

Limit[%,e->0]

Thank you

Initialization code

Get[ToFileName[{DirectoryName[ToFileName["FileName"/.NotebookInformation[EvaluationNotebook[]]]]},"DiffEqData.m"]]

DifferenceGraphics[pt_List]:=With[{xmin=Min[pt[[All,1]]]-1,xmax=Max[pt[[All,1]]]+1,ymin=Min[pt[[All,2]]]-1,ymax=Max[pt[[All,2]]]+1},
Graphics[{PointSize[0.1],Red,Point[pt]},GridLines->{Range[xmin,xmax],Range[ymin,ymax]},PlotRange->{{xmin,xmax},{ymin,ymax}}]]

DifferenceGraphics[pt_List,{{xmin_,xmax_},{ymin_,ymax_}}]:=
Graphics[{PointSize[0.1],Red,Point[pt]},GridLines->{Range[xmin,xmax],Range[ymin,ymax]},PlotRange->{{xmin,xmax},{ymin,ymax}}]

DifferenceGraphics[pt_List,{min_,max_}]:=
Graphics[{PointSize[0.1],Red,Point[pt]},GridLines->{Range[min,max],Range[min,max]},PlotRange->{{min,max},{min,max}}]

DifferenceGraphics3D[pt_List]:=
Deploy[Graphics3D[{PointSize[0.1],Red,Point[pt]},FaceGrids->All,FaceGridsStyle->Directive[Gray[0.4],Dotted],PlotRange->Automatic,BoxRatios->Automatic]]

Clear[KeepFactorsDependingOn];
KeepFactorsDependingOn[e_List,head_]:=KeepFactorsDependingOn[#,head]&/@e

KeepFactorsDependingOn[e_Times,head_]:=DeleteCases[e,_?(FreeQ[#,head]&)]

KeepFactorsDependingOn[e_,head__]:=e

RecurrenceAsymptote[eq_Equal,w__]:=With[{res=DSolve`RSolveAsymptoticSeries[eq,w]},res/;Head[res]===Rule]