# Matthew Russell # Experimental Math # Homework 16 # I give permission to post this read `public_html/em12/C14.txt`: read `public_html/em12/hw14.txt`: ##################################################################################### # Write a program, LargestAB(L) that inputs a list of complex numbers, and finds the # one with the largest absolute value. (Use the evalf and evalc Maple commands) # For example, # LargestAB([5,-3,-4,6,1]); # should give 6 and # LargestAB([5,-3,-4,6,-7,1]); # should give -7. LargestAB:=proc(L) local pos, currmax, i: pos:=1: currmax:=evalf(abs(evalc(L[1]))): for i from 2 to nops(L) do if evalf(abs(evalc(L[i])))>currmax then pos:=i: currmax:=evalf(abs(evalc(L[i]))): fi: od: return L[pos]: end: ##################################################################################### # Using the above program, write a program # FindMu(ope,n,N) # that inputs a linear recurrence operator with polynomial coefficients, ope(n,N), # and outputs the root of the largest absolute value of lcoeff(ope,n) (a certain polynomial in N), # if it is positive, or returns FAIL if it is negative, or if there are ties (i.e. both a and -a are roots), # or if the root(s) with maximal absolute value are complex. # For example, # FindMu((n+2)*N^2-(3*n+3)*N+(2*n+5),n,N); # should return: 2. FindMu:=proc(ope,n,N) local pol, L, cm: pol:=lcoeff(ope,n): L:=[solve(pol=0,N,AllSolutions)]: cm:=LargestAB(L): if Im(evalc(cm))<>0 then return `FAIL`: fi: if {op(L)} union {-cm} = {op(L)} then return `FAIL`: fi: if evalf(cm)<0 then return `FAIL`: fi: return cm: end: ##################################################################################### # Using the above, write a program # FindTheta(ope,n,N,Ini) # that inputs a linear recurrence operator ope, phrased in terms of n and the shift N, # and Ini, a list of length degree(ope,N) indicating the initial values of the sequence, # and outputs a guessed RATIONAL number such that the sequence defined by # a(m):=HtoSeq(ope,n,N,Ini,m)[m] # is asympotically given by # a(m) IS ASYMPOTIC TO C*mu^m n^theta # for some constant C, and the above mu that you found, using FindMu(ope,n,N). # [Hint: look at log(a(m+1)/a(m)/mu)/log(m/(m-1)) for very large m (you decide), and then use # convert(%,confrac,cvgts); # to get a good rational number approximation for theta .] FindTheta:=proc(ope,n,N,Ini) local l, cvgts, L, m, mu, c, i, flag, spot: L:=HtoSeq(ope,n,N,Ini,15001): if L=FAIL then return `FAIL`: fi: m:=15000: mu:=FindMu(ope,n,N): l:=evalf(log(L[m+1]/L[m]/mu)/log((m+1)/m)): c:=convert(l,confrac,cvgts): flag:=true: spot:=1: for i from 2 to nops(cvgts) while flag do if abs(c[i])>20*abs(c[i-1])+20 then flag:=false: spot:=i-1: fi: od: if flag then return cvgts[nops(cvgts)]: else return cvgts[spot]: fi: end: ##################################################################################### # Using all the above, write a program # FindC(ope,n,N,Ini), that estimates the C in: # a(m) IS ASYMPTOTIC TO C*mu^m m^theta # Use Maple's built-in command "identify" to see whether C can be expressed in terms of known constants. FindC:=proc(ope,n,N,Ini) local L,m,mu,theta,C: L:=HtoSeq(ope,n,N,Ini,15001): if L=FAIL then return `FAIL`: fi: m:=15000: mu:=FindMu(ope,n,N): theta:=FindTheta(ope,n,N,Ini): C:=evalf(L[m]/mu^m/m^theta): print(`Maple thinks that C is `, identify(C)): return identify(C): end: ##################################################################################### # Combine all the above, and write a procedure # Asy(ope,n,N,Ini) # that inputs ope and Ini as above and outputs the leading asymptotics. Asy:=proc(ope,n,N,Ini) local C, theta, mu: C:=FindC(ope,n,N,Ini): theta:=FindTheta(ope,n,N,Ini): mu:=FindMu(ope,n,N): return C*mu^n*n^theta: end: ##################################################################################### # By combining with # EmpiricalZeilberger(F,n,k,N,L); # Write a program: # AsyBS(F,n,k,L) # that inputs a binomial coefficient expression phrased in terms of n and k, and # a pos. integer L (for guessing, made large enough) and conjenctures asympotics for # a(n):=add(F(n,k),k=0..n) AsyBS:=proc(F,n,k,L) local ope,deg,Ini,n1,k1: ope:=EmpiricalZeilberger(F,n,k,N,L): deg:=degree(ope,N): Ini:=[seq(expand(add(subs({n=n1,k=k1},F),k1=0..n1)),n1=1..deg)]: return Asy(ope,n,N,Ini): end: # What are: # AsyBS(binomial(n,k)^2,n,k,50); = .5641848819*4^n/n^(1/2) # AsyBS(binomial(n,k)^3,n,k,50); = .3675444292*8^n/n # AsyBS(binomial(n,k)^4,n,k,50); = .2539650198*16^n/n^(3/2) # AsyBS(binomial(n+k,k)^2*binomial(n,k)^2,n,k,50); = .2200384509*(17+12*2^(1/2))^n/n^(3/2) ##################################################################################### # Challenge! [5 dollars to be divided among all correct answers] Write a program # BetterAsy(ope,n,N,Ini,r) # that empirically finds an asymptotic expression for the sequence, # let's call it a(n), annihilated by ope(n,N) with the # initial values given by Ini to order r, in other words # a(n)=C*mu^n*n^theta*(1+c[1]/n+c[2]/n^2+ ...+ c[r]/n^r +O(1/n^(r+1)) ), # for C,mu,theta as for Asy(ope,n,N,Ini), # and c[1],c[2], ..., c[r], "nice" rational numbers. If you can't find it, # it should return FAIL. Then write BetterAsyBS(F,n,k,L,r) and try out # BetterAsyBS(binomial(n,k)^2,n,k,50,2); # BetterAsyBS(binomial(n,k)^3,n,k,50,2); # BetterAsyBS(binomial(n,k)^4,n,k,50,2); # BetterAsyBS(binomial(n+k,k)^2*binomial(n,k)^2,n,k,50,2); #################################################################################### #################################################################################### #################################################################################### # STILL IN PROGRESS #################################################################################### #################################################################################### #################################################################################### BetterAsy:=proc(ope,n,N,Ini,r) local C, mu, theta, L, l, m, i, l1, m1: C:=FindC(ope,n,N,Ini): theta:=FindTheta(ope,n,N,Ini): mu:=FindMu(ope,n,N): #for i from 1 to r do #m:=15000+5000*i: m:=30000: L:=HtoSeq(ope,n,N,Ini,30000): l:=evalf((L[m]/C/mu^m/m^theta-1)*m): print(l): l:=RatApprox(l): print(l): m1:=30000: l1:=evalf(((L[m1]/C/mu^m1/m1^theta-1)*m1-l)*m1): print(l1): l1:=RatApprox(l1): #od: return C*mu^n*n^theta: end: BetterAsyBS:=proc(F,n,k,L,r) local ope,deg,Ini,n1,k1: ope:=EmpiricalZeilberger(F,n,k,N,L): deg:=degree(ope,N): Ini:=[seq(expand(add(subs({n=n1,k=k1},F),k1=0..n1)),n1=1..deg)]: return BetterAsy(ope,n,N,Ini,r): end: # RatApprox(r): Inputs a real number r and finds a reasonably nice # rational approximation to it, by calculating its continued fraction # representation RatApprox:=proc(r) local c, cvgts, flag, spot, i: c:=convert(r,confrac,cvgts): print(c): print(cvgts): flag:=true: spot:=1: for i from 2 to nops(cvgts) while flag do if abs(c[i])>20*abs(c[i-1])+50 then flag:=false: spot:=i-1: fi: od: if flag then print(cvgts[nops(cvgts)]): return cvgts[nops(cvgts)]: else print(cvgts[spot]): return cvgts[spot]: fi: end: