元胞自動(dòng)機(jī)與Matlab

1、元胞自動(dòng)機(jī)與 MATLAB引言元胞自動(dòng)機(jī)( CA )是一種用來仿真局部規(guī)則和局部聯(lián)系的方法。典型的元 胞自動(dòng)機(jī)是定義在網(wǎng)格上的，每一個(gè)點(diǎn)上的網(wǎng)格代表一個(gè)元胞與一種有限的狀 態(tài)。變化規(guī)則適用于每一個(gè)元胞并且同時(shí)進(jìn)行。 典型的變化規(guī)則， 決定于元胞的 狀態(tài)，以及其( 4 或 8 )鄰居的狀態(tài)。元胞自動(dòng)機(jī)已被應(yīng)用于物理模擬，生物 模擬等領(lǐng)域。本文就一些有趣的規(guī)則，考慮如何編寫有效的 MATLAB 的程序來 實(shí)現(xiàn)這些元胞自動(dòng)機(jī)。MATLAB 的編程考慮元胞自動(dòng)機(jī)需要考慮到下列因素， 下面分別說明如何用 MATLAB 實(shí)現(xiàn)這些部分。 并以 Conway 的生命游戲機(jī)的程序?yàn)槔?，說明怎樣實(shí)現(xiàn)一個(gè)元胞自動(dòng)機(jī)

