Simulating Experiences
 Excursions in Programming

 The Ising Model

by
Richard J. Gaylord
Department of Materials Science and Engineering
University of Illinois
Champaign, Illinois 61801
gaylord@ux1.cso.uiuc.edu

and

Kazume Nishidate
Department of Electrical and Electronic Engineering
Iwate University
Morioka 020
JAPAN
nisidate@wrion1.s.kanazawa-u.ac.jp

Mathematica in Education
Volume 3 Number 1
Winter 1994
(c) TELOS/Springer-Verlag

Introduction

In statistical physics, the "Ising" system serves as one of the simplest models of interacting bodies. The spin flip dynamics version of the model has been used to model ferromagnets, anti-ferromagnetism and phase separation in binary alloys. It has also been applied to spin glasses and neural networks. Elsewhere, it has been suggested [Callen and Shapero, 1974] that the Ising model is relevant to imitative behavior in general, including such disparate systems as flying birds, swimming fish, flashing fireflies, beating heart cells, spreading diseases and even fashion fads. 

We'll look at the two-dimensional version of the Ising model. Essentially, the model is comprised of an 'n by n' square lattice in which each lattice site has associated with it a value of 1 (known as an up spin) or -1 (known as a down spin). Spins on adjacent, nearest-neighbor sites (henceforth called neighbor sites) interact in a pair-wise manner with a strength J (known as the exchange energy or constant). When J is positive, the energy is lower when spins are in the same direction and when J is negative, the energy is lower when spins are in opposite directions. There may also be an external field of strength B (known as the magnetic field). The magnetization of the system is the difference between the number of up and down spins on the lattice. In spin flipping, lattice sites are selected and either flipped or not, based on the energy change in the system that would result from the flip. We will look at two implementations of the model.

The Probabilistic Algorithm

The probabilistic Ising model corresponds to a constant temperature condition (in statistical physics, this is known as the canonical ensemble formulation). A lattice site is selected at random and a decision is made on whether or not to flip the site using the Metropolis Algorithm. The implementation of this algorithm is described as follows:

(1) Start with some initial 'n by n' lattice of spins. We can point all of the spins in the same direction or use a random distribution spin directions.

lat = 2 Table[Random[Integer], {n}, {n}] - 1

(2) Randomly select a lattice site.

{i1,i2} = {Random[Integer, {1, n}], Random[Integer, {1, n}]} 

(note that we can also select lattice sites by going through the lattice in a so-called 'typewriter' fashion, sequentially selecting suceeding sites). 

(3) Determine the energy change that is involved in 'flipping' the spin at the lattice site. This calculation is done in a number of steps:

(3a) determine the neighbor sites to the selected sites.

The neighbor sites are lat[[dn,i2]], lat[[up,i2]], lat[[i1,rt]] and lat[[i1,lt]] where 

lat[[dn,i2]] is the lattice site beneath the site
lat[[up,i2]] is the lattice site above the site
lat[[i1,rt]] is the lattice site to the right of the site
lat[[i1,lt]] is the lattice site to the left of the site

and the values of dn, up, rt and lt are determined by the location of the selected lattice site. 

If[i1 == n, dn = 1, dn = i1+1]; 
If[i2 == n, rt = 1, rt = i2+1];
If[i1 == 1, up = n, up = i1-1]; 
If[i2 == 1, lt = n, lt = i2-1];

When the selected site is in the interior of the lattice, the neighbor sites are {{i1,i2+1},{i1,i2-1},{i1+1,i2},{i1-1,i2}}. When the selected site is along the lattice boundary we take neighbor sites from the opposing side of the lattice .This method is known as a reflecting or periodic boundary condition.

(3b) calculate the energy change that would result from flipping the spin of the selected lattice site. This is given by

 2 lat[[i1,i2]] (B + J nnvalsum) 

where nnvalsum is the total spin from the neighbor sites

nnvalsum = 
   lat[[dn,i2]] + lat[[up,i2]] + lat[[i1,rt]] + lat[[i1, lt]]
   
The value of lat[[i1,i2]] is either 1 or -1 and the value of nnvalsum is either  4, 2, 0, -2, or -4 (4 corresponds to all spins up, 2 to three spins up and one spin down, 0 to two spins up and two spins down, -2 to three spins down and one spin up and -4 to all spins down). Thus, there are ten possible energy changes as a result of spin flipping.

