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

-.z元胞自動(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)化，所以矩陣的顯示是可以真接實(shí)現(xiàn)的。如果矩陣cells的所有元素只包含兩種狀態(tài)且矩陣Z含有零，則用image函數(shù)來顯示cat命令建的RGB圖像，并且能夠返回句柄。imh=image(cat(3,cells,z,z));set(imh,'erasemode','none')a*isequala*istight矩陣和圖像可以相互轉(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:.75*n,n/2)=1; Matlab的代碼應(yīng)盡量簡(jiǎn)潔以減小運(yùn)算量。以下程序計(jì)算了最近鄰居總和，并按照CA規(guī)則進(jìn)展了計(jì)算。本段Matlab代碼非常靈活的表示了相鄰鄰居。*=2:n-1;y=2:n-1;sum(*,y)=cells(*,y-1)+cells(*,y+1)+...cells(*-1,y)+cells(*+1,y)+...cells(*-1,y-1)+cells(*-1,y+1)+...cells(*+1,y-1)+cells(*+1,y+1);cells=(sum==3)|(sum==2&cells);參加一個(gè)簡(jiǎn)單的圖形用戶界面是很容易的。在下面這個(gè)例子中，應(yīng)用了三個(gè)按鈕和一個(gè)文本框。三個(gè)按鈕，作用分別是運(yùn)行，停頓，程序退出按鈕。文框是用來顯示的仿真運(yùn)算的次數(shù)。%buildtheGUI%definetheplotbuttonplotbutton=uicontrol('style','pushbutton',...'string','Run',...'fontsize',12,...'position',[100,400,50,20],...'callback','run=1;');%definethestopbuttonerasebutton=uicontrol('style','pushbutton',...'string','Stop',...'fontsize',12,...'position',[200,400,50,20],...'callback','freeze=1;');%definetheQuitbuttonquitbutton=uicontrol('style','pushbutton',...'string','Quit',...'fontsize',12,...'position',[300,400,50,20],...'callback','stop=1;close;');number=uicontrol('style','te*t',...'string','1',...'fontsize',12,...'position',[20,400,50,20]);經(jīng)過對(duì)控件〔和CA〕初始化，程序進(jìn)入一個(gè)循環(huán)，該循環(huán)測(cè)試由回調(diào)函數(shù)的每個(gè)按鈕控制的變量。剛開場(chǎng)運(yùn)行時(shí)，只在嵌套的while循環(huán)和if語句中運(yùn)行。直到退出按鈕按下時(shí)，循環(huán)停頓。另外兩個(gè)按鈕按下時(shí)執(zhí)行相應(yīng)的if語句。stop=0;%waitforaquitbuttonpushrun=0;%waitforadrawfreeze=0;%waitforafreezewhile(stop==0)if(run==1) %nearestneighborsum sum(*,y)=cells(*,y-1)+cells(*,y+1)+... cells(*-1,y)+cells(*+1,y)+... cells(*-1,y-1)+cells(*-1,y+1)+...cells(3:n,y-1)+cells(*+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))endif(freeze==1) run=0; freeze=0;enddrawnow%needthisintheloopforcontrolstoworkend例子1.Conway的生命游戲機(jī)。規(guī)則是：對(duì)周圍的8個(gè)近鄰的元胞狀態(tài)求和如果總和為2的話，則下一時(shí)刻的狀態(tài)不改變?nèi)绻偤蜑?，則下一時(shí)刻的狀態(tài)為1否則狀態(tài)=0核心代碼：*=2:n-1;y=2:n-1;%nearestneighborsumsum(*,y)=cells(*,y-1)+cells(*,y+1)+... cells(*-1,y)+cells(*+1,y)+... cells(*-1,y-1)+cells(*-1,y+1)+... cells(3:n,y-1)+cells(*+1,y+1);%TheCArulecells=(sum==3)|(sum==2&cells);2.外表力規(guī)則是：對(duì)周圍的8近鄰的元胞以及自身的狀態(tài)求和如果總和<4或=5，下一時(shí)刻的狀態(tài)=0否則狀態(tài)=1核心代碼：*=2:n-1;y=2:n-1;%nearestneighborsumsum(*,y)=cells(*,y-1)+cells(*,y+1)+...cells(*-1,y)+cells(*+1,y)+...cells(*-1,y-1)+cells(*-1,y+1)+...cells(3:n,y-1)+cells(*+1,y+1)+...cells(*,y);%TheCArulecells=~((sum<4)|(sum==5));3.滲流集群規(guī)則：對(duì)周圍相鄰的8鄰居求和〔元胞只有兩種狀態(tài)，0或1〕。元胞也有一個(gè)單獨(dú)的狀態(tài)參量〔所謂'記錄'〕記錄它們之前是否有非零狀態(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)+...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)的容，并用getframe函數(shù)把它們寫入一個(gè)矩陣a*=a*es('units','pi*els','position',[11500400],'color','k');te*t('units','pi*els','position',[130,255,0],...'string','MCM','color','w','fontname','helvetica','fontsize',100)te*t('units','pi*els','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;經(jīng)過幾十個(gè)時(shí)間間隔〔從MCMCellularAutomata這個(gè)圖像開場(chǎng)〕，我們可以得到以下的圖像。4.激發(fā)介質(zhì)〔BZreactionorheart〕規(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。如果該時(shí)刻狀態(tài)=2，然后下一時(shí)刻狀態(tài)=3，對(duì)于其它的狀態(tài)依次類推，直到第9種狀態(tài)。如果狀態(tài)=9，然后下一狀態(tài)=0并且元胞回到休息狀態(tài)。核心代碼：*=[2:n-1];y=[2:n-1];sum(*,y)=((cells(*,y-1)>0)&(cells(*,y-1)<t))+((cells(*,y+1)>0)&(cells(*,y+1)<t))+...((cells(*-1,y)>0)&(cells(*-1,y)<t))+((cells(*+1,y)>0)&(cells(*+1,y)<t))+...((cells(*-1,y-1)>0)&(cells(*-1,y-1)<t))+((cells(*-1,y+1)>0)&(cells(*-1,y+1)<t))+...((cells(*+1,y-1)>0)&(cells(*+1,y-1)<t))+((cells(*+1,y+1)>0)&(cells(*+1,y+1)<t));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);一個(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〕開場(chǎng)燒〔因?yàn)殚W電〕。一個(gè)燃燒著的元胞〔狀態(tài)為1〕在下一時(shí)時(shí)刻變成空位的〔狀態(tài)為0〕。空元胞以一個(gè)低概率〔例如0.01〕變?yōu)樯忠阅M生長。出于矩陣邊界連接的考慮，如果左邊界開場(chǎng)著火，火勢(shì)將向右蔓延，右邊界同理。同樣適用于頂部和底部。核心代碼： 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)<Plightning)))+...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ì)頭，這是根底的動(dòng)力仿真。這種鄰居通常被稱為margolis鄰居并且這種鄰居通常由重疊的2*2塊的元胞構(gòu)成。在下面的表格中，偶數(shù)步長時(shí)左上方4元胞為鄰居關(guān)系，奇數(shù)步長時(shí)右下的4元胞為鄰居關(guān)系。*一特定元胞在每一個(gè)時(shí)間步長都有3個(gè)鄰居，但是具體的元胞構(gòu)成了鄰居的旋轉(zhuǎn)和反復(fù)。偶偶偶元胞奇奇奇規(guī)則：此規(guī)則叫作HPP-氣體規(guī)則。每個(gè)元胞有2種狀態(tài)。狀態(tài)=0是空的，狀態(tài)=1代表粒子。在任何一個(gè)時(shí)間步長，假設(shè)粒子是剛剛進(jìn)入2*2的網(wǎng)格塊。它將通過其網(wǎng)格塊的中心到達(dá)對(duì)角的網(wǎng)格中，所以在任何時(shí)間步長，每一個(gè)元胞與該元胞對(duì)角對(duì)元胞交換的容。如下所示，左邊顯示出來的元胞構(gòu)造經(jīng)過一個(gè)時(shí)間步長變?yōu)橛疫叺臉?gòu)造。以下是六種不同的情況，所有所有的元胞都遵循一樣的轉(zhuǎn)動(dòng)規(guī)則。下文還將考慮兩種特殊情況，即粒子-粒子碰撞和粒子-墻碰撞。000010000000000110101001010110011101111110111111為了實(shí)現(xiàn)粒子碰撞過程〔保證動(dòng)量和能量守恒〕，對(duì)于兩個(gè)處于對(duì)角線上的粒子，他們相互撞擊后偏轉(zhuǎn)90度。在一個(gè)時(shí)間步長里使其從一個(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);%margolisneighborhood,whereiisthetimestep%upperleftcellupdate*ind=[1+p:2:n*-2+p];yind=[1+p:2:ny-2+p];%Seeife*actlyonediagonalisones%only(atmost)oneofthefollowingcanbetrue!diag1(*ind,yind)=(sand(*ind,yind)==1)&(sand(*ind+1,yind+1)==1)&...(sand(*ind+1,yind)==0)&(sand(*ind,yind+1)==0);diag2(*ind,yind)=(sand(*ind+1,yind)==1)&(sand(*ind,yind+1)==1)&...(sand(*ind,yind)==0)&(sand(*ind+1,yind+1)==0);%Thediagonalsbothnotoccupiedbytwoparticlesand12(*ind,yind)=(diag1(*ind,yind)==0)&(diag2(*ind,yind)==0);%Onediagonalisoccupiedbytwoparticlesor12(*ind,yind)=diag1(*ind,yind)|diag2(*ind,yind);%foreverygasparticleseeifitneartheboundarysums(*ind,yind)=gnd(*ind,yind)|gnd(*ind+1,yind)|...gnd(*ind,yind+1)|gnd(*ind+1,yind+1);%celllayout:%*,y*+1,y%*,y+1*+1,y+1%If(nowalls)and(diagonalsarebothnotoccupied)%thenthereisnocollision,somoveoppositecelltocurrentcell%If(nowalls)and(onlyonediagonalisoccupied)%thenthereisacollisionsomoveccwcelltothecurrentcell%If(awall)%thendon'tchangethecell(causesareflection)sandNew(*ind,yind)=...(and12(*ind,yind)&~sums(*ind,yind)&sand(*ind+1,yind+1))+...(or12(*ind,yind)&~sums(*ind,yind)&sand(*ind,yind+1))+...(sums(*ind,yind)&sand(*ind,yind));sandNew(*ind+1,yind)=...(and12(*ind,yind)&~sums(*ind,yind)&sand(*ind,yind+1))+...(or12(*ind,yind)&~sums(*ind,yind)&sand(*ind,yind))+...(sums(*ind,yind)&sand(*ind+1,yind));sandNew(*ind,yind+1)=...(and12(*ind,yind)&~sums(*ind,yind)&sand(*ind+1,yind))+...(or12(*ind,yind)&~sums(*ind,yind)&sand(*ind+1,yind+1))+...(sums(*ind,yind)&sand(*ind,yind+1));sandNew(*ind+1,yind+1)=...(and12(*ind,yind)&~sums(*ind,yind)&sand(*ind,yind))+...(or12(*ind,yind)&~sums(*ind,yind)&sand(*ind+1,yind))+...(sums(*ind,yind)&sand(*ind+1,yind+1));sand=sandNew;8.擴(kuò)散限制聚集這個(gè)系統(tǒng)是模擬粘性顆粒的聚集，形成分形構(gòu)造。質(zhì)點(diǎn)以一個(gè)類似于例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);%margolisneighborhood%upperleftcellupdate*ind=[1+p:2:n*-2+p];yind=[1+p:2:ny-2+p];%randomvelocitychoicevary=rand(n*,ny)<.5;vary1=1-vary;%diffusionrule--margolusneighborhood%rotatethe4cellstorandomizevelocitysandNew(*ind,yind)=...vary(*ind,yind).*sand(*ind+1,yind)+...%cwvary1(*ind,yind).*sand(*ind,yind+1);%ccwsandNew(*ind+1,yind)=...vary(*ind,yind).*sand(*ind+1,yind+1)+...vary1(*ind,yind).*sand(*ind,yind);sandNew(*ind,yind+1)=...vary(*ind,yind).*sand(*ind,yind)+...vary1(*ind,yind).*sand(*ind+1,yind+1);sandNew(*ind+1,yind+1)=...vary(*ind,yind).*sand(*ind,yind+1)+...vary1(*ind,yind).*sand(*ind+1,yind);sand=sandNew;%foreverysandgrainseeifitnearthefi*ed,stickyclustersum(2:n*-1,2:ny-1)=gnd(2:n*-1,1:ny-2)+gnd(2:n*-1,3:ny)+...gnd(1:n*-2,2:ny-1)+gnd(3:n*,2:ny-1)+...gnd(1:n*-2,1:ny-2)+gnd(1:n*-2,3:ny)+...gnd(3:n*,1:ny-2)+gnd(3:n*,3:ny);%addtotheclustergnd=((sum>0)&(sand==1))|gnd;%andeliminatethemovingparticlesand(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è)粒子，可以在2*2塊中向著底部運(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)作都完成。000000001000010000100001100001001011011110000100011110111110110110110111核心代碼：p=mod(i,2);%margolisneighborhoodsand(n*/2,ny/2)=1;%addagrainatthetop%upperleftcellupdate*ind=[1+p:2:n*-2+p];yind=[1+p:2:ny-2+p];%randomizetheflow--10%ofthetimevary=rand(n*,ny)<.9;vary1=1-vary;sandNew(*ind,yind)=...gnd(*ind,yind).*sand(*ind,yind)+...(1-gnd(*ind,yind)).*sand(*ind,yind).*sand(*ind,yind+1).*...(sand(*ind+1,yind+1)+(1-sand(*ind+1,yind+1)).*sand(*ind+1,yind));sandNew(*ind+1,yind)=...gnd(*ind+1,yind).*sand(*ind+1,yind)+...(1-gnd(*