美国康奈尔大学BioNB441在Matlab中的元胞自动机

介绍

元胞自动机（CA）是用于计算计划利用当地的规则和本地通信。普遍CA定义一个网格，网格上的每个点代表一个有限数量的状态中的细胞。过渡规则同时应用到每一个细胞。 典型的转换规则依赖于细胞和它的（4个或8个）近邻的状态，虽然临近的细胞也同样使用。 CA的应用在并行计算研究，物理模拟和生物模拟。这个页面将考虑如何写出高效的MATLAB代码的CA的实施和看一些有趣的规则。
Matlab代码注意事项
以下注意事项将说明使用Matlab程序计算康威的生命。部分的代码的解释如下。
矩阵和图像可以被转换为一个另一个，所以显示是为straighforward的的。如果阵列细胞包含二进制的每一个电池单体的状态和数组Z含有零，那么cat命令建立一个RGB图像，显示图像命令。图像命令也返回一个句柄。

imh = image(cat(3,cells,z,z));
set(imh, 'erasemode', 'none')
axis equal
axis tight

矩阵和图像可以被转换为一个另一个，那么初始条件，可以计算矩阵或图形命令。下面的代码生成一个元素为零的数组，细胞状态初始化为零，然后使细胞的交叉状态=1的数组的中心。使用（渗流群集）下面的例子中的一个图形命令来初始化CA阵列。

z = zeros(n,n);
cells = z;
cells(n/2,.25*n:.75*n) = 1;
cells(.25*n:.75*n,n/2) = 1;

Matlab代码进行矢量化，以减少开销。下面的两个语句计算的总和近邻的计算CA规则。代码利用Matlab的非常灵活的索引指定的近邻。

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);  

添加一个简单的GUI是很容易的。在这个例子中，三个按钮和一个文本字段implmented。三个按钮允许用户运行，停止，或选择退出。文本字段中显示模拟执行的步骤数。

%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]);

后控制（和CA）初始化，程序进入一个循环的测试中设置的每个按钮的回调函数的变量的状态下降。就目前而言，只是看在while循环，if语句的嵌套结构。程序循环，直到按下退出按钮。另外两个按钮导致当推一个if语句来执行。

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

例子

1.康威的生命。(Conway’s life with GUI)

其规则是：

• 总结8近邻
• 如果总和=2 那么状态不改变
• 如果总和=3，然后在状态= 1
• 否则状态= 0

完整代码：

%Conway's life with GUI

clf
clear all
clc

%=============================================
%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]);


%=============================================
%CA setup

n=128;

%initialize the arrays
z = zeros(n,n);
cells = z;
sum = z;
%set a few cells to one
cells(n/2,.25*n:.75*n) = 1;
cells(.25*n:.75*n,n/2) = 1;

%cells(.5*n-1,.5*n-1)=1;
%cells(.5*n-2,.5*n-2)=1;
%cells(.5*n-3,.5*n-3)=1;
cells = (rand(n,n))<.5 ;
%how long for each case to stability or simple oscillators

%build an image and display it
imh = image(cat(3,cells,z,z));
set(imh, 'erasemode', 'none')
axis equal
axis tight

%index definition for cell update
x = 2:n-1;
y = 2:n-1;

%Main event loop
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

程序运行结果：

第一次测试初始状态

基本静止状态

第二次测试初始状态

运行中状态

运行中状态

基本静止状态

2.表面张力（Surface Tension ）

其规则是：

• 总结8近邻和细胞本身
• 如果总和<4或总和= 5，则状态= 0
• 否则状态= 1

完整代码：

%Surface Tension with GUI

clf
clear all
clc

%=============================================
%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]);


%=============================================
%CA setup

n=128;

%initialize the arrays
z = zeros(n,n);
cells = z;
sum = z;
%set a few cells to one
cells(n/2,.25*n:.75*n) = 1;
cells(.25*n:.75*n,n/2) = 1;

%cells(.5*n-1,.5*n-1)=1;
%cells(.5*n-2,.5*n-2)=1;
%cells(.5*n-3,.5*n-3)=1;
cells = (rand(n,n))<.5 ;
%how long for each case to stability or simple oscillators

%build an image and display it
imh = image(cat(3,cells,z,z));
set(imh, 'erasemode', 'none')
axis equal
axis tight

%index definition for cell update
x = 2:n-1;
y = 2:n-1;

%Main event loop
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)+...
            cells(x,y);
        % The CA rule
        cells = ~((sum<4) | (sum==5));       
        %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

程序运行结果

初始状态

第75次状态

第481次状态

第2601次状态，静止不在变化

3.渗流集群（Percolation Cluster）

规则：

