《元胞自动机与Matlab.doc》由会员分享，可在线阅读，更多相关《元胞自动机与Matlab.doc（14页珍藏版）》请在三一文库上搜索。
1、元胞自动机与MATLAB引言元胞自动机（CA）是一种用来仿真局部规则和局部联系的方法。典型的元 胞自动机是定义在网格上的，每一个点上的网格代表一个元胞与一种有限的状 态。变化规则适用于每一个元胞并且同时进行。 典型的变化规则，决定于元胞的 状态，以及其（4或8 ）邻居的状态。元胞自动机已被应用于物理模拟，生物 模拟等领域。本文就一些有趣的规则，考虑如何编写有效的MATLAB的程序来实现这些元胞自动机。MATLAB的编程考虑 元胞自动机需要考虑到下列因素，下面分别说明如何用MATLAB实现这些部分 并以Conway的生命游戏机的程序为例，说明怎样实现一个元胞自动机。矩阵和图像可以相互转化，所以矩
2、阵的显示是可以真接实现的。如果矩阵cells的所有元素只包含两种状态且矩阵 Z含有零，那么用image函数来显 示cat命令建的RGB图像，并且能够返回句柄。imh = image(cat(3,cells,z,z); set(imh, 'erasemode', 'non e') axis equalaxis tight矩阵和图像可以相互转化，所以初始条件可以是矩阵，也可以是图形。以下 代码生成一个零矩阵，初始化元胞状态为零，然后使得中心十字形的元胞状 态=1。z = zeros (n,n);cells = z;cells( n/2,.25* n:.75* n)
3、= 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 &
4、cells);加入一个简单的图形用户界面是很容易的。在下面这个例子中，应用了三个 按钮和一个文本框。三个按钮，作用分别是运行，停止，程序退出按钮。文 框是用来显示的仿真运算的次数。%build the GUI%defi ne the plot butt onplotbutt on=uic on trol('style','pushbutt on',.'stri ng','Ru n',.'fon tsize',12, .'positio n',100,400,50,20, .'callback
5、', 'run=1;');%defi ne the stop butt onerasebutt on=uic on trol('style','pushbutt on',.'stri ng','Stop', .'fon tsize',12, .'positio n',200,400,50,20, .'callback','freeze=1;');%defi ne the Quit butt onquitbutt on=uic on trol(&
6、#39;style','pushbutt on',.'stri ng','Quit', .'fon tsize',12, .'positio n',300,400,50,20, .'callback','stop=1;close;');nu mber = uic on trol('style','text', .'stri ng','1', .'fon tsize',12, .'posit
7、io n',20,400,50,20）;经过对控件（和CA）初始化，程序进入一个循环，该循环测试由回调函数的每 个按钮控制的变量。刚开始运行时，只在嵌套的 while循环和if语句中运行。 直到退出按钮按下时，循环停止。另外两个按钮按下时执行相应的if语句。stop= 0; %wait for a quit butt on pushrun = 0; %wait for a drawfreeze = 0; %wait for a freezewhile (stop=0)if (run=1)%n earest n eighbor sumsum(x,y) = cells(x,y-1) + c
8、ells(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 nu mber diaplay step nu mber = 1 + str2 nu m(get (nu
9、mber,'stri ng'); set (nu mber,'stri ng', nu m2str(step nu mber)endif (freeze=1)run = 0;freeze = 0;enddraw now%n eed this in the loop for con trols to workend例子1 .Conway的生命游戏机。规则是:? 对周围的8个近邻的元胞状态求和? 如果总和为2的话，则下一时刻的状态不改变? 如果总和为3，则下一时刻的状态为1? 否则状态=0核心代码:x = 2:n-1;y = 2:n-1;%n earest n ei
10、ghbor 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;%n earest
11、 n eighbor 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之间产生
12、一个随机数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 = ran d(a,b);cells = ce
13、lls | (sum>=1) & (pick>=threshold) & (visit=O)visit = (sum>=1);变量a和b是图像的尺寸。最初的图形是由图形操作决定的。以下程序设定坐 标系为一个固定的尺寸，在坐标系里写入文本，然后获得并返回坐标内的内容， 并用getframe函数把它们写入一个矩阵ax = axes('u nits','pixels','positio n',1 1 500 400,'color','k');text('u nits',
14、'pixels', 'positi on', 130,255,0,.'stri ng','MCM','color','w','fo ntn ame','helvetica','fo ntsize',100)text(' un its', 'pixels', 'positi on', 10,120,0,.'stri ng','Cellular Automata','c
15、olor','w','fo ntn ame','helvetica','fo ntsize',50)in itial = getframe(gca);a,b,c=size(i nitial.cdata);z=zeros(a,b);cells = double(i nitial.cdata(:,:,1)=255);visit = z ;sum = z;经过几十个时间间隔(从 MCM Cellular Automata这个图像开始)，我们可以得到以下的图像。50100150200250300350400501001502002
16、503003504004505005504 .激发介质(BZ reaction or heart )规则：? 元胞有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;
17、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)<0)&(cel
18、ls(x+1,y-1)<t)+(cells(x-1,y+1)>t)+(cells(x+1,y+1)>0)&(cells(x+1,y+1)< t);cells = (cells=O) & (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个不同的状态。
19、状态为0是空位，状态=1是燃烧着的树木，状态 =2是树木。? 如果4个邻居中有一个或一个以上的是燃烧着的并且自身是树木（状态为2），那么该元胞下一时刻的状态是燃烧（状态为1）。? 森林元胞（状态为2 ）以一个低概率（例如 0.000005 ）开始烧（因为闪 电）。? 一个燃烧着的元胞（状态为 1 ）在下一时时刻变成空位的（状态为 0 ）。? 空元胞以一个低概率（例如0.01 ）变为森林以模拟生长。? 出于矩阵边界连接的考虑，如果左边界开始着火，火势将向右蔓延，右边界同理。同样适用于顶部和底部核心代码：sum = (veg(1: n,n 1:n-1)=1) + (veg(1: n,2:n 1)=
20、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);注意环形连接是由序标实现的6 .气体动力学 这个CA （以及接下来的两个CA）是用来模拟粒子运动的。此元胞自动机需要 一种不同类型的元胞的邻居。此元胞的邻居时刻变化，因此某一个方向运动趋势， 将继续在同一个方向。换言之，此规则保存势头，这是基础的动力仿真。这种邻 居通
21、常被称为margolis邻居并且这种邻居通常由重叠的 2x2块的元胞构成。在 下面的表格中，偶数步长时左上方 4兀胞为邻居关系，奇数步长时右下的 4兀 胞为邻居关系。某一特定元胞在每一个时间步长都有 3个邻居，但是具体的元胞 构成了邻居的旋转和反复。偶偶偶元胞奇奇奇规则: ? 此规则叫作HPP-气体规则。? 每个元胞有2种状态。状态=0是空的，状态=1代表粒子。? 在任何一个时间步长，假设粒子是刚刚进入2x2的网格块。它将通过其网格 块的中心到达对角的网格中，所以在任何时间步长，每一个元胞与该元胞对 角对元胞交换的内容。如下所示，左边显示出来的元胞结构经过一个时间步 长变为右边的结构。以下是六
22、种不同的情况，所有所有的元胞都遵循相同的 转动规则。下文还将考虑两种特殊情况，即粒子 -粒子碰撞和粒子-墙碰撞。1101101111111111? 为了实现粒子碰撞过程（保证动量和能量守恒），对于两个处于对角线上的 粒子，他们相互撞击后偏转90度。在一个时间步长里使其从一个对角转成 另一个对角。你可以逆时针旋转这四个元胞来实现这个过程。则第三规则可 以表示为：10010110? 粒子撞击墙壁时，简单地使其离开且状态不变。这就引起反射现象。核心代码:p=mod(i,2); %margolis n eighborhood, where i is the time step %upper left
23、cell updatexi nd = 1+p:2: nx-2+p;yi nd = 1+p:2: ny-2+p;%See if exactly one diago nal is ones%only (at most) one of the following can be true!diag1(x ind,yind) = (sa nd(x ind,yin d)=1) & (sa nd(x in d+1, yin d+1)=1) & .(sa nd(x in d+1, yin d)=0) & (sa nd(x ind,yin d+1)=0);diag2(x ind,yind
24、) = (sa nd(x in d+1, yin d)=1) & (sa nd(x ind,yin d+1)=1) & .(sa nd(x ind,yin d)=0) & (sa nd(x in d+1, yin d+1)=0);%The diago nals both not occupied by two particlesand12(xind,yind) = (diag1(xind,yind)=0) & (diag2(xind,yind)=0);%One diago nal is occupied by two particlesor12(x ind,yi
25、nd)= diag1(x ind,yind) | diag2(x ind,yin d);%for every gas particle see if it n ear the boun dary sums(x ind,yind) = gn d(x ind,yind) | gn d(x in d+1, yind) | .gn d(x ind,yin d+1) | gn d(x in d+1, yin d+1);% cell layout:% x,y x+1,y% x,y+1x+1,y+1%lf (no walls) and (diag on als are both not occupied)%
26、the n there is no collisi on, so move opposite cell to curre nt cell%If (no walls) and (only one diago nal is occupied)%the n there is a collisi on so move ccw cell to the curre nt cell%If (a wall)%the n don't cha nge the cell (causes a reflecti on)san dNew(x ind,yind)=.(an d12(x ind,yind)&
27、sums(x ind,yind) & san d(x in d+1, yin d+1) + .(or12(x ind,yind) & sums(x ind,yind) & san d(x ind,yin d+1) + .(sums(x ind,yind) & san d(x ind,yin d);san dNew(x in d+1, yind)=.(an d12(x ind,yind)& sums(x ind,yind) & san d(x ind,yin d+1) + .(or12(x ind,yind) & sums(x ind,yi
28、nd) & san d(x ind,yin d)+ .(sums(x ind,yind) & san d(x in d+1, yin d);san dNew(x ind,yin d+1)=.(an d12(x ind,yind)& sums(x ind,yind) & san d(x in d+1, yin d) + .(or12(x ind,yind) & sums(x ind,yind) & san d(x in d+1, yin d+1)+ .(sums(x ind,yind) & san d(x ind,yin d+1);san
29、dNew(x in d+1, yin d+1)=.(an d12(x ind,yind)& sums(x ind,yind) & san d(x ind,yin d) + .(or12(x ind,yind) & sums(x ind,yind) & san d(x in d+1, yin d)+ . (sums(x ind,yind) & san d(x in d+1, yin d+1);sand = san dNew;8.扩散限制聚集这个系统是模拟粘性颗粒的聚集，形成分形结构。质点以一个类似于例6中的HPP-气体规则发生运动。不同的是粒子在一些高密度
30、（但看不见）的液体周围 被假定是弹跳的。效果是每一个粒子在每个时间步长在随机的方向上运动。换言 之，每一个时间步长是一个碰撞的过程。 这个仿真矩阵的中心确定了在一个固定 生长颗粒。任何弥散粒子触及它就会被它粘住， 并成为一个不能移动的，有粘性 颗粒。规则：? 使用Margolus型邻居。在每一个时间步，等概率地顺时针或逆时针旋转4个元胞。旋转使速度随机化。? 在移动后，如果八个最近的邻居有一个或一个以上元胞是固定的粘性颗粒，则下时刻该元胞将被冻结，并且使之有粘性。核心代码：p=mod(i,2); %margolis n eighborhood%upper left cell update xi
31、 nd = 1+p:2: nx-2+p;yind = 1+p:2: ny-2+p;%ran dom velocity choice vary = rand(nx,ny)< .5 ;vary1 = 1-vary;%diffusi on rule - margolus n eighborhood%rotate the 4 cells to ran domize velocitysan dNew(x ind,yind)=.vary(x ind,yin d).*sa nd(x in d+1, yind) + . %cwvary1(x ind,yin d).*sa nd(x ind,yin d+1
32、) ;%ccwsan dNew(x in d+1, yind)=.vary(x ind,yin d).*sa nd(x in d+1, yin d+1) + .vary1(x ind,yin d).*sa nd(x ind,yind);san dNew(x ind,yin d+1)=.vary(x ind,yin d).*sa nd(x ind,yind) + .vary1(x ind,yin d).*sa nd(x in d+1, yin d+1);san dNew(x in d+1, yin d+1)=.vary(x ind,yin d).*sa nd(x ind,yin d+1) + .
33、vary1(x ind,yin d).*sa nd(x in d+1, yind);sand = san dNew;%for every sand gra in see if it n ear the fixed, sticky clustersum(2: nx-1,2: ny-1) = gn d(2: nx-1,1: ny-2) + gn d(2: nx-1,3: ny) + .gnd(1:nx-2, 2:ny-1) + gnd(3:nx,2:ny-1) + .gn d(1: nx-2,1: ny-2) + gn d(1: nx-2,3: ny) + .gn d(3: nx,1: ny-2)
34、 + gn d(3:n x,3: ny);%add to the clustergnd = (sum> 0) & (sa nd=1) | gnd ;%and elim in ate the moving particlesan d(fi nd(g nd=1) = 0;以下经过很多时间步长后固定集聚后的图像显示。9 .砂堆规则一堆沙子的横截面，可以使用 Margolus型邻居仿真，但运动规则不同。 规则：? 元胞有2个状态。状态=0是空的，状态=1代表沙子。? 在任何时间步长，一个粒子，可以在 2x2块中向着底部运动。可能运动如下 所示。墙壁和底部将阻止粒子继续运动。? 为了让该运
35、动略有随机性，我亦补充说一项规则，有时处于下落状态的两个 元胞还旋转，直到所有的动作都完成。00000000核心代码：p=mod(i,2); %margolis n eighborhoodsand(nx/2,ny/2) = 1; %add a grain at the top%upper left cell updatexi nd =1+p:2: nx-2+p;yi nd =1+p:2: ny-2+p;%ran domize the flow - 10% of the timevary = rand(nx,ny)< .9 ;vary1 = 1-vary;san dNew(x ind,yi
36、nd)=.gn d(x ind,yin d).*sa nd(x ind,yind) + .(1-g nd(x ind,yin d).*sa nd(x ind,yin d).*sa nd(x ind,yin d+1) .* .(san d(x in d+1, yin d+1)+(1-sa nd(x in d+1, yin d+1).*sa nd(x in d+1, yin d);san dNew(x in d+1, yind)=.gn d(x in d+1, yin d).*sa nd(x in d+1, yind) + .(1-g nd(x in d+1,yi nd).*sa nd(xi nd
37、+1,yi nd).*sa nd(x in d+1,yi nd+1) .* . (sa nd(x ind,yin d+1)+(1-sa nd(x ind,yin d+1).*sa nd(x ind,yin d);san dNew(x ind,yin d+1)=.san d(x ind,yin d+1) + .(1-sa nd(xi nd,yi nd+1) .* .(san d(x ind,yin d).*(1-g nd(x ind,yin d) + .(1-sa nd(x in d,yi nd).*sa nd(xi nd+1,yi nd).*(1-g nd(xi nd+1,yi nd).*sa
38、 nd(x in d+1,yi nd+1) )；san dNew(x in d+1, yin d+1)=.san d(x in d+1, yin d+1) + .(1-sand(xind+1,yind+1) .* .(san d(x in d+1, yin d).*(1-g nd(x in d+1, yin d) + .(1-sa nd(x in d+1, yin d).*sa nd(x ind,yin d).*(1-g nd(x ind,yin d).*sa nd(x ind,yin d+1);%scramble the sites to make it look bettertemp1 =
39、 san dNew(x ind,yin d+1).*vary(x ind,yin d+1) + .san dNew(x in d+1, yin d+1).*vary1(x ind,yin d+1);temp2 = san dNew(x in d+1, yin d+1).*vary(x ind,yin d+1) + .san dNew(x ind,yin d+1).*vary1(x ind,yin d+1);san dNew(x ind,yin d+1) = temp1;san dNew(x in d+1, yin d+1) = temp2;sand = san dNew;参考文献1 Cellular Automata Modeli ng of Physical SystemsNorma n2 Cellular Automata Machines by Tommaso Toffoli and Margolus, MIT Press, 1987.3 Cellular Automata Modeli ng of Physical Systems by Bastie n Chopardand Michel Droz, Cambridge University Press, 1998.