#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: