#C27.txt, April 30, 2015, the Metropolis algorithm Help:=proc(): print(` raR() , Nick(L,i0,K) `): print(`CheckNic(L,i0,K1,K2), Energy(M), RandU(m,n) `): print(` Gibbs(M,z) `): end: raR:=proc(): evalf(rand(1..10^10)()/10^10): end: Nick:=proc(L,i0,K) local ra,n,Pa,i,j,roll: n:=nops(L): ra:=rand(1..n): Pa:=[i0]: i:=i0: while nops(Pa)<=K do j:=ra(): if L[j]/L[i]>=1 then Pa:=[op(Pa),j]: i:=j: else roll:=raR(): if L[j]/L[i]>roll then Pa:=[op(Pa),j]: i:=j: else Pa:=[op(Pa),i]: fi: fi: od: Pa: end: #CheckNic(L,i0,K1,K2): CheckNick:=proc(L,i0,K1,K2) local x,Pa,i1,Ka: Pa:=Nick(L,i0,K1+K2): Ka:=evalf(add(x[Pa[i1]],i1=K1+1..K1+K2)/K2): [seq(coeff(Ka,x[i1],1),i1=1..nops(L))],evalf(L) : end: #The "energy" of a (-1,1) matrix M given as a list of lists Energy:=proc(M) local i,j,n,m: n:=nops(M): m:=nops(M[1]): add(add(M[i][j]*M[i][j+1],j=1..m-1)+ M[i][m]*M[i][1],i=1..n) +add(add(M[i+1][j]*M[i][j],i=1..n-1)+M[n][j]*M[1][j],j=1..m): end: #RandU(m,n): a (uniformly) m by n matrix of (1,-1) RandU:=proc(m,n) local ra,i,j: ra:=rand(0..1): [seq([seq(2*ra()-1,i=1..n)],j=1..m)]: end: #The Boltzmann weight of M with parameter z #z=exp(1/k/T) (T abs. temperature, and k physical const.) Gibbs:=proc(M,z) z^Energy(M): end: