Physical Simulation with Cellular Automata | Details

Matlab code considerations

The following considerations will be illustrated using a Matlab program which computes Conway's life. Parts of the code are explained below.

Matrixes and images can be converted to one-another, so display is straighforward. If the array cells contains the binary state of each cell and the array z contains zeros, then the cat command builds an RGB image, which is displayed by the image command. The image command also returns a handle.
```
imh = image(cat(3,cells,z,z));
set(imh, 'erasemode', 'none')
axis equal
axis tight
  
```
Matrixes and images can be converted to one-another, so initial conditions may be computed by matrix or graphics commands. The following code generates an array of zeros, initializes the cell state to zero, then makes a cross of cells in state=1 in the center of the array. One of the examples below (percolation cluster) uses graphics commands to initialize the CA array.
```
z = zeros(n,n);
cells = z;
cells(n/2,.25*n:.75*n) = 1;
cells(.25*n:.75*n,n/2) = 1;
	
```

The Matlab code should be vectorized to minimize overhead. The following two statements compute the sum of nearest neighbors and compute the CA rule. The code makes use of Matlab's very flexible indexing to specify the nearest neighbors.

x = 2:n-1;
y = 2:n-1;
sum(x,y) = cells(x,y-1) + cells(x,y+1) + ...
            cells(x-1, y) + cells(x+1,y) + ...
            cells(x-1,y-1) + cells(x-1,y+1) + ...
            cells(x+1,y-1) + cells(x+1,y+1);
cells = (sum==3) | (sum==2 & cells);

Adding a simple GUI is easy. In this example three buttons and a text field were implmented. The three buttons allow a user to Run, Stop, or Quit. The text field displays the number of simulation steps executed.

%build the GUI
%define the plot button
plotbutton=uicontrol('style','pushbutton',...
   'string','Run', ...
   'fontsize',12, ...
   'position',[100,400,50,20], ...
   'callback', 'run=1;');

%define the stop button
erasebutton=uicontrol('style','pushbutton',...
   'string','Stop', ...
   'fontsize',12, ...
   'position',[200,400,50,20], ...
   'callback','freeze=1;');

%define the Quit button
quitbutton=uicontrol('style','pushbutton',...
   'string','Quit', ...
   'fontsize',12, ...
   'position',[300,400,50,20], ...
   'callback','stop=1;close;');

number = uicontrol('style','text', ...
    'string','1', ...
   'fontsize',12, ...
   'position',[20,400,50,20]);

After the controls (and CA) are initialized, the program drops into a loop which tests the state of the variables which are set in the callback functions of each button. For the moment, just look at the nested structure of the while loop and if statements. The program loops until the Quit button is pushed. The other two buttons cause an if statement to execute when pushed.

stop= 0; %wait for a quit button push
run = 0; %wait for a draw 
freeze = 0; %wait for a freeze
while (stop==0) 
       if (run==1)
	      %nearest neighbor sum
	      sum(x,y) = cells(x,y-1) + cells(x,y+1) + ...
		            cells(x-1, y) + cells(x+1,y) + ...
		            cells(x-1,y-1) + cells(x-1,y+1) + ...
		           cells(3:n,y-1) + cells(x+1,y+1);
	      % The CA rule
	      cells = (sum==3) | (sum==2 & cells);       
	      %draw the new image
	      set(imh, 'cdata', cat(3,cells,z,z) )
	      %update the step number diaplay
	      stepnumber = 1 + str2num(get(number,'string'));
	      set(number,'string',num2str(stepnumber))
      end
      if (freeze==1)
	     run = 0;
	     freeze = 0;
      end
      drawnow  %need this in the loop for controls to work  
end

Examples

Conway's life.
The rule is:

Sum the 8 nearest neighbors
If the sum=2 then the state does not change
If the sum=3 then the state=1
Otherwise the state=0

The code:

x = 2:n-1;
y = 2:n-1;
  %nearest neighbor sum
sum(x,y) = cells(x,y-1) + cells(x,y+1) + ...
		   cells(x-1, y) + cells(x+1,y) + ...
		   cells(x-1,y-1) + cells(x-1,y+1) + ...
		   cells(3:n,y-1) + cells(x+1,y+1);
% The CA rule
cells = (sum==3) | (sum==2 & cells);

Surface Tension
The rule is:

Sum the 8 nearest neighbors and the cell itself
If the sum< 4 or sum=5 then the state=0
Otherwise the state=1

The code:

x = 2:n-1;
y = 2:n-1;
  sum(x,y) = cells(x,y-1) + cells(x,y+1) + ...
            cells(x-1, y) + cells(x+1,y) + ...
            cells(x-1,y-1) + cells(x-1,y+1) + ...
            cells(3:n,y-1) + cells(x+1,y+1)+...
            cells(x,y);
        % The CA rule
        cells = ~((sum< 4) | (sum==5));

