元胞自動機與Matlab

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

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

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

4、是很容易的。在下面這個例子中，應(yīng)用了三個 按鈕和一個文本框。三個按鈕，作用分別是運行，停止，程序退出按鈕。文 框是用來顯示的仿真運算的次數(shù)。%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, run=1;);%defi ne the stop butt onerasebutt on=uic on trol(style,pushbutt on,.stri n

5、g,Stop, .fon tsize,12, .positio n,200,400,50,20, .callback,freeze=1;);%defi ne the Quit butt onquitbutt on=uic on trol(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, .positio n,20,400,50,

6、20）;經(jīng)過對控件（和CA）初始化，程序進(jìn)入一個循環(huán)，該循環(huán)測試由回調(diào)函數(shù)的每 個按鈕控制的變量。剛開始運行時，只在嵌套的 while循環(huán)和if語句中運行。 直到退出按鈕按下時，循環(huán)停止。另外兩個按鈕按下時執(zhí)行相應(yīng)的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) + cells(x,y+1) + .cells

7、(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 mber,stri ng); set (nu mber,stri n

8、g, 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的生命游戲機。規(guī)則是:? 對周圍的8個近鄰的元胞狀態(tài)求和? 如果總和為2的話，則下一時刻的狀態(tài)不改變? 如果總和為3，則下一時刻的狀態(tài)為1? 否則狀態(tài)=0核心代碼:x = 2:n-1;y = 2:n-1;%n earest n eighbor sumsum(x,y) = cells(x,y-1) + cells(x,y+1) + . ce

9、lls(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ī)則是：? 對周圍的8近鄰的元胞以及自身的狀態(tài)求和? 如果總和4或=5，下一時刻的狀態(tài)=0? 否則狀態(tài)=1核心代碼:x = 2:n-1;y = 2:n-1;%n earest n eighbor sumsum(x,y) = cells(x,y-1) + cells(x,y+1) + . c

10、ells(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=1) & (pick=threshold) & (visit=O)visit = (sum=1);變量a和b是圖像的尺寸。最初的圖形是由圖形操作決定的。以下程序設(shè)定坐 標(biāo)系為一個固定的尺寸，在坐標(biāo)系里寫入文本，然后獲得并返回坐標(biāo)內(nèi)的內(nèi)容， 并用getframe函數(shù)把它們寫入一個矩陣ax = axes(u nits,p

11、ixels,positio n,1 1 500 400,color,k);text(u nits, 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,color,w,fo ntn ame,helvetica,fo ntsize,50)in itial = getframe(gca);a,b,c=size(i nitial.cdata);z=

12、zeros(a,b);cells = double(i nitial.cdata(:,:,1)=255);visit = z ;sum = z;經(jīng)過幾十個時間間隔(從 MCM Cellular Automata這個圖像開始)，我們可以得到以下的圖像。50100150200250300350400501001502002503003504004505005504 .激發(fā)介質(zhì)(BZ reaction or heart )規(guī)則：? 元胞有10個不同的狀態(tài)。狀態(tài)0是體眠。1-5為活躍狀態(tài)，、6-9為是極活 躍狀態(tài)。? 計算每一個處于活躍的狀態(tài)的元胞近鄰的 8個元胞。? 如果和大于或等于3 （至少有三個

13、活躍的鄰居），則下一時刻該元胞=1 0)&(cells(x,y-1)文案大全0)&(cells(x,y+1) 0)&(cells(x-1,y)0)&(cells(x+1,y)0)&(cells(x-1,y+1)0)&(cells(x-1,y-1)0)&(cells(x+1,y-1)t)+(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) +

14、.0*(cells=9);一個CA初始圖形經(jīng)過螺旋的變化，得到下圖5 .森林火災(zāi)規(guī)則：? 元胞有3個不同的狀態(tài)。狀態(tài)為0是空位，狀態(tài)=1是燃燒著的樹木，狀態(tài) =2是樹木。? 如果4個鄰居中有一個或一個以上的是燃燒著的并且自身是樹木（狀態(tài)為2），那么該元胞下一時刻的狀態(tài)是燃燒（狀態(tài)為1）。? 森林元胞（狀態(tài)為2 ）以一個低概率（例如 0.000005 ）開始燒（因為閃 電）。? 一個燃燒著的元胞（狀態(tài)為 1 ）在下一時時刻變成空位的（狀態(tài)為 0 ）。? 空元胞以一個低概率（例如0.01 ）變?yōu)樯忠阅M生長。? 出于矩陣邊界連接的考慮，如果左邊界開始著火，火勢將向右蔓延，右邊界同理。同樣適用于頂

15、部和底部核心代碼：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) & rand(n,n ) Pgrowth);注意環(huán)形連接是由序標(biāo)實現(xiàn)的6 .氣體動力學(xué) 這個CA （以及接下來的兩個CA）是用來模擬粒子運動的。此元胞自動機需要 一種不同類型的元胞的鄰居。此元胞的鄰居時刻變化，因此某一個方向運動趨勢， 將繼

16、續(xù)在同一個方向。換言之，此規(guī)則保存勢頭，這是基礎(chǔ)的動力仿真。這種鄰 居通常被稱為margolis鄰居并且這種鄰居通常由重疊的 2x2塊的元胞構(gòu)成。在 下面的表格中，偶數(shù)步長時左上方 4兀胞為鄰居關(guān)系，奇數(shù)步長時右下的 4兀 胞為鄰居關(guān)系。某一特定元胞在每一個時間步長都有 3個鄰居，但是具體的元胞 構(gòu)成了鄰居的旋轉(zhuǎn)和反復(fù)。偶偶偶元胞奇奇奇規(guī)則: ? 此規(guī)則叫作HPP-氣體規(guī)則。? 每個元胞有2種狀態(tài)。狀態(tài)=0是空的，狀態(tài)=1代表粒子。? 在任何一個時間步長，假設(shè)粒子是剛剛進(jìn)入2x2的網(wǎng)格塊。它將通過其網(wǎng)格 塊的中心到達(dá)對角的網(wǎng)格中，所以在任何時間步長，每一個元胞與該元胞對 角對元胞交換的內(nèi)容。如

17、下所示，左邊顯示出來的元胞結(jié)構(gòu)經(jīng)過一個時間步 長變?yōu)橛疫叺慕Y(jié)構(gòu)。以下是六種不同的情況，所有所有的元胞都遵循相同的 轉(zhuǎn)動規(guī)則。下文還將考慮兩種特殊情況，即粒子 -粒子碰撞和粒子-墻碰撞。1101101111111111? 為了實現(xiàn)粒子碰撞過程（保證動量和能量守恒），對于兩個處于對角線上的 粒子，他們相互撞擊后偏轉(zhuǎn)90度。在一個時間步長里使其從一個對角轉(zhuǎn)成 另一個對角。你可以逆時針旋轉(zhuǎn)這四個元胞來實現(xiàn)這個過程。則第三規(guī)則可 以表示為：10010110? 粒子撞擊墻壁時，簡單地使其離開且狀態(tài)不變。這就引起反射現(xiàn)象。核心代碼:p=mod(i,2); %margolis n eighborhood, w

18、here i is the time step %upper left 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

19、d+1)=0);diag2(x ind,yind) = (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(

20、x ind,yind)= 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 oc

21、cupied)%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 dont cha nge the cell (causes a reflecti on)san dNew(x ind,yind)=.(an d12(x ind,yind)& s

22、ums(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,yind) & san d(x ind,yin d)+ .(s

23、ums(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 dNew(x in d+1, yin d+1)=.(an d12(x ind,yind)& sums(x ind,

24、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.擴散限制聚集這個系統(tǒng)是模擬粘性顆粒的聚集，形成分形結(jié)構(gòu)。質(zhì)點以一個類似于例6中的HPP-氣體規(guī)則發(fā)生運動。不同的是粒子在一些高密度（但看不見）的液體周圍 被假定是彈跳的。效果是每一個粒子在每個時間步長在隨機的方向上運動。換言 之，每一個時間步長是一個碰撞的過程。 這個仿真矩陣的中心確

25、定了在一個固定 生長顆粒。任何彌散粒子觸及它就會被它粘住， 并成為一個不能移動的，有粘性 顆粒。規(guī)則：? 使用Margolus型鄰居。在每一個時間步，等概率地順時針或逆時針旋轉(zhuǎn)4個元胞。旋轉(zhuǎn)使速度隨機化。? 在移動后，如果八個最近的鄰居有一個或一個以上元胞是固定的粘性顆粒，則下時刻該元胞將被凍結(jié)，并且使之有粘性。核心代碼：p=mod(i,2); %margolis n eighborhood%upper left cell update xi nd = 1+p:2: nx-2+p;yind = 1+p:2: ny-2+p;%ran dom velocity choice vary = rand

26、(nx,ny) 0) & (sa nd=1) | gnd ;%and elim in ate the moving particlesan d(fi nd(g nd=1) = 0;以下經(jīng)過很多時間步長后固定集聚后的圖像顯示。9 .砂堆規(guī)則一堆沙子的橫截面，可以使用 Margolus型鄰居仿真，但運動規(guī)則不同。 規(guī)則：? 元胞有2個狀態(tài)。狀態(tài)=0是空的，狀態(tài)=1代表沙子。? 在任何時間步長，一個粒子，可以在 2x2塊中向著底部運動?？赡苓\動如下 所示。墻壁和底部將阻止粒子繼續(xù)運動。? 為了讓該運動略有隨機性，我亦補充說一項規(guī)則，有時處于下落狀態(tài)的兩個 元胞還旋轉(zhuǎn)，直到所有的動作都完成。00000

27、000核心代碼：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,yind)=.gn d(x ind,yin d).*sa nd(x ind,yind) + .(1-g nd(x

28、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+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,y