Help:=proc(): print(` Wt(A,z,w) , GenM(m,n,z,w,K), OneStep(m,n,A,z,w)`): end: HelpOld:=proc(): print(`En(A), PF(m,n,z) , RealAve1(m,n,z,w,p,q)`): print(`RealAve2(m,n,z,w,p1,q1,p2,q2), Cor(m,n,z,w,p1,q1,p2,q2) `): end: #En(A): the sum of the products of all nearest-neighbors of #the matrix A given as a list of lists En:=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-1][j]*A[i][j],i=2..nops(A)),j=1..nops(A[1])): end: #PF(m,n,z) PF:=proc(m,n,z) local s,c,ab: s:=0: for c from 0 to 2^(m*n)-1 do ab:=convert(c,base,2): ab:=[op(ab),0\$(m*n-nops(ab))]: ab:=subs(0=-1,ab): ab:=[seq([op((i-1)*n+1..i*n,ab)],i=1..m)]: s:=s+z^En(ab): od: s: end: #Ex: evaluate f(z)=lim (m,n->infinity) PF(m,n,z)^(1/(m*n)) #(equivalenty: limit of log(PF(m,n,z))/(m*n) ) #gPF(m,n,z,w) gPF:=proc(m,n,z,w) local s,c,ab,ab1: s:=0: for c from 0 to 2^(m*n)-1 do ab:=convert(c,base,2): ab:=[op(ab),0\$(m*n-nops(ab))]: ab:=subs(0=-1,ab): ab1:=[seq([op((i-1)*n+1..i*n,ab)],i=1..m)]: s:=s+z^En(ab1)*w^convert(ab,`+`): od: s: end: #Ex: evaluate g(z,w)=lim (m,n->infinity) PF(m,n,z,w)^(1/(m*n)) #(equivalenty: limit of log(PF(m,n,z,w))/(m*n) ) #RealAve1(m,n,z,w,p,q): the (weighted) average of A[i][j] #amongsts all m by n matrices with {-1,1} with Boltzmann weight RealAve1:=proc(m,n,z,w,p,q) local t,s,ab,ab1,i,j,c: s:=0: for c from 0 to 2^(m*n)-1 do ab:=convert(c,base,2): ab:=[op(ab),0\$(m*n-nops(ab))]: ab:=subs(0=-1,ab): ab1:=[seq([op((i-1)*n+1..i*n,ab)],i=1..m)]: s:=s+z^En(ab1)*w^convert(ab,`+`)*ab1[p][q]: od: s/gPF(m,n,z,w): end: #RealAve2(m,n,z,w,p1,q1,p2,q2): the (weighted) average of A[p1][q1]*A[p2][q2] #amongsts all m by n matrices with {-1,1} with Boltzmann weight RealAve2:=proc(m,n,z,w,p1,q1,p2,q2) local t,s,ab,ab1,i,j,c: s:=0: for c from 0 to 2^(m*n)-1 do ab:=convert(c,base,2): ab:=[op(ab),0\$(m*n-nops(ab))]: ab:=subs(0=-1,ab): ab1:=[seq([op((i-1)*n+1..i*n,ab)],i=1..m)]: s:=s+z^En(ab1)*w^convert(ab,`+`)*ab1[p1][q1]*ab1[p2][q2]: od: s/gPF(m,n,z,w): end: #E(XY)-E(X)E(Y)=E((X-E(X))(Y-E(Y)) Cor:=proc(m,n,z,w,p1,q1,p2,q2): RealAve2(m,n,z,w,p1,q1,p2,q2)- RealAve1(m,n,z,w,p1,q1)*RealAve1(m,n,z,w,p2,q2): end: #Wt(A,z,w) Wt:=proc(A,z,w) local s,c,ab,ab1: z^En(A)*w^convert(convert(A,`+`),`+`) end: #OneStep(m,n,A,z,w): performs one step in the Metropolis arlgorithm #applied to the Ising model with parameters z and w (importatnce given #by z^En(A)*w^convert(convert(A,`+`),`+`) OneStep:=proc(m,n,A,z,w) local i,j,A1,rat,dennis: i:=rand(1..m)(): j:=rand(1..n)(): A1:=A: A1[i][j]:=-A[i][j]: rat:=Wt(A1,z,w)/Wt(A,z,w): if rat>=1 then RETURN(A1): fi: dennis:=stats[random,uniform[0,1]](1): if dennis