Percolation Cluster
The rule:
- Sum the 8 nearest neighbors of each cell (cells are binary valued, 0/1). Cells also have a separate state variable (called 'visited') which records if they have ever had a nonzero neighbor before.
- Compute a random number r between 0 and 1.
- If the sum> 0 (at least one neighbor) and the r>threshold, and the cell has never had a neighbor then cell=1.
- If the sum> 0 set the "visited" flag for the cell to record that the cell has a nonzero neighbor.
The update code:
```
	sum(2:a-1,2:b-1) = cells(2:a-1,1:b-2) + cells(2:a-1,3:b) + ...
                        cells(1:a-2, 2:b-1) + cells(3:a,2:b-1) + ...
                        cells(1:a-2,1:b-2) + cells(1:a-2,3:b) + ...
                        cells(3:a,1:b-2) + cells(3:a,3:b);
    
    pick =  rand(a,b); 
    cells = cells  | ((sum>=1) & (pick>=threshold) & (visit==0))  ;
    visit = (sum>=1) ;
	
```
The variables a and b are the size of the image. The inital image is determined by graphics operations. The following statements set up axes of a fixed size, write text into the axes, then grab the contents of the axes and put them back into an array using the getframe function..
```
ax = axes('units','pixels','position',[1 1 500 400],'color','k');
text('units', 'pixels', 'position', [50,255,0],...
    'string','BioNB','color','w','fontname','helvetica','fontsize',100)
text('units', 'pixels', 'position', [120,120,0],...
    'string','441','color','w','fontname','helvetica','fontsize',100)
initial = getframe(gca);

[a,b,c]=size(initial.cdata);
z=zeros(a,b);
cells = double(initial.cdata(:,:,1)==255);
visit = z ;
sum = z;
	
```
After a few dozen time steps (starting with the BioNB 441 image) we get the following image. Click on it to see a full size image.

Excitable media (BZ reaction or heart)
The rule:

Cells can be in 10 different states. State 0 is resting. States 1-5 are active, and states 6-9 are refractory.
Count the 8 nearest neighbors of each cell which are in one of the active states.
If the sum is greater or equal to 3 (at least three active neighbors) then cell=1.
States 1 to 9 occur stepwise with no more input. If state=1 then the next state=2. If state=2 then the next state=3, and similarly of all the states up to 9. If state=9 then the next state=0 and the cell is back at rest.

The update code:

x = [2:n-1];
y = [2:n-1];

    sum(x,y) = ((cells(x,y-1)> 0)&(cells(x,y-1)< t)) + ((cells(x,y+1)> 0)&(cells(x,y+1)< t)) + ...
        ((cells(x-1, y)> 0)&(cells(x-1, y)< t)) + ((cells(x+1,y)> 0)&(cells(x+1,y)< t)) + ...
        ((cells(x-1,y-1)> 0)&(cells(x-1,y-1)< t)) + ((cells(x-1,y+1)> 0)&(cells(x-1,y+1)< t)) + ...
        ((cells(x+1,y-1)> 0)&(cells(x+1,y-1)< t)) + ((cells(x+1,y+1)> 0)&(cells(x+1,y+1)< t));
       
    cells = ((cells==0) & (sum>=t1)) + ...
            2*(cells==1) + ...
            3*(cells==2) + ...
            4*(cells==3) + ...
            5*(cells==4) + ...
            6*(cells==5) +...
            7*(cells==6) +...
            8*(cells==7) +...
            9*(cells==8) +...
            0*(cells==9);

An image from this CA after spiral waves develop is shown below.

Forest Fire
The rule:
- Cells can be in 3 different states. State=0 is empty, state=1 is burning and state=2 is forest.
- If one or more of the 4 neighbors if a cell is burning and it is forest (state=2) then the new state is burning (state=1).
- There is a low probablity (say 0.000005) of a forest cell (state=2) starting to burn on its own (from lightning).
- A cell which is burning (state=1) becomes empty (state=0).
- There is a low probability (say, 0.01) of an empty cell becoming forest to simulate growth.
- The array is considered to be toroidly connected, so that fire which burns to left side will start fires on the right. The top and bottom are similarly connected.
The update code:
```
 
	sum = (veg(1:n,[n 1:n-1])==1) + (veg(1:n,[2:n 1])==1) + ...
           (veg([n 1:n-1], 1:n)==1) + (veg([2:n 1],1:n)==1) ;
 
    veg = ...
         2*(veg==2) - ((veg==2) & (sum> 0 | (rand(n,n)< Plightning))) + ...
         2*((veg==0) & rand(n,n)< Pgrowth) ;
	
```
Note that the toroidal connection is implemented by the ordering of subscripts.

