Source Code for Program "Ising" (Mathematica version)

For use with An Introduction to Thermal Physics by Daniel V. Schroeder.

This code is written in Mathematica (version 6 or later). You can copy and paste each code block into a Mathematica notebook, or download the whole thing in notebook format here.

Mathematica code is marvelously concise and easy to modify for interactive explorations. On the other hand, it cannot do continuous animation, and the software is expensive to purchase and keep up to date. If you don't like (or don't have) Mathematica, click here for a list of versions of the Ising program in other languages.

The first block of code sets the size of the lattice, initializes it to a random array, defines a function to plot the array in pretty colors, and then plots the initial array:

size =50;
Clear[s];
s = Table[If[Random[Integer]==0, 1, -1], {i,1,size}, {j,1,size}];
plotit := ListDensityPlot[s, Mesh->False, InterpolationOrder->0, PlotRange->{-1,1},
		ColorFunction->(RGBColor[1-#,1-#,#]&)];
plotit

Next, define a function to compute the change in energy caused by a hypothetical flip of the dipole at location (i,j), using periodic boundary conditions:

deltaU[i_, j_] := 2*s[[i, j]]*(
        If[i == 1, s[[size, j]], s[[i-1, j]]]
      + If[i == size, s[[1, j]], s[[i+1, j]]]
      + If[j == 1, s[[i, size]], s[[i, j-1]]]
      + If[j == size, s[[i, 1]], s[[i, j+1]]])

The following function executes the Metropolis algorithm at the specified temperature for any desired number of steps, then plots the array to show the results.

runawhile[T_, stepcount_] := (
  For[step = 0, step < stepcount, step++,
    i = RandomInteger[{1, size}];
    j = RandomInteger[{1, size}];
    Ediff = deltaU[i, j];
    If[Ediff <= 0, s[[i, j]] = -s[[i, j]],
        If[ RandomReal[] < Exp[-Ediff/T], s[[i, j]] = -s[[i, j]] ]
        ]
    ];
  plotit)

For example, this invocation executes 100 steps per lattice site, at the critical temperature:

runawhile[2.27, 100*size^2]

The following instruction produces an animated sequence of plots, with one step per lattice site between successive plots:

ListAnimate[Table[runawhile[2.27, size^2], {plotnum, 1, 50}]]

Last modified on January 20, 2013.