<< HolonomicFunctions.m

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

(* Generalizations of Theorem 2.9 to Jacobi polynomials *)

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

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

FindRelation[jacann,Support->{S[a]S[b],S[a]S[b]S[n],S[a]S[b]S[n]^2,S[n]^2}]

(* Linearization coefficients for products of Legendre polynomials *)

time = Timing[ann=Factor[Annihilator[LegendreP[k,x] LegendreP[m,x] LegendreP[n,x] (2k+1)/2,{S[k],S[m],S[n],Der[x]}]];][[1]]

time = Timing[conn = Factor[CreativeTelescoping[ann,Der[x],{S[k],S[m],S[n]}]];][[1]]

(* principal part *)

conn[[1]]

(* delta part: one needs to verify that this vanishes indeed (factors: (1-x)(1+x)) *)

conn[[2]]

(* determine a recurrence solely in shifts in k *)

krec = Factor[FindRelation[conn[[1]],Support->{1,S[k],S[k]^2}]]

krec1=ApplyOreOperator[krec,a[k]]

FullSimplify[RSolve[First[krec1]==0,a[k],k]]

(* product recurrence *)

ann = Factor[Annihilator[JacobiP[i,2,0,x] JacobiP[j,1,1,x],{S[i],S[j]}]]

rec = FindRelation[ann,Eliminate->{x}]//Factor

Support[rec]

(* extract the recurrence from the annihilator *)

rec1 = ApplyOreOperator[rec,A[i,j]]

(* graphical interpretation of the support *)

Show[Graphics[{PointSize[0.1],Point[#]&/@(Cases[rec1,A[_,_],Infinity]/.A[i+a_.,j+b_.]->{a,b})}]]

Collect[Solve[(rec1/.{i->i-1,j->j-2})==0,A[i,j]][[1,1,2]],A[_,_],Factor]

(* compute the matrix entries recursively *)

ClearAll[A];
A[-1,j_Integer]:=0;
A[i_Integer,-1]:=0;
A[i_Integer,0]:=2;
A[i_Integer,j_Integer]:=A[i,j]=((i (2+i) (1+j) (1+2 j) A[-1+i,-1+j])/((1+i) (3+2 i) j (2+j)))-(((1+j) A[i,-2+j])/(2+j))-(((1+j) (1+2 j) A[i,-1+j])/((1+i) (2+i) j (2+j)))+(((1+i) (3+i) (1+j) (1+2 j) A[1+i,-1+j])/((2+i) (3+2 i) j (2+j)));

Timing[mat = Table[A[i,j],{i,0,50},{j,0,50}];]

Table[A[i,j],{i,0,15},{j,0,15}]//MatrixForm