Friday, September 28, 2007

Ising model simulation

This posting shows how to use Mathematica (version 6) to create an interactive simulation of a 2-dimensional Ising model, which is based on some work on Boltzmann machines that I did back in the mid 1980s. This nicely illustrates the use of the Manipulate function in Mathematica, which is by far the easiest way of creating interactive graphical demonstrations that I have ever used.

Define a function for randomly initialising the state of the 2-dimensional Ising model.

init[narray_] := array = RandomInteger[{0,1}, {narray,narray}];

Define a function for doing a single Monte Carlo update of the 2-dimensional Ising model:

update[{pe_,pne_,pn_,pnw_}, narray_] :=
Module[
{probpatch={{pnw,pn,pne}, {pe,0.5,pe}, {pne,pn,pnw}}, i, j, di, dj, arraypatch, freq1, freq0},
{i,j} = RandomInteger[{1,narray}, {2}];
arraypatch = Table[array[[Mod[i+di,narray,1],Mod[j+dj,narray,1]]], {di,-1,1}, {dj,-1,1}];
freq1 = Apply[Times, probpatch arraypatch+(1-probpatch)(1-arraypatch), {0,1}];
freq0 = Apply[Times, (1-probpatch)arraypatch+probpatch(1-arraypatch), {0,1}];
array[[i,j]] = If[Random[] < freq1/(freq0+freq1), 1, 0];
];


Define a function for implementing the interactive graphical demonstration.

Manipulate[
If[narray!=Length[array], init[narray], ];
If[u, Do[update[{p2en[[1]],p2nenw[[1]],p2en[[2]],p2nenw[[2]]}, narray], {nupdate}], ];
ArrayPlot[array],
{{p2en, {0.5,0.5}, "probability(\[RightArrow],\[UpArrow])"}, {e, e}, {1-e, 1-e}, ControlPlacement->Left},
{{p2nenw, {0.5,0.5}, "probability(\[UpperRightArrow],\[UpperLeftArrow])"}, {e, e}, {1-e, 1-e}, ControlPlacement->Left},
{{u,"False","update"}, {True,False}, ControlPlacement->Top},
{{narray,32,"size"}, 3, 64, 1, Appearance->"Labeled", ControlPlacement->Top},
FrameLabel->"Ising Model",
Initialization->(e=0.001; array={}; nupdate=50)
]


When Manipulate is evaluated the initial output typically looks like this:


Use the check box at the top of the output window to switch the simulation on and off. The number of Monte Carlo updates that is applied between the display of each "frame" of the simulation is hardwired to be 50.

The initial size of the 2-dimensional array of "spins" is 32 by 32, which is controlled by the 1D slider at the top of the output window. This slider can be moved at any time, even whilst a simulation is being run.

Each Monte Carlo update is implemented by selecting a spin at random, and calculating 8 probability factors due to its interaction with each of its 8 surrounding spins. Each such factor specifies a multiplicative contribution to the relative likelihood that the central spin has the same/different value compared to its neighbouring spin (this is not the most general interaction that is possible). There are 4 independent probability factors corresponding to the following directions from the central spin: east, north, north-east, and north-west. These probability factors are controlled by the two 2D sliders on the left of the output window. The top 2D slider controls the east (left/right slider movement) and north (up/down slider movement) factors, and analogously the bottom 2D slider controls the north-east and north-west factors.

Denote the probability factor for a pair of spins being the same as p, and for being different as 1-p. The central position of each 2D slider corresponds to a probability factor p=1/2 (i.e. the same contribution whether the spins are the same or different), the left and bottom edges of each 2D slider correspond to p=0 (i.e. force the spins to be different), and the right and top edges correspond to p=1 (i.e. force the spins to be the same). These sliders can be moved at any time, even whilst a simulation is being run.

So, using a 1-dimensional example, the spin configuration 0?0 gives a product of probability factors (1-p)^2 for 010 and p^2 for 000, whereas 1?0 gives p(1-p) for 110 and (1-p)p for 100 (i.e. the same in each case). A Monte Carlo update of the central spin selects ? as being 1 or 0 with a relative probability given by the ratio of the corresponding products of probability factors. This ratio is ((1-p)/p)^2 and 1 in the two examples above, respectively. The 2-dimensional case is a straightforward generalisation of the 1-dimensional case.

Each snapshot of the output below shows the typical result of running a randomly initialised simulation for a few seconds with the 2D slider settings as shown in the snapshot.





Note that extreme values are used for the probability factors in each of these simulations, and because the simulations each have a short duration they do not actually reach equilibrium. For instance, the first simulation enforces such strong positive correlations between adjacent spins that it should eventually degenerate to all the spins having the same colour.

Enjoy playing with this Ising model simulation. It is quite interesting to move the 2D sliders to vary the probability factors as the simulation is running, because the speed of the simulation is sufficiently fast that you get an almost real-time response as the Ising model dynamically adjusts its equilibrium state.

2 comments:

Anonymous said...

Great post. I just started learning simulation in economics and came across this post. Interesting and helpful. Thanks.

Stephen Luttrell said...

Thank you for your support. I am still experimenting with the best way to publish useful Mathematica snippets. At the moment you have to cut & paste my code into your own copy of Mathematica, but in due course I hope to be able to publish standalone programs. This will have to wait for Wolfram Research to start selling the means to do this sort of thing.