#March 12, 2020 "Virtual Class". Corrected version of March 13, 2020 #C15.txt: Using Metropolis to estimate macroscopic quantities (i.e. expectations) of physical quantities #and comparing them for the actual values for small lattices #The new procedures for C15.txt are listed in Help(). To see the procedures done in Classes 13 and 14 #type: Help13(); and Help14(); Help:=proc(): print(` NickAv(M,N,x,y,K1,K2,G,Eng,Mag), ExactAv(M,N,x,y,G,En,Ma) , ExactAvBrute(M,N,x,y,G,En,Ma) `): end: #ExactAvBrute(M,N,x,y,G,En,Ma): Like ExactAv(M,N,x,y,G,En,Ma) but by brute force ExactAvBrute:=proc(M,N,x,y,G,En,Ma) local S,tot,su, A: S:=Mats(M,N): tot:=0: su:=0: for A in S do tot:=evalf(tot+x^Ene(A)*y^Mag(A)): su:=evalf(su+x^Ene(A)*y^Mag(A)*subs({En=Ene(A),Ma=Mag(A)},G)): od: su/tot: end: #ExactAv(M,N,x,y,G,En,Ma): Corrected version of March 13, 2020 #Inputs positive integers M and N and real numbers x and y, positive integers K1 and K2 #and an expression G that is a MONOMIAL in the symbols E and M (denoting energy and magnetization, respectively), outputs the #the EXACT average of G(En,Ma) over all 2^(M*N) possible M by N {1,-1} matrices with the Ising prob. distribution with #temperature parameter x and external field parameter y. It uses "generating-functionlogy". If the #Prob. gen. fun of a r.v. is f(z) then its expectation is (z*d/dz f(z))|z=1 and its r-th moment is #(zd/z)^r) f(z) |z=1 , similarly for more than variable #For example, to the exact average of the energy with x=1.1 and y=1 on the 4 by 4 rectangular torus do: #ExactAv(4,4, 1.1, 1.0, En,En,Ma). #To get the exact average of En^2*Ma^2 do #ExactAv(4,4, 1.1, 1.0, En^2*Ma^2,En,Ma). ExactAv:=proc(M,N,x,y,G,En,Ma) local d1,d2,C,i,f,g,X,Y: #We check that G is indeed a monomial in En and Ma, i.e. is of the form C*En^d1*Ma^d2 for some number C d1:=degree(G,En): d2:=degree(G,Ma): C:=normal(G/(En^d1*Ma^d2)): if not type(C,numeric) then print(G, `is not a monomial in`, En, Ma ): RETURN(FAIL): fi: #Introduce FORMAL variables X and Y to represent the energy and magnetization #this is needed for the "generating-functionlogy" part, we need x and y to be numeric #to reflect the probabilty distribution #f is the exact expression in X and Y for the so-called (unnormalized) partition function for the M by N torodial rectangle f:=ParF(M,N,X,Y): #Let's normalize f f:=f/subs({X=1,Y=1},f): #Now let's use the numbers x and y to make it a numeric prob. generating function where the coeff. of the monomial X^a*Y^b #is the exact probability that the energy of the spin configuration is a and its magnetization is b #f:=subs({X=x*X,y=y*Y},f): THIS WAS A BUG SPOTTED BY ROBERT DOUGHERTY-BLISS y=y*Y should be, of course, Y=y*Y. It is now corrected f:=subs({X=x*X,Y=y*Y},f): THIS WAS A BUG SPOTTED BY ROBERT DOUGHERTY-BLISS #Now let's have a copy of f, let's call it g g:=f: #Now let's apply (X*d/dX) d1 times followed by (Y*d/dY) d2 times for i from 1 to d1 do g:=expand(X*diff(g,X)): od: for i from 1 to d2 do g:=expand(Y*diff(g,Y)): od: #now let's plug-in X=1, Y=1 to get a number that only depends on the physical parameters x and y. #Note that x and y can be left symbolic, if desired (but for comparing it to the output of NickAv, they have to be numeric) #we also have to multiply by the coefficient of the monomial #and divide by the numerator C*subs({X=1,Y=1},g)/subs({X=1,Y=1},f): end: #####Old erroneous code #ExactAvNoGood(M,N,x,y,G,En,Ma): Inputs positive integers M and N and real numbers x and y, positive integers K1 and K2 #and an expression G that is a MONOMIAL in the symbols E and M (denoting energy and magnetization, respectively), outputs the #the EXACT average of G(En,Ma) over all 2^(M*N) possible M by N {1,-1} matrices with the Ising prob. distribution with #temperature parameter x and external field parameter y. It uses "generating-functionlogy". If the #Prob. gen. fun of a r.v. is f(z) then its expectation is (z*d/dz f(z))|z=1 and its r-th moment is #(zd/z)^r) f(z) |z=1 , similarly for more than variable #For example, to the exact average of the energy with x=1.1 and y=1 on the 4 by 4 rectangular torus do: #ExactAv(4,4, 1.1, 1.0, En,En,Ma). #To get the exact average of En^2*Ma^2 do #ExactAvNoGood(4,4, 1.1, 1.0, En^2*Ma^2,En,Ma). ExactAvNoGood:=proc(M,N,x,y,G,En,Ma) local d1,d2,C,i,f,X,Y: #We check that G is indeed a monomial in En and Ma, i.e. is of the form C*En^d1*Ma^d2 for some number C d1:=degree(G,En): d2:=degree(G,Ma): C:=normal(G/(En^d1*Ma^d2)): if not type(C,numeric) then print(G, `is not a monomial in`, En, Ma ): RETURN(FAIL): fi: #Introduce FORMAL variables X and Y to represent the energy and magnetization #this is needed for the "generating-functionlogy" part, we need x and y to be numeric #to reflect the probabilty distribution #f is the exact expression in X and Y for the so-called (unnormalized) partition function for the M by N torodial rectangle f:=ParF(M,N,X,Y): #Let's normalize f f:=f/subs({X=1,Y=1},f): #Now let's use the numbers x and y to make it a numeric prob. generating function where the coeff. of the monomial X^a*Y^b #is the exact probability that the energy of the spin configuration is a and its magnetization is b f:=subs({X=x*X,y=y*Y},f): #Now let's apply (X*d/dX) d1 times followed by (Y*d/dY) d2 times for i from 1 to d1 do f:=expand(X*diff(f,X)): od: for i from 1 to d2 do f:=expand(Y*diff(f,Y)): od: #now let's plug-in X=1, Y=1 to get a number that only depends on the physical parameters x and y. #Note that x and y can be left symbolic, if desired (but for comparing it to the output of NickAv, they have to be numeric) #we also have to multiply by the coefficient of the monomial C*expand(subs({X=1,Y=1},f)): end: #NickAv(M,N,x,y,K1,K2,G,En,Ma): Inputs positive integers M and N and real numbers x and y, positive integers K1 and K2 #and an expression G in the symbols E and M (denoting energy and magnetization, respectively), outputs the #estimate produced by running the Metropolis algorithm K1 in the "warm-up" stage (w/o taking records) #followed by K2 iterations where a record of G(E,M) is taken, and added to the average. #For example, to get an estimate for the average of the energy-squared, try: #NickAv(4,4,2. ,1.1, 100,1000, E^2, E,M); NickAv:=proc(M,N,x,y,K1,K2,G,En,Ma) local A,i,su: A:=RM(M,N): for i from 1 to K1 do A:=OSM(A,M,N,x,y): od: su:=0: for i from 1 to K2 do A:=OSM(A,M,N,x,y): su:=su+subs({En=Ene(A),Ma=Mag(A)},G): od: evalf(su/K2): end: with(Statistics): ###stuff from C13.txt Help13:=proc(): print(`RM(M,N), PlotMat(A), Ene(A) , Mag(A)`): print(`Prob(A,x,y) , Vecs(N) , Mats(M,N) , ParF(M,N,x,y)`): end: RM:=proc(M,N) local i,j,ra: ra:=rand(0..1): subs(0=-1,[seq([seq(ra(),j=1..N)],i=1..M)]): end: with(plots): #PlotMat(A): plots the matrix A where 1 is black and -1 is white PlotMat:=proc(A) local i,j,d: if A[1][1]=1 then d:=pointplot([[1,1]],axes=none,color=black,symbol=circle): else d:=pointplot([[1,1]],axes=none,color=red,symbol=circle): fi: for j from 2 to nops(A[1]) do if A[1][j]=1 then d:=d,pointplot([[1,j]],color=black,symbol=circle): else d:=d,pointplot([[1,j]],color=red,symbol=circle): fi: od: for i from 2 to nops(A) do for j from 1 to nops(A[i]) do if A[i][j]=1 then d:=d,pointplot([[i,j]],color=black,symbol=circle): else d:=d,pointplot([[i,j]],color=red,symbol=circle): fi: od: od: display(d); end: #Ene(A): The energy of the matrix A of {1,-1} Ene:=proc(A) local i,j: add(add( A[i][j]*A[i][j+1] , j=1..nops(A[i])-1) ,i=1..nops(A))+ add(add( A[i][j]*A[i+1][j] , j=1..nops(A[1]) ) ,i=1..nops(A)-1)+ add(A[1][i]*A[-1][i],i=1..nops(A[1]))+ add(A[i][1]*A[i][-1],i=1..nops(A)): end: #Mag(A): the magnetization of the spin-config. A Mag:=proc(A): add(add(A)): end: #Vecs(N): the set of all 2^N {1,-1} lists of length N Vecs:=proc(N) local S,v: option remember: if N=0 then RETURN({[]}): fi: S:=Vecs(N-1): {seq([op(v),1], v in S),seq([op(v),-1], v in S)}: end: #Mats(M,N): the set of all M by N {1,-1} matrices Mats:=proc(M,N) local S,S1,A,v: if M=0 then RETURN({[]}): fi: S:=Mats(M-1,N): S1:=Vecs(N): {seq(seq([op(A),v],v in S1),A in S)}: end: #x=exp(k/T) (T=abs. temp. (in Kelvin)) y=exp(External Magnetication) #x^Ene*y^Mag #ProbA(A,x,y):Unnormalied prob. ProbA:=proc(A,x,y): x^Ene(A)*y^Mag(A):end: #ParF(M,N,x,y): the Gibbs partition function of M by N Ising model ParF:=proc(M,N,x,y) local S,A: S:=Mats(M,N): add(ProbA(A,x,y),A in S): end: ###end from C13.txt ###start stuff from C14.txt Help14:=proc(): print(` Neis(c,M,N) , NewA(A,c), OSM(A,M,N,x,y), Nick(M,N,x,y,K1)`): end: #Neis(c,M,N): Inputs a location c=[i,j] and pos. integers M and N outputs the SET #of four neighbors of the cell [i,j], {[i \pm 1,j ], [i, j \pm 1] but with toral convention #CODE by Alison Bu and Yukun Yao, with improvements by Charles Kenney Neis:=proc(c,M,N) local i,j: i:=c[1]: j:=c[2]: {[(i-2 mod M)+1,j ],[(i mod M)+1, j ],[i, (j-2 mod N)+1 ] ,[i,(j mod N)+1]}: end: #NewA(A,c): The result of flipping the spin at location c in the matrix A #followed by the change in Energy, and change in Magnetization NewA:=proc(A,c) local i,j,s, A1,Ne,DiffE, DiffM,a: i:=c[1]: j:=c[2]: s:=A[i][j]: A1:=[op(1..i-1,A), [op(1..j-1,A[i]), -s, op(j+1..nops(A[i]), A[i])], op(i+1..nops(A),A)]: Ne:=Neis(c,nops(A),nops(A[1])): DiffE:=-2*s*add(A[a[1]][a[2]], a in Ne): DiffM:=-2*s: [A1,DiffE,DiffM]: end: #OSM(A,M,N,x,y): one step in the Metropolis algorithm inputs an M by N matrix of {1,-1} #and pos. numbers x,y outputs one step in the Metropolis algorithm according to #the importance x^Ene(A)*y^Mag(A) OSM:=proc(A,M,N,x,y) local c,Zidong, A1,Rat,r,U: c:=[rand(1..M)(), rand(1..N)()]: Zidong:=NewA(A,c): A1:=Zidong[1]: Rat:=x^Zidong[2]*y^Zidong[3]: r:=Sample(RandomVariable(Uniform(0,1)),1)[1]: if r