元胞自動(dòng)機(jī)(CA)代碼及應(yīng)用(DOC)

1、元胞自動(dòng)機(jī)（CA）代碼及應(yīng)用引言元胞自動(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ī)。矩陣和圖像可以相互轉(zhuǎn)化

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

3、5*n,n/2) = 1;Matlab的代碼應(yīng)盡量簡(jiǎn)潔以減小運(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 = (sum=3) | (sum=2 & cells);加入一個(gè)簡(jiǎn)單的圖形用戶界面是

4、很容易的。在下面這個(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, .callback, run=1;);%define the stop button erasebutton=uicontrol(style,pushbutton,.string,Stop, .font

5、size,12, . position,200,400,50,20, .callback,freeze=1;);%define the Quit button 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);經(jīng)過對(duì)控件(和CA)初始化，程序進(jìn)入一個(gè)循環(huán)

6、，該循環(huán)測(cè)試由回調(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) + cells(x,y+1) + .cells(x-1, y) + cells(x+1,y) + . c

7、ells(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 numberdiaplay stepnumber = 1 + str2num(get(number,string);set(number,string,num2str(stepnumber) end if (freeze=

8、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)不改變?nèi)绻偤蜑?，則下一時(shí)刻的狀態(tài)為1否則狀態(tài)= 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) + ce

9、lls(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 = 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) + cell

10、s(x-1,y+1) + . cells(3:n,y-1) +cells(x+1,y+1)+.cells(x,y);% The CA rulecells = (sum 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) + .

11、 cells(3:a,1:b-2) + cells(3:a,3:b);pick = rand(a,b); cells = cells | (sum=1) & (pick=threshold) & (visit=O) visit = (sum=1);變量a和b是圖像的尺寸。最初的圖形是由圖形操作決定的。 以下程序設(shè)定坐標(biāo) 系為一個(gè)固定的尺寸，在坐標(biāo)系里寫入文本，然后獲得并返回坐標(biāo)內(nèi)的內(nèi)容，并 用getframe函數(shù)把它們寫入一個(gè)矩陣ax = axes(u nits,pixels,positio n,1 1 500 400,color,k); text(u nits, pixels

12、, 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,color,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); vi

13、sit = z ;sum = z;經(jīng)過幾十個(gè)時(shí)間間隔（從MCM Cellular Automata這個(gè)圖像開始），我們可以得 到以下的圖像。4 .激發(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 0)&(cells(x,y-1) 0)&(cells(x,y+1) 0)&(cells(x-1, y) 0)&(cells(x+1,y) 0)&(cells(x

14、-1,y-1) 0)&(cells(x-1,y+1) 0)&(cells(x+1,y-1) 0)&(cells(x+1,y+1)=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);一個(gè)CA初始圖形經(jīng)過螺旋的變化，得到下圖5 .森林火災(zāi)規(guī)則：元胞有3個(gè)不同的狀態(tài)。狀態(tài)為0是空位，狀態(tài)=1是燃燒著的樹木，狀態(tài)=2是樹木。如果4個(gè)鄰居中有一個(gè)

15、或一個(gè)以上的是燃燒著的并且自身是樹木(狀態(tài)為2)，那么該元胞下一時(shí)刻的狀態(tài)是燃燒(狀態(tài)為1)。森林元胞(狀態(tài)為2)以一個(gè)低概率(例如0.000005)開始燒(因?yàn)殚W電) 一個(gè)燃燒著的元胞(狀態(tài)為1)在下一時(shí)時(shí)刻變成空位的(狀態(tài)為0)空元胞以一個(gè)低概率(例如0.01)變?yōu)樯忠阅M生長(zhǎng)。出于矩陣邊界連接的考慮，如果左邊界開始著火，火勢(shì)將向右蔓延，右邊界 同理。同樣適用于頂部和底部。核心代碼：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 =.

