元胞自动机与MATLAB引言元胞自动机（CA）是一种用来仿真局部规则和局部联系的方法。
典型的元胞自动机是定义在网格上的，每一个点上的网格代表一个元胞与一种有限的状态。
变化规则适用于每一个元胞并且同时进行。
典型的变化规则，决定于元胞的状态，以及其（4或8 ）邻居的状态。
元胞自动机已被应用于物理模拟，生物模拟等领域。
本文就一些有趣的规则，考虑如何编写有效的MATLAB的程序来实现这些元胞自动机。
MATLAB的编程考虑元胞自动机需要考虑到下列因素，下面分别说明如何用MATLAB实现这些部分。
并以Conway的生命游戏机的程序为例，说明怎样实现一个元胞自动机。
●矩阵和图像可以相互转化，所以矩阵的显示是可以真接实现的。
如果矩阵cells的所有元素只包含两种状态且矩阵Z含有零，那么用image函数来显示cat命令建的RGB图像，并且能够返回句柄。
imh = image(cat(3,cells,z,z));set(imh, 'erasemode', 'none')axis equalaxis tight●矩阵和图像可以相互转化，所以初始条件可以是矩阵，也可以是图形。
以下代码生成一个零矩阵，初始化元胞状态为零，然后使得中心十字形的元胞状态= 1。
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);●加入一个简单的图形用户界面是很容易的。
在下面这个例子中，应用了三个按钮和一个文本框。
三个按钮，作用分别是运行，停止，程序退出按钮。
文框是用来显示的仿真运算的次数。
%build the GUI%define the plot buttonplotbutton=uicontrol('style','pushbutton',...'string','Run', ...'fontsize',12, ...'position',[100,400,50,20], ...'callback', 'run=1;');%define the stop buttonerasebutton=uicontrol('style','pushbutton',...'string','Stop', ...'fontsize',12, ...'position',[200,400,50,20], ...'callback','freeze=1;');%define the Quit buttonquitbutton=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 pushrun = 0; %wait for a drawfreeze = 0; %wait for a freezewhile (stop==0)if (run==1)%nearest neighbor sumsum(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 rulecells = (sum==3) | (sum==2 & cells);%draw the new imageset(imh, 'cdata', cat(3,cells,z,z) )%update the step number diaplaystepnumber = 1 + str2num(get(number,'string'));set(number,'string',num2str(stepnumber))endif (freeze==1)run = 0;freeze = 0;enddrawnow %need this in the loop for controls to workend例子1 .Conway的生命游戏机。
规则是：➢对周围的8个近邻的元胞状态求和➢如果总和为2的话，则下一时刻的状态不改变➢如果总和为3 ，则下一时刻的状态为1➢否则状态= 0核心代码：x = 2:n-1;y = 2:n-1;%nearest neighbor sumsum(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 rulecells = (sum==3) | (sum==2 & cells);2 .表面张力规则是：➢对周围的8近邻的元胞以及自身的状态求和➢如果总和< 4或= 5 ，下一时刻的状态= 0➢否则状态= 1核心代码：x = 2:n-1;y = 2:n-1;%nearest neighbor sumsum(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 rulecells = ~((sum< 4) | (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',[1 1 500 400],'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','Cellular Automata','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;经过几十个时间间隔（从MCM Cellular Automata 这个图像开始） ，我们可以得到以下的图像。