###################################################################### # Guessdiff: Save this file as Guessdiff. # # To use it, stay in the # # same directory, get into Maple (by typing: maple ) # # then type: read Guessdiff: # # Then follow the instructions given there # # # # Written by Arvind Ayyer, Rutgers University. # ###################################################################### with(SolveTools): with(linalg): with(gfun): read TILINGS: print(` This is Guessdiff: A Maple package to use the package TILINGS`): print(` to guess various properties about tilings in the plane`): print(` Written for the course: Experimental Mathematics, Spring, 2006`): print(` Uses the package TILINGS by Prof. Zeilberger `): Help:=proc() if args=NULL then print(` Guessdiff: A Maple package for empirically guessing ODEs`): print(): print(`For help with a specific procedure, type "Help(procedure_name);"`): print(`Contains procedures: finddiff, finddiff2, splitdiff2, splittile, tilingpattern `): elif nargs=1 and args[1]=finddiff then print(`finddiff(f,DEGREE,ORDER,x,D): guesses a differential operator annihilating`): print(`the taylor series f of degree DEGREE and order ORDER.`): print(`For example, try: finddiff([seq(1/i!,i=0..10)],0,1,x,D);`): print(`This is the series for the exponential function and so should return D-1`): elif nargs=1 and args[1]=finddiff2 then print(`finddiff2(f,DEGREE,ORDER,x,D): guesses a differential operator annihilating`): print(`the function f(x) of degree DEGREE and order ORDER.`): print(`Please keep the function variable as x.`); print(`For example, try: finddiff2(exp(x),0,1,x,D);`): print(`This should return D-1`): elif nargs=1 and args[1]=splitdiff2 then print(`splitdiff2(eq,y(z),DEGREE,ORDER) tries to find if the differential equation eq`): print(` with indep variable z and dep variable y can be split into two other differential equations`): print(` eq1 and eq2 such that eq1(f1)=0 and eq2(f2)=0 => eq(f1+f2)=0`): elif nargs=1 and args[1]=splittile then print(` splittile(tiles,m,n,DEGREE,ORDER) tries very hard to see if the `): print(` differential equation of deg DEGREE and ord ORDER satisfied by `): print(` the enumerator of an m*n board with tiles given by the set tiles `): print(` can be split as a sum of two differential equations. `): print(` Uses the package TILINGS`): elif nargs=1 and args[1]=tilingpattern then print(`tilingpattern(tiles,m,n,z) returns the roots of the polynomials `): print(`corresponding to tilings of mx1,...,mxn for the set of tiles `): print(`given by "tiles". Uses the package TILINGS.`): fi: end: Yafe:=proc(ope,N) local i,ope1,coe1,L: if ope=0 then RETURN(1,0): fi: ope1:=expand(ope): L:=degree(ope1,N): coe1:=coeff(ope1,N,L): #ope1:=normal(ope1/coe1): ope1:=normal(ope1): ope1:= convert( [seq(factor(coeff(ope1,N,i))*N^i,i=ldegree(ope1,N)..degree(ope1,N))],`+`): factor(coe1),ope1: end: #finddiff(f,DEGREE,ORDER,x,D): guesses a differential operator annihilating #the sequence f of degree DEGREE and order ORDER #For example, try: finddiff([seq(1/i!,i=0..10)],0,1,x,D); finddiff:=proc(f,DEGREE,ORDER,x,D) local ope,var,eq,i,j,k,n0,kv,var1,eq1,mu,a, sol: if (1+DEGREE)*(1+ORDER)+5+ORDER>nops(f) then ERROR(`Insufficient data for a recurrence of order`,ORDER, `degree`,DEGREE): fi: ope:=0: var:={}: for i from 0 to ORDER do for j from 0 to DEGREE do ope:=ope+a[i,j]*x^j*D^i: var:=var union {a[i,j]}: od: od: eq:={}: #for n0 from 2+DEGREE to (1+DEGREE)*(1+ORDER)+5+ORDER do #replaced 1+DEGREE in the starting point of n0 to DEGREE and it seems to work for n0 from 1 to nops(f)-ORDER do eq1:=0: for i from 0 to ORDER do for j from 0 to DEGREE do if n0+i-j>0 then eq1 := eq1 + coeff(coeff(ope,D,i),x,j)*mul(n0-1+i-j-k,k=0..i-1)*op(n0+i-j,f): fi: od: od: eq:= eq union {eq1}: od: #print(eq): var1:=solve(eq,var): kv:={}: for i from 1 to nops(var1) do mu:=op(i,var1): if op(1,mu)=op(2,mu) then kv:= kv union {op(1,mu)}: fi: od: ope:=subs(var1,ope): if ope=0 then RETURN(FAIL): fi: ope:={seq(coeff(expand(ope),kv[i],1),i=1..nops(kv))} minus {0}: sol := {}: if nops(ope) = 0 then return FAIL: else for i from 1 to nops(ope) do sol := sol union {Yafe(ope[i],N)[2]}: od: fi: return op(sol): end: #finddiff2(f,DEGREE,ORDER,x,D): guesses a differential operator annihilating #the function f(x) of degree DEGREE and order ORDER #For example, try: finddiff2(exp(x),0,1,x,D); finddiff2 := proc(f,DEGREE,ORDER,x,D): return finddiff([seq(coeftayl(f,x=0,n),n=0..(1+DEGREE)*(1+ORDER)+5+ORDER)],DEGREE,ORDER,x,D): end: #logfn(g,n) takes as input a function g of x and an integer n and returns another function of x whose k'th term in its taylor series is log of the kth term in the taylor series of g(x) for k=0..n logfn := proc(g,n) local i: return log( mul( ((D@@i)(g)(0)/i!)^(x^i), i=0..n)): end: #polysolve(eqns,x) returns the solution of a group of polynomial equations equated to zero. #example polysolve({ax+b-cx+d,cx+d},x) should return {a=c=0,b=d=0} polysolve := proc(eqns,x) local i,j: return solve({seq(seq(coeff(eqns[i],x,j),j=0..nops(eqns[i])),i=1..nops(eqns))}): end: test1 := proc(set) local i,j: for i from 1 to nops(set) do for j from 1 to nops(set[i]) do print(i,j,set[i][j]): od: od: end: #splitdiff2(eq,yofz,DEGREE,ORDER) tries to find if the differential equation eq with indep variable z and dep variable y can be split into two other differential equations eq1 and eq2 such that eq1(f1)=0 and eq2(f2)=0 => eq(f1+f2)=0 splitdiff2 := proc(eq, yofz,DEGREE,ORDER) local y,z,neweq,eqn,eq1,eq2,i,i2,j,j2,k,l,soln,p0,q0,p,q,ch,S: y := op(0,yofz): z := op(yofz): eqn := convert(eq,D): S := {}: for i from 1 to ORDER-1 do eq1 := add(add( p[l][k]*z^k, k=0..DEGREE/2)*(D@@(l-1))(y)(z),l=1..i+1); eq2 := add(add( q[l][k]*z^k, k=0..DEGREE/2)*(D@@(l-1))(y)(z),l=1..ORDER-i+1); #print(eq1,eq2): neweq := `diffeq+diffeq`(eq1, eq2, y(z)): #return neweq: if nops(neweq)=1 then neweq := convert(neweq,D): else for k from 1 to nops(neweq) do if type(neweq[k],algebraic) = true then neweq := convert(neweq[k],D): break; fi: od: fi: #print(neweq): #The main solve command soln := polysolve( { seq( coeff(eqn, (D@@n)(y)(z)) - coeff(neweq, (D@@n)(y)(z)), n=0..ORDER) },z ): #print(`neweq=`,neweq,`soln=`,soln); #Checking if it is a nontrivial solution if soln <> {} then for k from 1 to nops({soln}) do for l from 1 to nops({soln}[k]) do if lhs({soln}[k][l]) <> rhs({soln}[k][l]) and rhs({soln}[k][l]) <> 0 then S := S union {factor({subs({soln}[k],eq1),subs({soln}[k],eq2)})}: fi: od: od: return S: fi: p := 'p': q := 'q': od: return FAIL: end: #splittile(tiles,m,n,DEGREE,ORDER) tries very hard to see if the differential equation of deg DEGREE and ord ORDER satisfied by the enumerator of an m*n board with tiles given by the set tiles can be split as a sum of two differential equations. Uses the package TILINGS from Dr. Z's page. splittile := proc(tiles,m,n,DEGREE,ORDER) local i,j,y,gf,diffgf,diffgf2,h,t,H,yel,D,splits: gf := GFtH(n,tiles,t,H): yel := BuildTree1(n,tiles)[1][2]: h := [seq(H[yel[i][2]],i=1..nops(yel))]: diffgf := [seq(finddiff2(coeftayl(gf,t=0,m),DEGREE,ORDER,h[i],D),i=1..nops(h))]: diffgf := expand(diffgf): #for i from 1 to nops(diffgf) do # if diffgf[i] = FAIL then # print(`No differential equation of given degree and order found`): # return FAIL: # fi: #od: for i from 1 to nops(h) do diffgf2[i] := subs( {seq(D^j = (D@@j)(y)(h[i]),j=1..ORDER)}, diffgf[i]): #splits[i] := splitdiff2(diffgf2[i],y(h[i]),DEGREE,ORDER): od: return diffgf: end: #eq1 := D(y)(x) - x*y(x); eq2 := D(y)(x)-2*y(x); eq := `diffeq+diffeq`(eq1,eq2,y(x)); #eq1 := D(y)(x) - 3*y(x); eq2 := D(y)(x)-(x^2+1)*y(x); eq :=`diffeq+diffeq`(eq1,eq2,y(x)); #tilingpattern(tiles,m,n,z) returns the roots of the polynomials corresponding to tilings of mx1,...,mxn for the set of tiles given by "tiles". Uses the package TILINGS. tilingpattern := proc(tiles,m,n,z) local gf,yel,h,H,t,i,pols,roots,mins: Digits := 10: gf := GFtH(m,tiles,t,H): yel := BuildTree1(m,tiles)[1][2]: h := [seq(H[yel[i][2]],i=1..nops(yel))]: pols := {}: for i from 0 to n do if coeftayl(gf,t=0,i) <> FAIL then pols := pols union {subs({h[1]^2=z, seq(h[i]=1,i=2..nops(h))},coeftayl(gf,t=0,i))}: fi: od: #return pols: roots := [ seq([fsolve(pols[i])],i=1..nops(pols))]: #mins := [seq( min(op(roots[i])),i=1..nops(roots))]: return roots: end: