ca.htm

交通流研究中的元胞自动机模型源码

        <TD width="18%">1</TD>
        <TD width="18%">0</TD>
        <TD width="32%">to</TD>
        <TD width="15%">1</TD>
        <TD width="17%">1</TD></TR></TBODY></TABLE>
    <P>
    <TABLE width="23%" border=2>
      <TBODY>
      <TR>
        <TD width="18%">1</TD>
        <TD width="18%">1</TD>
        <TD width="32%">goes</TD>
        <TD width="15%">1</TD>
        <TD width="17%">1</TD></TR>
      <TR>
        <TD width="18%">1</TD>
        <TD width="18%">1</TD>
        <TD width="32%">to</TD>
        <TD width="15%">1</TD>
        <TD width="17%">1</TD></TR></TBODY></TABLE>
    <P></P>
    <LI>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. The third rule above becomes: 
    <P>
    <TABLE width="23%" border=2>
      <TBODY>
      <TR>
        <TD width="18%">1</TD>
        <TD width="18%">0</TD>
        <TD width="32%">goes</TD>
        <TD width="15%">0</TD>
        <TD width="17%">1</TD></TR>
      <TR>
        <TD width="18%">0</TD>
        <TD width="18%">1</TD>
        <TD width="32%">to</TD>
        <TD width="15%">1</TD>
        <TD width="17%">0</TD></TR></TBODY></TABLE>
    <P></P>
    <LI>To make a particle collide with a wall, simply leave its state 
    unchanged. This causes a reflection. 
    <P></P></LI></UL>The update code: <BR><PRE>	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) &amp; (sand(xind+1,yind+1)==1) &amp; ...
        (sand(xind+1,yind)==0) &amp; (sand(xind,yind+1)==0);
    
    diag2(xind,yind) = (sand(xind+1,yind)==1) &amp; (sand(xind,yind+1)==1) &amp; ...
        (sand(xind,yind)==0) &amp; (sand(xind+1,yind+1)==0);
    
    %The diagonals both not occupied by two particles
    and12(xind,yind) = (diag1(xind,yind)==0) &amp; (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)  &amp; ~sums(xind,yind) &amp; sand(xind+1,yind+1)) + ... 
        (or12(xind,yind) &amp; ~sums(xind,yind) &amp; sand(xind,yind+1)) + ...
        (sums(xind,yind) &amp; sand(xind,yind)); 
        
    sandNew(xind+1,yind) = ...
        (and12(xind,yind)  &amp; ~sums(xind,yind) &amp; sand(xind,yind+1)) + ... 
        (or12(xind,yind) &amp; ~sums(xind,yind) &amp; sand(xind,yind))+ ...
        (sums(xind,yind) &amp; sand(xind+1,yind));  
        
    sandNew(xind,yind+1) = ...    
        (and12(xind,yind)  &amp; ~sums(xind,yind) &amp; sand(xind+1,yind)) + ... 
        (or12(xind,yind) &amp; ~sums(xind,yind) &amp; sand(xind+1,yind+1))+ ...
        (sums(xind,yind) &amp; sand(xind,yind+1)); 
        
     sandNew(xind+1,yind+1) = ...    
        (and12(xind,yind)  &amp; ~sums(xind,yind) &amp; sand(xind,yind)) + ... 
        (or12(xind,yind) &amp; ~sums(xind,yind) &amp; sand(xind+1,yind))+ ...
        (sums(xind,yind) &amp; sand(xind+1,yind+1)); 
    
    sand = sandNew;
	</PRE>
  <LI><A 
  href="http://instruct1.cit.cornell.edu/courses/bionb441/CA/diffusion.m">Diffusion 
  limited aggregation</A> <BR>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.<BR>
  <P>The rule: <BR>
  <UL>
    <LI>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. 
    <LI>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. 
  </LI></UL>The update code: <BR><PRE>	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)&lt; .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&gt; 0) &amp; (sand==1)) | gnd ;
    %and eliminate the moving particle
    sand(find(gnd==1)) = 0;
	</PRE>The following image shows the fixed cluster after many time steps. 
  <P><IMG src="CA.files/Diffusion1.gif"> 
  <P></P>
  <LI><A href="http://instruct1.cit.cornell.edu/courses/bionb441/CA/sand.m">Sand 
  pile</A> <BR>The cross-section of a pile of sand can be modeled using a 
  Margolus neighborhood to propagate cells, but with a different motion rule<BR>
  <P>The rule: <BR>
  <UL>
    <LI>See reference [2] Chapter 2.2.6.. 
    <LI>Cells have 2 states. State=0 is empty, state=1 represents agrain of 
    sand. 
    <LI>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. 
    <LI>To make the motion slightly random, I added a rule which sometimes 
    reverses the state of the two lower cells after all the moves are 
    complete.<BR>
    <TABLE width="23%" border=2>
      <TBODY>
      <TR>
        <TD width="18%">0</TD>
        <TD width="18%">0</TD>
        <TD width="32%">goes</TD>
        <TD width="15%">0</TD>
        <TD width="17%">0</TD></TR>
      <TR>
        <TD width="18%">0</TD>
        <TD width="18%">0</TD>
        <TD width="32%">to</TD>
        <TD width="15%">0</TD>
        <TD width="17%">0</TD></TR></TBODY></TABLE>
    <P>
    <TABLE width="23%" border=2>
      <TBODY>
      <TR>
        <TD width="18%">1</TD>
        <TD width="18%">0</TD>
        <TD width="32%">goes</TD>
        <TD width="15%">0</TD>
        <TD width="17%">0</TD></TR>
      <TR>
        <TD width="18%">0</TD>
        <TD width="18%">0</TD>
        <TD width="32%">to</TD>
        <TD width="15%">1</TD>
        <TD width="17%">0</TD></TR></TBODY></TABLE>
    <P>
    <TABLE width="23%" border=2>
      <TBODY>
      <TR>
        <TD width="18%">0</TD>
        <TD width="18%">1</TD>
        <TD width="32%">goes</TD>
        <TD width="15%">0</TD>
        <TD width="17%">0</TD></TR>
      <TR>
        <TD width="18%">0</TD>
        <TD width="18%">0</TD>
        <TD width="32%">to</TD>
        <TD width="15%">0</TD>
        <TD width="17%">1</TD></TR></TBODY></TABLE>
    <P>
    <TABLE width="23%" border=2>
      <TBODY>
      <TR>
        <TD width="18%">1</TD>
        <TD width="18%">0</TD>
        <TD width="32%">goes</TD>
        <TD width="15%">0</TD>
        <TD width="17%">0</TD></TR>
      <TR>
        <TD width="18%">1</TD>
        <TD width="18%">0</TD>
        <TD width="32%">to</TD>
        <TD width="15%">1</TD>
        <TD width="17%">1</TD></TR></TBODY></TABLE>
    <P>
    <TABLE width="23%" border=2>
      <TBODY>
      <TR>
        <TD width="18%">1</TD>
        <TD width="18%">0</TD>
        <TD width="32%">goes</TD>
        <TD width="15%">0</TD>
        <TD width="17%">0</TD></TR>
      <TR>
        <TD width="18%">0</TD>
        <TD width="18%">1</TD>
        <TD width="32%">to</TD>
        <TD width="15%">1</TD>
        <TD width="17%">1</TD></TR></TBODY></TABLE>
    <P>
    <TABLE width="23%" border=2>
      <TBODY>
      <TR>
        <TD width="18%">0</TD>
        <TD width="18%">1</TD>
        <TD width="32%">goes</TD>
        <TD width="15%">0</TD>
        <TD width="17%">0</TD></TR>
      <TR>
        <TD width="18%">1</TD>
        <TD width="18%">0</TD>
        <TD width="32%">to</TD>
        <TD width="15%">1</TD>
        <TD width="17%">1</TD></TR></TBODY></TABLE>
    <P>
    <TABLE width="23%" border=2>
      <TBODY>
      <TR>
        <TD width="18%">0</TD>
        <TD width="18%">1</TD>
        <TD width="32%">goes</TD>
        <TD width="15%">0</TD>
        <TD width="17%">0</TD></TR>
      <TR>
        <TD width="18%">0</TD>
        <TD width="18%">1</TD>
        <TD width="32%">to</TD>
        <TD width="15%">1</TD>
        <TD width="17%">1</TD></TR></TBODY></TABLE>
    <P>
    <TABLE width="23%" border=2>
      <TBODY>
      <TR>
        <TD width="18%">1</TD>
        <TD width="18%">1</TD>
        <TD width="32%">goes</TD>
        <TD width="15%">1</TD>
        <TD width="17%">0</TD></TR>
      <TR>
        <TD width="18%">1</TD>
        <TD width="18%">0</TD>
        <TD width="32%">to</TD>
        <TD width="15%">1</TD>
        <TD width="17%">1</TD></TR></TBODY></TABLE>
    <P>
    <TABLE width="23%" border=2>
      <TBODY>
      <TR>
        <TD width="18%">1</TD>
        <TD width="18%">1</TD>
        <TD width="32%">goes</TD>
        <TD width="15%">0</TD>
        <TD width="17%">1</TD></TR>
      <TR>
        <TD width="18%">0</TD>
        <TD width="18%">1</TD>
        <TD width="32%">to</TD>
        <TD width="15%">1</TD>
        <TD width="17%">1</TD></TR></TBODY></TABLE>
    <P>The update code: <BR><PRE>	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)&lt; .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;
	</PRE></LI></UL></LI></OL>
<HR>

<P><B>References</B></P>
<P>[1] <I>Cellular Automata Machines</I> by Tommaso Toffoli and Norman Margolus, 
MIT Press, 1987.</P>
<P>[2] <I>Cellular Automata Modeling of Physical Systems</I> by Bastien Chopard 
and Michel Droz, Cambridge University Press, 1998.</P></BODY></HTML>