16、2*(veg=2) - (veg=2) & (sum 0 | (ra nd( n,n) Plight nin g) + . 2*(veg=0) &rand(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)趨勢(shì)， 將繼續(xù)在同一個(gè)方向。換言之，此規(guī)則保存勢(shì)頭，這是基礎(chǔ)的動(dòng)力仿真。這種鄰 居通常被稱為margolis鄰居并且這種鄰居通常由重疊的2x2塊的元胞構(gòu)成。在下 面的表格中，偶數(shù)步長(zhǎng)時(shí)左上方4元胞為鄰居關(guān)系，奇數(shù)步長(zhǎng)時(shí)

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

18、對(duì)于兩個(gè)處于對(duì)角線上的 粒子，他們相互撞擊后偏轉(zhuǎn)90度。在一個(gè)時(shí)間步長(zhǎng)里使其從一個(gè)對(duì)角轉(zhuǎn)成 另一個(gè)對(duì)角。你可以逆時(shí)針旋轉(zhuǎn)這四個(gè)元胞來實(shí)現(xiàn)這個(gè)過程。則第三規(guī)則可 以表示為：10010110粒子撞擊墻壁時(shí)，簡(jiǎn)單地使其離開且狀態(tài)不變。這就引起反射現(xiàn)象。核心代碼：p=mod(i,2); %margolis n eighborhood, 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 m

19、ost) one of the following can be true!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 diago

20、nals both not occupied by two particles 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,y

21、ind) | .gnd(xind,yind+1) | gnd(xind+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 t

22、o the current cell%If (a wall)%then dont 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) = .

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

24、) & 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.擴(kuò)散限制聚集 這個(gè)系統(tǒng)

25、是模擬粘性顆粒的聚集， 形成分形結(jié)構(gòu)。 質(zhì)點(diǎn)以一個(gè)類似于例6中的HPP-氣體規(guī)則發(fā)生運(yùn)動(dòng)。不同的是粒子在一些高密度（但看不見）的液體周圍 被假定是彈跳的。 效果是每一個(gè)粒子在每個(gè)時(shí)間步長(zhǎng)在隨機(jī)的方向上運(yùn)動(dòng)。 換言 之，每一個(gè)時(shí)間步長(zhǎng)是一個(gè)碰撞的過程。 這個(gè)仿真矩陣的中心確定了在一個(gè)固定 生長(zhǎng)顆粒。 任何彌散粒子觸及它就會(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=m

26、od(i,2); %margolis neighborhood%upper left cell update xind = 1+p:2:nx-2+p; yind = 1+p:2:ny-2+p; %random velocitychoice vary = rand(nx,ny) 0) & (sand=1) | gnd ;%and eliminate the moving particle sand(find(gnd=1) = 0;以下經(jīng)過很多時(shí)間步長(zhǎng)后固定集聚后的圖像顯示。9 .砂堆規(guī)則一堆沙子的橫截面，可以使用Margolus型鄰居仿真，但運(yùn)動(dòng)規(guī)則不同。 規(guī)則：元胞有2個(gè)狀態(tài)。狀態(tài)=0

27、是空的，狀態(tài)=1代表沙子。在任何時(shí)間步長(zhǎng)，一個(gè)粒子，可以在2x2塊中向著底部運(yùn)動(dòng)?？赡苓\(yùn)動(dòng)如下 所示。墻壁和底部將阻止粒子繼續(xù)運(yùn)動(dòng)。為了讓該運(yùn)動(dòng)略有隨機(jī)性，我亦補(bǔ)充說一項(xiàng)規(guī)則，有時(shí)處于下落狀態(tài)的兩個(gè) 元胞還旋轉(zhuǎn)，直到所有的動(dòng)作都完成。00000000%upper left cell updatexi nd = 1+p:2: nx-2+p;yind = 1+p:2: ny-2+p;%ran domize the flow - 10% of the time vary = rand(nx,ny) .9 ;vary1 = 1-vary;san dNew(x ind,yind)=.gn d(x ind

28、,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) .* .(sa nd(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)=.gnd(xind+1,yind).*sand(xind+1,yind) + .(1-gnd(xind+1,yind).*sand(xind+1,yind).*sand(xind+1,yind+1) .* .sand(nx/2,ny/2) = 1; %add a grain at the top(sand(xind,yind+1)+(1-sand(xind,yind+1).*sand(xind,yind)