For use with An Introduction to Thermal Physics by Daniel V. Schroeder.
This code is written in Mathematica. You can copy and paste it into a Mathematica notebook, or click here to download a copy in notebook format.
Besides being slower than other languages, Mathematica has the disadvantage of not providing an easy way to change the color of one square without redrawing the whole plot. On the other hand, it is an interactive language in which you can easily try executing any number of steps between successive plots.
If you don't like Mathematica, here are versions of this program written in True Basic and Java. There's also a much more elaborate version written in REALbasic, which you can download (as source code or a compiled application) from here.
size =50;
Clear[s];
s = Table[If[Random[Integer]==0, 1, -1], {i,1,size}, {j,1,size}];
plotit := ListDensityPlot[s, Mesh->False, 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 = Random[Integer, {1, size}];
j = Random[Integer, {1, size}];
Ediff = deltaU[i, j];
If[Ediff <= 0, s[[i, j]] = -s[[i, j]],
If[ Random[] < Exp[-Ediff/T], s[[i, j]] = -s[[i, j]] ]
]
];
plotit;)
For example, this invokation executes 100 steps per lattice site,
at the critical temperature. (Slow, you say? Then why are
you using Mathematica?)
runawhile[2.27, 100*size^2]The following instruction produces a sequence of plots, with one step per lattice site between successive plots. To animate the sequence, click on the bar at the right that groups them together, then choose "Animate Selected Graphics" from the Cell menu.
For[plotnum = 1, plotnum <= 20, plotnum++, runawhile[2.27, size^2]]
Last modified on January 22, 2003.