• 每个单元格的8个最近邻的求和（细胞是二进制值，0/1）。细胞也有一个独立的状态变量（称为’参观’）记录，如果他们曾经有过一个非零的前邻居。
• 计算0和1之间的随机数r。
• 如果总和> 0（至少一个相邻）的r>阈值，和小区从未有过的邻居小区= 1。
• 如果总和>0记录细胞，该细胞具有非零邻居“访问”的标志。
• 变量a和b的图像的大小。 >初始的图像是由图形操作。设立轴大小是固定的，下面的语句写文字转换成轴，然后抢轴的内容，并把它们放回一个数组，用的getFrame功能.

完整代码：

%Percolation Cluster with GUI
clf;clc
clear all
threshold = .63;
%
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;

imh = image(cat(3,z,cells,z));
set(imh, 'erasemode', 'none')

%return

for i=1:100
    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); 
    %edges only
    %cells = (cells & (sum<8)) | ((sum>=1) & (pick>=threshold) & (visit==0))  ;
    cells = cells  | ((sum>=1) & (pick>=threshold) & (visit==0))  ;
    visit = (sum>=1) ;%& (pick<threshold)  ;

    set(imh, 'cdata', cat(3,z,cells,z) )
    drawnow
end

return
figure(2)
image(cat(3,z,cells,z))

经过几十个时间间隔（从BioNB 441这个图像开始） ，我们可以得到以下的图像。

4.激发介质（BZ反应或心脏）（Excitable media (BZ reaction or heart) ）

规则：

• 细胞可以在10个不同的状态。休息状态0。国家1-5活跃，状态6-9难治。
• 计数8个最近邻的每个单元格是在一个激活状态。
• 如果总和大于或等于3（至少三个活性邻居的），那么细胞=1。
• 状态1〜9中，没有更多的输入逐步发生。如果状态=1，那么下一个状态=2。如果状态= 2那么下一状态= 3，同样的所有的州多达9个。如果状态= 9，那么下一个状态=0和细胞是在休息。

完整代码：

%
%CA driver
%
%excitable media
clc
clf
clear all

n=128;

z=zeros(n,n);
cells=z;

cells = (rand(n,n))<.1 ;
%cells(n/2,n*.25:n*.75) = 1;
%cells(n*.25:n*.75,n/2) = 1;
sum=z;

imh = image(cat(3,cells,z,z));
set(imh, 'erasemode', 'none')
axis equal
axis tight

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

t = 6; % center value=6; 7 makes fast pattern; 5 analiating waves
t1 = 3; % center value=3
for i=1:1200


    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);

    set(imh, 'cdata', cat(3,z,cells/10,z) )
    drawnow
end

一个CA初始图形经过螺旋的变化，得到下图

再次运行得到下图

5.森林火灾(Forest Fire )

规则：

• 细胞可以在3个不同的状态。状态= 0是空的，state = 1时燃烧，状态=2是森林。
• 如果一个或更多的4个邻居，如果细胞在燃烧，这是森林（状态= 2），那么新的状态（状态 = 1时）燃烧。
• 有森林细胞（状态=2），开始自身的燃烧（闪电）（0.000005）是一种低概率。
• 燃烧（state = 1时）后一个细胞变成空的（状态= 0）。
• 有一个空单元格成为森林模拟增长是一种低概率（比如说，0.01）。
• 该阵列被视为toroidly连接的，使燃烧的火焰，以左侧将启动火灾在右边。类似地连接的顶部和底部。

完整代码：

%
%CA driver
%
%forest fire

clf
clear all

n=100;

Plightning = .000005;
Pgrowth = .01; %.01

z=zeros(n,n);
o=ones(n,n);
veg=z;
sum=z;


imh = image(cat(3,z,veg*.02,z));
set(imh, 'erasemode', 'none')
axis equal
axis tight

% burning -> empty
% green -> burning if one neigbor burning or with prob=f (lightning)
% empty -> green with prob=p (growth)
% veg = {empty=0 burning=1 green=2}
for i=1:3000
    %nearby fires?

     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) ;

    set(imh, 'cdata', cat(3,(veg==1),(veg==2),z) )
    drawnow
end

变化图像

6.气体动力学（Gas dynamics）

此CA（及之后的两个）被用来计算粒子的运动。此应用程序需要不同种类的邻里。瓦里斯的小区附近的每一个时间步长，使运动在一定的方向上，将继续在同一方向。换句话说，规则节省势头，这是根据动力计算。所用的邻域被称为一个的马戈利斯附近，由重叠的2×2块的细胞。在下面的图中，左上角的4细胞是邻近地区的偶数时间步长期间，在奇数时间步长，右下4个细胞分数。一个给定的单元有3个邻居在每个时间步长，但具体的细胞构成的邻居来回翻转。

规则：

• 这就是所谓的HPP-GAS规则。见参考文献[1]第12章..
• 细胞有2个状态。状态= 0是空的，状态=1表示运动中的粒子。
• 在任何一个步骤中，粒子被假设刚才输入的2x2块对角线路径。通过它的中心，因此，它必须跨块。因此，在任何时间步长，每个单元的内容交换与细胞斜对面。对块如下所示。左边的一个是前一段时间后步的权利。六例示，但他们每个人的所有旋转版本相同的方式处理。两种特殊的情况，颗粒与颗粒的碰撞和粒子壁被认为是如下。