2、。矩陣和圖像可以相互轉(zhuǎn)化，所以矩陣的顯示是可以真接實(shí)現(xiàn)的。如果矩陣 cells的所有元素只包含兩種狀態(tài)且矩陣 Z 含有零，那么用 image函數(shù)來顯示 cat命令建的 RGB 圖像，并且能夠返回句柄。imh = image(cat(3,cells,z,z);set(imh, 'erasemode', 'none')axis equalaxis tight矩陣和圖像可以相互轉(zhuǎn)化，所以初始條件可以是矩陣，也可以是圖形。以下 代碼生成一個(gè)零矩陣，初始化元胞狀態(tài)為零，然后使得中心十字形的元胞狀 態(tài)= 1。z = zeros(n,n); cells = z; cells(

3、n/2,.25*n:.75*n) = 1; cells(.25*n:.75*n,n/2) = 1;Matlab 的代碼應(yīng)盡量簡潔以減小運(yùn)算量。以下程序計(jì)算了最近鄰居總和，并 按照 CA 規(guī)則進(jìn)行了計(jì)算。本段 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 = (s

4、um=3) | (sum=2 & cells);加入一個(gè)簡單的圖形用戶界面是很容易的。在下面這個(gè)例子中，應(yīng)用了三個(gè) 按鈕和一個(gè)文本框。三個(gè)按鈕，作用分別是運(yùn)行，停止，程序退出按鈕。文 框是用來顯示的仿真運(yùn)算的次數(shù)。%build the GUI%define the plot button plotbutton=uicontrol('style','pushbutton',.'string','Run', .'fontsize',12, . 'position',100,400,50,20,

5、.'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=uicontro

6、l('style','pushbutton',.'string','Quit', .'fontsize',12, . 'position',300,400,50,20, .'callback','stop=1;close;');number = uicontrol('style','text', .'string','1', . 'fontsize',12, .'position&

7、#39;,20,400,50,20); 經(jīng)過對(duì)控件(和 CA)初始化，程序進(jìn)入一個(gè)循環(huán)，該循環(huán)測試由回調(diào)函數(shù)的每 個(gè)按鈕控制的變量。剛開始運(yùn)行時(shí)，只在嵌套的 while 循環(huán)和 if 語句中運(yùn)行。 直到退出按鈕按下時(shí)，循環(huán)停止。另外兩個(gè)按鈕按下時(shí)執(zhí)行相應(yīng)的 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 sum sum(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 image set(imh, 'cdata', cat(3,cells,z,z) ) %update the step number diaplay stepnumber = 1 + str2num(get(number,

9、'string'); set(number,'string',num2str(stepnumber) end if (freeze=1)run = 0;freeze = 0;enddrawnow %need this in the loop for controls to work end例子1 .Conway 的生命游戲機(jī)。 規(guī)則是：? 對(duì)周圍的 8 個(gè)近鄰的元胞狀態(tài)求和? 如果總和為 2 的話，則下一時(shí)刻的狀態(tài)不改變? 如果總和為 3 ，則下一時(shí)刻的狀態(tài)為 1? 否則狀態(tài) = 0 核心代碼：x = 2:n-1;y = 2:n-1;%nearest neigh

10、bor 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 rulecells = (sum=3) | (sum=2 & cells);2 . 表面張力 規(guī)則是： ? 對(duì)周圍的 8 近鄰的元胞以及自身的狀態(tài)求和 ? 如果總和 < 4 或 = 5 ，下一時(shí)刻的狀態(tài) = 0? 否則狀態(tài) = 1 核心代碼：x = 2:n-1;y

11、= 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 . 滲流集群規(guī)則：? 對(duì)周圍相鄰的 8 鄰居求和(元胞只有兩種狀態(tài)， 0 或 1 )。元胞也有一個(gè)單 獨(dú)的狀態(tài)參量(所謂 '記錄 &#

12、39; )記錄它們之前是否有非零狀態(tài)的鄰居。? 在 0 與 1 之間產(chǎn)生一個(gè)隨機(jī)數(shù) r 。? 如果總和 > 0 (至少一個(gè)鄰居)并且 r > 閾值，或者元胞從未有過一個(gè)鄰居， 則元胞 = 1 。? 如果總和 > 0 則設(shè)置"記錄 "的標(biāo)志，記錄這些元胞有一個(gè)非零的鄰居。 核心代碼：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)

13、+ .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是圖像的尺寸。 最初的圖形是由圖形操作決定的。 以下程序設(shè)定坐標(biāo) 系為一個(gè)固定的尺寸， 在坐標(biāo)系里寫入文本， 然后獲得并返回坐標(biāo)內(nèi)的內(nèi)容， 并 用 getframe 函數(shù)把它們寫入一個(gè)矩陣ax = axes('units','pixels','posit

14、ion',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'

15、;, 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;經(jīng)過幾十個(gè)時(shí)間間隔 （從 MCM Cellular Au

16、tomata 這個(gè)圖像開始） ，我們可以得 到以下的圖像。100 150 200 250 300 350 40050 100 150 200 250 300 350400 450 500 5504 . 激發(fā)介質(zhì)( BZ reaction or heart ) 規(guī)則：? 元胞有 10個(gè)不同的狀態(tài)。狀態(tài) 0是體眠。 1-5為活躍狀態(tài)，、6-9為是極活躍 狀態(tài)。? 計(jì)算每一個(gè)處于活躍的狀態(tài)的元胞近鄰的 8 個(gè)元胞。? 如果和大于或等于 3 （至少有三個(gè)活躍的鄰居） ，則下一時(shí)刻該元胞 = 1? 不需要其它輸入， 1 至 9 種狀態(tài)依次出現(xiàn)。 如果該時(shí)刻狀態(tài) = 1 那么下一時(shí)刻 狀態(tài) = 2 。如果

17、該時(shí)刻狀態(tài) = 2 ，然后下一時(shí)刻狀態(tài) = 3 ，對(duì)于其它的狀態(tài) 依次類推，直到第 9 種狀態(tài)。如果狀態(tài) = 9 ，然后下一狀態(tài) = 0 并且元胞回 到休息狀態(tài)。核心代碼：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

18、,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*(cell

19、s=3) + .5*(cells=4) + .6*(cells=5) +.7*(cells=6) +.8*(cells=7) +.9*(cells=8) +.0*(cells=9);一個(gè) CA 初始圖形經(jīng)過螺旋的變化，得到下圖。5 . 森林火災(zāi)規(guī)則：? 元胞有 3 個(gè)不同的狀態(tài)。狀態(tài)為 0 是空位，狀態(tài) = 1 是燃燒著的樹木，狀態(tài) = 2 是樹木。? 如果 4 個(gè)鄰居中有一個(gè)或一個(gè)以上的是燃燒著的并且自身是樹木（狀態(tài)為2 ） ，那么該元胞下一時(shí)刻的狀態(tài)是燃燒（狀態(tài)為 1 ） 。? 森林元胞（狀態(tài)為 2 ）以一個(gè)低概率 （例如 0.000005 ）開始燒（因?yàn)殚W電）? 一個(gè)燃燒著的元胞（狀態(tài)為

20、 1 ）在下一時(shí)時(shí)刻變成空位的（狀態(tài)為 0 ）? 空元胞以一個(gè)低概率（例如 0.01 ）變?yōu)樯忠阅M生長。? 出于矩陣邊界連接的考慮，如果左邊界開始著火，火勢將向右蔓延，右邊界 同理。同樣適用于頂部和底部。核心代碼：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) & ran

21、d(n,n)< Pgrowth) ;注意環(huán)形連接是由序標(biāo)實(shí)現(xiàn)的。6 . 氣體動(dòng)力學(xué)這個(gè) CA （以及接下來的兩個(gè) CA ）是用來模擬粒子運(yùn)動(dòng)的。此元胞自動(dòng)機(jī)需要 一種不同類型的元胞的鄰居。 此元胞的鄰居時(shí)刻變化， 因此某一個(gè)方向運(yùn)動(dòng)趨勢， 將繼續(xù)在同一個(gè)方向。換言之，此規(guī)則保存勢頭，這是基礎(chǔ)的動(dòng)力仿真。這種鄰 居通常被稱為 margolis 鄰居并且這種鄰居通常由重疊的 2x2 塊的元胞構(gòu)成。 在下 面的表格中，偶數(shù)步長時(shí)左上方 4元胞為鄰居關(guān)系， 奇數(shù)步長時(shí)右下的 4元胞為 鄰居關(guān)系。某一特定元胞在每一個(gè)時(shí)間步長都有 3 個(gè)鄰居，但是具體的元胞構(gòu)成 了鄰居的旋轉(zhuǎn)和反復(fù)。偶偶偶元 胞奇奇