Rather than spending time calculating the energy change for the selected site, we create a 'look-up' table of the possible energy changes 

energydiff[ 1, 4] =  2 (B + 4J);
energydiff[ 1, 2] =  2 (B + 2J);
energydiff[ 1, 0] =  2 (B + 0J);
energydiff[ 1,-2] =  2 (B - 2J);
energydiff[ 1,-4] =  2 (B - 4J);
energydiff[-1, 4] = -2 (B + 4J);
energydiff[-1, 2] = -2 (B + 2J);
energydiff[-1, 0] = -2 (B + 0J);
energydiff[-1,-2] = -2 (B - 2J);
energydiff[-1,-4] = -2 (B - 4J);
   
(note: by convention, a negative energydiff corresponds to a lowering of the lattice energy and a positive energydiff corresponds to a raising of the lattice energy).

(4) Decide if the spin of the selected lattice site should be flipped and if so, perform the flip. 

The flip operation, given by  

lat[[i1,i2]] = -lat[[i1,i2]] 

is performed if either of two conditions is met:

(4a) if the energy of the lattice is lowered as a result of the flip (ie., if energydiff is negative).

energydiff[lat[[i1,i2]],nnvalsum] < 0

(4b) if the energy of the lattice is raised as a result of the flip, but the exponential of (-energydiff) is greater than a random number between 0 and 1 

