Page 111

```Matlab: Dice Rolls via Markov Chain Monte Carlo

f=[0, 1, 2, 3, 4, 5, 6, 5, 4, 3, 2, 1];
d=zeros(1,20000); % for histogramming
x=5;
for i = 1:20000
U = rand;
if x == 2
if U  <  0.5
y = 3;
else
y = 2;
end
elseif x == 12
if U  <  0.5
y = 11;
else
y = 12;
end
else
if U  <  0.5
y = x-1;
else
y = x+1;
end
end
h = min(1,f(y)/f(x));
U = rand;
if U  <  h
x = y;
end
d(i) = x; % record the state
end
a = 1:1:12;
hist(d,a)
[N,h] = hist(d,a)
```

Page 122

```Java: Ising Simulation

import java.io.*;
import java.util.*;

public class isingModel
{

static String[] helptext =
{
"Usage:\n\t < invoke program>",
""};

int nTrials = 9900;
int startUp = 3000;
int width = 16, hghth = 16;
double beta = .01; //beta=1/(kT), beta=1 low temp, .06 high temp

static Random r;
static long seed;

//ccc===
public
isingModel()
{
double magMoment = 0;

//create a random starting +1/-1 vector
isingCritter wCrit = new isingCritter(width,hghth,1,0,'r');
isingCritter tCrit = new isingCritter(width,hghth,1,0,'r');
isingCritter x, y; //pointers

int k, row, col;
double a=0;

//main loop
x = tCrit;  //random start
for(int i=0; i <  nTrials+startUp; ++i)
{
//x.printSpins();
//x = current, y = proposed
if( x == tCrit )
y = wCrit;
else
y = tCrit;

//select site to flip
k = (int)(x.sigma.length*r.nextDouble());
for( int j=0; j <  x.sigma.length; ++j ) y.sigma[j] = x.sigma[j];
y.sigma[k] *= -1;  //flip bit k
y.calcHamiltonian();

//System.out.format("Delta E= %.1f\n",y.E - x.E);
a = Math.min(1,Math.exp(-beta*(y.E - x.E)));

//decide to accept or not
if( r.nextDouble()  <  a ) x = y;

//update the magnetic moment
if( i > startUp )
{
magMoment += (double)x.calcMagneticMoment();
}
}
x.printSpins();

magMoment /= nTrials;
System.out.format("avg mag. moment per spin= %.3f\n",
magMoment/(width*hghth));
}
/*=====mmm=====*/
public static void
main(String[] argv)
{

seed = System.currentTimeMillis();
//seed = 17;  // remove // for testing
r = new Random(seed);

isingModel im = new isingModel();
}
}
//kkk===
class isingCritter
{
int width, hghth;
byte[] sigma; //up = +1, down = -1
double J, H;  //J=spin-spin interaction, H=external field
double E; //hamiltonian (= energy)
char boundaryType = 'p'; //p=periodic, f=free
int M; //magnetic moment

//ccc===
public
isingCritter(int w, int h)
{
hghth = h; width = w;
sigma = new byte[w*h];
J = 1; H = 0;

randomInit();
calcHamiltonian();
}
//ccc===
public
isingCritter(int w, int h, double J, double H, char initType)
{
hghth = h; width = w;
sigma = new byte[w*h];
this.J = J; this.H = H;

if( initType == 'u' )
spinUpInit();
else
randomInit();
calcHamiltonian();
}
//===
void
randomInit()
{
for( int i=0; i < sigma.length; ++i)
{
if( isingModel.r.nextDouble()  <  .5 )
sigma[i] = 1;
else
sigma[i] = -1;
}
}
//===
void
spinUpInit()
{
for( int i=0; i < sigma.length; ++i)
{
sigma[i] = 1;
}
}
//===
/**
* Here we use the free boundary condition, i.e. the spins
* on the boundary are coupled only to spins of the interior
* and not to spins on the opposite end of the lattice
*/
void
calcHamiltonian()
{
int r,c;
//printSpins();
E=0;
/**
* free boundary calculation
*/
//go along the columns first between each row
for( r=0; r < hghth-1; ++r)
{
for( c=0; c < width; ++c)
{
E += -J*(sigma[c+r*width]*sigma[c+(r+1)*width]);
}
}

//now go along the rows between each column
for( c=0; c < width-1; ++c)
{
for( r=0; r < hghth; ++r)
{
E += -J*(sigma[c+r*width]*sigma[c+1+r*width]);
}
}

/**
* if boundary condition is periodic, do this too
*/
if( boundaryType == 'p' )
{
r = (hghth-1)*width;
for( c=0; c < width; ++c)
E += -J*sigma[c]*sigma[c+r];

c= width-1;
for( r=0; r < width*hghth; r += width)
E += -J*sigma[r]*sigma[c+r];
}

/**
* now figure in the external field
* assume aligned with 'up'
*/
if( H != 0)
{
for( int i=0; i < sigma.length; ++i)
E += -H*sigma[i];
}
//System.out.println("total: E= "+E+"\n");
}
//===
/**
* M(x) = net up vs down spins in absence of external field
*/
int
calcMagneticMoment()
{
M = 0;
for( int i=0; i < sigma.length; ++i)
M += sigma[i];
//System.out.println("magnet moment: M= "+M+"\n");
return M;
}
//===
void
printSpins()
{
int j=0;
for( int k=0; k < sigma.length; ++k)
{
System.out.format("%c ",(sigma[k]==1)?'+':'-');
//System.out.format("%2d ",sigma[k]);
if( ++j==width )
{
System.out.println("");
j=0;
}
}
System.out.println("");
}
}
```

