## Kristen Lew, Homework 16 ## ExpMath, Spring 2012 # Post if you wish. ### Problem 1 ## 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 largest, k, current, returnThis: largest:=0: current:=0: returnThis:=0: if nops(L)=0 then RETURN(FAIL): fi: for k in L do current:=evalf(abs(evalc(k))): if largest < current then largest:=current: returnThis:=k: fi: od: returnThis: end: ### Problem 2 ## 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 lcoeff1, roots, Mu: lcoeff1:=lcoeff(ope,n): roots:=[solve(lcoeff1,N)]: Mu:=LargestAB(roots): if Mu=FAIL or evalf(Mu)<0 then RETURN(FAIL): fi: Mu: end: ### Problem 3 ## 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 given by: # a(m) IS ASYMPOTIC TO C*mu^m * m^theta # for some constant C, and the above mu that you found, using FindMu(ope,n,N). FindTheta:=proc(ope,n,N,Ini) local Mu, a, ApproxT,Theta: Mu:=FindMu(ope,n,N): if Mu=FAIL then RETURN(FAIL): fi: a:=HtoSeq(ope,n,N,Ini,1000): ApproxT:=evalf(log(a[1000]/(mu*a[999]))/log(1000/999)): convert(GuessT, confrac, `cvgts`): Theta:=cvgts[3]: end: HtoSeq:=proc(ope,n,N,Ini,M) local L,n1,ORD,BenMiles,kirsten,i: ORD:=degree(ope,N): if nops(Ini)<>ORD then RETURN(FAIL): fi: BenMiles:=expand(subs(n=n-ORD,ope)/N^ORD): L:=Ini: for n1 from nops(Ini)+1 to M do kirsten:=subs(n=n1,coeff(BenMiles,N,0)): if kirsten=0 then print(`blows up at`, n1): print(`so far we have`, L): RETURN(FAIL): fi: L:=[op(L), -add(subs(n=n1,coeff(BenMiles,N,-i))*L[n1-i],i=1..ORD)/kirsten]: od: L: end: ### Problem 4 ## Using all the above, write a program FindC(ope,n,N,Ini), that estimates the C in: # a(m) IS ASYMPTOTIC TO C*mum* mtheta # 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 a,Mu,Theta: a:=HtoSeq(ope,n,N,Ini,1000): Mu:=FindMu(ope,n,N): Theta:=FindTheta(ope,n,N,Ini): identify(evalf(a[1000]/Mu^1000/1000^Theta)): end: ### Problem 5 ## 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 Mu, Theta, C: Mu:=FindMu(ope,n,N): Theta:=FindTheta(ope,n,N,Ini): C:=FindC(ope,n,N,Ini): C*(Mu^n)*(n^Theta): end: ### Problem 6 ## 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 conjectures asymptotics for # a(n):=add(F(n,k),k=0..n) AsyBS:=proc(F,n,k,L) local ope, Order1, i , j, Ini: ope:=EmpericalZeilberger(F,n,k,N,L): Order1:=degree(ope,N): Ini:=[seq(evalf(add(subs({n=i,k=j},F),j=0..i)), i=1..Order)]: Asy(ope,n,N,Ini): end: EmpiricalZeilberger:=proc(F,n,k,N,M) local L: with(combinat): L:=eval([seq(eval(add(subs({n=n1,k=k1},F),k1=0..n1)),n1=1..M)]): GuessH(L,n,N): end: GuessH:=proc(L,n,N) local K,ope,DEG,ORD: ORD:=1: DEG:=0: for K from 1 while nops(L)>=(1+ORD)*(1+DEG)+ORD+6 do for ORD from 1 to K do DEG:=K-ORD: ope:=GuessH1(L,ORD,DEG,n,N): if ope<>FAIL then RETURN(ope): fi: od: od: FAIL: end: GuessH1:=proc(L,ORD,DEG,n,N) local ope,a,i,j,var,eq,n1,var1: if nops(L)<(1+ORD)*(1+DEG)+ORD+6 then print(`We need at least`, (1+ORD)*(1+DEG)+ORD+6, `terms in L`): RETURN(FAIL): fi: ope:=add(add(a[i,j]*n^i*N^j,i=0..DEG),j=0..ORD): var:={seq(seq(a[i,j],i=0..DEG),j=0..ORD)}: eq:={seq(add(subs(n=n1, coeff(ope,N,j))*L[n1+j],j=0..ORD), n1=1..(1+DEG)*(1+ORD)+6)}: var1:=solve(eq,var): ope:=subs(var1,ope): if ope=0 then RETURN(FAIL): fi: ope:=subs({seq(v=1, v in var)}, ope): ope:=add(factor(coeff(ope,N,j))*N^j, j=0..degree(ope,N)): eq:=expand({seq(add(subs(n=n1, coeff(ope,N,j))*L[n1+j],j=0..ORD), n1=1..nops(L)-ORD)}): if eq<>{0} then print(ope, `didn't work out`): RETURN(FAIL): else RETURN(ope): fi: end: # What are: # AsyBS(binomial(n,k)^2,n,k,50); # AsyBS(binomial(n,k)^3,n,k,50); # AsyBS(binomial(n,k)^4,n,k,50); # AsyBS(binomial(n+k,k)^2*binomial(n,k)^2,n,k,50); ### Problem 7 Challenge ## 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);