Gas dynamics
This CA (and the next two) are used to compute motion of particles. This application requires a different kind of neighborhood. The neighborhood of a cell varys on each time step so that motion in a certain direction will continue in the same direction. In other words, the rule conserves momentum, which is the basis of dynamical calculations. The neighborhood used is called a Margolis neighborhood and consists of overlapping 2x2 blocks of cells. In the following diagram, the upper-left 4 cells are the neighborhood during an even time-step, the bottom-right 4 cells are the neighborhood during an odd time-step. A given cell has 3 neighbors at each time-step, but the specific cells which constitute the neighbors flips back and forth.

even	even
even	cell	odd
odd	odd

The rule:

This is called a HPP-GAS rule. See reference [1] Chapter 12..
Cells have 2 states. State=0 is empty, state=1 represents a particle in motion.
To make a particle collision (which conserves momentum and energy), treat the case of exactly two particles on a diagonal as if they hit and deflected each other 90 degrees. This is done by converting one diagonal to the other on the time-step. You can implement this by rotating the 4 cells counterclockwise by one cell.
To make a particle collide with a wall, simply leave its state unchanged. This causes a reflection.

The update code:

	p=mod(i,2); %margolis neighborhood, where i is the time step
    %upper left cell update
    xind = [1+p:2:nx-2+p];
    yind = [1+p:2:ny-2+p];
    
    %See if exactly one diagonal is ones
    %only (at most) one of the following can be true!
    diag1(xind,yind) = (sand(xind,yind)==1) & (sand(xind+1,yind+1)==1) & ...
        (sand(xind+1,yind)==0) & (sand(xind,yind+1)==0);
    
    diag2(xind,yind) = (sand(xind+1,yind)==1) & (sand(xind,yind+1)==1) & ...
        (sand(xind,yind)==0) & (sand(xind+1,yind+1)==0);
    
    %The diagonals both not occupied by two particles
    and12(xind,yind) = (diag1(xind,yind)==0) & (diag2(xind,yind)==0);
    
    %One diagonal is occupied by two particles
    or12(xind,yind)  = diag1(xind,yind) | diag2(xind,yind);
    
    %for every gas particle see if it near the boundary
    sums(xind,yind) = gnd(xind,yind) | gnd(xind+1,yind) | ...
                        gnd(xind,yind+1) | gnd(xind+1,yind+1) ;
    
    % cell layout:
    % x,y    x+1,y
    % x,y+1  x+1,y+1
    %If (no walls) and (diagonals are both not occupied)  
    %then there is no collision, so move opposite cell to current cell
    %If (no walls) and (only one diagonal is occupied) 
    %then there is a collision so move ccw cell to the current cell
    %If (a wall) 
    %then don't change the cell (causes a reflection)
    sandNew(xind,yind) = ...
        (and12(xind,yind)  & ~sums(xind,yind) & sand(xind+1,yind+1)) + ... 
        (or12(xind,yind) & ~sums(xind,yind) & sand(xind,yind+1)) + ...
        (sums(xind,yind) & sand(xind,yind)); 
        
    sandNew(xind+1,yind) = ...
        (and12(xind,yind)  & ~sums(xind,yind) & sand(xind,yind+1)) + ... 
        (or12(xind,yind) & ~sums(xind,yind) & sand(xind,yind))+ ...
        (sums(xind,yind) & sand(xind+1,yind));  
        
    sandNew(xind,yind+1) = ...    
        (and12(xind,yind)  & ~sums(xind,yind) & sand(xind+1,yind)) + ... 
        (or12(xind,yind) & ~sums(xind,yind) & sand(xind+1,yind+1))+ ...
        (sums(xind,yind) & sand(xind,yind+1)); 
        
     sandNew(xind+1,yind+1) = ...    
        (and12(xind,yind)  & ~sums(xind,yind) & sand(xind,yind)) + ... 
        (or12(xind,yind) & ~sums(xind,yind) & sand(xind+1,yind))+ ...
        (sums(xind,yind) & sand(xind+1,yind+1)); 
    
    sand = sandNew;

Diffusion limited aggregation
This system simulates sticky particles aggregating to form a fractal structure. particle motion takes place with a rule similar to the HPP-GAS rule in example 6. The difference is that particles are assumed to be bouncing around in some dense (but invisible) liquid. The effect is to randomize the direction of motion of every particle at every time step. Put differently, every time-step is a collision step. The simulation is also seeded with one fixed particle in the center of the array. Any diffusing particle which touches it sticks to it, and itself becomes a non-moving, sticky particle.

The rule:

Use a Margolus neighborhood. At every time step, rotate the 4 cells either clockwise or counterclockwise by one cell with equal probability. The rotation randomizes the velocities.
After the move, if one or more of the eight nearest neighboors is a fixed, sticky particle, then freeze the particle and make it sticky.

The update code:

	p=mod(i,2); %margolis neighborhood
   
    %upper left cell update
    xind = [1+p:2:nx-2+p];
    yind = [1+p:2:ny-2+p];
    %random velocity choice
    vary = rand(nx,ny)< .5 ;
    vary1 = 1-vary;
    
    %diffusion rule -- margolus neighborhood
    %rotate the 4 cells to randomize velocity
    sandNew(xind,yind) = ...
        vary(xind,yind).*sand(xind+1,yind) + ... %cw
        vary1(xind,yind).*sand(xind,yind+1) ;    %ccw
        
    sandNew(xind+1,yind) = ...
        vary(xind,yind).*sand(xind+1,yind+1) + ...
        vary1(xind,yind).*sand(xind,yind) ;
        
    sandNew(xind,yind+1) = ...    
        vary(xind,yind).*sand(xind,yind) + ...
        vary1(xind,yind).*sand(xind+1,yind+1) ;
        
     sandNew(xind+1,yind+1) = ...    
        vary(xind,yind).*sand(xind,yind+1) + ...
        vary1(xind,yind).*sand(xind+1,yind) ;
    
    sand = sandNew;
    
    %for every sand grain see if it near the fixed, sticky cluster
    sum(2:nx-1,2:ny-1) = gnd(2:nx-1,1:ny-2) + gnd(2:nx-1,3:ny) + ...
                        gnd(1:nx-2, 2:ny-1) + gnd(3:nx,2:ny-1) + ...
                        gnd(1:nx-2,1:ny-2) + gnd(1:nx-2,3:ny) + ...
                        gnd(3:nx,1:ny-2) + gnd(3:nx,3:ny);
    
    %add to the cluster
    gnd = ((sum> 0) & (sand==1)) | gnd ;
    %and eliminate the moving particle
    sand(find(gnd==1)) = 0;

The following image shows the fixed cluster after many time steps.

Sand pile
The cross-section of a pile of sand can be modeled using a Margolus neighborhood to propagate cells, but with a different motion rule

The rule:

See reference [2] Chapter 2.2.6..
Cells have 2 states. State=0 is empty, state=1 represents agrain of sand.
On any step, a particle can fall toward the bottom of the the 2x2 block. The possible transitions are shown below. Walls and floors stop motion.

		see Full doc link for the matrix state transitions

The update code:

	p=mod(i,2); %margolis neighborhood
    sand(nx/2,ny/2) = 1; %add a grain at the top
   
    %upper left cell update
    xind = [1+p:2:nx-2+p];
    yind = [1+p:2:ny-2+p];
    %randomize the flow -- 10% of the time 
    vary = rand(nx,ny)< .9 ;
    vary1 = 1-vary;
    
    sandNew(xind,yind) = ...
        gnd(xind,yind).*sand(xind,yind) + ...
        (1-gnd(xind,yind)).*sand(xind,yind).*sand(xind,yind+1) .* ...
            (sand(xind+1,yind+1)+(1-sand(xind+1,yind+1)).*sand(xind+1,yind));
    
    sandNew(xind+1,yind) = ...
        gnd(xind+1,yind).*sand(xind+1,yind) + ...
        (1-gnd(xind+1,yind)).*sand(xind+1,yind).*sand(xind+1,yind+1) .* ...
            (sand(xind,yind+1)+(1-sand(xind,yind+1)).*sand(xind,yind));
        
    sandNew(xind,yind+1) = ...    
        sand(xind,yind+1) + ...
        (1-sand(xind,yind+1)) .* ...
        (  sand(xind,yind).*(1-gnd(xind,yind)) + ...
           (1-sand(xind,yind)).*sand(xind+1,yind).*(1-gnd(xind+1,yind)).*sand(xind+1,yind+1));
        
     sandNew(xind+1,yind+1) = ...    
        sand(xind+1,yind+1) + ...
        (1-sand(xind+1,yind+1)) .* ...
        (  sand(xind+1,yind).*(1-gnd(xind+1,yind)) + ...
           (1-sand(xind+1,yind)).*sand(xind,yind).*(1-gnd(xind,yind)).*sand(xind,yind+1));
    
    %scramble the sites to make it look better
    temp1 = sandNew(xind,yind+1).*vary(xind,yind+1) + ...
        sandNew(xind+1,yind+1).*vary1(xind,yind+1);
    
    temp2 = sandNew(xind+1,yind+1).*vary(xind,yind+1) + ...
        sandNew(xind,yind+1).*vary1(xind,yind+1);
    sandNew(xind,yind+1) = temp1;
    sandNew(xind+1,yind+1) = temp2;
    
    sand = sandNew;

References

[1] Cellular Automata Machines by Tommaso Toffoli and Norman Margolus, MIT Press, 1987.

[2] Cellular Automata Modeling of Physical Systems by Bastien Chopard and Michel Droz, Cambridge University Press, 1998.

Project Details