• 为了使粒子碰撞（节省动量和能量），治疗的情况下，如果他们正好一个对角线上的两个粒子击中对方90度偏转。这是通过将一条对角线上的其他时间步长。您可以实现由一个细胞的细胞逆时针旋转。上述第三条规则变为：

• 为了使粒子碰撞与墙，只需将其状态不变。这导致的一种体现。

完整代码：

%CA driver
%HPP-gas

clc
clear all
clf

nx=52; %must be divisible by 4
ny=100;

z=zeros(nx,ny);
o=ones(nx,ny);
sand = z ;
sandNew = z;
gnd = z ;
diag1 = z;
diag2 = z;
and12 = z;
or12 = z;
sums = z;
orsum = z;

gnd(1:nx,ny-3)=1 ; % right ground line
gnd(1:nx,3)=1 ; % left ground line
gnd(nx/4:nx/2-2,ny/2)=1; %the hole line
gnd(nx/2+2:nx,ny/2)=1; %the hole line
gnd(nx/4, 1:ny) = 1; %top line
gnd(3*nx/4, 1:ny) = 1 ;%bottom line

%fill the left side
r = rand(nx,ny);
sand(nx/4+1:3*nx/4-1, 4:ny/2-1) = r(nx/4+1:3*nx/4-1, 4:ny/2-1)<0.3;
%sand(nx/4+1:3*nx/4-1, ny*.75:ny-4) = r(nx/4+1:3*nx/4-1, ny*.75:ny-4)<0.75;
%sand(nx/2,ny/2) = 1;
%sand(nx/2+1,ny/2+1) = 1;

imh = image(cat(3,z,sand,gnd));
set(imh, 'erasemode', 'none')
axis equal
axis tight


for i=1:1000
    p=mod(i,2); %margolis neighborhood

    %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;

    set(imh, 'cdata', cat(3,z,sand,gnd) )
    drawnow
end

程序运行结果

7.扩散限制聚集

该系统模拟粘性颗粒聚集形成的分形结构的。粒子的运动发生，例如6 HPP-气体规则类似的规则。所不同的是该颗粒被认为是在某些致密的（但不可见）的液体蹦跳着。其效果是随机的每一个粒子的运动方向在每个时间步长。换句话说，每一个时间步长的碰撞的步骤。仿真也接种与固定在阵列中心的粒子。任何接触它的散射粒子粘到它，本身成为一个不可移动，粘性颗粒。
规则：

• 使用附近Margolus。在每一个时间步长，旋转4顺时针或逆时针以相等的概率由一个细胞的细胞。旋转随机速度。
• 移动后，如果一个或多个的八个最接近的邻居的是一个固定的，粘性颗粒，然后冻结粒子，并使其粘。

完整代码：

%diffusion + dla
clc
clear all
clf

nx=200; %must be divisible by 4
ny=200;

z=zeros(nx,ny);
o=ones(nx,ny);
sand = z ;
sandNew = z;
sum = z;
gnd = z;
gnd(nx/2,ny/2) = 1 ;

sand = rand(nx,ny)<.1;

imh = image(cat(3,z,sand,gnd));
set(imh, 'erasemode', 'none')
axis equal
axis tight



for i=1:10000
    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;

    set(imh, 'cdata', cat(3,gnd,gnd,(sand==1)) )
    drawnow
end

程序运行结果：

8.砂桩

一堆沙子的横截面可以仿照使用一个Margolus的邻里传播细胞，但具有不同的运动规律
规则：

• 见参考文献[2]第2.2.6章..
• 细胞有2个状态。状态= 0是空的，状态=1表示沙agrain。
• 在任何步骤中，粒子可以落在朝下方的2×2的块。可能的转换如下所示。墙壁和地板停止运动。
• 为了使议案略有随机，我添加了一个规则，有时反转的两个较低的细胞状态，所有的动作都完成后。

完整代码：

%sand pile
clear all
clf

nx=52; %must be divisible by 4
ny=100;

Pbridge = .05; 

z=zeros(nx,ny);
o=ones(nx,ny);
sand = z ;
sandNew = z;
gnd = z ;
gnd(1:nx,ny-3)=1 ; % the ground line
gnd(nx/4:nx/2+4,ny-15)=1; %the hole line
gnd(nx/2+6:nx,ny-15)=1; %the hole line
gnd(nx/4, ny-15:ny) = 1; %side line
gnd(3*nx/4, 1:ny) = 1 ;

imh = image(cat(3,z',sand',gnd'));
set(imh, 'erasemode', 'none')
axis equal
axis tight


for i=1:1000
    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];
    vary = rand(nx,ny)<.95 ;
    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;
    set(imh, 'cdata', cat(3,z',sand',gnd') )
    drawnow
end

程序运行结果：

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.