Page 125

```Matlab: Ising Simulation

nTrials=2000; startup=3000;
d=zeros(1,nTrials); %to histogram
beta=.01; M=0;
x=2*round(rand(1,4))-1; %random 4-vector of +/-1
Hx = -(x(1)*x(2)+x(3)*x(4)+x(1)*x(3)+x(2)*x(4));
for t = 1:startup+nTrials
k = fix(1+4*rand); %index to flip
y = x; y(k)= -x(k);
Hy = -(y(1)*y(2)+y(3)*y(4)+y(1)*y(3)+y(2)*y(4));
h = min(1,exp(-beta*(Hy-Hx)));
if rand  <  h
x = y;
end
if t  >  startup
Mx = x(1)+x(2)+x(3)+x(4);
M = M + Mx; % to obtain avg mag. mom.
d(t-startup) = Mx;
end
end
r = -4:1:4;
hist(d,r)
% compute avg moment per spin
M = (M/nTrials)/4;
```

Page 128

```Matlab: Metropolis Algorithm for Counting, Monomer/Dimer Problem

% this pertains to figure 3.8; in the figure the vertices are "a"
% through "f" on the left and "1" through "6" on the right but here
% we use 1 through 6 for the left vertices and 7 through 12 for the
% right ones.  E(i,j)=k means this: j designates the edge so j runs
% from 1 to 14, i=1 is for the left vertex of the edge and i=2 is
% for the right vertex and k is which vertex in the figure. Edge 1
% runs from "a" to "1" so E(1,1)=1 and E(2,1)=7. Edge 2 runs from "a"
% to "3" so E(1,2)=1 and E(2,2)=9 and so on.
% initialize the first two rows of E to the given graph, a 3x14 matrix
% initialize third row to 0, start with an empty match set
E=zeros(3,14);
E(1,1)=1; E(2,1)=7; E(1,2)=1; E(2,2)=9;
E(1,3)=1; E(2,3)=11; E(1,4)=2; E(2,4)=8;
E(1,5)=2; E(2,5)=11; E(1,6)=3; E(2,6)=10;
E(1,7)=3; E(2,7)=12; E(1,8)=4; E(2,8)=7;
E(1,9)=4; E(2,9)=9; E(1,10)=4; E(2,10)=12;
E(1,11)=5; E(2,11)=11; E(1,12)=5; E(2,12)=12;
E(1,13)=6; E(2,13)=8; E(1,14)=6; E(2,14)=10;
s=zeros(1,14); % array to count # matches of size k
s0=0; % need to get size 0 too
k=0; % the initial match set is empty
for t=1:990000 % 990000 trials
i=1+floor(14*rand); % select edge unif. at random
if E(3,i) == 1 %check if i is in the current match set
E(3,i)=0; % remove it
k=k - 1; % decrement matching size
else
l=E(1,i); r=E(2,i); % left, right ends of i
kl=0; kr=0; %search current match set for l and r
for j=1:14 % loop over all edges
if E(3,j) == 1 %this edge is in the match set
if l == E(1,j) | l == E(2,j) % l is covered
kl=j;
end
if r == E(1,j) | r == E(2,j) % r is covered
kr == j;
end
end % if j in match set
end %loop over edges
if kl == 0 & kr == 0 % neither l or r covered
E(3,i)=1; %add edge i to match
k=k+1; % increment match set size
elseif kl ~= 0 & kr == 0 % only left end covered
E(3,kl)=0; E(3,i)=1; % swap edges
elseif kl == 0 & kr ~= 0 % only right end covered
E(3,kr)=0; E(3,i)=1; % swap edges
end
end % check i in current match set
% increment the count for the current size
if k == 0
s0=s0+1;
else
s(k)=s(k)+1;
end
end % loop over trials
s0 % print results
s
```