元胞自动机与Matlab.docx
《元胞自动机与Matlab.docx》由会员分享，可在线阅读，更多相关《元胞自动机与Matlab.docx（18页珍藏版）》请在冰豆网上搜索。
元胞自动机与Matlab
元胞自动机与MATLAB
引言
元胞自动机（CA）是一种用来仿真局部规则和局部联系的方法。
典型的元胞自动机是定义在网格上的，每一个点上的网格代表一个元胞与一种有限的状态。
变化规则适用于每一个元胞并且同时进行。
典型的变化规则，决定于元胞的状态，以及其（4或8）邻居的状态。
元胞自动机已被应用于物理模拟，生物模拟等领域。
本文就一些有趣的规则，考虑如何编写有效的MATLAB的程序来实现这些元胞自动机。
MATLAB的编程考虑
元胞自动机需要考虑到下列因素，下面分别说明如何用MATLAB实现这些部分。
并以Conway的生命游戏机的程序为例，说明怎样实现一个元胞自动机。
●矩阵和图像可以相互转化，所以矩阵的显示是可以真接实现的。
如果矩阵cells的所有元素只包含两种状态且矩阵Z含有零，那么用image函数来显示cat命令建的RGB图像，并且能够返回句柄。
imh=image（cat（3,cells,z,z））;
set（imh,'erasemode','none'）
axisequal
axistight
●矩阵和图像可以相互转化，所以初始条件可以是矩阵，也可以是图形。
以下代码生成一个零矩阵，初始化元胞状态为零，然后使得中心十字形的元胞状态=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）;
●加入一个简单的图形用户界面是很容易的。
在下面这个例子中，应用了三个按钮和一个文本框。
三个按钮，作用分别是运行，停止，程序退出按钮。
文框是用来显示的仿真运算的次数。
%buildtheGUI
%definetheplotbutton
plotbutton=uicontrol（'style','pushbutton',...
'string','Run',...
'fontsize',12,...
'position',[100,400,50,20],...
'callback','run=1;'）;
%definethestopbutton
erasebutton=uicontrol（'style','pushbutton',...
'string','Stop',...
'fontsize',12,...
'position',[200,400,50,20],...
'callback','freeze=1;'）;
%definetheQuitbutton
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;%waitforaquitbuttonpush
run=0;%waitforadraw
freeze=0;%waitforafreeze
while（stop==0）
if（run==1）
%nearestneighborsum
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）;
%TheCArule
cells=（sum==3）|（sum==2&cells）;
%drawthenewimage
set（imh,'cdata',cat（3,cells,z,z））
%updatethestepnumberdiaplay
stepnumber=1+str2num（get（number,'string'））;
set（number,'string',num2str（stepnumber））
end
if（freeze==1）
run=0;
freeze=0;
end
drawnow%needthisintheloopforcontrolstowork
end
例子
1.Conway的生命游戏机。
规则是：
Ø对周围的8个近邻的元胞状态求和
Ø如果总和为2的话，则下一时刻的状态不改变
Ø如果总和为3，则下一时刻的状态为1
Ø否则状态=0
核心代码：
x=2:
n-1;
y=2:
n-1;
%nearestneighborsum
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）;
%TheCArule
cells=（sum==3）|（sum==2&cells）;
2.表面张力
规则是：
Ø对周围的8近邻的元胞以及自身的状态求和
Ø如果总和<4或=5，下一时刻的状态=0
Ø否则状态=1
核心代码：
x=2:
n-1;
y=2:
n-1;
%nearestneighborsum
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）;
%TheCArule
cells=~（（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',[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）0）&（cells（x,y+1）（（cells（x-1,y）>0）&（cells（x-1,y）0）&（cells（x+1,y）（（cells（x-1,y-1）>0）&（cells（x-1,y-1）0）&（cells（x-1,y+1）（（cells（x+1,y-1）>0）&（cells（x+1,y-1）0）&（cells（x+1,y+1）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）&（sum>0|（rand（n,n）2*（（veg==0）&rand（n,n）注意环形连接是由序标实现的。
6.气体动力学
这个CA（以及接下来的两个CA）是用来模拟粒子运动的。
此元胞自动机需要一种不同类型的元胞的邻居。
此元胞的邻居时刻变化，因此某一个方向运动趋势，将继续在同一个方向。
换言之，此规则保存势头，这是基础的动力仿真。
这种邻居通常被称为margolis邻居并且这种邻居通常由重叠的2x2块的元胞构成。
在下面的表格中，偶数步长时左上方4元胞为邻居关系，奇数步长时右下的4元胞为邻居关系。
某一特定元胞在每一个时间步长都有3个邻居，但是具体的元胞构成了邻居的旋转和反复。
偶
偶
偶
元
胞
奇
奇
奇
规则：
Ø此规则叫作HPP-气体规则。
Ø每个元胞有2种状态。
状态=0是空的，状态=1代表粒子。
Ø在任何一个时间步长，假设粒子是刚刚进入2x2的网格块。
它将通过其网格块的中心到达对角的网格中，所以在任何时间步长，每一个元胞与该元胞对角对元胞交换的内容。
如下所示，左边显示出来的元胞结构经过一个时间步长变为右边的结构。
以下是六种不同的情况，所有所有的元胞都遵循相同的转动规则。
下文还将考虑两种特殊情况，即粒子-粒子碰撞和粒子-墙碰撞。
0
0
0
0
1
0
0
0
0
0
0
0
0
0
0
1
1
0
1
0
1
0
0
1
0
1
0
1
1
0
0
1
1
1
0
1
1
1
1
1
1
0
1
1
1
1
1
1
Ø为了实现粒子碰撞过程（保证动量和能量守恒），对于两个处于对角线上的粒子，他们相互撞击后偏转90度。
在一个时间步长里使其从一个对角转成另一个对角。
你可以逆时针旋转这四个元胞来实现这个过程。
则第三规则可以表示为：
1
0
0
1
0
1
1
0
Ø粒子撞击墙壁时，简单地使其离开且状态不变。
这就引起反射现象。
核心代码：
p=mod（i,2）;%margolisneighborhood,whereiisthetimestep
%upperleftcellupdate
xind=[1+p:
2:
nx-2+p];
yind=[1+p:
2:
ny-2+p];
%Seeifexactlyonediagonalisones
%only（atmost）oneofthefollowingcanbetrue!
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）;
%Thediagonalsbothnotoccupiedbytwoparticles
and12（xind,yind）=（diag1（xind,yind）==0）&（diag2（xind,yind）==0）;
%Onediagonalisoccupiedbytwoparticles
or12（xind,yind）=diag1（xind,yind）|diag2（xind,yind）;
%foreverygasparticleseeifitneartheboundary
sums（xind,yind）=gnd（xind,yind）|gnd（xind+1,yind）|...
gnd（xind,yind+1）|gnd（xind+1,yind+1）;
%celllayout:
%x,yx+1,y
%x,y+1x+1,y+1
%If（nowalls）and（diagonalsarebothnotoccupied）
%thenthereisnocollision,somoveoppositecelltocurrentcell
%If（nowalls）and（onlyonediagonalisoccupied）
%thenthereisacollisionsomoveccwcelltothecurrentcell
%If（awall）
%thendon'tchangethecell（causesareflection）
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;
8.扩散限制聚集
这个系统是模拟粘性颗粒的聚集，形成分形结构。
质点以一个类似于例6中的HPP-气体规则发生运动。
不同的是粒子在一些高密度（但看不见）的液体周围被假定是弹跳的。
效果是每一个粒子在每个时间步长在随机的方向上运动。
换言之，每一个时间步长是一个碰撞的过程。
这个仿真矩阵的中心确定了在一个固定生长颗粒。
任何弥散粒子触及它就会被它粘住，并成为一个不能移动的，有粘性颗粒。
规则：
Ø使用Margolus型邻居。
在每一个时间步，等概率地顺时针或逆时针旋转4个元胞。
旋转使速度随机化。
Ø在移动后，如果八个最近的邻居有一个或一个以上元胞是固定的粘性颗粒，则下时刻该元胞将被冻结，并且使之有粘性。
核心代码：
p=mod（i,2）;%margolisneighborhood
%upperleftcellupdate
xind=[1+p:
2:
nx-2+p];
yind=[1+p:
2:
ny-2+p];
%randomvelocitychoice
vary=rand（nx,ny）<.5;
vary1=1-vary;
%diffusionrule--margolusneighborhood
%rotatethe4cellstorandomizevelocity
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;
%foreverysandgrainseeifitnearthefixed,stickycluster
sum
升级会员 