Matlab编写元胞自动机程序

元胞自动机与MATLAB引言元胞自动机（CA）是一种用来仿真局部规则和局部联系的方法。典型的元胞自动机是定义在网格上的，每一个点上的网格代表一个元胞与一种有限的状态。变化规则适用于每一个元胞并且同时进行。典型的变化规则，决定于元胞的状态，以及其（4或8）邻居的状态。元胞自动机已被应用于物理模拟，生物模拟等领域。本文就一些有趣的规则，考虑如何编写有效的MATLAB的程序来实现这些元胞自动机。MATLAB的编程考虑元胞自动机需要考虑到下列因素，下面分别说明如何用MATLAB实现这些部分。并以Conway的生命游戏机的程序为例，说明怎样实现一个元胞自动机。z矩阵和图像可以相互转化，所以矩阵的显示是可以真接实现的。如果矩阵cells的所有元素只包含两种状态且矩阵Z含有零，那么用image函数来显示cat命令建的RGB图像，并且能够返回句柄。imh=image(cat(3,cells,z,z));set(imh,'erasemode','none')axisequalaxistightz矩阵和图像可以相互转化，所以初始条件可以是矩阵，也可以是图形。以下代码生成一个零矩阵，初始化元胞状态为零，然后使得中心十字形的元胞状态=1。z=zeros(n,n);cells=z;cells(n/2,.25*n:.75*n)=1;cells(.25*n:.75*n,n/2)=1;zMatlab的代码应尽量简洁以减小运算量。以下程序计算了最近邻居总和，并按照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);z加入一个简单的图形用户界面是很容易的。在下面这个例子中，应用了三个按钮和一个文本框。三个按钮，作用分别是运行，停止，程序退出按钮。文框是用来显示的仿真运算的次数。%buildtheGUI%definetheplotbuttonplotbutton=uicontrol('style','pushbutton',...'string','Run',...'fontsize',12,...'position',[100,400,50,20],...'callback','run=1;');%definethestopbuttonerasebutton=uicontrol('style','pushbutton',...'string','Stop',...'fontsize',12,...'position',[200,400,50,20],...'callback','freeze=1;');%definetheQuitbuttonquitbutton=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;%waitforaquitbuttonpushrun=0;%waitforadrawfreeze=0;%waitforafreezewhile(stop==0)if(run==1)%nearestneighborsumsum(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);%TheCArulecells=(sum==3)|(sum==2&cells);%drawthenewimageset(imh,'cdata',cat(3,cells,z,z))%updatethestepnumberdiaplaystepnumber=1+str2num(get(number,'string'));set(number,'string',num2str(stepnumber))endif(freeze==1)run=0;freeze=0;enddrawnow%needthisintheloopforcontrolstoworkend例子1.Conway的生命游戏机。规则是：¾对周围的8个近邻的元胞状态求和¾如果总和为2的话，则下一时刻的状态不改变¾如果总和为3，则下一时刻的状态为1¾否则状态=0核心代码：x=2:n-1;y=2:n-1;%nearestneighborsumsum(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);%TheCArulecells=(sum==3)|(sum==2&cells);2.表面张力规则是：¾对周围的8近邻的元胞以及自身的状态求和¾如果总和4或=5，下一时刻的状态=0¾否则状态=1核心代码：x=2:n-1;y=2:n-1;%nearestneighborsumsum(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);%TheCArulecells=~((sum4)|(sum==5));3.渗流集群规则：¾对周围相邻的8邻居求和（元胞只有两种状态，0或1）。元胞也有一个单独的状态参量（所谓'记录'）记录它们之前是否有非零状态的邻居。¾在0与1之间产生一个随机数r。¾如果总和0（至少一个邻居）并且r阈值，或者元胞从未有过一个邻居，则元胞=1。¾如果总和0则设置记录的标志，记录这些元胞有一个非零的邻居。核心代码：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);变量a和b是图像的尺寸。最初的图形是由图形操作决定的。以下程序设定坐标系为一个固定的尺寸，在坐标系里写入文本，然后获得并返回坐标内的内容，并用getframe函数把它们写入一个矩阵ax=axes('units','pixels','position',[11500400],'color','k');text('units','pixels','position',[130,255,0],...'string','MCM','color','w','fontname','helvetica','fontsize',100)text('units','pixels','position',[10,120,0],...'string','CellularAutomata','color','w','fontname','helvetica','fontsize',50)initial=getframe(gca);[a,b,c]=size(initial.cdata);z=zeros(a,b);cells=double(initial.cdata(:,:,1)==255);visit=z;sum=z;经过几十个时间间隔（从MCMCellularAutomata这个图像开始），我们可以得到以下的图像。4.激发介质（BZreactionorheart）规则：¾元胞有10个不同的状态。状态0是体眠。1-5为活跃状态，、6-9为是极活跃状态。¾计算每一个处于活跃的状态的元胞近邻的8个元胞。¾如果和大于或等于3（至少有三个活跃的邻居），则下一时刻该元胞=1。¾不需要其它输入，1至9种状态依次出现。如果该时刻状态=1那么下一时刻状态=2。如果该时刻状态=2，然后下一时刻状态=3，对于其它的状态依次类推，直到第9种状态。如果状态=9，然后下一状态=0并且元胞回到休息状态。核心代码：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);一个CA初始图形经过螺旋的变化，得到下图。5.森林火灾规则：¾元胞有3个不同的状态。状态为0是空位，状态=1是燃烧着的树木，状态=2是树木。¾如果4个邻居中有一个或一个以上的是燃烧着的并且自身是树木（状态为2），那么该元胞下一时刻的状态是燃烧（状态为1）。¾森林元胞（状态为2）以一个低概率（例如0.000005）开始烧（因为闪电）。¾一个燃烧着的元胞（状态为1）在下一时时刻变成空位的（状态为0）。¾空元胞以一个低概率（例如0.01）变为森林以模拟生长。¾出于矩阵边界连接的考虑，如果左边界开始着火，火势将向右蔓延，右边界同理。同样适用于顶部和底部。核心代码：sum=(veg(1:n,[n1:n-1])==1)+(veg(1:n,[2:n1])==1)+...(veg([n1:n-1],1:n)==1)+(veg([2:n1],1:n)==1);veg=...2*(veg==2)-((veg==2)&(sum0|(rand(n,n)Plightning)))+...2*((veg==0)&rand(n,n)Pgrowth);注意环形连接是由序标实现的。6.气体动力学这个CA（以及接下来的两个CA）是用来模拟粒子运动的。此元胞自动机需要一种不同类型的元胞的邻居。此元胞的邻居时刻变化，因此某一个方向运动趋势，将继续在同一个方向。换言之，此规则保存势头，这是基础