### Homework 12 - Patrick Devlin - March 4 2012 ### ## I have no preference about keeping my work private ## Help:=proc(): print(`GuessDE(L, x, Di, Max), AlgToDE(F, P, x, Di, K)`): print(`From class 13:`): print(`LIF(p,u,N)`): print(`GuessHGivenDegAndOrd(L,DEG,ORD,n,N,printErrors:=false)`): print(`-- Type HelpAll() to see all functions --`): end: HelpAll:=proc(): print(`From class 12:`): print(`AlgToSeq(F, P, x, K)\n`): print(`From class 10:`): print(`GuessAlg(L, P, x, MaxDegP:=(nops(L)/2 - 10)), GuessAlgHelp(L, P, x, degP, maxDegX:=infinity)`): print(`GuessAlgExactDegrees(L, P, x, degP, degX)\n`): print(`From class 9:`): print(`GenHomPolG(X, d, a, DegList:=[infinity], co:=0), GPg(X, d, a, DegList:=[infinity], co:=0)`): print(`GuessPR1(L, x, ORDER, DEGREE)`): print(`From previous classes:`): print(`GenPol(X, d, a, co), GenHOMPol(X, d, a, co)`): print(`GP(X, d, a), GuessPR(L, x, ORDER, DEGREE)`): end: ################### # Problem 1 and 2 # ################### # There's no good reason that the iterative method used ought to converge to anything at all. # Moreover, there are (usually) several power series that satisfy a given expression. # Also! If a method DOES converge to a power series, there's no good reason # each coefficient ought to converge, and there's kind of "super" no good # reason for a method to converge exactly after just a handful of iterations. # There are some necessary and sufficient conditions I can think of in terms # of abstract algebra, but they're not enlightening, and they're hard to # check. # # If we're going for something, I'd probably just write out a power series # with undetermined coefficients and try to solve a bunch of nonlinear # equations. That method at least has a good reason to work, and it's clear # why we should expect multiple solutions and such. ############# # Problem 3 # ############# GuessDE:=proc(L,x,Di, Max) local temp, i: for i from 0 to Max do temp:=GuessDEGivenOrderPlusDegree(L,x,Di,i): if(temp <> FAIL) then return temp: fi: od: return FAIL: end proc: GuessDEGivenOrderPlusDegree:=proc(L, x, Di, M) local temp, ord: for ord from 0 to M do temp:=GuessDEGivenOrderAndDegree(L, x, Di, ord, M-ord): if(temp <> FAIL) then return temp: fi: od: return FAIL: end proc: GuessDEGivenOrderAndDegree:=proc(L,x,Di, ord, deg) local ope, a, vars, eqns, f, i, j, stuffAppliedToF, solvedStuff: f:=sum(L[i] *x^i, i=1..nops(L)): ope:=sum(sum(a[i,j] *x^i * Di^j, i=0..deg), j=0..ord): vars:={seq(seq(a[i,j],i=0..deg),j=0..ord)}: stuffAppliedToF:= sum(a[i,0] *x^i, i=0..deg)*f + sum( sum(a[i,j]*x^i, i=0..deg) *diff(f,x$j), j=1..ord): stuffAppliedToF:= expand(stuffAppliedToF): eqns:={seq( coeff(stuffAppliedToF, x, i) = 0, i=1..(nops(L)- ord-2))}: solvedStuff:=solve(eqns, vars): ope:=eval(ope, solvedStuff): if(testeq(ope,0) or (solvedStuff=null) or (solvedStuff=[])) then return FAIL: fi: return eval(ope, [seq(v=1, v in vars)]): end proc: ############# # Problem 4 # ############# AlgToDE:=proc(F, P, x, Di, K, Max:=5): return GuessDE(AlgToSeq(F, P, x, K+1), x, Di, Max): end proc: #################### # From other stuff # #################### LIF:=proc(p, u, N) local n, s: s:=taylor(p,u,N+2): return [seq(coeftayl(s^n, u=0, n-1)/n,n=1..N)]: # not even checking if p is in fact a function in u holomorphic in the neighborhood u=0. end proc: ## This takes a list, L, and tries to output a linear recurrence operator with polynomial coefficients # of order exactly ORD and maximum degree exactly DEG that annihilates L. # It's in terms of the shift operator N and polynomials n. # ## It will output something like (a_0 (n) + a_1 (n) N + a_2 (n) N^2 + ... + a_ORD (n) N^ORD) L = 0 # or FAIL if it fails. # If the polynomials a_0 (n) = a_1 (n) = ... = a_ORD (n) = 0 are the only solution, then it will return FAIL GuessHGivenDegAndOrd:=proc(L, ORD, DEG, n, N, printErrors:=false) local ope, a, i, j, vars, eqns, n1, solvedVars: if(nops(L) < (1+ORD)*(1+DEG)+ORD+6) then if(printErrors) then print(`List is too short! Need at least `, (1+ORD)*(1+DEG)+ORD+6, ` terms in L`): fi: return FAIL: fi: ope:=add(add(a[i,j]*n^i * N^j, i=0..DEG),j=0..ORD): vars:={seq(seq(a[i,j], i=0..DEG),j=0..ORD)}: eqns:=[seq( add(eval(coeff(ope,N,i)*L[n1+i],{n=n1,N=1}), i=0..ORD) = 0, n1=1..min((1+DEG)*(1+ORD)+6, (nops(L)-ORD)))]: #this "min" is actually just for good measure. We already checked this at first. # We don't solve for ALL the equations cause that would be solving perhaps a LOT of stuff if L is large. solvedVars:=solve(eqns, vars): ope:=subs(solvedVars, ope): if(testeq(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)): return ope: end proc: ############ # Class 12 # ############ ### AlgToSeq(F,P,x,K) ## Inputs: # (i) an algebraic expression F in variables P and x, where P is a the gen func. of a mystery sequence in the variable x # (ii) P, the generating function # (iii) an indeterminant x that is the dummy variable in the gen function P # (iv) a positive integer K ## Outputs: # A sequence of length K+1 denoting the coefficient of x^0 , x^1 , ... , x^K in the Maclaurin expansion of P # - If this fails, it outputs FAIL. - # ## Examples: # AlgToSeq(P*(1-x)-1, P, x, 4) --outputs--> [1,1,1,1,1] # AlgToSeq(P-1-x*P^2, P, x, 4) --outputs--> [1,1,2,5,14] AlgToSeq:=proc(F,P,x,K,start:=1) local H, G, i, F1, f: H:=coeff(F,P,1): if(testeq(0, eval(H,x=0))) then return FAIL: fi: # Need to check this or the method will fail # Think...! The expression F says... P = (P*H-F)/H. So set this F1:=(P*H-F)/H: # Then F1 = -(a0 (x) + a2(x)*P^2 + ...) / H(x). # So if we iterate this, the P's ought to go further and further away, and we'll be left with F1(x) = STUFF(x) + P^BIG * OTHERSTUFF(P,x) f:=start: for i from 1 to K do f:=expand(eval(F1, P=f)): f:=series(f,x=0,2*K+10): # truncate it for speed od: f:=series(f,x=0,K+1): return [seq(coeff(f,x,i), i=0..K)]: end proc: ############ # Class 10 # ############ # GuessAlg(L, P, x, MaxDegP:=(nops(L)/2 - 10)): Takes an input sequence of numbers L and # tries to guess and algebraic equation of the form: # F(P,x), where P is the symbol for the generating function of L, [as best as we can tell!], # and x is an indeterminant (corresponding to generating function for L). GuessAlg:=proc(L, P, x, MaxDegP:=(nops(L)/2-10), MaxDegX:=infinity) local g, DegP: for DegP from 1 to MaxDegP do g:=GuessAlgHelp(L, P, x, DegP, MaxDegX): if(g <> FAIL) then return g: fi: od: return FAIL: end proc: GuessAlgHelp:=proc(L, P, x, DegP, MaxDegX:=infinity) local g, DegX: for DegX from 0 while( ((1+DegX)*(1+DegP) < nops(L) - 6) and (DegX <= MaxDegX)) do g:=GuessAlgExactDegrees(L, P, x, DegP, DegX): if(g <> FAIL) then return g: fi: od: return FAIL: end proc: ### ## Given a list L, formal symbols, P, and x, and two integers, DegP, DegX, ## it tries to return a polynomial of the form ## F(P, x) [not equal 0], with max degree of P equal to DegP, and max degree of x equal to DegX ## such that if genP(x) is the generating function of L, then F(genP, x) is close enough to 0. GuessAlgExactDegrees:=proc(L, P, x, DegP, DegX) local F, a, i, j, vars, genP, FofP, eqns, solvedVars, v: F:=add(add(a[i,j]*x^i * P^j, i=0..DegX), j=0..DegP): vars:={seq(seq(a[i,j],i=0..DegX), j=0..DegP)}: genP:=add(L[i] * x^(i-1), i=1..nops(L)): FofP:=series(eval(F, P=genP), x, nops(L)+10): eqns:={ seq( coeff(FofP, x, i) = 0, i=0..nops(L)-1)}: solvedVars:=solve(eqns, vars): F:=eval(F, solvedVars): if (testeq(F,0)) then return FAIL: fi: # This next line is "optional". Dr. Z. put it in to get rid of "stupid looking" free parameters F:=eval(F, {seq(v = 1, v in vars)}): return add(factor(coeff(F,P, i))*P^i, i=0..degree(F,P)): # this is just F cleaned up a little end proc: ########### # Class 9 # ########### # GenHomPolG(X,d,a,DegList, co): # ## Inputs: # (i) a list of variables X e.g. [x,y,z] (let m=nops(X), the number of variables) # (ii) a pos. integer d, for the degree # (iii) a symbol a by which the coeff. (undetermined) of the pol. are expressed # (iv) a list of degrees for the individual variables (corresponding to X) # (v) co, starting counter # ## Outputs: # (i) A generic HOMOG. polynomial in the variables X of degree d # with coeff. expressed as a[co],a[co+1], ..., a[co+AsNeeded] # and with the restriction that the degree of X[i] <=DegList[i] # (ii) the set of coefficients # (iii) the new value of the counter # # For example, GenHomPolG([x,y],3,a,[2,2],0) --outputs--> a[0]*x^2*y + a[1]*x*y^2, {a[0],a[1]}, 2 GenHomPolG:=proc(X,d,a,DegList:=[infinity],co:=0) local P, var, co1, m, i, x, X1, P1, newDegList: option remember: m:=nops(X): newDegList:=DegList: if(m>nops(DegList)) then newDegList:=[op(DegList), infinity$m]: fi: # or maybe we wanna return FAIL instead? if(convert(newDegList,`+`) < d) then return (0, {}, co): fi: x:=X[1]: if m=1 then RETURN(a[co]*x^d, {a[co]},co+1): # We don't need to wonder if d > DegList[1] by that above check fi: co1:=co: X1:=[op(2..m,X)]: newDegList:=[op(2..m, newDegList)]: P:=0: var:={}: for i from 0 to min(d,DegList[1]) do P1:=GenHomPolG(X1, d-i, a, newDegList, co1): P:=expand(P+P1[1]*x^i): var:=var union P1[2]: co1:=P1[3]: od: return (P,var,co1): end proc: #GPg(X,d,a, DegList:=[infinity]): # ## Inputs: # (i) a list of variables X e.g. [x,y,z] (let m=nops(X), the number of variables) # (ii) a pos. integer d, for the degree # (iii) a symbol a by which the coeff. (undetermined) of the pol. are expressed # (iv) an optional parameter for DegList, which is a list of maximum degrees for the variables # (v) an optional parameter for where the index of the coefficients should start # ## Outputs: # (i) A generic polynomial in the variables X of total degree d with coeff. expressed as a[co],ac[co+1], ..., a[co+AsNeeded] # (ii) the set of coefficients # (iii) the new value of the counter # # For example, GenPol([x],1,a,0) --outputs--> a[0]+a[1]*x,{a[0],a[1]},2 GPg:=proc(X,d,a, DegList:=[infinity], co:=0) local P, var,co1,d1,P1: co1:=co: P:=0: var:={}: for d1 from 0 to d do P1:=GenHomPolG(X,d1,a,DegList,co1): var:=var union P1[2]: co1:=P1[3]: P:=P+P1[1]: od: return (P,var): end proc: # GuessPR1(L,x,ORDER,DEGREE) ## Inputs # (i) a sequene L of numbers # (ii) a SYMBOL x (for the indexed variables x[1],x[2], ...) # (iii) a pos. integer, ORDER # (iv) a pos. integer Degree ## Output # A polynomial, let's call it, P(x[0],x[1], ..., x[ORDER]) such # that for all n P(L[n],L[n+1], ..., L[n+ORDER])=0 # and the degree of X[ORDER] is 1 GuessPR1:=proc(L,x,ORDER,DEGREE) local X,i,P,var,eq,a: X:=[seq(x[i],i=0..ORDER)]: P:=GPg(X,DEGREE,a,[DEGREE$ORDER, 1]): # Note, this is the only line that was changed from GuessPR [below] var:=P[2]: if nops(var)+6>nops(L) then print(`Please make the list longer, we need`, nops(var)+6 , `terms`): RETURN(FAIL): fi: P:=P[1]: eq:={ seq( subs( {seq(x[i]=L[n+i],i=0..ORDER)}, P)=0, n=1.. nops(var)+6 ) }: var:=solve(eq,var): factor(subs(var,P)): end: #GuessPR(L,x,ORDER,DEGREE) #Inputs #(i)a sequene L of numbers #(ii) a SYMBOL x (for the indexed variables x[1],x[2], ...) #(iii) a pos. integer, ORDER #(iv) a pos. integer Degree #Output #A polynomial, let's call it, P(x[0],x[1], ..., x[ORDER]) such #that for all n P(L[n],L[n+1], ..., L[n+ORDER])=0 GuessPR:=proc(L,x,ORDER,DEGREE) local X,i,P,var,eq,a: X:=[seq(x[i],i=0..ORDER)]: P:=GP(X,DEGREE,a): var:=P[2]: if nops(var)+6>nops(L) then print(`Please make the list longer, we need`, nops(var)+6 , `terms`): RETURN(FAIL): fi: P:=P[1]: #eq:={ seq( subs( {x[0]=L[n],x[1]=L[n+1], ..., x[ORDER]=L[n+ORDER]}, # P)=0, n=1..min(nops(L)-ORDER, nops(var)+6) }: eq:={ seq( subs( {seq(x[i]=L[n+i],i=0..ORDER)}, P)=0, n=1.. nops(var)+6 ) }: var:=solve(eq,var): factor(subs(var,P)): end: