元胞自动机简单例子

1、1. sierpinski 直角垫片 function sierpinski3_by_CA(n);% 使用元胞自动机生成 sierpinski 直角垫片% Example:% sierpinski3_by_CA(256);% %算法见 : 孙博文 , 分形算法与程序设计：用 Visual C+ 实现 if nargin=0;n=256;endX=zeros(n);X(1,round(n/2)=1;H=imshow(X,); set(gcf,'doublebuffer','on');k=1;while k<round(n/2);X(k+1,2:end-1)=

2、and(xor(X(k,1:end-2),X(k,3:end),. X(k,2:end-1);set(H,'CData',1-X);pause(0.05); k=k+1;endnm=round(n/2);k=1;while k<nm;X(nm+k,1:end)=X(nm-k,1:end);set(H,'CData',1-X);pause(0.05);k=k+1;end2. sierpinski 直角垫片function sierpinski(n);% 使用元胞自动机生成 sierpinski 直角垫片% Example:% sierpinski(256);

3、% %算法见 : 孙博文 , 分形算法与程序设计：用 Visual C+ 实现 if nargin=0;n=256;endX=ones(n);X(1,n-1)=0;H=imshow(X,); set(gcf,'doublebuffer','on');k=1;while k<n;X(k+1,1:end-1)=xor(X(k,1:end-1),X(k,2:end);X(k+1,n)=1;set(H,'CData',X);pause(0.1);k=k+1;end3. 扩散限制凝聚 clc;clear;close all;S=ones(40,100

4、);% state matrixS(end,:)=0; % initial sttaeSs=zeros(size(S)+1,0); % top line is origin of particleSs(2:end,:)=S; % showing matrixN=size(S,2);II=imagesc(Ss);axis equal;colormap(gray) set(gcf,'DoubleBuffer','on');while sum(1-S(1,:)<0.5;y=1;x=round(rand*N-1+1); % random positionD=0;w

5、hile D<0.5; % random travelr=rand;if abs(x-1)<0.1;SL=1;elseSL=S(y,x-1);endif abs(x-N)<0.1;SR=1;elseSR=S(y,x+1);endif SL+SR+S(y+1,x)<2.5; % check the neighbor: left, right, under D=1;S(y,x)=0; % stop in the positionendif r<=1/3; % travel randomlyx=x-1;elseif r<=2/3;x=x+1;elsey=y+1;e

6、ndSs(2:end,:)=S;if x<0.5|x>N+0.5;D=1; % out of the rangeelseSs(y,x)=0; % to show the moving particle end set(II,'CData',Ss); % to show pause(0.1);endend模拟卫星云图function CA_sim_cloud;% 使用元胞自动机模拟地球卫星的云图% reference:% Piazza, E.; Cuccoli, F.;% Cellular Automata Simulation of Clouds in Satell

7、ite Images,% Geoscience and Remote Sensing Symposium, 2001. IGARSS '01.% IEEE 2001 International Volume 4, 9-13 July 2001 Page(s):% 1722 - 1724 vol.4 Digital Object Identifier 10.1109/IGARSS.% 2001.977050time=888; % 程序执行步数 M=240;N=320;S=round(rand(M,N)*15); p=1,2,1,6,6,1,2,1;p=sum(tril(meshgrid(

8、p),2)/20;rand('state',0);SS=S;R=rand(M,N);G=R;B=R;C=cat(3,R,G,B);fig=figure; set(fig,'DoubleBuffer','on'); mov = avifile('example2.avi'); cc=imshow(C,);set(gcf,'Position',13 355 157 194) x1=(1:3)+round(M/2);y1=(1:3)+round(N/3); x2=(1:3)+round(M/3);y2=(1:3)+rou

9、nd(N/2); x3=(1:3)+round(M/1.5);y3=(1:3)+round(N/2);q=0;qq=15/4; while q<time;SS=zeros(M,N);for k=1:15;r=rand(M,N); %生成几率 rK=zeros(M+2,N+2);T=(S-k>=0); %粒子数矩阵K(2:end-1,2:end-1)=T;SS=K(1:end-2,1:end-2).*(r<p(1)+.K(1:end-2,2:end-1).*(r<p(2) & r>=p(1)+. K(1:end-2,3:end).*(r<p(3) &a

10、mp; r>=p(2)+.K(2:end-1,1:end-2).*(r<p(4) & r>=p(3)+. K(2:end-1,3:end).*(r<p(5) & r>=p(4)+. K(3:end,1:end-2).*(r<p(6) & r>=p(5)+. K(3:end,2:end-1).*(r<p(7) & r>=p(6)+. K(3:end,3:end).*(r>=p(7)+SS;endS=SS; %SS 是粒子扩散后的分布 S(S>15)=15;S(x1,y1)=15;S(x2,y2)=1

11、5;S(x3,y3)=15; % 粒子源赋值 G=(S<=7.5);B=(S>qq);R=(S>qq & S<=7.5); C=double(cat(3,R,G,B); set(cc,'CData',C);q=q+1; pause(0.2);title('q=',num2str(q); Nu(q)=sum(S(1:end); F = getframe(gca); mov = addframe(mov,F);endmov = close(mov); figure;plot(Nu)奇偶规则function edwards(N)% 简单

12、元胞自动机奇偶规则 ( 模式 3) 同或运算 % N is the size of calculational matrix% Examples:% figure% edwards(200) warning off M=ones(N);M(fix(29*N/59):fix(30*N/59),fix(29*N/59):fix(30*N/59)=0; close all imshow(M,) for t=1:187;M,Nu=jisuan(M);pause(0.1) imshow(M) HH(t)=Nu; end figure; plot(HH) function Y,Nu=jisuan(M);

13、x,y=find(M=0); Nu=prod(size(x); Xmax=max(max(x); Xmin=min(min(x); Ymax=max(max(y); Ymin=min(min(y); T=ones(Xmax-Xmin+3,Ymax-Ymin+3); T(2:end-1,1:end-2)=M(Xmin:Xmax,Ymin:Ymax); Su=T;T=ones(Xmax-Xmin+3,Ymax-Ymin+3); T(2:end-1,3:end)=M(Xmin:Xmax,Ymin:Ymax); Su=xor(Su,T);Su=not(Su);T=ones(Xmax-Xmin+3,Ym

14、ax-Ymin+3); T(1:end-2,2:end-1)=M(Xmin:Xmax,Ymin:Ymax); Su=xor(Su,T);Su=not(Su);T=ones(Xmax-Xmin+3,Ymax-Ymin+3);T(3:end,2:end-1)=M(Xmin:Xmax,Ymin:Ymax); Su=xor(Su,T);Su=not(Su);M(Xmin-1:Xmax+1,Ymin-1:Ymax+1)=Su;Y=M;森林火灾模拟 close all;clc;clear; figure;p=0.3; % 概率 p f=6e-5; % 概率 f axes;rand('state&#

15、39;,0); set(gcf,'DoubleBuffer','on');% S=round(rand(300)/2+0.5)*2); S=round(rand(300)*2);% copyright: zjliu% Author's email: zjliu2001Sk=zeros(302);Sk(2:301,2:301)=S;% 红色表示正在燃烧 (S 中等于 2 的位置 )% 绿色表示绿树 (S 中等于 1 的位置 )% 黑色表示空格位 (S 中等于 0 的位置 )C=zeros(302,302,3);R=zeros(300);G=zeros(30

16、0);R(S=2)=1;G(S=1)=1;C(2:301,2:301,1)=R;C(2:301,2:301,2)=G;Ci=imshow(C);ti=0;tp=title('T = ',num2str(ti);while 1;ti=ti+1;St=Sk;St(Sk=2)=0; % for rule (1)Su=zeros(302);Sf=Sk;Sf(Sf<1.5)=0;Sf=Sf/2;Su(2:301,2:301)=Sf(1:300,1:300)+Sf(1:300,2:301)+Sf(1:300,3:302)+.Sf(2:301,1:300)+Sf(2:301,3:302

17、)+Sf(3:302,1:300)+.Sf(3:302,2:301)+Sf(3:302,3:302);St(Sf>0.5)=2; % for rule (2)Se=Sk(2:301,2:301);Se(Se<0.5)=4;Se(Se<3)=0;Se(Se>3)=1;St(2:301,2:301)=St(2:301,2:301)+Se.*(rand(300)<p); %for rule (3)Ss=zeros(302);Ss(Sk=1)=1;Ss(2:301,2:301)=Ss(1:300,1:300)+Ss(1:300,2:301)+Ss(1:300,3:302)+.Ss(2:301,1:300)+Ss(2:301,3:302)+Ss(3:302,1:300)+.Ss(3:302,2:301)+Ss(3:302,3:302);Ss(Ss<7.5)=0;Ss(Ss&g