#C14.txt: March 11, 2019, The P-recursive (holonomic) "ansatz"
Help:=proc(): print(` SeqToP1(L,n,k,d), SeqToP1([seq(i!,i=1..20)],n,1,1); , SeqToP(L,n,MaxC), `):
print(`EstAsy(L,n), PtoSeq(P,N) , Asy(rec) `): end:
##from C12.txt
#C12.txt: March 4, 2019, (virtual class), Experimental Mathematics class, Math640 (Rutgers University)(Dr. Z.)
#Recall that a sequence {a(n), n=0..infinity} is C-finite (or order k) if there exist NUMBERS c1, ..., ck such that
#a(n)=c1*a(n-1)+ ...+ ck*a(n-k)
#One you know c1, c2, ..., ck, and the k-1 initial values a(0), ..., a(k-1), the sequence is uniquely determined
#Our data structure for C-finite sequence is
#C=[INI,Rec] where
#INI=[a(0),..., a(k-1)] and Rec=[c1,.., ck]. See
#http://sites.math.rutgers.edu/~zeilberg/EM19/C10.txt
#The imporant procedures were:
#CtoSeq(C,N) inputs a C-finite sequence (in our data-structure) and outputs the first N+1 terms of the sequence
#SeqToC(L): that inputs a list L and tries to find a linear-recurrence with CONSTANT coefficients satisfied by it
#(with enough `safey data'). Of course, if the list is not long enough, or the sequence is not C-finite it would return FAIL.
#The NATURAL generalization of C-finite sequences is P-recursive sequences (aka discrete holonomic sequences)
#that satisfy a linear recurrence with POLYNOMIAL coefficients. If its order is k, then for some polynomials
#P_0(n),..., P_k(n), we have
#P_0(n)*a(n)+P_1(n)*a(n-1)+ ... + P_k(n)*a(n-k)=0
#Of course we also need initial conditions, so our data structure for them would be [INI,Rec]
#but now Rec would be written as the triple, [n,[P_0,P_1,..., P_k]]. the first entry ,n, just denotes
#the symbol used to describe the polynomial coefficients of the recurrence.
#For example the sequence a(n)=n! satisfies (by definition of factorial)
#a(n)-n*a(n-1)=0, a(1)=1, so its representation is going to be
#[[1],[n,[1,-n]]. Also C-finite sequences can be represented in this data structure (now the polynomials P0,P1,..., Pk
#have degree 0, so the representation of the Fibonaii sequence is
#[[1,1],[n,[1,-1,-1]] ] or [[1,1],[n, [-1,1,1]] (why?)
#The representation of the sequence a(n)=2^n is
#[[2],[n,[1,-2]] ], etc.
#Note that our convention is that INI starts at n=1, rather than n=0, since it is more convenient, because in Maple
#L[i] is the i-th member of the list
#Our first goal is to write a procedure SeqToP(L,n), that inputs a sequence (usually of numbers, but of course, they may
#be symbolic), and a discrete variable name, n, and guesses a P-recursive description. Of course, we must start
#with an auxiliary procedure SeqToP1(L,n,k,d) that looks for such a description with order k and degree of the
#k+1 polynomials P0(n), ..., Pk(n) equal to d. This requires (k+1)*(d+1) degrees of freedom, and to be safe
#against `over-fitting' let's require 4 extra data points,
#and since we can only start collecting data for n>=k, we need the `training data', the input list L to be at least
#of length (k+1)*(d+1)+k+4. If it is less, we would get an ERROR message.
#SeqToP1(L,n,k,d): Inputs a list of length >= (k+1)*(d+1)+k+4, outputs a P-recursive representation (if it exists) that
#fits the list. For example, try:
#SeqToP1([seq(i!,i=1..20)],n,1,1);
SeqToP1:=proc(L,n,k,d) local eq,var,P,a,i,j,n1,Rec,c:
#Some checking for the input
if not (type(L,list) and type(n,symbol) and type(k,integer) and type(d,integer) and k>=1 and d>=0) then
print(`Bad input`):
RETURN(FAIL):
fi:
if nops(L)< (k+1)*(d+1)+k+4 then
print(`List too short`):
RETURN(FAIL):
fi:
#end checking the input
#Setting up the k+1 generic polynomials P[i], i=0..,k, each of degree d using a[i,j] for the coefficient of n^j in P[i]
for i from 0 to k do
P[i]:=add(a[i,j]*n^j,j=0..d):
od:
#Since the set of unknowns (undetermined coefficients) is {a[i,j], i=0..k,j=0..d)} and we are going use linear algebra (the solve command)
#we need to define it
var:={seq(seq(a[i,j],j=0..d),i=0..k)}:
#We now need to set up the set of equations, let's call it eq, by plugging-in into the tentative recurrence the values of the sequence
#we use n1 to access the terms of the sequence
eq:={seq(add(subs(n=n1,P[i])*L[n1-i],i=0..k),n1=k+1..nops(L))}:
#we now solve the system, renaming var as the solution
var:=solve(eq,var):
#We now plug-in into the generting Rec
Rec:=[seq(subs(var,P[i]),i=0..k)]:
#Since the 0 recurrence is always a recurrence, we have to check that it is not 0
if Rec=[0$(k+1)] then
RETURN(FAIL):
fi:
#since our system is homogeneous, we will get an annoying aribitrary factor, so let's divide by it
c:=lcoeff(Rec[1],n):
Rec:=[seq(normal(Rec[i]/c),i=1..nops(Rec))]:
#We now return the sequence [INI,[n,Rec]] in our data structure
[[op(1..k,L)],[n,Rec]]:
end:
#Now we are ready to keep trying SeqToP1(L,n,k,d) for the `smallest' k and d until it is no longer true that
#nops(L)>=(k+1)*(d+1)+k+4 and then we have to admit failure. Let's define the "complexity" of the P-recurrence
#sequence as k+d, but let's prefer low order to lower degree.
#SeqToP(L,n,MaxC): Inputs a list outputs a P-recursive representation (if it exists), and maximum complexity that
#fits the list. The output is the lowest Degree+Order that fits it, and we prefer lower order. For example, try:
#SeqToP([seq(i!,i=1..20)],n,2);
SeqToP:=proc(L,n,MaxC) local k,d,C,tem,L1,Rec,i,n1:
for C from 1 to MaxC do
for k from 1 to C do
for d from 0 to C-k while nops(L)>=(k+1)*(d+1)+k+4 do
#For the sake of efficiency, we truncate the list to shorter one
L1:=[op(1..(k+1)*(d+1)+k+4,L)]:
tem:=SeqToP1(L1,n,k,d):
if tem<>FAIL then
#To make sure, we check that the guessed recurrence holds for the full sequence
Rec:=tem[2][2]:
if normal({seq(add(subs(n=n1,Rec[i+1])*L[n1-i],i=0..k),n1=k+1..nops(L))})<>{0} then
print(tem, `did not work out`):
RETURN(FAIL):
else
RETURN(tem):
fi:
fi:
od:
od:
od:
#If nothing worked, we return FAIL
FAIL:
end:
##end from C12.txt
#PtoSeq(P,N): inputs a P-recursive sequence in the format [INI,[n,[P0,.., Pk]] where INI is the list of initial values for n=1, ...,k and
#P0,..., Pk are polynomials in n, representing the coefficients of the linear recurrence with positive coefficients
#P0(n)*x(n)+ ...+ Pk(n)*x(n-k)=0 with those initial conditions, and outputs the first N terms of the sequence
#(stolen from Victoria Chayes' hw12.txt)
PtoSeq:=proc(P,N) local L,n, nt, i, l,k,P0:
L:=P[1]:
n:=P[2][1]:
l:=nops(P[2][2])-1:
k:=nops(L)+1:
P0:=-1/P[2][2][1]:
while k<=N do
nt:=subs(n=k,P0)*add(subs(n=k, P[2][2][i+1])*L[-i],i=1..l):
L:=[op(L),nt]:
k:=k+1:
od:
L:
end:
#EstAsy(L,n): inputs a long list of pos. numbers L and a variable n, estimates
#C,mu,theta, such that L[n] is asympt. to C*mu^n*n^theta
#log(L[n])= logC+n*logmu+theta*log(n). Outputs [mu,theta]
EstAsy:=proc(L,n) local logC,logmu,theta,eq,var,n1,JS,mu,C:
n1:=nops(L):
var:={logmu,logC,theta}:
JS:=logC+n*logmu+theta*log(n):
eq:={seq(evalf(log(L[i]))=subs(n=i,JS),i=n1-2..n1)}:
var:=solve(eq,var):
mu:=exp(subs(var,logmu)):
theta:=subs(var,theta):
C:=exp(subs(var,logC)):
[C,mu,theta]:
#C*mu^n*n^theta:
end:
#P0(n)*x(n)+P1(n)*x(n-1)+ ...+Pk(n)*x(n-k)=0
#Asy(rec): inputs a P-recursive sequence rec=[INI, [n, [P0,...,Pk]]]
#and outputs mu and theta such L[n]=C*mu^n*n^theta for some C
Asy:=proc(rec) local n,rec0,d,INI,i,rec00,mu,p,mu0,LS,x,theta:
INI:=rec[1]:
n:=rec[2][1]:
rec0:=rec[2][2]:
d:=max(seq(degree(rec0[i],n),i=1..nops(rec0))):
rec00:=[seq(coeff(rec0[i],n,d),i=1..nops(rec0))]:
#add(rec00[i+1]*mu^(n-i),i=0..nops(rec00)-1)=0
p:=numer(add(rec00[i+1]/mu^i,i=0..nops(rec00)-1)):
mu0:=max(solve(p,mu)):
if not evalf(mu0)> 0 then
RETURN(FAIL):
fi:
#y(n)=x(n)/mu0^n:
#rec0:=[rec0[1],rec0[1]/mu, ..., rec0[k]/mu^k]
rec0:=[seq(rec0[i+1]/mu0^i,i=0..nops(rec0)-1)]:
#x(n)=n^theta:
LS:=add(rec0[i+1]*(1-i/n)^theta/n^d,i=0..nops(rec0)-1):
LS:=normal(subs(n=1/x,LS)):
LS:=taylor(LS,x=0,3);
if coeff(LS,x,0)<>0 then
RETURN(FAIL):
fi:
LS:=coeff(LS,x,1):
[mu0,solve(LS,theta)]:
end:
#BetterAsy(rec,n,K): inputs a P-recursive sequence rec=[INI, [n, [P0,...,Pk]]]
#and outputs the asympotic expansion of thesequence to order K
#L[n]=C*mu^n*n^theta*(1+a1/n+a2/n^2+...+aK/n^K)
# BetterAsy:=proc(rec,n,K) local n,rec0,d,INI,i,rec00,mu,p,mu0,LS,x:
#to be hw14.txt
# end: