A Simple Rule-Based Programming Approach to 
Modeling Complex Random Walk Systems with Cellular Automata

by
Richard J. Gaylord and Kazume Nishidate

Introduction

Many, apparently disparate, processes in nature can be modeled as random walks in which one or more objects move about in space, taking successive steps in randomly chosen directions. While a great deal of analytical (equational) work has been done on the properties of a lone random walker, the multiple random walker system is less well analyzed and exact analytical treatments soon become intractable when interactions between the walkers are introduced. Computer simulation provides an invaluable tool for studying these complex systems. We show here, in detail, how cellular automata and a rule-based programming approach, can be used to model a system of interacting random walkers.

Definition of a cellular automaton

A cellular automaton (CA) is a dynamic system of discrete lattice sites. Starting with various initial values, these sites evolve in discrete time steps as each site assumes a new value based on the values of some neighborhood of sites and a finite number of previous time steps.

The lattice

We will use a two dimensional rectangular lattice for our random walkers CA. The lattice sites on a small, empty square lattice are shown below.

Table[0, {5}, {5}]//MatrixForm

The neighborhood

The 'neighborhood' of a lattice site is comprised of the site and the sites lying in its vicinity. The figure below shows some of the neighboring sites to a site lying in the interior of a two-dimensional rectangular lattice.

The set of five sites consisting of the X site and its four nearest neighbor sites lying north (above), east (right), south (below), and west (left) of X are known as the ``von Neumann'' neighborhood of X.

The set of nine sites consisting of the X site, the four nearest neighbor sites in its von Neumann neighborhood, and the four nearest neighbor sites lying northeast, southeast,  southwest, and northwest of X are known as the ``Moore'' neighborhood of X.

The four outermost sites shown in the figure, lying north of the north nearest neighbor site, east of the east nearest neighbor site, south of the south nearest neighbor site, and west of the west nearest neighbor site of X are next nearest neighbors of X.

The neighbors of a site lying along a side (border) of a rectangular lattice depend on the CA boundary conditions. Under periodic boundary conditions the nearest neighbor left of a site on the left border is the site in the same row on the right border, the nearest neighbor right of a site on the right border is the site in the same row on the left border, the nearest neighbor above a site on the top border is the site in the same column on the bottom border, and the nearest neighbor below a site on the bottom border is the site in the same column on the top border.

The matrices shown below illustrate the Moore neighborhoods of the sites on a small rectangular  lattice, Lat, under periodic boundary conditions.

(Lat = Partition[Range[9], 3]) // MatrixForm

Wrap = Join[{Last[#]}, #, {First[#]}]&;

(LatNeighborhoods = 
Partition[Wrap[Map[Wrap, Lat]], {3, 3}, {1, 1}]) // MatrixForm

We begin our demonstration of the creation of a cellular automaton model of a system of interacting multiple random walkers by looking at the case of a single random walker on a two-dimensional rectangular lattice.

A single (lone) random walker

The system

The  system consists of a (2n + 1) by (2n + 1) square lattice with periodic boundary conditions. The lattice sites all have value 0 except for one site, initially located at the center of the lattice, whose value is an integer between 1 and 4.

initConf =
ReplacePart[Table[0, {2n+1}, {2n+1}],
           Random[Integer, {1,4}], {n+1, n+1}]

The site values represent the following quantities:

0 - an empty site
1 - a site occupied by a north-facing walker
2 - a site occupied by a east-facing walker  
3 - a site occupied by a south-facing walker  
4 - a site occupied by a west-facing walker 

The walker always moves in the direction it is facing in a given time step and hence, a value of 1, 2, 3, 4 indicates that the walker will move in the north, east, south, and west direction, respectively.

Note: Symbols can be used instead of numbers to represent directions, by replacing 1, 2, 3, and 4 with N, E, S and W respectively, and replacing Random[Integer, {1, 4}] with {N, E, S, W}[[Random[Integer, {1, 4}]]] in the program that follow.

The update rules

In a given time step, the walker moves from the site it is occupying to one of the four adjacent sites. Accordingly, the rules for updating the lattice sites in the lone walker CA take five arguments. The arguments are, in order: the value of a site, the value of the nearest neighbor site above the site, the value of the nearest neighbor site to the right of the site, the value of the nearest neighbor site below the site, and the value of the nearest neighbor site to the left of the site. The update rules can be stated both in terms of what happens to the walker and in terms of what happens to the sites of the lattice:

The walker moves from the site it is occupying. 
The value of an occupied site (ie. a site with value 1, 2, 3, or 4) beomes 0.

update[_?Positive, _, _, _, _] := 0

The walker steps to the empty adjacent site it is facing and randomly chooses a direction to face. The value of an empty site (ie. a site with value 0) which is above, right, beneath, or left of a site with value 1, 2, 3, or 4 respectively, is given a randomly selected integer value 1, 2, 3, or 4.

update[_, 3, _, _, _] := RND
update[_, _, 4, _, _] := RND
update[_, _, _, 1, _] := RND
update[_, _, _, _, 2] := RND

where

RND := Random[Integer, {1, 4}]

A empty site that is not faced by a walker on the adjacent site,  remains empty. 
The values of any lattice site, except the occupied site and one adjacent site, remains unchanged.

update[x_, _, _, _, _] := x

Note: It is possible to be more specific in writing these update rules by using 0s rather than underspaces (which will pattern-match any value, including 0). However, it is not necessary, and the use of underspaces will allow us to use these rules in the gladiator walker model that will be discussed later.

Applying the rules

It is not necessary to break apart (decompose) the CA lattice (as we did in the Lat example earlier) in order to form the vonNeumann neighborhoods of the lattice sites which the update rules can then be applied to. Instead, we can use the MapThread function to work with the CA lattice and various rotated versions of the lattice in their entirety.

To illustrate how the MapThread function works, we can use the general, undefined, function g, and the previously defined Lat lattice.

MapThread[g, 
         {#, RotateRight[#, {0,1}], RotateRight[#, {1,0}],
             RotateRight[#, {0,-1}], RotateRight[#, {-1,0}]},
          2]&[Lat] // MatrixForm

Comparing the five arguments of g in each of the nine quantities shown above with the matrices of the neighborhoods of the sites in Lat that were shown earlier, we can see that the arguments of the various g's above are the values of a site in Lat and the four nearest neighbor sites in its von Neumann neighborhood. Thus, given an update rule (or rule set) and a CA lattice, we can apply the rule to each lattice site and its neighbor sites by simply mapping the rule over the second level of the lattice and the appropriate rotated versions of the lattice.

Note: As shown above, for CA's whose update rules depend on the von Neumann neighborhood of a site, four rotated lattices are used. For CA's  whose update rules depend on the Moore neighborhood of a site, eight rotated lattices are used.

To apply the update rules defined above for the lone walker CA, to all of the lattice sites, we apply the following anonymous function to the CA lattice.

updateVonNeumann =
MapThread[update, {#, RotateRight[#, {0,1}],
                      RotateRight[#, {1,0}], 
                      RotateRight[#, {0,-1}],
                      RotateRight[#, {-1,0}]}, 
         2]&

Note: In the lone random walker CA, we could give the update rules the Listable attribute and then update the CA lattice sites by applying the anonymous function update[#, RotateRight[#, {0, 1}], RotateRight[#, {1, 0}], RotateRight[#, {0, -1}],  RotateRight[#, {-1, 0}]]& to the lattice [Gaylord and Nishidate, 1995a]. The reason that we use the MapThread operation instead, is that in more complex CA's (such as we will discuss shortly), the site values are tuplets (lists) rather than simple numbers or symbols and MapThread applies the update rules to the tuplet values of the neighborhood sites.

The evolution of the CA over t time steps is determined using the nest operation

NestList[updateVonNeumann, initConf, t]

The program

Meander[n_, t_] :=
Module[{RND, update, updateVonNeumann, initConf},
      RND := Random[Integer, {1, 4}];
      initConf = ReplacePart[Table[0, {2n+1}, {2n+1}],
                 Random[Integer, {1,4}], {n+1, n+1}];
        
      update[_?Positive, _, _, _, _] := 0;
      update[_, 3, _, _, _] := RND;
      update[_, _, 4, _, _] := RND;
      update[_, _, _, 1, _] := RND;
      update[_, _, _, _, 2] := RND;
      update[x_, _, _, _, _] := x;
        
      updateVonNeumann =
      MapThread[update,
                {#, RotateRight[#, {0,1}],
                    RotateRight[#, {1,0}],
                    RotateRight[#, {0,-1}],
                    RotateRight[#, {-1,0}]}, 
                2]&;
         
      NestList[updateVonNeumann,initConf,t]]

Multiple random walkers

The gladiator walkers CA

The simplest extension of the CA for a single walker to the case of a number of walkers  is straightforward. The Warriors program is identical to the Meander program with one modification. The value of initConfig is given by

initConfig =
Table[Random[Integer, {1, 4}] Floor[1 + p - Random[]],
     {2 n + 1}, {2 n + 1}]

where p is the probability of a site being occupied. This initConfig randomly places a number of walkers on the lattice (in the limit of a large lattice, p approaches the density of walkers in the system).

When the Warriors program is run, the walkers randomly move about the lattice and whenever two or more walkers simultaneously move to the same lattice site, only one walker survives. This evolutionary strategy might represent a society of ``mean-spirited'' individuals practicing a brutish sort of social darwinism. We next turn to the modeling of other, less aggressive, kinds of interactions between the walkers.

The excluded volume constraint

In a system of objects moving about in space, the most basic kind of interaction between the objects is that of volume exclusion, which prohibits having more than one objects in a given place at a given time.

The gladiator walker program satisfies the excluded volume constraint by eliminating all but one walker on any lattice site. This results in a dwindling population of walkers over time.  If we also want to conserve the number of walkers in the system over time, we need to prevent more than one walker from moving to a given lattice site during a time step.

Note: We don't mean to imply that all multiple random walkers CA's must satisfy the constraints of volume exclusion and walker conservation. There are models in which  one or both of these constraints are inappropriate (eg., a chemical reaction model or a model of the migration of ethnic populations to cities). These other multiple random walkers models can also be handled using the rule-based cellular automaton approach presented here and we intend to look at these systems in the future.

We will now show two strategies for conserving the number of walkers while satisfying the volume exclusion constraint.

The shy walkers CA

The system

Initially, the sites of a n by n square lattice with periodic boundary conditions are randomly occupied by random walkers with probability p, and the values of the occupied sites are randomly selected, using

initConf = 
   Table[Floor[p + Random[]], {n}, {n}] * 
   Table[Random[Integer, {1, 4}], {n}, {n}]

The site values represent the following quantities:
 
0 - an empty site
1 - a site occupied by a north-facing walker
2 - a site occupied by a east-facing walker  
3 - a site occupied by a south-facing walker  
4 - a site occupied by a west-facing walker

Note: The density of walkers in initConf (which approaches p as the system size increases) is conserved in the system over time, and is given by

walkerDensity = N[Count[Flatten[initConf], _?Positive]/n^2]

The update rules

The volume exclusion constraint is met by having the walkers avoid any confrontation with one another using the following procedure:

In a given time step, each walker in the system moves to the nearest neighbor site lying in the direction the walker is facing if that site is empty, with one exception: when two or more walkers face the same empty site, they remain where they are (eg., when a right-facing walker lies to the left of an empty site and a south-facing walker lies above the same empty site, both walkers stay put).

The movement of a particle under this avoidance strategy depends on the values of twelve neighboring sites (these sites were shown earlier in the figure of the neighbors of a site X) and therefore, the update rules take thirteen arguments.

update[site, N, E, S, W, NE, SE, SW, NW, Nn, Ee, Ss, Ww]

where the thirteen arguments represent the value of the site, the values of the four nearest neighbors in the N, E, S, W directions, the values of the four nearest neighbors in the NE, SE, SW, NW directions, and the values of four next nearest neighbors in the N, E, S, W directions.

The following twenty-eight (28) update rules are used:

At each time step, a walker, regardless of whether it moves or remains in place, randomly chooses a direction to face, using

RND := Random[Integer, {1, 4}]

A walker facing an empty site moves from the site it is occupying unless another walker faces the same empty site in which case it remains in place and randomly chooses a direction to face.

update[1, 0, _, _, _, 4, _, _, _, _, _, _, _] := RND
update[1, 0, _, _, _, _, _, _, 2, _, _, _, _] := RND
update[1, 0, _, _, _, _, _, _, _, 3, _, _, _] := RND
update[1, 0, _, _, _, _, _, _, _, _, _, _, _] := 0

update[2, _, 0, _, _, 3, _, _, _, _, _, _, _] := RND
update[2, _, 0, _, _, _, 1, _, _, _, _, _, _] := RND
update[2, _, 0, _, _, _, _, _, _, _, 4, _, _] := RND
update[2, _, 0, _, _, _, _, _, _, _, _, _, _] := 0

update[3, _, _, 0, _, _, 4, _, _, _, _, _, _] := RND
update[3, _, _, 0, _, _, _, 2, _, _, _, _, _] := RND
update[3, _, _, 0, _, _, _, _, _, _, _, 1, _] := RND
update[3, _, _, 0, _, _, _, _, _, _, _, _, _] := 0

update[4, _, _, _, 0, _, _, 1, _, _, _, _, _] := RND
update[4, _, _, _, 0, _, _, _, 3, _, _, _, _] := RND
update[4, _, _, _, 0, _, _, _, _, _, _, _, 2] := RND
update[4, _, _, _, 0, _, _, _, _, _, _, _, _] := 0

Any other walker remains in place and randomly selects a direction to face.

update[_?Positive, _, _, _, _, _, _, _, _, _, _, _, _] := RND

An empty site remains empty if two or more walkers face it.

update[0, 3, 4, _, _, _, _, _, _, _, _, _, _] := 0
update[0, 3, _, 1, _, _, _, _, _, _, _, _, _] := 0
update[0, 3, _, _, 2, _, _, _, _, _, _, _, _] := 0
update[0, _, 4, 1, _, _, _, _, _, _, _, _, _] := 0
update[0, _, 4, _, 2, _, _, _, _, _, _, _, _] := 0
update[0, _, _, 1, 2, _, _, _, _, _, _, _, _] := 0

An empty site becomes occupied if it is faced by exactly one walker.

update[0, 3, _, _, _, _, _, _, _, _, _, _, _] := RND
update[0, _, 4, _, _, _, _, _, _, _, _, _, _] := RND
update[0, _, _, 1, _, _, _, _, _, _, _, _, _] := RND
update[0, _, _, _, 2, _, _, _, _, _, _, _, _] := RND

Any other empty site remains empty

update[0, _, _, _, _, _, _, _, _, _, _, _, _] := 0

Applying the rules

The sites in the lattice are updated at each time step by applying the following anonymous function to the CA lattice.

updateMooreNNN =  
     MapThread[update, {#, RotateRight[#, {0, 1}],
                           RotateRight[#, {1, 0}],
                           RotateRight[#, {0, -1}],
                           RotateRight[#, {-1, 0}],
                           RotateRight[#, {1, 1}],
                           RotateRight[#, {1, -1}],
                           RotateRight[#, {-1, -1}],
                           RotateRight[#, {-1, 1}],
                           RotateRight[#, {0, 2}],
                           RotateRight[#, {2, 0}],
                           RotateRight[#, {0, -2}],
                           RotateRight[#, {-2, 0}]}, 
              2]&

The positions of the random walkers evolve over t time steps, starting with the initial lattice configuration, using the following nest operation.

NestList[updateMooreNNN, initConf, t]

The program

Hermits[n_, p_, t_] := Module[{update, initConf, RND},
   
   initConf = 
      Table[Floor[p + Random[]], {n}, {n}] *
      Table[Random[Integer, {1, 4}], {n}, {n}];
   
   RND := Random[Integer, {1, 4}];
   
   update[1, 0, _, _, _, 4, _, _, _, _, _, _, _] := RND;
   update[1, 0, _, _, _, _, _, _, 2, _, _, _, _] := RND;
   update[1, 0, _, _, _, _, _, _, _, 3, _, _, _] := RND;
   update[1, 0, _, _, _, _, _, _, _, _, _, _, _] := 0;
   update[2, _, 0, _, _, 3, _, _, _, _, _, _, _] := RND;
   update[2, _, 0, _, _, _, 1, _, _, _, _, _, _] := RND;
   update[2, _, 0, _, _, _, _, _, _, _, 4, _, _] := RND;
   update[2, _, 0, _, _, _, _, _, _, _, _, _, _] := 0;
   update[3, _, _, 0, _, _, 4, _, _, _, _, _, _] := RND;
   update[3, _, _, 0, _, _, _, 2, _, _, _, _, _] := RND;
   update[3, _, _, 0, _, _, _, _, _, _, _, 1, _] := RND;
   update[3, _, _, 0, _, _, _, _, _, _, _, _, _] := 0;
   update[4, _, _, _, 0, _, _, 1, _, _, _, _, _] := RND;
   update[4, _, _, _, 0, _, _, _, 3, _, _, _, _] := RND;
   update[4, _, _, _, 0, _, _, _, _, _, _, _, 2] := RND;
   update[4, _, _, _, 0, _, _, _, _, _, _, _, _] := 0;
   update[_?Positive, _, _, _, _, _, _, _, _, _, _, _, _] := RND;
   update[0, 3, 4, _, _, _, _, _, _, _, _, _, _] := 0;
   update[0, 3, _, 1, _, _, _, _, _, _, _, _, _] := 0;
   update[0, 3, _, _, 2, _, _, _, _, _, _, _, _] := 0;
   update[0, _, 4, 1, _, _, _, _, _, _, _, _, _] := 0;
   update[0, _, 4, _, 2, _, _, _, _, _, _, _, _] := 0;
   update[0, _, _, 1, 2, _, _, _, _, _, _, _, _] := 0;
   update[0, 3, _, _, _, _, _, _, _, _, _, _, _] := RND;
   update[0, _, 4, _, _, _, _, _, _, _, _, _, _] := RND;
   update[0, _, _, 1, _, _, _, _, _, _, _, _, _] := RND;
   update[0, _, _, _, 2, _, _, _, _, _, _, _, _] := RND;
   update[0, _, _, _, _, _, _, _, _, _, _, _, _] := 0;
   
  updateMooreNNN =  
     MapThread[update, {#,  
                                      RotateRight[#, {0, 1}],
                                      RotateRight[#, {1, 0}],  
                                      RotateRight[#, {0, -1}],  
                                      RotateRight[#, {-1, 0}],
                                      RotateRight[#, {1, 1}],
                                      RotateRight[#, {1, -1}], 
                                      RotateRight[#, {-1, -1}],
                                      RotateRight[#, {-1, 1}],
                                      RotateRight[#, {0, 2}],  
                                      RotateRight[#, {2, 0}], 
                                      RotateRight[#, {0, -2}],    
                                      RotateRight[#, {-2, 0}]}, 
                        2]&;
  
  NestList[updateMooreNNN, initConf, t] 
  ]

The courteous walkers CA

The system

The system consists of a n by n square lattice with periodic boundary conditions whose sites have values which are ordered pairs, {x, y}.  The first component of the ordered pair is an integer ranging from 1 to 4 where the values have the following meanings:

0  - an empty site
1  - a site occupied by a north-facing walker 
2  - a site occupied by a east-facing walker 
3  - a site occupied by a south-facing walker 
4  - a site occupied by a west-facing walker

The second component of the ordered pair is an random real number between 0 and 1 where a zero value is associated with an empty site and a non-zero value is associated with a site occupied by a walker. This random number determines if the walker will move or stay put in a time step when it faces an empty site that is also faced by another walker.

We start by randomly placing walkers on the sites of a n by n lattice, with a probability p.

initConfig = 
    Table[{{Random[Integer, {1, 4}], Random[] }* Floor[Random[] + p], 
    {n}, {n}]

Note: The density of walkers in initConfig (which approaches p as the system size increases) is conserved in the system over time, and is given by

walkerDensity = N[Count[Flatten[initConfig, 1], {_?Positive, _}]/n^2]

The update rules

The volume exclusion constraint was dealt with in the Hermits CA by freezing in place, any particles facing the same vacant lattice site on a given time step. An alternative way to deal with this constraint is to use the following procedure:

When two or more walkers face the same empty site, one of the random walkers is randomly selected to move into the empty site while the other walkers remain in place.

This behavior is implemented using the random numbers  associated with the walkers.

The rules for updating the lattice sites during a time step, take thirteen arguments in the following order:

update[site, N, E, S, W, NE, SE, SW, NW, Nn, Ee, Ss, Ww]

where the thirteen arguments represent the value of the site, the values of the four nearest neighbors in the N, E, S, W directions, the values of the four nearest neighbors in the NE, SE, SW, NW directions, and the values of four next nearest neighbors in the N, E, S, W directions.

At  each time step, both the direction faced by a walker and the random number associated with the walker change randomly, regardless of whether the walker moves or remains in place. This is done, using

rnd:= {Random[Integer, {1, 4}], Random[]}

The behavior of the lattice sites during a time step can be described by the following thirty-eight (38) rules:

A site occupied by a walker who is facing an empty adjacent site that is also faced by zero or more particles (eg., a right-facing walker that lies to the left of an empty site that is faced by an left-facing walker lying to the right of the empty site) remains occupied and randomly selects a new associated random number and direction to face, if the walker's associated random number is less than the random numbers associated with the other walkers. If the walker's associated random number is greater than the random numbers associated with the other walkers, the site becomes empty as the walker moves.

update[{1, a_}, {0, 0}, _, _, _, {4, b_},  _,  _, {2, c_}, {3, d_}, _, _, _] := rnd /; a != Max[a, b, c, d]
update[{1, a_}, {0, 0}, _, _, _, {4, b_},  _,  _, {2, c_}, _, _, _, _] := rnd /; a != Max[a, b, c]                                                                                              
update[{1, a_}, {0, 0}, _, _, _, {4, b_},  _,  _, _, {3, d_}, _, _, _] := rnd /; a != Max[a, b, d]
update[{1, a_}, {0, 0}, _, _, _, _,  _,  _, {2, c_}, {3, d_}, _, _, _] := rnd /; a != Max[a, c, d]
update[{1, a_}, {0, 0}, _, _, _, {4, b_}, _,  _, _, _, _, _, _] := rnd /; a != Max[a, b]  
update[{1, a_}, {0, 0}, _, _, _, _,  _, _, {2, c_}, _, _, _, _] := rnd /; a != Max[a, c]
update[{1, a_}, {0, 0}, _, _, _, _,  _, _, _, {3, d_}, _, _, _] := rnd /; a != Max[a, d]
update[{1, a_}, {0, 0}, _, _, _, _, _, _, _, _, _, _, _] := {0, 0}

update[{2, a_}, _, {0, 0}, _, _, {3, b_}, {1, c_},  _, _, _, {4, d_}, _, _] := rnd /; a != Max[a, b, c, d]
update[{2, a_}, _, {0, 0}, _, _, {3, b_}, {1, c_}, _,  _, _, _, _, _] := rnd /; a != Max[a, b, c]                                                                                              
update[{2, a_}, _, {0, 0}, _, _, {3, b_}, _, _, _, _, {4, d_}, _, _] := rnd /; a != Max[a, b, d]
update[{2, a_}, _, {0, 0}, _, _, _, {1, c_}, _,  _, _, {4, d_}, _, _] := rnd /; a != Max[a, c, d]
update[{2, a_}, _, {0, 0}, _, _, {3, b_}, _,  _,  _, _, _, _, _] := rnd /; a != Max[a, b]  
update[{2, a_}, _, {0, 0}, _, _, _, {1, c_},  _,  _, _, _, _, _] := rnd /; a != Max[a, c]
update[{2, a_}, _, {0, 0}, _, _, _, _,  _,  _, _, {4, d_}, _, _] := rnd /; a != Max[a, d]
update[{2, a_}, _, {0, 0}, _, _, _, _,  _,  _, _, _, _, _] := {0, 0}

update[{3, a_}, _, _, {0, 0}, _, _, {4, b_}, {2, c_}, _, _, _, {1, d_}, _] := rnd /; a != Max[a, b, c, d]
update[{3, a_}, _, _, {0, 0}, _, _, {4, b_}, {2, c_}, _, _, _, _, _] := rnd /; a != Max[a, b, c]
update[{3, a_}, _, _, {0, 0}, _, _, {4, b_}, _, _, _, _, {1, d_}, _] := rnd /; a != Max[a, b, d]
update[{3, a_}, _, _, {0, 0}, _, _,_, {2, c_}, _, _, _, {1, d_}, _] := rnd /; a != Max[a, c, d]
update[{3, a_}, _, _, {0, 0}, _, _, {4, b_}, _, _, _, _, _, _] := rnd /; a != Max[a, b]
update[{3, a_}, _, _, {0, 0}, _, _, _, {2, c_}, _, _, _, _, _] := rnd /; a != Max[a, c]
update[{3, a_}, _, _, {0, 0}, _, _, _, _, _, _, _, {1, d_}, _] := rnd /; a != Max[a, d]
update[{3, a_}, _, _, {0, 0}, _, _, _, _, _, _, _, _, _] := {0, 0}

update[{4, a_}, _, _, _, {0, 0}, _, _, {1, b_}, {3, c_}, _, _, _, {2, d_}] := rnd /; a != Max[a, b, c, d]
update[{4, a_}, _, _, _, {0, 0}, _, _, {1, b_}, {3, c_}, _, _, _, _] := rnd /; a != Max[a, b, c]
update[{4, a_}, _, _, _, {0, 0}, _, _, {1, b_}, _, _, _, _, {2, d_}] := rnd /; a != Max[a, b, d]
update[{4, a_}, _, _, _, {0, 0}, _, _, _, {3, c_}, _, _, _, {2, d_}] := rnd /; a != Max[a, c, d]
update[{4, a_}, _, _, _, {0, 0}, _, _, {1, b_}, _, _, _, _, _] := rnd /; a != Max[a, b]
update[{4, a_}, _, _, _, {0, 0}, _, _, _, {3, c_}, _, _, _, _] := rnd /; a != Max[a, c]
update[{4, a_}, _, _, _, {0, 0}, _, _, _, _, _, _, _, {2, d_}] := rnd /; a != Max[a, d]
update[{4, a_}, _, _, _, {0, 0}, _, _, _, _, _, _, _, _] := {0, 0}

A site occupied by a walker who is facing an occupied site remains occupied by the walker who randomly selects a new associated random number and direction to face.

update[{_?Positive, _}, _, _, _, _, _,_, _, _, _, _, _, _] := rnd

An empty site faced by one or more walkers becomes occupied by a walker who randomly selects a new associated random number and direction to face.

update[{0, 0}, {3, _}, _, _, _, _, _, _, _, _, _, _, _] := rnd
update[{0, 0}, _, {4, _}, _, _, _, _, _, _, _, _, _, _] := rnd
update[{0, 0}, _, _, {1, _}, _, _, _, _, _, _, _, _, _] := rnd
update[{0, 0}, _, _, _, {2, _}, _, _, _, _, _, _, _, _] := rnd

An empty site that is not faced by a walker remains empty.

update[{0, 0}, _, _, _, _, _, _, _, _, _, _, _, _] := {0, 0}

Applying the rules

The sites in the lattice are updated at each time step by applying the anonymous updateMooreNNN function defined earlier, to the lattice.

The positions of the random walkers evolve over t time steps, starting with the initial lattice configuration, using the following nest operation.

NestList[updateMooreNNN, initConfig, t]

The program

Gentlemen[n_, p_, t_]:=
   Module[{initConfig, update, rnd, evolve},

rnd := {Random[Integer, {1, 4}], Random[]};

initConfig = 
    Table[rnd* Floor[Random[] + p], {n}, {n}];      

update[{1, a_}, {0, 0}, _, _, _, {4, b_},  _,  _, {2, c_}, {3, d_}, _, _, _] := rnd /; a != Max[a, b, c, d];
update[{1, a_}, {0, 0}, _, _, _, {4, b_},  _,  _, {2, c_}, _, _, _, _] := rnd /; a != Max[a, b, c];                                                                                             
update[{1, a_}, {0, 0}, _, _, _, {4, b_},  _,  _, _, {3, d_}, _, _, _] := rnd /; a != Max[a, b, d];
update[{1, a_}, {0, 0}, _, _, _, _,  _,  _, {2, c_}, {3, d_}, _, _, _] := rnd /; a != Max[a, c, d];
update[{1, a_}, {0, 0}, _, _, _, {4, b_}, _,  _, _, _, _, _, _] := rnd /; a != Max[a, b];  
update[{1, a_}, {0, 0}, _, _, _, _,  _, _, {2, c_}, _, _, _, _] := rnd /; a != Max[a, c];
update[{1, a_}, {0, 0}, _, _, _, _,  _, _, _, {3, d_}, _, _, _] := rnd /; a != Max[a, d];
update[{1, a_}, {0, 0}, _, _, _, _, _, _, _, _, _, _, _] := {0, 0};

update[{2, a_}, _, {0, 0}, _, _, {3, b_}, {1, c_},  _, _, _, {4, d_}, _, _] := rnd /; a != Max[a, b, c, d];
update[{2, a_}, _, {0, 0}, _, _, {3, b_}, {1, c_}, _, _, _, _, _, _] := rnd /; a != Max[a, b, c];                                                                                              
update[{2, a_}, _, {0, 0}, _, _, {3, b_}, _, _, _, _, {4, d_}, _, _] := rnd /; a != Max[a, b, d];
update[{2, a_}, _, {0, 0}, _, _, _, {1, c_}, _, _, _, {4, d_}, _, _] := rnd /; a != Max[a, c, d];
update[{2, a_}, _, {0, 0}, _, _, {3, b_}, _, _, _, _, _, _, _] := rnd /; a != Max[a, b];  
update[{2, a_}, _, {0, 0}, _, _, _, {1, c_}, _, _, _, _, _, _] := rnd /; a != Max[a, c];
update[{2, a_}, _, {0, 0}, _, _, _, _, _, _, _, {4, d_}, _, _] := rnd /; a != Max[a, d];
update[{2, a_}, _, {0, 0}, _, _, _, _, _, _, _, _, _, _] := {0, 0};

update[{3, a_}, _, _, {0, 0}, _, _, {4, b_}, {2, c_}, _, _, _, {1, d_}, _] := rnd /; a != Max[a, b, c, d];
update[{3, a_}, _, _, {0, 0}, _, _, {4, b_}, {2, c_}, _, _, _, _, _] := rnd /; a != Max[a, b, c];
update[{3, a_}, _, _, {0, 0}, _, _, {4, b_}, _, _, _, _, {1, d_}, _] := rnd /; a != Max[a, b, d];
update[{3, a_}, _, _, {0, 0}, _, _,_, {2, c_}, _, _, _, {1, d_}, _] := rnd /; a != Max[a, c, d];
update[{3, a_}, _, _, {0, 0}, _, _, {4, b_}, _, _, _, _, _, _] := rnd /; a != Max[a, b];
update[{3, a_}, _, _, {0, 0}, _, _, _, {2, c_}, _, _, _, _, _] := rnd /; a != Max[a, c];
update[{3, a_}, _, _, {0, 0}, _, _, _, _, _, _, _, {1, d_}, _] := rnd /; a != Max[a, d];
update[{3, a_}, _, _, {0, 0}, _, _, _, _, _, _, _, _, _] := {0, 0};

update[{4, a_}, _, _, _, {0, 0}, _, _, {1, b_}, {3, c_}, _, _, _, {2, d_}] := rnd /; a != Max[a, b, c, d];
update[{4, a_}, _, _, _, {0, 0}, _, _, {1, b_}, {3, c_}, _, _, _, _] := rnd /; a != Max[a, b, c];
update[{4, a_}, _, _, _, {0, 0}, _, _, {1, b_}, _, _, _, _, {2, d_}] := rnd /; a != Max[a, b, d];
update[{4, a_}, _, _, _, {0, 0}, _, _, _, {3, c_}, _, _, _, {2, d_}] := rnd /; a != Max[a, c, d];
update[{4, a_}, _, _, _, {0, 0}, _, _, {1, b_}, _, _, _, _, _] := rnd /; a != Max[a, b];
update[{4, a_}, _, _, _, {0, 0}, _, _, _, {3, c_}, _, _, _, _] := rnd /; a != Max[a, c];
update[{4, a_}, _, _, _, {0, 0}, _, _, _, _, _, _, _, {2, d_}] := rnd /; a != Max[a, d];
update[{4, a_}, _, _, _, {0, 0}, _, _, _, _, _, _, _, _] := {0, 0};

update[{_?Positive, _}, _, _, _, _, _,_, _, _, _, _, _, _] := rnd;

update[{0, 0}, {3, _}, _, _, _, _, _, _, _, _, _, _, _] := rnd;
update[{0, 0}, _, {4, _}, _, _, _, _, _, _, _, _, _, _] := rnd;
update[{0, 0}, _, _, {1, _}, _, _, _, _, _, _, _, _, _] := rnd;
update[{0, 0}, _, _, _, {2, _}, _, _, _, _, _, _, _, _] := rnd;

update[{0, 0}, _, _, _, _, _, _, _, _, _, _, _, _] := {0, 0};

NestList[MapThread[update, {#, RotateRight[#, {1, 0}], RotateRight[#, {0, -1}],
                                                 RotateRight[#, {-1, 0}], RotateRight[#, {0, 1}],
                                                 RotateRight[#, {1, -1}], RotateRight[#, {-1, -1}], 
                                                 RotateRight[#, {-1, 1}], RotateRight[#, {1, 1}], 
                                                 RotateRight[#, {2, 0}], RotateRight[#, {0, -2}], 
                                                 RotateRight[#, {-2, 0}], RotateRight[#, {0, 2}] }, 
                                    2]&, 
                  initConfig, t]
]

Applications of multiple random walker systems

We can illustrate the use of our approach to the CA modeling of a system of mutiple random walkers using examples taken from the physical, materials, and biological sciences.  We will look at the expansion of the universe [Gaylord and Nishidate, 1995a], the formation of banding in a solid electrolyte [Gaylord and Nishidate, 1995b], and the activity of an ant colony [Gaylordand Nishidate, 1995c]. The focus here will be on the development of the programs (the graphical output andthe physical interpretation of these models is given in the referenced articles).

Note: While we could use variants of either the Hermits CA or the Gentlemen CA in these examples (see project #3 in the computer simulation projects section), we will only work with the Hermits CA so that we can use ordered pairs, rather than triplets, for the lattice site values and fewer rules, which will make the development of the code a bit easier to follow.

An example from physical science

The expansion of the universe

The universe is generally thought to have originated in a ``big bang'' explosion. We can model the expansion of a system of initially clustered particles using a CA model of random walkers who move randomly in any direction subject to the excluded volume constraint (see project #4 in the computer simulation projects section).

The system

The only modifications that need to be made to the Hermits program are to: (1) decorate the initial lattice configuration with a border of empty sites, and (2) apply the following anonymous conditional function to the lattice that results from the application of the update function to the lattice in each time step.

increaseLat = 
  If[Max[Map[{Max[First[#]], Max[Last[#]]}&, {#, Transpose[#]}&[#]]] != 0, zeros[#], #]&

where

zeros =
 (Prepend[Append[Map[Prepend[Append[#, 0], 0]&, #], 
                              Table[0, {Length[#] + 2}]], Table[0, {Length[#] + 2}]])&

The zeros function adds a border of zeros to the lattice and the increaseLat function adds this borders of 0s whenever a particle moves to a site along the border. The reason for maintaining a border of zeros around the lattice is that in its absense, the periodic boundary conditions might result in particles disappearing off one side of the lattice and reappearing on the opposite side of the lattice (which would be a poor representation of an expanding universe). The reason for adding a new border of zeros only when a particle appears along the border is to minimize the number of lattice sites that have to be updated in a given time step.

The program

Boom[n_, p_, t_] :=
 Module[{update, initConf, RND},

RND := Random[Integer, {1, 4}];

update[1, 0, _, _, _, 4, _, _, _, _, _, _, _] := RND;
update[1, 0, _, _, _, _, _, _, 2, _, _, _, _] := RND;
update[1, 0, _, _, _, _, _, _, _, 3, _, _, _] := RND;
update[1, 0, _, _, _, _, _, _, _, _, _, _, _] := 0;
update[2, _, 0, _, _, 3, _, _, _, _, _, _, _] := RND;
update[2, _, 0, _, _, _, 1, _, _, _, _, _, _] := RND;
update[2, _, 0, _, _, _, _, _, _, _, 4, _, _] := RND;
update[2, _, 0, _, _, _, _, _, _, _, _, _, _] := 0;
update[3, _, _, 0, _, _, 4, _, _, _, _, _, _] := RND;
update[3, _, _, 0, _, _, _, 2, _, _, _, _, _] := RND;
update[3, _, _, 0, _, _, _, _, _, _, _, 1, _] := RND;
update[3, _, _, 0, _, _, _, _, _, _, _, _, _] := 0;
update[4, _, _, _, 0, _, _, 1, _, _, _, _, _] := RND;
update[4, _, _, _, 0, _, _, _, 3, _, _, _, _] := RND;
update[4, _, _, _, 0, _, _, _, _, _, _, _, 2] := RND;
update[4, _, _, _, 0, _, _, _, _, _, _, _, _] := 0;
update[_?Positive,___] := RND;
update[0, 3, 4, _, _, _, _, _, _, _, _, _, _] := 0;
update[0, 3, _, 1, _, _, _, _, _, _, _, _, _] := 0;
update[0, 3, _, _, 2, _, _, _, _, _, _, _, _] := 0;
update[0, _, 4, 1, _, _, _, _, _, _, _, _, _] := 0;
update[0, _, 4, _, 2, _, _, _, _, _, _, _, _] := 0;
update[0, _, _, 1, 2, _, _, _, _, _, _, _, _] := 0;
update[0, 3, _, _, _, _, _, _, _, _, _, _, _] := RND;
update[0, _, 4, _, _, _, _, _, _, _, _, _, _] := RND;
update[0, _, _, 1, _, _, _, _, _, _, _, _, _] := RND;
update[0, _, _, _, 2, _, _, _, _, _, _, _, _] := RND;
update[0, _, _, _, _, _, _, _, _, _, _, _, _] := 0;
Attributes[update]=Listable;

zeros = (Prepend[Append[Map[Prepend[Append[#, 0], 0]&, #],
              Table[0, {Length[#] + 2}]], Table[0, {Length[#] + 2}]])&;

increaseLat = 
  If[Max[Map[{Max[First[#]], Max[Last[#]]}&, {#, Transpose[#]}&[#]]] != 0, zeros[#], #]&;


initConf = zeros[Table[Floor[p + Random[]], {n}, {n}] *
                           Table[Random[Integer, {1, 4}], {n}, {n}]];
                  
NestList[increaseLat[update[#,
                                         RotateRight[#,{0,1}], RotateRight[#,{1,0}],
                                         RotateRight[#,{0,-1}], RotateRight[#,{-1,0}],
                                         RotateRight[#,{1,1}], RotateRight[#,{1,-1}],
                                         RotateRight[#,{-1,-1}], RotateRight[#,{-1,1}],
                                         RotateRight[#,{0,2}], RotateRight[#,{2,0}],
                                         RotateRight[#,{0,-2}], RotateRight[#,{-2,0}]]]&, 
                initConf, t]
]

An example from materials science

Materials known as superionic conductors (aka fast ion conductors, super ionic solids, and solid electrolytes) exhibit high electrical conductivity and have potential for use in batteries, fuel cells, etc. One simulation model of these materials is based on the biased diffusion of a system of two species of particles (eg.,  oppositely charged species of particles) in the presense of an external field which drives each type of particle in opposite directions.

The driven diffusion of two species

The system

There are two types of particles, A and B, on a n by n square lattice with periodic boundary conditions. 

Initially, A and B type particles, facing randomly chosen directions, are randomly placed on the sites of the lattice, with probability p, using

initConf = 
   Table[{{1, 2, 3, 4}[[Random[Integer, {1, 4}] ]], {A, B}[[Random[Integer, {1, 2}]]] }
             * Floor[Random[] + p], 
             {n}, {n}];

The values of the sites in initConf are ordered pairs. The first component of the ordered pair is an integer ranging from 0 to 4 where the values have the following meanings:

0  - an empty site
1  - a site occupied by a north-facing particle 
2  - a site occupied by a east-facing particle
3  - a site occupied by a south-facing particle
4  - a site occupied by a west-facing particle

The second component of the ordered pair is either  0, A or B where 0 indicates an empty site, A indicates a site occupied by an A type particle, and B indicates a site occupied by an B type particle. Hence, each lattice site has a tuplet of the form {1 | 2 | 3 | 4, A | B} (indicating a site occupied by a particle) or {0, 0} (indicating an empty site) associated with it.

The update rules

The particles in the system have hardcore volume exclusion, so the rules for the behavior of a particle during a time step, take thirteen arguments

update[site, N, E, S, W, NE, SE, SW, NW, Nn, Ee, Ss, Ww]

where the arguments represent the value of the site the values of the nearest neighbors in the N, E, S, W directions, the values of the nearest neighbors in the NE, SE, SW, NW directions, and thevalues of the next nearest neighbors in the N, E, S, W directions.

In a given time step, regardless of whether or not it moves, a particle randomly picks a direction to face and this direction will determine the direction in which the walker will try to move on the next time step. The presense of an external field will cause the particles of one species to prefer to move in a given direction and causes the particles of the other species to prefer to move in the opposite direction.

Specifically, an A particle will have a preference to move to the north and the east while a B particle will prefer to move to the south and the west. This preference can be implemented in the CA by biasing the liklihood that a particle will choose a particular direction to face on a given time step. We will use rndA to choose the direction faced by A particles and rndB to choose the direction faced by B particles. These quantities are defined as follows:

rndA := 
{Random[Integer, {1, 2}], 
 Random[Integer, {3, 4}]}[[2 - Floor[Random[] + a]]];

rndA chooses Random[Integer, {1, 2}] a% of the time and Random[Integer, {3, 4}] (1-a)% of the time so that  1, 2, 3, and 4 are chosen (a/2)%, (a/2)%, (1 - a)/2%, and 
(1 - a)/2% of the time, respectively.

rndB := {Random[Integer, {1, 2}], Random[Integer, {3, 4}]}[[1 + Floor[Random[] + a]]];

rndB chooses Random[Integer, {1, 2}] (1 - a)% of the time and Random[Integer, {3, 4}] a% of the time so that  1, 2, 3, and 4 are chosen (1 - a)/2%, (1 - a)/2%, (a/2)%, and (a/2)% of the time, respectively

Using these quantities, the behavior of the lattice sites during a time step can be described by the following fifty-one rules:

A particle facing a site occupied by another particle (of either type) remains in place and randomly chooses a direction to face.

update[{1, A}, {x_?Positive, _}, _, _, _, _, _, _, _, _, _, _, _] := {rndA, A}
update[{2, A}, _, {x_?Positive, _}, _, _, _, _, _, _, _, _, _, _] := {rndA, A}
update[{3, A}, _, _, {x_?Positive, _}, _, _, _, _, _, _, _, _, _] := {rndA, A}
update[{4, A}, _, _, _, {x_?Positive, _}, _, _, _, _, _, _, _, _] := {rndA, A}
update[{1, B}, {x_?Positive, _}, _, _, _, _, _, _, _, _, _, _, _] := {rndB, B}
update[{2, B}, _, {x_?Positive, _}, _, _, _, _, _, _, _, _, _, _] := {rndB, B}
update[{3, B}, _, _, {x_?Positive, _}, _, _, _, _, _, _, _, _, _] := {rndB, B}
update[{4, B}, _, _, _, {x_?Positive, _}, _, _, _, _, _, _, _, _] := {rndB, B}

A particle facing an empty adjacent site that is also faced by another particle of either type (e.g., an right-facing particle that lies to the left of an empty site that is faced by an left-facing particle lying to the right of the empty site) stays put and randomly chooses a direction to face.

update[{1, A}, {0, 0}, _, _, _, {4, _}, _, _, _, _, _, _, _] := {rndA, A}
update[{1, A}, {0, 0}, _, _, _, _, _, _, {2, _}, _, _, _, _] := {rndA, A}
update[{1, A}, {0, 0}, _, _,_ ,_ ,_ ,_ , _, {3, _}, _, _, _] := {rndA, A}
update[{2, A}, _, {0, 0}, _, _, {3, _}, _, _, _, _, _, _, _] := {rndA, A}
update[{2, A}, _, {0, 0}, _, _, _, {1, _}, _, _, _, _, _, _] := {rndA, A}
update[{2, A}, _, {0, 0}, _, _, _, _, _, _, _, {4, _}, _, _] := {rndA, A}
update[{3, A}, _, _, {0, 0}, _, _, {4, _}, _, _, _, _, _, _] := {rndA, A}
update[{3, A}, _, _, {0, 0}, _, _, _, {2, _}, _, _, _, _, _] := {rndA, A}
update[{3, A}, _, _, {0, 0}, _, _, _, _, _, _, _, {1, _}, _] := {rndA, A}
update[{4, A}, _, _, _, {0, 0}, _, _, {1, _}, _, _, _, _, _] := {rndA, A}
update[{4, A}, _, _, _, {0, 0}, _, _, _, {3, _}, _, _, _, _] := {rndA, A}
update[{4, A}, _, _, _, {0, 0}, _, _, _, _, _, _, _, {2, _}] := {rndA, A}
update[{1, B}, {0, 0}, _, _, _, {4, _}, _, _, _, _, _, _, _] := {rndB, B}
update[{1, B}, {0, 0}, _, _, _, _, _, _, {2, _}, _, _, _, _] := {rndB, B}
update[{1, B}, {0, 0}, _, _,_ ,_ ,_ ,_ , _, {3, _}, _, _, _] := {rndB, B}
update[{2, B}, _, {0, 0}, _, _, {3, _}, _, _, _, _, _, _, _] := {rndB, B}
update[{2, B}, _, {0, 0}, _, _, _, {1, _}, _, _, _, _, _, _] := {rndB, B}
update[{2, B}, _, {0, 0}, _, _, _, _, _, _, _, {4, _}, _, _] := {rndB, B}
update[{3, B}, _, _, {0, 0}, _, _, {4, _}, _, _, _, _, _, _] := {rndB, B}
update[{3, B}, _, _, {0, 0}, _, _, _, {2, _}, _, _, _, _, _] := {rndB, B}
update[{3, B}, _, _, {0, 0}, _, _, _, _, _, _, _, {1, _}, _] := {rndB, B}
update[{4, B}, _, _, _, {0, 0}, _, _, {1, _}, _, _, _, _, _] := {rndB, B}
update[{4, B}, _, _, _, {0, 0}, _, _, _, {3, _}, _, _, _, _] := {rndB, B}
update[{4, B}, _, _, _, {0, 0}, _, _, _, _, _, _, _, {2, _}] := {rndB, B}

A site occupied by a particle facing an empty adjacent site that is not faced by another particle is vacated.

update[{1, _}, {0, 0}, _, _, _, _, _, _, _, _, _, _, _] := {0, 0}
update[{2, _}, _, {0, 0}, _, _, _, _, _, _, _, _, _, _] := {0, 0}
update[{3, _}, _, _, {0, 0}, _, _, _, _, _, _, _, _, _] := {0, 0}
update[{4, _}, _, _, _, {0, 0}, _, _, _, _, _, _, _, _] := {0, 0}

An empty site faced by two or more particles remains empty.

update[{0, 0}, {3, _}, {4, _}, _, _, _, _, _, _, _, _, _, _] := {0, 0}
update[{0, 0}, {3, _}, _, {1, _}, _, _, _, _, _, _, _, _, _] := {0, 0}
update[{0, 0}, {3, _}, _, _, {2, _}, _, _, _, _, _, _, _, _] := {0, 0}
update[{0, 0}, _, {4, _}, {1, _}, _, _, _, _, _, _, _, _, _] := {0, 0}
update[{0, 0}, _, {4, _}, _, {2, _}, _, _, _, _, _, _, _, _] := {0, 0}
update[{0, 0}, _, _, {1, _}, {2, _}, _, _, _, _, _, _, _, _] := {0, 0}

An empty site faced by exactly one particle beomes occupied.

update[{0, 0}, {3, A}, _, _, _, _, _, _, _, _, _, _, _] := {rndA, A}
update[{0, 0}, _, {4, A}, _, _, _, _, _, _, _, _, _, _] := {rndA, A}
update[{0, 0}, _, _, {1, A}, _, _, _, _, _, _, _, _, _] := {rndA, A}
update[{0, 0}, _, _, _, {2, A}, _, _, _, _, _, _, _, _] := {rndA, A}
update[{0, 0}, {3, B}, _, _, _, _, _, _, _, _, _, _, _] := {rndB, B}
update[{0, 0}, _, {4, B}, _, _, _, _, _, _, _, _, _, _] := {rndB, B}
update[{0, 0}, _, _, {1, B}, _, _, _, _, _, _, _, _, _] := {rndB, B}
update[{0, 0}, _, _, _, {2, B}, _, _, _, _, _, _, _, _] := {rndB, B}

An empty site remains unchanged.

update[{0, 0}, _, _, _, _, _, _, _, _, _, _, _, _] := {0, 0}

Note: We used fifty-one update rules above in order to make the rules as clear as possible. If greater brevity is desired, the number of update rules can be reduced to thirty-one (with a resulting small loss in runtime efficiency) by making the following three changes:

(a) using the values 1 and 2, rather than A and B, to indicate the type of particle 

(b) replacing rndA and rndB with

rnd := {
        {Random[Integer, {1, 2}], Random[Integer, {3, 4}]}[[2 - Floor[Random[] + a]]],
        {Random[Integer, {1, 2}], Random[Integer, {3, 4}]}[[1 + Floor[Random[] + a]]]
       }

(c) combining the rules that are the same, except for the type of particle they apply to. For example, the two rules

update[{1, A}, {x_?Positive, _}, _, _, _, _, _, _, _, _, _, _, _] := {rndA, A}
update[{1, B}, {x_?Positive, _}, _, _, _, _, _, _, _, _, _, _, _] := {rndB, B}

can be merged into the single rule

update[{1, y_}, {x_?Positive, _}, _, _, _, _, _, _, _, _, _, _, _] := {rnd[[y]], y}

Applying the rules

The rules are applied to the lattice using the updateMooreNNN function defined earlier.

NestList[updateMooreNNN, initConf, t]

The program

DrivenDiffusionTwoSpecies[n_, p_, a_, t_]:=
   Module[{initPhase, update, rndA, rndB, evolve},
   
initPhase = 
   Table[{{1, 2, 3, 4}[[Random[Integer, {1, 4}] ]], {A, B}[[Random[Integer, {1, 2}]]] }
             * Floor[Random[] + p], 
             {n}, {n}]; 

rndA := {Random[Integer, {1, 2}], Random[Integer, {3, 4}]}[[2 - Floor[Random[] + a]]];
      
rndB := {Random[Integer, {1, 2}], Random[Integer, {3, 4}]}[[1 + Floor[Random[] + a]]];
      
update[{1, A}, {x_?Positive, _}, _, _, _, _, _, _, _, _, _, _, _] := {rndA, A};
update[{2, A}, _, {x_?Positive, _}, _, _, _, _, _, _, _, _, _, _] := {rndA, A};
update[{3, A}, _, _, {x_?Positive, _}, _, _, _, _, _, _, _, _, _] := {rndA, A};
update[{4, A}, _, _, _, {x_?Positive, _}, _, _, _, _, _, _, _, _] := {rndA, A};
update[{1, B}, {x_?Positive, _}, _, _, _, _, _, _, _, _, _, _, _] := {rndB, B};
update[{2, B}, _, {x_?Positive, _}, _, _, _, _, _, _, _, _, _, _] := {rndB, B};
update[{3, B}, _, _, {x_?Positive, _}, _, _, _, _, _, _, _, _, _] := {rndB, B};
update[{4, B}, _, _, _, {x_?Positive, _}, _, _, _, _, _, _, _, _] := {rndB, B};
update[{1, A}, {0, 0}, _, _, _, {4, _}, _, _, _, _, _, _, _] := {rndA, A};
update[{1, A}, {0, 0}, _, _, _, _, _, _, {2, _}, _, _, _, _] := {rndA, A};
update[{1, A}, {0, 0}, _, _,_ ,_ ,_ ,_ , _, {3, _}, _, _, _] := {rndA, A};
update[{2, A}, _, {0, 0}, _, _, {3, _}, _, _, _, _, _, _, _] := {rndA, A};
update[{2, A}, _, {0, 0}, _, _, _, {1, _}, _, _, _, _, _, _] := {rndA, A};
update[{2, A}, _, {0, 0}, _, _, _, _, _, _, _, {4, _}, _, _] := {rndA, A};
update[{3, A}, _, _, {0, 0}, _, _, {4, _}, _, _, _, _, _, _] := {rndA, A};
update[{3, A}, _, _, {0, 0}, _, _, _, {2, _}, _, _, _, _, _] := {rndA, A};
update[{3, A}, _, _, {0, 0}, _, _, _, _, _, _, _, {1, _}, _] := {rndA, A};
update[{4, A}, _, _, _, {0, 0}, _, _, {1, _}, _, _, _, _, _] := {rndA, A};
update[{4, A}, _, _, _, {0, 0}, _, _, _, {3, _}, _, _, _, _] := {rndA, A};
update[{4, A}, _, _, _, {0, 0}, _, _, _, _, _, _, _, {2, _}] := {rndA, A};
update[{1, B}, {0, 0}, _, _, _, {4, _}, _, _, _, _, _, _, _] := {rndB, B};
update[{1, B}, {0, 0}, _, _, _, _, _, _, {2, _}, _, _, _, _] := {rndB, B};
update[{1, B}, {0, 0}, _, _,_ ,_ ,_ ,_ , _, {3, _}, _, _, _] := {rndB, B};
update[{2, B}, _, {0, 0}, _, _, {3, _}, _, _, _, _, _, _, _] := {rndB, B};
update[{2, B}, _, {0, 0}, _, _, _, {1, _}, _, _, _, _, _, _] := {rndB, B};
update[{2, B}, _, {0, 0}, _, _, _, _, _, _, _, {4, _}, _, _] := {rndB, B};
update[{3, B}, _, _, {0, 0}, _, _, {4, _}, _, _, _, _, _, _] := {rndB, B};
update[{3, B}, _, _, {0, 0}, _, _, _, {2, _}, _, _, _, _, _] := {rndB, B};
update[{3, B}, _, _, {0, 0}, _, _, _, _, _, _, _, {1, _}, _] := {rndB, B};
update[{4, B}, _, _, _, {0, 0}, _, _, {1, _}, _, _, _, _, _] := {rndB, B};
update[{4, B}, _, _, _, {0, 0}, _, _, _, {3, _}, _, _, _, _] := {rndB, B};
update[{4, B}, _, _, _, {0, 0}, _, _, _, _, _, _, _, {2, _}] := {rndB, B};
update[{1, _}, {0, 0}, _, _, _, _, _, _, _, _, _, _, _] := {0, 0};
update[{2, _}, _, {0, 0}, _, _, _, _, _, _, _, _, _, _] := {0, 0};
update[{3, _}, _, _, {0, 0}, _, _, _, _, _, _, _, _, _] := {0, 0};
update[{4, _}, _, _, _, {0, 0}, _, _, _, _, _, _, _, _] := {0, 0};
update[{0, 0}, {3, _}, {4, _}, _, _, _, _, _, _, _, _, _, _] := {0, 0};
update[{0, 0}, {3, _}, _, {1, _}, _, _, _, _, _, _, _, _, _] := {0, 0};
update[{0, 0}, {3, _}, _, _, {2, _}, _, _, _, _, _, _, _, _] := {0, 0};
update[{0, 0}, _, {4, _}, {1, _}, _, _, _, _, _, _, _, _, _] := {0, 0};
update[{0, 0}, _, {4, _}, _, {2, _}, _, _, _, _, _, _, _, _] := {0, 0};
update[{0, 0}, _, _, {1, _}, {2, _}, _, _, _, _, _, _, _, _] := {0, 0};
update[{0, 0}, {3, A}, _, _, _, _, _, _, _, _, _, _, _] := {rndA, A};
update[{0, 0}, _, {4, A}, _, _, _, _, _, _, _, _, _, _] := {rndA, A};
update[{0, 0}, _, _, {1, A}, _, _, _, _, _, _, _, _, _] := {rndA, A};
update[{0, 0}, _, _, _, {2, A}, _, _, _, _, _, _, _, _] := {rndA, A};
update[{0, 0}, {3, B}, _, _, _, _, _, _, _, _, _, _, _] := {rndB, B};
update[{0, 0}, _, {4, B}, _, _, _, _, _, _, _, _, _, _] := {rndB, B};
update[{0, 0}, _, _, {1, B}, _, _, _, _, _, _, _, _, _] := {rndB, B};
update[{0, 0}, _, _, _, {2, B}, _, _, _, _, _, _, _, _] := {rndB, B};
update[{0, 0}, _, _, _, _, _, _, _, _, _, _, _, _] := {0, 0};

NestList[MapThread[update, 
                                  {#, RotateRight[#, {1, 0}], RotateRight[#, {0, -1}],
                                    RotateRight[#, {-1, 0}], RotateRight[#, {0, 1}],
                                    RotateRight[#, {1, -1}], RotateRight[#, {-1, -1}], 
                                    RotateRight[#, {-1, 1}], RotateRight[#, {1, 1}], 
                                    RotateRight[#, {2, 0}], RotateRight[#, {0, -2}], 
                                    RotateRight[#, {-2, 0}], RotateRight[#, {0, 2}] }, 
                                   2]&, 
                initPhase, t]
]

An example from biological science

The activity of an ant colony

A society of ants shows rather fascinating complex behavior, even though an individual ant has no apparent native intelligence (an ant cannot be trained to perform even the simplest of tasks) and a group of ants has no apparent social hierarchical order or chain of command (no ``boss'' ant directs the action of other ants). It has been suggested by biologists that the behaviors of an ant colony might arise from local interactions of the sort found in a CA.

The system

An ant colony is comprised of ants who, at any given time, are either active or inactive. Inactive ants are immobile while active ants move about within the territory of the colony. It appears that active ants become inactive after a period of time, while inactive ants become active either spontaneously or as a result of coming into contact with an active ant.

The ant colony CA consists of a n by n square lattice whose sites have values which are ordered pairs, {x, y}.  The first component of the ordered pair is either the symbol, b, or an integer ranging from -1 to 4 where the values have the following meanings:

b  - a border site
0  - an empty site
-1 - a site occupied by an inactive ant
1  - a site occupied by a north-facing active ant 
2  - a site occupied by a east-facing active ant 
3  - a site occupied by a south-facing active ant 
4  - a site occupied by a west-facing active ant

The second component of the ordered pair, which we will call the 'activity time clock', is an integer ranging from 0 to s where the values have the following meanings:

0 - an empty site, a border site or a site occupied by an inactive ant
1 - a site occupied by an ant that has just become active on that time step
2 - a site occupied by an ant that became active on the previous time step
...
s - a site occupied by an ant that became active on on the sth previous time step

We start by randomly placing inactive and active ants on the sites of a (n - 2) by (n - 2) lattice, with a probability p. If the ant is active, it faces a randomly chosen direction and its activity time clock is randomly set between 1 and s.

antPopulation = 
    Table[{{-1, 1, 2, 3, 4}[[Random[Integer, {1, 5}] ]], 
       Random[Integer, {1, s}] }* Floor[Random[] + p], {n - 2}, {n - 2}] /. {-1, _} -> {-1, 0}

The ant density on the ant farm is conserved in the system over time. It can be calculated using

antDensity = N[Apply[Plus, Abs[Sign[Transpose[Flatten[antPopulation, 1]][[1]]]]]/(n - 2)^2]

When the system is large,  the ant density approaches the value of p.

The antPopulation lattice is decorated with a border of sites having value {b, 0} using

antFarm = border[antPopulation]

where

border =
 (Prepend[Append[Map[Prepend[Append[#, {b, 0}], {b, 0}]&, #], 
                              Table[{b, 0}, {Length[#] + 2}]], Table[{b, 0}, {Length[#] + 2}]])&

The update rules

The rules for updating the values of the CA lattice sites during a time step, take thirteen arguments

update[site, N, E, S, W, NE, SE, SW, NW, Nn, Ee, Ss, Ww]

where the arguments represent the value of the site and the values of twelve neighbor sites; the nearest neighbors in the N, E, S, W directions, the nearest neighbors in the NE, SE, SW, NW directions, and the next nearest neighbors in the N, E, S, W directions.

The following thirty-eight (38) update rules describe the behavior of the lattice sites during a time step:

The function RND is a  random number generator used to determine the direction faced by an ant.

RND := Random[Integer, {1, 4}]

An inactive ant (a site whose value is {-1, 0}) adjacent to an active ant becomes active, randomly chooses a direction to face, and resets its activity time clock to 1.

update[{-1,  0}, {x_?Positive, _}, _, _, _, _, _, _, _, _, _, _, _] := {RND, 1}
update[{-1,  0},  _, {x_?Positive,_}, _, _, _, _, _, _, _, _, _, _] := {RND, 1}
update[{-1,  0},  _, _, {x_?Positive,_}, _, _, _, _, _, _, _, _, _] := {RND, 1}
update[{-1,  0},  _, _, _, {x_?Positive,_}, _, _, _, _, _, _, _, _] := {RND, 1}

An inactive ant becomes active, randomly chooses a direction to face, and resets its its activity time clock to 1, with probability r.

update[{-1, 0}, _, _, _, _, _, _, _, _, _, _, _, _] := 
   { -1 + # * Random[Integer, {2, 5}], #}&[Floor[Random[] + r] ]

An active ant (a site whose value is an ordered pair whose first component is 1, 2, 3, or 4) become inactive after s time steps and resets its activity time clock to 0.

update[{x_?Positive, s}, _, _, _, _, _, _, _, _, _, _, _, _] :=
{-1, 0}

An active ant facing an adjacent site that is either on the border or is occupied by another ant, stays put, randomly chooses a direction to face, and increments its activity time clock by 1.

update[{1, y_}, {-1 | x_?Positive | b, _}, _, _, _, _, _, _, _, _, _, _, _] := {RND, y + 1}
update[{2, y_}, _, {-1 | x_?Positive | b, _}, _, _, _, _, _, _, _, _, _, _] := {RND, y + 1}
update[{3, y_}, _, _, {-1 | x_?Positive | b, _}, _, _, _, _, _, _, _, _, _] := {RND, y + 1}
update[{4, y_}, _, _, _, {-1 | x_?Positive | b, _}, _, _, _, _, _, _, _, _] := {RND, y + 1}

An active ant facing an empty adjacent site that is also faced by another active ant stays put, randomly chooses a direction to face, and increments its activity time clock by 1.

update[{1, y_}, {0, 0}, _, _, _, {4, _}, _, _, _, _, _, _, _] := {RND, y + 1}
update[{1, y_}, {0, 0}, _, _, _, _, _, _, {2, _}, _, _, _, _] := {RND, y + 1}
update[{1, y_}, {0, 0}, _, _,_ ,_ ,_ ,_ , _, {3, _}, _, _, _] := {RND, y + 1}
update[{2, y_}, _, {0, 0}, _, _, {3, _}, _, _, _, _, _, _, _] := {RND, y + 1}
update[{2, y_}, _, {0, 0}, _, _, _, {1, _}, _, _, _, _, _, _] := {RND, y + 1}
update[{2, y_}, _, {0, 0}, _, _, _, _, _, _, _, {4, _}, _, _] := {RND, y + 1}
update[{3, y_}, _, _, {0, 0}, _, _, {4, _}, _, _, _, _, _, _] := {RND, y + 1}
update[{3, y_}, _, _, {0, 0}, _, _, _, {2, _}, _, _, _, _, _] := {RND, y + 1}
update[{3, y_}, _, _, {0, 0}, _, _, _, _, _, _, _, {1, _}, _] := {RND, y + 1}
update[{4, y_}, _, _, _, {0, 0}, _, _, {1, _}, _, _, _, _, _] := {RND, y + 1}
update[{4, y_}, _, _, _, {0, 0}, _, _, _, {3, _}, _, _, _, _] := {RND, y + 1}
update[{4, y_}, _, _, _, {0, 0}, _, _, _, _, _, _, _, {2, _}] := {RND, y + 1}

An active ant facing an empty adjacent site that is not faced by another active ant vacates the site it is occupying.

update[{1, _}, {0, 0}, _, _, _, _, _, _, _, _, _, _, _] := {0, 0}
update[{2, _}, _, {0, 0}, _, _, _, _, _, _, _, _, _, _] := {0, 0}
update[{3, _}, _, _, {0, 0}, _, _, _, _, _, _, _, _, _] := {0, 0}
update[{4, _}, _, _, _, {0, 0}, _, _, _, _, _, _, _, _] := {0, 0}

An empty site (a site having an ordered pair {0, 0}) faced by two or more active ants stays empty.

update[{0, 0}, {3, _}, {4, _}, _, _, _, _, _, _, _, _, _, _] := {0, 0}
update[{0, 0}, {3, _}, _, {1, _}, _, _, _, _, _, _, _, _, _] := {0, 0}
update[{0, 0}, {3, _}, _, _, {2, _}, _, _, _, _, _, _, _, _] := {0, 0}
update[{0, 0}, _, {4, _}, {1, _}, _, _, _, _, _, _, _, _, _] := {0, 0}
update[{0, 0}, _, {4, _}, _, {2, _}, _, _, _, _, _, _, _, _] := {0, 0}
update[{0, 0}, _, _, {1, _}, {2, _}, _, _, _, _, _, _, _, _] := {0, 0}

An empty site faced by one ant that has been active less than s time steps beomes occupied by the ant which randomly chooses a direction to face and increments its activity time clock by 1.

update[{0, 0}, {3, y_?(# < s &)}, _, _, _, _, _, _, _, _, _, _, _] := {RND, y + 1} 
update[{0, 0}, _, {4, y_?(# < s &)}, _, _, _, _, _, _, _, _, _, _] := {RND, y + 1}
update[{0, 0}, _, _, {1, y_?(# < s &)}, _, _, _, _, _, _, _, _, _] := {RND, y + 1}
update[{0, 0}, _, _, _, {2, y_?(# < s &)}, _, _, _, _, _, _, _, _] := {RND, y + 1}

An empty site remains unchanged.

update[{0, 0}, _, _, _, _, _, _, _, _, _, _, _, _] := {0, 0}

A border site (a site having an ordered pair {b, 0}) remains unchanged.

update[{b, 0}, _, _, _, _, _, _, _, _, _, _, _, _] := {b, 0}

Applying the rules

The rules are applied to the lattice using the updateMooreNNN function defined previously.

NestList[updateMooreNNN, antPopulation, t]

The program

MaggiesFarm[n_, p_, s_, r_, t_]:=
   Module[{antPopulation, border, antFarm, RND, update},

antPopulation = 
    Table[{{-1, 1, 2, 3, 4}[[Random[Integer, {1, 5}] ]], 
     Random[Integer, {1, s}] }* Floor[Random[] + p], {n - 2}, {n - 2}] /.  {-1, _} -> {-1, 0}; 
                       
border =
 (Prepend[Append[Map[Prepend[Append[#, {b, 0}], {b, 0}]&, #], 
                              Table[{b, 0}, {Length[#] + 2}]], Table[{b, 0}, {Length[#] + 2}]])&;                                                                                      
                                                                                               
antFarm = border[antPopulation];

RND := Random[Integer, {1, 4}];

update[{-1,  0}, {x_?Positive, _}, _, _, _, _, _, _, _, _, _, _, _] := {RND, 1};
update[{-1,  0},  _, {x_?Positive,_}, _, _, _, _, _, _, _, _, _, _] := {RND, 1};
update[{-1,  0},  _, _, {x_?Positive,_}, _, _, _, _, _, _, _, _, _] := {RND, 1};
update[{-1,  0},  _, _, _, {x_?Positive,_}, _, _, _, _, _, _, _, _] := {RND, 1};
update[{-1, 0}, _, _, _, _, _, _, _, _, _, _, _, _] := 
                            { -1 + # * Random[Integer, {2, 5}], #}&[Floor[Random[] + r] ];
update[{x_?Positive, s}, _, _, _, _, _, _, _, _, _, _, _, _] := {-1, 0};
update[{1, y_}, {-1 | x_?Positive | b, _}, _, _, _, _, _, _, _, _, _, _, _] := {RND, y + 1};
update[{2, y_}, _, {-1 | x_?Positive | b, _}, _, _, _, _, _, _, _, _, _, _] := {RND, y + 1};
update[{3, y_}, _, _, {-1 | x_?Positive | b, _}, _, _, _, _, _, _, _, _, _] := {RND, y + 1};
update[{4, y_}, _, _, _, {-1 | x_?Positive | b, _}, _, _, _, _, _, _, _, _] := {RND, y + 1};
update[{1, y_}, {0, 0}, _, _, _, {4, _}, _, _, _, _, _, _, _] := {RND, y + 1};
update[{1, y_}, {0, 0}, _, _, _, _, _, _, {2, _}, _, _, _, _] := {RND, y + 1};
update[{1, y_}, {0, 0}, _, _,_ ,_ ,_ ,_ , _, {3, _}, _, _, _] := {RND, y + 1};
update[{2, y_}, _, {0, 0}, _, _, {3, _}, _, _, _, _, _, _, _] := {RND, y + 1};
update[{2, y_}, _, {0, 0}, _, _, _, {1, _}, _, _, _, _, _, _] := {RND, y + 1};
update[{2, y_}, _, {0, 0}, _, _, _, _, _, _, _, {4, _}, _, _] := {RND, y + 1};
update[{3, y_}, _, _, {0, 0}, _, _, {4, _}, _, _, _, _, _, _] := {RND, y + 1};
update[{3, y_}, _, _, {0, 0}, _, _, _, {2, _}, _, _, _, _, _] := {RND, y + 1};
update[{3, y_}, _, _, {0, 0}, _, _, _, _, _, _, _, {1, _}, _] := {RND, y + 1};
update[{4, y_}, _, _, _, {0, 0}, _, _, {1, _}, _, _, _, _, _] := {RND, y + 1};
update[{4, y_}, _, _, _, {0, 0}, _, _, _, {3, _}, _, _, _, _] := {RND, y + 1};
update[{4, y_}, _, _, _, {0, 0}, _, _, _, _, _, _, _, {2, _}] := {RND, y + 1};
update[{1, _}, {0, 0}, _, _, _, _, _, _, _, _, _, _, _] := {0, 0};
update[{2, _}, _, {0, 0}, _, _, _, _, _, _, _, _, _, _] := {0, 0};
update[{3, _}, _, _, {0, 0}, _, _, _, _, _, _, _, _, _] := {0, 0};
update[{4, _}, _, _, _, {0, 0}, _, _, _, _, _, _, _, _] := {0, 0};
update[{0, 0}, {3, _}, {4, _}, _, _, _, _, _, _, _, _, _, _] := {0, 0};
update[{0, 0}, {3, _}, _, {1, _}, _, _, _, _, _, _, _, _, _] := {0, 0};
update[{0, 0}, {3, _}, _, _, {2, _}, _, _, _, _, _, _, _, _] := {0, 0};
update[{0, 0}, _, {4, _}, {1, _}, _, _, _, _, _, _, _, _, _] := {0, 0};
update[{0, 0}, _, {4, _}, _, {2, _}, _, _, _, _, _, _, _, _] := {0, 0};
update[{0, 0}, _, _, {1, _}, {2, _}, _, _, _, _, _, _, _, _] := {0, 0};
update[{0, 0}, {3, y_?(# < s &)}, _, _, _, _, _, _, _, _, _, _, _] := {RND, y + 1}; 
update[{0, 0}, _, {4, y_?(# < s &)}, _, _, _, _, _, _, _, _, _, _] := {RND, y + 1};
update[{0, 0}, _, _, {1, y_?(# < s &)}, _, _, _, _, _, _, _, _, _] := {RND, y + 1};
update[{0, 0}, _, _, _, {2, y_?(# < s &)}, _, _, _, _, _, _, _, _] := {RND, y + 1};
update[{0, 0}, _, _, _, _, _, _, _, _, _, _, _, _] := {0, 0};
update[{b, 0}, _, _, _, _, _, _, _, _, _, _, _, _] := {b, 0};

NestList[MapThread[update, 
       {#, RotateRight[#, {1, 0}], RotateRight[#, {0, -1}], RotateRight[#, {-1, 0}],
        RotateRight[#, {0, 1}], RotateRight[#, {1, -1}], RotateRight[#, {-1, -1}], 
        RotateRight[#, {-1, 1}], RotateRight[#, {1, 1}], RotateRight[#, {2, 0}], 
        RotateRight[#, {0, -2}], RotateRight[#, {-2, 0}], RotateRight[#, {0, 2}] }, 2]&, 
       antFarm, t]
]

Conclusion

The rule-based programming approach to the cellular automaton modeling of complex random walker systems is both powerful and easy to use. Essentially what one has to do is:

(1) determine the rules for the behavior of the walkers at each time step.

(2) figure out how to convert the rules for the walkers into rules for updating the lattice sites.

(3) implement the update rules as Mathematica code.

The first step requires an in-depth knowledge of the particular system being modeled.

Note: A word of caution to random walk specialists, who have a tendency, in common with algorithmic and theoretical physicists in general, to fall prey to The Law of the Hammer, which states that a child who is given a hammer finds that everything it encounters requires pounding:

It is important to establish a more than superficial resemblance of a natural phenomenon to a random walk process before modeling it as such (an animal should waddle and quack, as well as look, like a duck before it is cooked and served a la orange).

The second step is the hardest part of the process of creating a CA model because it requires one to think most deeply. One of the great advantages of using the rule-based programming approach is it can be carried out away from the computer. All that is required is paper and pencil (and a coffeehouse).

The third and final step is where the actual programming of the CA rules occurs. The very sophisticated pattern-matching facilities of the Mathematica programming language (and the ability to create anonymous predicate functions to make the pattern-matching conditional) allows you to express these rules in a very direct manner and the interpreted nature of the language allows you to debug the rules one-by-one, making Mathematica ideal for the task.

Computer Simulation Projects

1. Run timing comparisions of the Hermits and Gentlemen CAs.

2. An interesting system of interacting random walkers, known as the  ``vicious'' walkers model, consists of walkers who ``shoot each other dead on sight'' so that when more than one walker reaches the same site at the same time, they are all eliminated (this behavior is reminiscent of the final scene in the movie ``Reservoir Dogs''). Modify the Hermits program to model a system of N vicious walkers.

3. Implement the driven diffusion of two species model using the courteous walker CA approach.

Note: you will need to use triplets, {direction faced, associated random number, species identification}, for the lattice site values.

4. Modify the Boom program to incorporate the tendency of a particle to move in an outward direction, such as occurs in an explosion.

Note: you'll have to identify which quadrant a particle is in, and then use an appropriate random integer generator function (e.g., a particle in the NE quadrant might use rndA from the DrivenDiffusionTwoSpecies program).

General References

Gaylord, Richard J. and Wellin, Paul R. Computer Simulations with Mathematica: Explorations in Complex Physical and Biological Systems. TELOS/Springer-Verlag (1995).

Gaylord, Richard J. and Nishidate, Kazume. The lattice gas model of interacting random walkers. Mathematica World, February (1995a).

Gaylord, Richard J. and Nishidate, Kazume. Ant colony activity.  Mathematica World, March (1995b).

Gaylord, Richard J. and Nishidate, Kazume. The driven diffusion of a two-species system, Mathematica World, April (1995c).