22、奇規(guī)則：? 此規(guī)則叫作 HPP-氣體規(guī)則。? 每個(gè)元胞有 2 種狀態(tài)。狀態(tài) = 0 是空的，狀態(tài) = 1 代表粒子。? 在任何一個(gè)時(shí)間步長，假設(shè)粒子是剛剛進(jìn)入 2x2 的網(wǎng)格塊。它將通過其網(wǎng)格 塊的中心到達(dá)對(duì)角的網(wǎng)格中，所以在任何時(shí)間步長，每一個(gè)元胞與該元胞對(duì) 角對(duì)元胞交換的內(nèi)容。如下所示，左邊顯示出來的元胞結(jié)構(gòu)經(jīng)過一個(gè)時(shí)間步 長變?yōu)橛疫叺慕Y(jié)構(gòu)。以下是六種不同的情況，所有所有的元胞都遵循相同的 轉(zhuǎn)動(dòng)規(guī)則。下文還將考慮兩種特殊情況，即粒子 -粒子碰撞和粒子 -墻碰撞。? 為了實(shí)現(xiàn)粒子碰撞過程（保證動(dòng)量和能量守恒） ，對(duì)于兩個(gè)處于對(duì)角線上的 粒子，他們相互撞擊后偏轉(zhuǎn) 90 度。在一個(gè)時(shí)間步長里使其

23、從一個(gè)對(duì)角轉(zhuǎn)成 另一個(gè)對(duì)角。你可以逆時(shí)針旋轉(zhuǎn)這四個(gè)元胞來實(shí)現(xiàn)這個(gè)過程。則第三規(guī)則可 以表示為：10010110? 粒子撞擊墻壁時(shí)，簡單地使其離開且狀態(tài)不變。這就引起反射現(xiàn)象。核心代碼：p=mod(i,2); %margolis neighborhood, where i is the time step%upper left cell updatexind = 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!

24、 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 particl

25、es 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(xin

26、d+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)%the

27、n 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

28、(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

29、)+ . (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.擴(kuò)散限制聚集 這個(gè)系統(tǒng)是模擬粘性顆粒的聚集，形成分形結(jié)構(gòu)。質(zhì)點(diǎn)以一個(gè)類似于

30、例 6 中的 HPP-氣體規(guī)則發(fā)生運(yùn)動(dòng) 。不同的是粒子在一些高密度（但看不見）的液體周圍 被假定是彈跳的。 效果是每一個(gè)粒子在每個(gè)時(shí)間步長在隨機(jī)的方向上運(yùn)動(dòng)。 換言 之，每一個(gè)時(shí)間步長是一個(gè)碰撞的過程。 這個(gè)仿真矩陣的中心確定了在一個(gè)固定 生長顆粒。 任何彌散粒子觸及它就會(huì)被它粘住， 并成為一個(gè)不能移動(dòng)的， 有粘性 顆粒。規(guī)則：? 使用 Margolus 型鄰居。在每一個(gè)時(shí)間步，等概率地順時(shí)針或逆時(shí)針旋轉(zhuǎn) 4 個(gè)元胞。旋轉(zhuǎn)使速度隨機(jī)化。? 在移動(dòng)后，如果八個(gè)最近的鄰居有一個(gè)或一個(gè)以上元胞是固定的粘性顆粒， 則下時(shí)刻該元胞將被凍結(jié)，并且使之有粘性。核心代碼：p=mod(i,2); %margol

31、is 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(

32、xind,yind).*sand(xind,yind+1) ; %ccwsandNew(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)

33、+ . 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);

34、%add to the cluster gnd = (sum> 0) & (sand=1) | gnd ;%and eliminate the moving particle sand(find(gnd=1) = 0;以下經(jīng)過很多時(shí)間步長后固定集聚后的圖像顯示。9 . 砂堆規(guī)則 一堆沙子的橫截面，可以使用 Margolus 型鄰居仿真，但運(yùn)動(dòng)規(guī)則不同。 規(guī)則：? 元胞有 2 個(gè)狀態(tài)。狀態(tài) = 0 是空的，狀態(tài) = 1 代表沙子。? 在任何時(shí)間步長，一個(gè)粒子，可以在 2x2 塊中向著底部運(yùn)動(dòng)?？赡苓\(yùn)動(dòng)如下 所示。墻壁和底部將阻止粒子繼續(xù)運(yùn)動(dòng)。? 為了讓該運(yùn)動(dòng)略有隨機(jī)性，我亦補(bǔ)充說一

35、項(xiàng)規(guī)則，有時(shí)處于下落狀態(tài)的兩個(gè) 元胞還旋轉(zhuǎn)，直到所有的動(dòng)作都完成。00000000核心代碼：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;%randomize the flow - 10% of the time vary = rand(nx,ny)< .9 ; 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-sa