## 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.*