#!/usr/local/bin/maple # -*- maplev -*- # Nathaniel Shar # HW 16 # Experimental Mathematics # It is okay to link to this assignment on the course webpage. Help := proc(): print(`LargestAB(L), FindMu(ope, n, N), HtoSeq(ope,n,N,Ini,M), FindTheta(ope, n, N, ini), FindC(ope, n, N, ini, M:=1000), Asy(ope, n, N, ini), EmpiricalZeilberger(F,n,k,N,M), AsyBS(F, n, k, L)`): end: ############# # Problem 1 # ############# LargestAB := proc(L) local i, mc: mc := 0: for i in L do: if evalf(abs(i)) >= evalf(abs(mc)) then: mc := i: fi: od: return mc: end: ############# # Problem 2 # ############# FindMu := proc(ope, n, N) local roots, largest: roots := convert({solve(lcoeff(ope, n) = 0)}, list): largest := LargestAB(roots): if Im(largest) <> 0 then: return FAIL: fi: if LargestAB(convert(convert(roots, set) minus {largest}, list)) = largest then: return FAIL: fi: if evalf(largest) < 0 then: return FAIL: fi: return largest: end: ############# # Problem 3 # ############# # We need HtoSeq from class 14: 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: ## FindTheta := proc(ope, n, N, ini) local L, mu, absr2, m, cvgts: mu := FindMu(ope, n, N): if mu = FAIL then: return FAIL: fi: # To decide how many terms to compute, let r2 be the # second-largest root of lcoeff(ope, n, N). We'll continue until # (abs(r2)/abs(mu))^m < 0.001, where 0.001 is just a constant # chosen arbitrarily but is hopefully good enough. absr2 := max(seq(abs(i), i in {solve(lcoeff(ope, n))} minus {mu})): if mu = 0 then: # We're in trouble! return FAIL: fi: if evalf(absr2) <= 0 then: m := 100 # aribtrary else: m := ceil(log(0.001)*(log(absr2)-log(abs(mu)))): fi: L := HtoSeq(ope, n, N, ini, m): if L = FAIL then: return FAIL: fi: convert(log(L[m]/L[m-1]/mu)/log((m+1)/m), confrac, cvgts): if nops(cvgts) = 1 then: return cvgts[1]: else: return cvgts[2]: fi: end: ############# # Problem 4 # ############# # M = number of terms to compute in approximating C. # The call to identify is very slow and the user can provide it # manually if desired. FindC := proc(ope, n, N, ini, M:=1000) local L, mu, theta, L2: L := HtoSeq(ope, n, N, ini, M): mu := FindMu(ope, n, N): theta := FindTheta(ope, n, N, ini): return evalf(L[M]/mu^M/M^theta): end: ############# # Problem 5 # ############# Asy := proc(ope, n, N, ini): return [FindMu(ope, n, N), FindTheta(ope, n, N, ini), FindC(ope, n, N, ini)]; end: ############################################################################### # From HW 14: EmpiricalZeilberger := proc(F,n,k,N,M) local L, n1, k1: L := [seq(sum(subs({n=n1, k=k1}, F), k1=0..n1), n1=1..M)]; return GuessH(L,n,N): end: # From Class 14: #GuessH1(L,ORD,DEG,n,N): Inputs a sequence of numbers L #and pos. integer ORD and non-neg. integer DEG and letters #n and N and outputs an linear recurrence operator with #poly. coeffs. (conj.) annihilating the sequence L #nops(L)>=(1+ORD)*(1+DEG)+ORD+6 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: #end old stuff #start new stuff #GuessH(L,n,N): Inputs a sequence of numbers L # and symbols n and N #and finds a minimal (in the sense of ord+deg) # linear recurrence operator with #poly. coeffs. (conj.) annihilating the sequence L with #nops(L)>=(1+ORD)*(1+DEG)+ORD+6 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: ############################################################################### ############# # Problem 6 # ############# AsyBS := proc(F, n, k, L) local N, ope, ord: ope := EmpiricalZeilberger(F,n,k,N,L): printf("Found operator %s...\n", convert(ope, string)): ord := degree(ope, N): return Asy(ope, n, N, [seq(sum(subs({n=i,k=j}, F), j=0..i), i=1..ord)]): end: # AsyBS(binomial(n,k)^2,n,k,50) = [4, -1/2, 0.2820595321] # (that is, 0.28206*4^n/sqrt(n) # identify is not able to determine a formula for 0.2820595321, but it # is easy to get (from Stirling's formula) that the constant should be # 1/(2sqrt(pi)), which is 0.2820947918. # AsyBS(binomial(n,k)^3,n,k,50) = [8, -1, 0.3674300930] # AsyBS(binomial(n,k)^4,n,k,50) = [16, -3/2, 0.2538317153] # AsyBS(binomial(n+k,k)^2*binomial(n,k)^2,n,k,50) : I had to time it # out after 10 minutes. (I can't figure out what is so slow!)