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