### Homework 16 - Patrick Devlin - March 23 2012 ### ## I have no preference about keeping my work private ## Help:=proc(): print(`LargestAB(L), FindMu(ope, n, N), FindTheta(ope, n, N, Ini), FindC(ope, n, N, Ini)`): print(`Asy(ope, n, N, Ini), AsyBS(F, n, k, L), BetterAsy(ope, n, N, Ini, r)`): print(`-- Type HelpAll() to see all functions --`): end: HelpAll:=proc(): print(`HtoSeq(ope, n, N, Ini, M), EmpiricalZeilberger(F, n, k, N, L), GuessH(L, n, N)`): end: ############# # Problem 1 # ############# LargestAB:=proc(L) local biggestModulus, biggestIndex, i: biggestIndex:=1: biggestModulus:=evalf(abs(L[1])): for i from 1 to nops(L) do if (evalf(abs(L[i])) > biggestModulus) then biggestModulus:=evalf(abs(L[i])): biggestIndex:=i: fi: od: return L[biggestIndex]: end proc: ############# # Problem 2 # ############# FindMu:=proc(ope, n, N) local L, largest: L:=[fsolve(lcoeff(ope,n)=0, N)]: try: largest:=LargestAB(L): if(largest>0) then return largest: fi: catch: return FAIL: end try: return FAIL: end proc: ############# # Problem 3 # ############# FindTheta:=proc(ope, n, N, Ini, big:=10000): return FindThetaGivenMu(ope, n, N, Ini, FindMu(ope, n, N), big): end proc: FindThetaGivenMu:=proc(ope, n, N, Ini, mu, big:=10000) local L, th: if(mu = FAIL) then return FAIL: fi: L:=HtoSeq(ope, n, N, Ini, big): try: th:=convert( log(L[big]/(L[big-1] *mu))/log(big/(big-1)), confrac, cvgsts): return th[1] + 1/th[2]: catch: return FAIL: end try: end proc: ############# # Problem 4 # ############# FindC:=proc(ope, n, N, Ini, big:=10000) local mu: mu:=FindMu(ope, n, N): FindCGivenMuAndTheta(ope, n, N, Ini, mu, FindThetaGivenMu(ope, n, N, Ini, mu, big)): end proc: FindCGivenMuAndTheta:=proc(ope, n, N, Ini, mu, theta, big:=10000) local L, C: if( (mu = FAIL) or (theta = FAIL)) then return FAIL: fi: L:=HtoSeq(ope, n, N, Ini, big): try: return identify(L[big]/( (mu^big) * (big ^ theta))): catch: return FAIL: end try: end proc: ############# # Problem 5 # ############# Asy:=proc(ope, n, N, Ini, big:=100000) local mu, theta, C, m: mu:=FindMu(ope, n, N): theta:=FindThetaGivenMu(ope, n, N, mu, big): C:=FindCGivenMuAndTheta(ope, n, N, Ini, mu, theta, big): if((mu = FAIL) or (theta = FAIL) or (C = FAIL)) then return FAIL: fi: return (C* (mu^m) * (m ^theta)): end proc: ############# # Problem 6 # ############# AsyBS:=proc(F, n, k, L) local Ini, n1, k1, ope: Ini:=[seq(expand(add(subs({n=n1,k=k1},F),k1=0..n1)),n1=1..L)]: ope:=EmpiricalZeilberger(F, n, k, N, L): if(ope = FAIL) then return FAIL: fi: return Ays(ope, n, N, Ini): end proc: #################### # From other stuff # #################### EmpiricalZeilberger:=proc(F,n,k,N,M) local L: L:=[seq(expand(add(subs({n=n1,k=k1},F),k1=0..n1)),n1=1..M)]: return GuessH(L,n,N): end: #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: #HtoSeq(ope,n,N,Ini,M): Inputs a linear recurrence operator #ope phrased in terms of the discrete variable n and #its corresponding shift operator N, and ORD=degree(ope,N) #values for the sequence outputs the first M terms of #the sequence. For example, #HtoSeq(n+1-N,n,N,[1],4); should yield #[1,2,6,24] #N^(-ORD)*ope(N,n) 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: return L: end: