Help:=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: