#hw3.txt, Feb. 2 2014, Anthony Zaleski #Packages we're using #(including C3.txt from class,which should be put #in the same folder as this file): with(linalg): with(combinat): read("C3.txt"): #Help() Help:=proc() print('Help(),HilDiff(n,d),MinDigits(N,eps,digits),inv1(p),gf(n,q),gfGuess(n,q)'); end: ##########Problem 1: Garvan pp. 60-90########## #Note: to see the output--which includes plots #& animation, you must execute this. #Differential Equations f:='f': y:=f(x): dsolve(diff(y,x)=y,y); #2D Plots plot(sin(1/x),x=0..Pi); plot([t*cos(t),t*sin(t),t=0..8*Pi]); with(plots): polarplot(cos(10*t),t=0..2*Pi,numpoints=10000); implicitplot(sin(x)+sin(y)=.5,x=-4*Pi..4*Pi,y=-4*Pi..4*Pi,numpoints=10000); #3D Plots x:='x': y:='y': plot3d(sin(x^2+y^2),x=-Pi..Pi,y=-Pi..Pi); plot3d([r*sin(t),r*cos(t),t],r=0..1,t=0..3*Pi); spacecurve([cos(t),sin(t),t],t=0..6*Pi); implicitplot3d(y^2-x^2=z,x=-2..2,y=-2..2,z=-2..2); animate3d([r*sin(t+a),r*cos(t+a),t+a],r=0..1,t=0..3*Pi,a=0..6*Pi,title="Rotating helix animation.",frames=40); #Linear Algebra with(linalg): A:=matrix([[1,2],[3,4]]); B:=matrix(2,2,[2,1,3,4]); multiply(A,B); A+B; evalm(%); ##########Problem 2: Floating Point Error########## #Here we record the number of digits d required to get #evalf(det(Hil(n))) and det(evalf(Hil(n))) to agree. #The following function outputs the descrepancy #in terms of n and d. HilDiff:=proc(n,d): Digits:=d: abs(evalf(det(Hil(n)))-det(evalf(Hil(n)))): end: #Next, MinDigits takes a list N of n-values #and outputs a list of the min number of #digits to get HilDiffeps do d:=d+1: od: L:=[op(L),d]: print(L): d:=1: od: L: end: #Here is what we find: #MinDigits(L,10^(-100),[seq(i,i=1..10)]); #[10, 6, 4, 3, 1, 2, 1, 1, 1, 1] #these are the d-values #It doesn't seem to make much sense. There may #be a problem since keeping a small number of digits #misleadingly makes two numbers agree; for example, #if we set Digits:=1, then evalf(Pi-3)=0. #I think we need better criteria for determining #when the two computed values "agree." ##########Problem 3: Counting Inversions########## #inv1(p) counts inversions in the list p. inv1:=proc(p) local n,i,j,N: n:=nops(p): i:=1: j:=1: N:=0: for i from 1 to n-1 do for j from i+1 to n do if p[i]>p[j] then N:=N+1: fi: od: od: N: end: #gf(n,q) sums q^inv1(p) over all perms p of length n. gf:=proc(n,q) local P,p: P:=permute(n): add(q^inv1(p),p in P): end: #If we factor gf(n,x) for n from 1 to 8, we get #the following: (* for n from 1 to 8 do factor(gf(n,x)); od; 1 1+x (1+x)*(x^2+x+1) (x^2+1)*(x^2+x+1)*(1+x)^2 (x^2+x+1)*(x^4+x^3+x^2+x+1)*(x^2+1)*(1+x)^2 (x^2-x+1)*(x^4+x^3+x^2+x+1)*(x^2+1)*(x^2+x+1)^2*(1+x)^3 (x^2-x+1)*(x^4+x^3+x^2+x+1)*(x^6+x^5+x^4+x^3+x^2+x+1)*(x^2+1)*(x^2+x+1)^2*(1+x)^3 (x^2-x+1)*(x^4+x^3+x^2+x+1)*(x^4+1)*(x^6+x^5+x^4+x^3+x^2+x+1)*(x^2+1)^2*(x^2+x+1)^2*(1+x)^4 *) #There seems to be a pattern in the first 3 lines, #but then things get weird; perhaps we are overfactoring #and destroying the pattern. Let's make a function #gfGuess(n,q) that outputs the expected formula #and compare this with the actual gf: gfGuess:=proc(n,q) local i: mul(sum(q^j,j=0..i),i=0..n-1): end: (* for n from 1 to 8 do expand(gf(n,x)-gfGuess(n,x)); od; 0 0 0 0 0 0 0 0 *) #It works! Apparently our guess is indeed correct. #As for rigorously proving this..... #.....left to the reader!