Random[] < Exp[-energydiff[lat[[i1,i2]],nnvalsum]

(note: the more the energy is raised by the flip, the more positive energydiff is and hence the smaller the exponential of (-energydiff) is and the less likely it is that it will be greater than Random[]).

The two-step process in 4a,b is known as the Metropolis Algorithm and it is expressed using 

If[energydiff[lat[[i1,i2]],nnvalsum] < 0 || 
       Random[] < Exp[-energydiff[lat[[i1,i2]],nnvalsum],
   lat[[i1,i2]] = -lat[[i1,i2]]] 


The Metroplis method allows the system to reach a 'global' energy minimum rather than getting stuck in a 'local' minimum. 

(5) Perform the sequence of steps 1-4 many times, building a list of the lattice configuration after each flip attempt.

flipLis = NestList[flip, initconfig, m];

(6) Create a sublist,monteCarloStepLis, containing every (n^2)th element from the list of lattice configurations. Each element in the sublist is the lattice configuration due to one 'monte carlo step' (the use of every (n^2)th element corresponds to giving each lattice site an equal chance to be selected in a single monte carlo step).

monteCarloStepLis =  flipLis[[Range[1,m,n^2]]]; 
     
(7) Calculate for each element in monteCarloStepLis, some global property of the lattice such as the long range order, which is the absolute value of the magnetization of the lattice

Abs[Apply[Plus, Flatten[lat]]/n^2]

and create a ListPlot of the long range order as a function of the number of monte carlo steps.

(8) Create an animation of the evolution of lattice configurations from monteCarloStepLis.

The Program

(* 
n = lattice size 
m = number of steps 
B = external field in 1/kT units 
J = interaction energy in 1/kT units 
*)

isingMetropolis[n_, m_, B_, J_] :=
 Module[{},

  energydiff[ 1, 4] =  2 (B + 4J);
  energydiff[ 1, 2] =  2 (B + 2J);
  energydiff[ 1, 0] =  2 (B + 0J);
  energydiff[ 1,-2] =  2 (B - 2J);
  energydiff[ 1,-4] =  2 (B - 4J);
  energydiff[-1, 4] = -2 (B + 4J);
  energydiff[-1, 2] = -2 (B + 2J);
  energydiff[-1, 0] = -2 (B + 0J);
  energydiff[-1,-2] = -2 (B - 2J);
  energydiff[-1,-4] = -2 (B - 4J);
  
  initconfig = 2 Table[Random[Integer],{n},{n}] - 1;
 
  
  flip = (
          lat = #;
          {i1, i2} = {Random[Integer,{1,n}],
                      Random[Integer,{1,n}]};
          If[i1 == n, dn = 1, dn = i1+1]; 
          If[i2 == n, rt = 1, rt = i2+1];
          If[i1 == 1, up = n, up = i1-1]; 
          If[i2 == 1, lt = n, lt = i2-1];
          nnvalsum = lat[[dn,i2]] + lat[[up,i2]] + 
                     lat[[i1,rt]] + lat[[i1,lt]];
          If[energydiff[lat[[i1,i2]],nnvalsum] < 0 || 
            Random[] < Exp[-energydiff[#[[i1,i2]],nnvalsum]],
             lat[[i1,i2]] =-lat[[i1,i2]]; lat, 
             lat]
          )&;
             
  flipLis = NestList[flip, initconfig, m];
  
  monteCarloStepLis =  flipLis[[Range[1,m,n^2]]]; 
  
  longRangeOrderList = 
    Map[Abs[Apply[Plus,Flatten[#]]/n^2]&, monteCarloStepLis];
                                 
  ListPlot[
     longRangeOrderList,
     PlotJoined -> True,
     PlotLabel -> 
     FontForm["Long-range order parameter",{"Times",12}]
          ]; 
          
  (*
  Map[(Show[Graphics[
      RasterArray[#]/.{1->RGBColor[0,1,0],
                       -1->RGBColor[1,0,1]}]])&,
                                        monteCarloStepLis];
  *)                             
  Map[
    (DensityGraphics[ListDensityPlot[#,Mesh->False]])&, 
                                        monteCarloStepLis];                                                                             
       ]

In the isingMetropolis program,we have included as a comment (a comment is text enclosed in (* and *) that is not evalauted when entered), code for two possible graphics displays; one in color and the other in black and white.

Running the Program

Performing the 'flip' operation 20000 times on a '100 by 100' lattice system for B/kT = 0 and J/kT = -0.1, the relaxation of the magnetization is shown below.

SeedRandom[4] 
isingMetropolis[100, 20000, 0, -0.1]

The Determinstic Algorithm

The nature of the spin-flipping process brings up the question of the possibility of a cellular automaton representation of the model. In fact, there is such a model, known as the "Q2R" model. This model is entirely deterministic (except for the starting lattice configuration) and hence runs faster than the probabilistic model.

The Q2R model corresponds to a constant energy condition (in statistical physics, this is known as the microcanonical ensemble formulation). Multiple sites on the lattice are simultaneously selected and each of these sites is flipped only if there is no resulting energy change in the total interaction energy of the site with its neighbors. This condition requires that a site have an equal number of up spin and down spin neighbor sites. 

Because we require that a flip does not alter the energy of the system, we cannot update all of the lattice sites simultaneously. Instead we can simultaneously flip only those spins that do not directly interact. The way we do this is to separate the lattice into two sublattices arranged like a checkerboard. Half of the sites are 'red' and the other half are 'black' such that no sites of the same color are nearest neighbors. On a lattice having an odd number of sites in each row, all of the odd-numbered sites form one sublattice and all of the even-numbered sites belong to the other sublattice.

The update of the entire lattice in one time step is carried out in two stages:
(1) at time t+0.5, all of the red sites are simultaneously updated.
(2) at time t+1, all of the black sites are simultaneously updated, using the updated sites of the first sublattice as the neighbors of the sites in the second sublattice.

The implementation of this algorithm is described as follows:

(1) Start with some initial '(n+1) by n' lattice of spins.

lat = Table[2 Random[Integer]-1,{n+1},{n}]

(2) Divide the lattice into two sublattices.

subLatOdd = Flatten[Partition[Flatten[#],1,2]]

subLatEven = Flatten[Partition[Rest[Flatten[#]],1,2]]

(3) Create the neighborhoods of the lattice sites. This is done using two anonymous functions.

bc = Join[{Last[#]},#,{First[#]}]&
 
nbrhd = Partition[
           bc[Partition[bc[Flatten[#]],{n+2},{n}]],
                                              {3,3},{1,1}]&

The first argument of the outermost Partition operation in the nbrhd function creates an 'n+3 by n+2' matrix from a 'n+1 by n' matrix in two steps: 
(1) to each row, the last element in the preceding row is prepended and the first element in the following row is appended (in the case of the first row, the prepended element is the last element in the last row and in the case of the last row, the appended element is the first element in the first row). 
(2) a copy of the last row is placed before the first row and a copy of the first row is placed after the last row. 

These two steps can be illustrated using a simple lattice example.

Partition[Range[12],{3}]

Applying the first step to this lattice gives

Partition[bc[Flatten[#]],{5},{3}]&[%,{3}]

and applying the second step to this result gives

bc[%]

The outermost Partition operation in nbrhd 'tesselates' the expanded matrix into '3 by 3' matrices with overlapping elements. Applying this to our simple lattice example gives

Partition[%,{3,3},{1,1}]//MatrixForm

Comparing this to the original lattice,

Partition[Range[12],{3}]//MatrixForm

we see that neighbors of sites along the border are determined as follows: 

The left neighbor of the first lattice site is the last lattice site.
The left neighbor of any other site on the left border is taken from the site on the right border in the preceding row.
The right neighbor of the last lattice site is the first site.
The right neighbor of any other site on the right border is taken from the site on the left border in the following row.
The neighbor above any site on the top border is taken from the site in the same column in the last row.
The neighbor below any site on the bottom border is taken from the site in the same column in the first row.    

This way of assigning to border sites, neighbor sites from opposing borders, is known as a spiral or skewed reflecting boundary condition.

(4)  A list of the neighborhoods of the sites in the first sublattice is created.

oddNbrhd = Flatten[nbrhd[#],1][[Range[1,(n (n+1)),2]]]

(5) Using a rule set 

rule[{{_, 1, _}, {1, x_, -1}, {_, -1, _}}] = -x ;
rule[{{_, 1, _}, {-1, x_, 1}, {_, -1, _}}] = -x ;
rule[{{_, 1, _}, {-1, x_, -1}, {_, 1, _}}] = -x ;
rule[{{_, -1, _}, {1, x_, 1}, {_, -1, _}}] = -x ;
rule[{{_, -1, _}, {1, x_, -1}, {_, 1, _}}] = -x ;
rule[{{_, -1, _}, {-1, x_, 1}, {_, 1, _}}] = -x ;
rule[{{_, _, _}, {_, x_, _}, {_, _, _}}] =  x ;

consisting of functions which flip a lattice site if it has an equal number of up-spins and down-spins and leave it unchanged otherwise, the "rule" function is applied to each of the neighborhoods of the sites in the first sublattice to determine the updated values of the sites in the first sublattice.

oddNewVals = Map[rule,oddNbrhd]

(6) The updated values of the sites in the first sublattice are combined with the un-updated values of the sites in the second sublattice.

 newLat =
      Partition[Flatten[
             MapThread[List,{oddNewVals,subLatEven}]],n]
             
(7) The second sublattice is treated in the same way as the first sublattice. 
 
(7a) A list of the neighborhoods of the sites of the second sublattice is created.
 
nbrhdNewLat = nbrhd[newLat]
 
evenNbrhd = 
        Flatten[nbrhdNewLat,1][[Range[2,(n (n+1)),2]]]       

(7b) The rule function is applied to the neighborhoods of each of the sites in the second sublattice (the values of the neighbors to these sites are given by the updated values of the sites in the first sublattice) to determine the updated values of the sites in the second sublattice.

evenNewVals = Map[rule,evenNbrhd]

(7c) The new lattice is created by combining the updated values of the first sublattice and the updated values of the second sublattice. 

Partition[Flatten[
              MapThread[List,{oddNewVals,evenNewVals}]],n] 
              
(8) The updating process given in steps 2-7 is applied many times.

NestList[update,lat,m]

The Program

Putting the pieces of code together, we have the following program for the CA version of the Ising model (the graphics  commands in the isingMetropolis program can be also be used in the isingCA1 program).

isingCA1[n_?OddQ,m_] :=
 Module[{rule,lat,bc,nbrhd,update,subLatOdd,
         subLatEven,oddNbrhd,oddNewVals,newLat,
         evenNbrhd,evenNewVals},
         
 rule[{{_, 1, _}, {1, x_, -1}, {_, -1, _}}] = -x ;
 rule[{{_, 1, _}, {-1, x_, 1}, {_, -1, _}}] = -x ;
 rule[{{_, 1, _}, {-1, x_, -1}, {_, 1, _}}] = -x ;
 rule[{{_, -1, _}, {1, x_, 1}, {_, -1, _}}] = -x ;
 rule[{{_, -1, _}, {1, x_, -1}, {_, 1, _}}] = -x ;
 rule[{{_, -1, _}, {-1, x_, 1}, {_, 1, _}}] = -x ;
 rule[{{_, _, _}, {_, x_, _}, {_, _, _}}] =  x ;
 
 lat = Table[2 Random[Integer]-1,{n+1},{n}];
 
 bc = Join[{Last[#]},#,{First[#]}]&;
 
 nbrhd = Partition[
           bc[Partition[bc[Flatten[#]],{n+2},{n}]],
                                              {3,3},{1,1}]&;
 update = (
    subLatOdd = Flatten[Partition[Flatten[#],1,2]];
    subLatEven = Flatten[Partition[Rest[Flatten[#]],1,2]];
    oddNbrhd = Flatten[nbrhd[#],1][[Range[1,(n (n+1)),2]]];
    oddNewVals = Map[rule,oddNbrhd];
    newLat =
      Partition[Flatten[
             MapThread[List,{oddNewVals,subLatEven}]],n];
    evenNbrhd = 
        Flatten[nbrhd[newLat],1][[Range[2,(n (n+1)),2]]];
    evenNewVals = Map[rule,evenNbrhd];
    Partition[Flatten[
              MapThread[List,{oddNewVals,evenNewVals}]],n]            
    )&;
    
  NestList[update,lat,m]
]

An Alternative Program

We can speed up a CA program by not creating neighborhoods. The program below does this for the deterministic Ising model.

isingCA2[n_?OddQ, m_] :=
  Module[{},
     lat = Table[2 Random[Integer] - 1,{n+1},{n}];
	    subLatOdd = Flatten[Partition[Flatten[lat],1,2]];
	    subLatEven = Flatten[Partition[Rest[Flatten[lat]],1,2]];
	    
	    rule[x_, -1, -1, 1, 1] := -x;
	    rule[x_, -1, 1, -1, 1] := -x;
	    rule[x_, -1, 1, 1, -1] := -x;
	    rule[x_, 1, -1, -1, 1] := -x;
	    rule[x_, 1, -1, 1, -1] := -x;
	    rule[x_, 1, 1, -1, -1] := -x;
	    rule[x_, _, _, _, _] := x;
	    Attributes[rule] = Listable;
	    
	    update[{oldOdd_, oldEven_}] :=
	      Module[{newOdd, newEven},
	      
	       newOdd =
	         rule[oldOdd, RotateLeft[oldEven,(-n-1)/2],
	              RotateLeft[oldEven,-1],RotateLeft[oldEven,0],
	              RotateLeft[oldEven,(n-1)/2]];
	              
	       newEven =
	         rule[oldEven, RotateLeft[newOdd,(-n+1)/2],
	              RotateLeft[newOdd,0],RotateLeft[newOdd,1],
	              RotateLeft[newOdd,(n+1)/2]];
	              
	       {newOdd, newEven}
	            ];
	                    
	 Map[Partition[Flatten[Transpose[#]], n]&,
	            NestList[update, {subLatOdd, subLatEven}, m]] 
]

This program works in essentially the same way as the same isingCA1, except that the nn's to a sublattice site are determined by rotating the entire list of the other sublattice, rather than by creating neighborhoods. 

The rule function takes as its arguments the value of a lattice site and the values of its four nearest-neighbor sites and it changes the sign of the lattice site if it has the same number of up-spin and down-spin nn sites. It is a listable function so that it can be applied directly to the appropriately rotated sublattice lists.

The update function first updates the values of the odd sites and then uses the updated odd site values to update the values of the old even sites.

Note: the NestList operation in isingCA2 takes as arguments and returns as results, ordered pairs whose components are the odd and even sublattices. Hence, it is not necessary to begin the program by creating a lattice and then breaking it apart to form sublattices. This is done here simply to facilitate a comparision of isingCA1 and isingCA2. For the same reason, the last step in the program uses the Map operation to combine the pairs of sublattices into lattices. A timing comparision indicates that isingCA2 is twice as fast as isingCA1.

Comparing the Programs

We should first note that isingCA is not 'exactly' equivalent to isingMetropolis since the cellular automaton is not ergodic. Nonetheless, the two models give many comparable results.   

In order to compare the efficiency of the two programs, only the steps in the isingMetropolis program down to the flipLis operation were retained; the construction of the monteCarloStepLis and longRangeOrderLis lists were deleted and no graphics were created. 

One step in isingCA consists of applying the update function one time. This updates all (n (n+1)) lattice sites and is therefore equal to one monte carlo step. One monte carlo step in isingMetropolis requires n^2 applications of the flip function (in general, m monte Carlo steps consists of m calls of update in isingCA and ((n^2) m) calls of flip in isingMetropolis).

Running a timing comparision by taking one monte carlo step on a '50 by 50' lattice, we find that isingCA2 is about eight times faster than isingMetropolis (were we to create graphics in the two programs, the speed difference would increase since the isingMetropolis program requires us to first create the monteCarloStepLis and longRangeOrderLis lists before creating the graphics).

Conclusion

The reader might wish to carry out some computer 'experiments' based on these programs. One thing that can be done is to replace the spin flip dynamics used in the isingMetropolis program with spin exchange dynamics (this version of the model has been used to model a lattice gas).  This means that rather than selecting a single site and calculating the energy change resulting from performing a flip of the site, a pair of nearest-neighbor sites is selected, the energy change resulting from interchanging their spins is determined and the interchange is determined using the Metropolis algorithm (note that the use of spin exchange dynamics conserves the number of spins in the lattice). Another interesting area to explore is the clustering of sites with the same spin. This might be done using the cluster labelling algorithm [Gaylord, 1992]. Modifying that labelling scheme to determine clusters of alternating spins requires a bit of thought.    

As a final comment, consider the implications of the equivalence between the probabilistic and deterministic algorithms. If any natural (physical, chemical or biological) process can be simulated by a cellular automaton, does that tell us anything about the underlying mechanism by which the universe works? I think that the answer to this question may relevant to the long-standing question as to why mathematics works so unreasonably well in describing natural phenomena.

References

Callen, Earl and Shapero, Don. 1974. "A theory of social imitation", Physics Today, 27(7), pp.23-28.

Gaylord, Richard, 1993. "Percolation Clustering", Mathematica in Education, 2(2), pp.21-24.

What's Coming Next

I'll continue my exploration of cellular automata modelling. A model of neuron excitation and relaxation will be described. This model has relevance to various physiological operations, such as the contraction of cardiac tissue (a subject of special interest to me). I'll also show a related CA which produces some really neat graphics output.

Addendum

This column is the result of a collaboration carried out entirely by e-mail and it may be worthwhile to give its chronology. 

The genesis of the Ising model column was in a Mathematica notebook I came across during my weekly browsing through MathSource, the central archive for Mathematica-based material (including this journal). There was a program by Kazume for the probabilistic Ising model based on a straight-forward implementation of a published FORTRAN program. While I had long been interested in writing a program for the Ising model since it is an exemplar for the use of computers in statistical physics and hence useful for my computer simulation course, I had never gotten around to writing the program because the only Ising model code I could find was written in FORTRAN and whenever I tried to read it, I became nauseous. 

Initially, I sent Kazume some suggestions for improving his procedural program in minor ways (eg., creating the initial lattice configuration more efficiently). After a few exchanges, I suggested to Kazume that we work together to completely re-write the program to be as functional as possible. Kazume agreed to the project, saying that he wanted to become more proficient in his use of the functional style of Mathematica programming. After we wrote the program, Kazume it on a fast machine in Japan to produce the graphical output.

I then asked Kazume if he'd like to work with me on implementing the cellular automaton version of the Ising model. Kazume eagerly signed on to this project, saying that he was not familiar with CA's (a subject that is woefully  neglected in science curricula) and would like to learn about them. The principal problem we had in writing the program was that the published algorithm for the CA Ising model didn't indicate what boundary conditions were used (this is characterstic of both the "real programmers don't document" syndrome and the prevalent attitude amongst physicists that computational details are unimportant (despite the GIGO (garbage in - garbage out) phenomena). In order to figure out what boundary conditions were being used, Kazume executed the FORTRAN program by hand (he has a stronger stomach for FORTRAN code than me).


Working with Kazume has been a very enjoyable experience for both of us and we now talk often on all sorts of subjects. Meanwhile our collaboration continues; we are presently working on a CA model of traffic flow and the laminar-turbulant flow transition which will eventually appear in this column.

One of the reasons that I am giving this history is that I would like to encourage any reader who has a program that they think might be appropriate for this column to get in touch with me (my e-mail address is gaylord@ux1.cso.uiuc.edu). I'll be happy to look at whatever you send me and if it's sufficiently interesting to get me out of the movie theatre (my primary non-Mathematica vocation), I will work with you on optimizing the program and writing it up for this column. 

The best way to make sure that this column (and the MiE journal) is relevant and interesting to you is for you to contribute to it!