元胞自动机(CA)代码及应用
元胞自动机简介
元胞自动机基础元胞自动机（cellular automaton, CA）是最近一个比较热门的研究课题，其是物理、数学、计算机和生物等学科的交叉产物。
在计算机领域中，CA在人工智能、计算复杂性分析以及加密等多个领域中有着较大的用途。
特别是在大约十年前，密码学家H. Gutowitz根据CA的基本原理，提出了分块加密算法CA-1.1，使得CA在密码学中真正的迈出了第一步，也使得越来越多的密码学家开始了对CA的研究。
最近，我也开始对这个方面产生了浓厚的兴趣，并开始了一些学习，就先来简单的说说什么是CA吧！简单的说，元胞自动机是一个空间、时间和状态上都离散的动态系统。
构成CA的基本单位成为元胞（cellular），规则的分布在元胞空间（spatial lattice）的格点上，且各自的状态随着时间按照一定的局部规则变化。
也就是说，元胞的状态只能从一个有限的状态集中取值，每个时刻元胞的状态仅与其自身和邻居在上一时刻的状态有关，并且，所有的元胞在每个时刻均是同时更新的。
以上即是对CA的一个定性的描述，下面给出一个基于集合论的定量描述（L. Hurd等）：设d为CA空间的维数，k代表元胞的状态，集合S表示CA的整体状态，r表示元胞的邻居半径。
为了简单起见，我们在d=1，即一维空间上对CA进行讨论。
CA的动态性可以由一个全局函数F: St→St+1决定，并且，每个元胞的状态可以由一个局部函数f:kt→kt+1决定。
由于多维空间的CA具有很强的复杂性，故目前对CA的研究主要集中在一维和二维空间。
就一维空间而言，CA的结构显然只有可能是线性结构。
在二维空间，CA的结构可能有三角、四边或多边等构成方式。
显然，结构上的差异会对其在计算机表示及其他部分特性上带来一定的差异。
而CA 的邻居结构也通常包括Von. Neumann、Moore、扩展Moore和Margolus等多种形态，不同的邻居结构带来的特性和复杂度也不尽相同。
7月18日元胞自动机
3. 元胞空间 元胞所分布在的空间网点集合就是 元胞空间 A、元胞空间的几何划分 任意维数的欧几里 德空间规则划分。对于一维元胞自动机，元胞 空间划分只有一种。而高维的元胞自动机，元 胞空间的划分则可能有多种形式。对于常见的 二维自动机，元胞空间通常可按三角形、四边 形或六边形三种网格排列。
Triangle
基于个体的自底向上的研究方法： 程序的行为完全由它的内部机制决定， 通常将个体与程序相连，所模拟的复杂 现象包括许多个体。在计算机里生成一 个与真实世界对等的虚拟的人工世界， 通常这个虚拟的世界包括许多个体，而 这许多个体的行为呈现为复杂性。以此 来探讨微观的个体行为和宏观复杂性之 间的关系。
二.常用元胞自动机
元胞自动机特性
• 把一个空间划分成网络，每一个点表示一个元胞，它们的状态赋值， 在网格中用颜色的变化来表示，在事先设定的规则下，元胞的演化就用 网格颜色的变化来描述，这样的模型就是元胞自动机。 • • 通过对元胞自动机这些网络中的格点的不同定义，以及初始条件的 元胞自动机的基本特征： 离散性：元胞自动机是高度离散的。它不仅仅空间离散时间 离散，而且在函数值，即元胞的状态值也是离散的。 不同，可以模拟出不同的现象和过程。
元胞自动机
元胞自动机是探索复杂系统中局部 — 整体 互动关系的最简单模式，视为演化分析的 基本计算模型；
物理学：一个离散的无穷维的动力学系统
数学：一个时空离散的数学模型，描述连续对 象的偏微分方程的对立体 计算机科学：新兴的人工智能、人工生命的分 支 生物学：生命现象的一种描述
Square
Hexagon
三角网格拥有较少的邻居数目，这在某些时候很 有用。缺点是计算机的表达与显示不方便。 四边形网格直观简单，特别适合于计算机环境下 进行表达显示。 六边形网格能较好的模拟各向同性的现象，因此， 模型能更加自然而真实。其缺点同正三角网格 一样，在表达显示上较为困难和复杂。
元胞自动机在城市扩展方面的应用综述
元胞自动机在城市扩展方面的应用综述摘要本文在介绍元胞自动机各要素的基础上，综述了元胞自动机用于城市扩展模拟的历史、元胞自动机用于城市扩展模拟的具体研究方向，即在具体的模型中如何确定模型的结构和参数，并对其未来的发展趋势进行了展望，并指出CA中的转换规则的扩展是在将来的研究中的一个首要问题。
关键字：元胞自动机；城市扩展模拟；转换规则一引言元胞自动机（CA）是一种时间、空间、状态都离散，空间的相互作用及时间上的因果关系皆局部的网格动力学模型，其“自下而上”的研究思路，强大的复杂计算功能、固有的平行计算能力、高度动态以及具有空间概念等特征,使得它在模拟空间复杂系统的时空动态演变方面具有很强的能力。
在城市空间动态变化的模拟研究方面，CA模型已应用到除非洲、南极洲的所有大洲的城市模拟研究当中。
CA模型和GIS的集成,一方面增强GIS的空间模型运算及分析能力，另一方面，GIS提供的强大空间处理能力可以为CA模型准备数据和定义有效的元胞转换规则以及对模拟结果进行可视化。
同时CA模型还可以与神经网络、主成分分析、遗传算法、模糊逻辑以及其他研究方法相结合，以增强其在城市空间变化模拟研究方面的能力。
将CA与MAS技术相结合，建立一个能够模拟多个不同参与因子（自然系统）、不同决策者（人文系统）共同影响下的城市发展模型，以此来模拟与预测城市发展的真实状况，将是CA模型在城市空间变化模拟与预测研究中的未来发展趋势。
国内元胞自动机应用研究起步较晚，受国际研究的推动，20世纪90年代末地理学界才开始类似的尝试研究，主要集中在基于元胞自动机的LUC（和城市增长模拟，罗平从经典地理过程分析的基本理论人手，分析和阐述了CA寸于经典，地理过程分析概念的表达程度的局限性，综合地理系统的几何属性和非几何属性提出了基于地理特征概念的元胞自动机（GeoFeature 一CA）,周成虎等人在Batty和Xie的DUE模型的基础上，构建了面向对象的、随机的、不同构的和两个CA莫型耦合的GeoCA- Urban模型，并成功模拟了深圳特区土地利用动态演化过程。
三维元胞自动机及其在计算机上的实现 2007
第07卷 第08期 中 国 水 运 Vol.7 No.08 2007年 08月 China Water Transport August 2007收稿日期：2007-5-15作者简介：许 林 男（1980—） 蓝天学院制造系 教师 （330098）徐 明 蓝天学院机械系 （330098）研究方向：数值模拟三维元胞自动机及其在计算机上的实现许 林 徐 明摘 要：元胞自动机是一种时空、状态离散的数学模型，适于计算机模拟实施。
特别适合于处理那些难以用数学定量描述的复杂动态体系问题，如材料微观组织的演变。
本文阐述了元胞自动机法的产生和发展。
并以Visual C++为编译平台，运用OpenGL 图形函数库建立了一种三维元胞自动机模型。
该模型具备了经典元胞自动机的基本特征，因此可以根据需要进行扩展。
文中运用该模型进行了简化的枝晶生长模拟，并与二维的模拟结果进行比较，验证了该模拟的正确性。
关键词：元胞自动机 三维建模 OpenGL中图分类号：TP211 文献标识码：A 文章编号：1006-7973（2007）08-0171-02 一、引言元胞自动机(CA) 又称细胞自动机是建立于细胞发育演化基础上的时空离散、状态离散的并行数学模型[1]。
元胞自动机最早是由数学家、物理学家John Von Neumann 和Stanislaw Ulam 在1940年提出的[2]。
但从应用角度看，直到1960年John Horton Conway 运用元胞自动机建立了一种“生命游戏”后[3,4,5]，元胞自动机才得到广泛的运用。
80年代，由于元胞自动机这类简单模型能十分方便地复制出复杂的现象或动态演变过程中的吸引、自组织和混沌现象，从而引起了物理学家、计算机科学家的极大兴趣，并在许多领域得到了应用，如混沌、分形的产生[6]，模式分类[7]，图像处理[8]，智能材料[9]，复杂现象[1]等，并提出了许多变形的元胞自动机，如以凝固理论为演化规则的元胞自动机[10]，模糊元胞自动机[11]，神经元胞自动机[12]等。
元胞自动机模拟害虫防治python代码
元胞自动机模拟害虫防治python代码元胞自动机（Cellular Automaton）是一种离散模型，由一组按照某种规则进行演化的格子组成。
我们可以使用元胞自动机来模拟害虫的防治过程。
下面是一个简单的Python代码示例，用于模拟害虫的扩散和防治过程。
注意：这个代码只是一个基础的示例，并未考虑许多真实世界中的复杂因素，例如环境变量、害虫种类、防治策略等。
pythonimport numpy as npimport matplotlib.pyplot as plt# 初始化元胞自动机def initialize_grid(size):grid = np.zeros((size, size))# 假设害虫在初始状态位于中心grid[size//2, size//2] = 1return grid# 更新元胞自动机状态def update_grid(grid, size, threshold, kill_prob):new_grid = grid.copy()for i in range(1, size-1):for j in range(1, size-1):# 计算邻居中的害虫数量neighbors = np.sum(grid[(i-1:i+2, j-1:j+2)]) - grid[i, j] if neighbors > threshold:# 如果邻居中的害虫数量超过阈值，当前格子可能会变为害虫new_grid[i, j] = 1 if np.random.rand() < 0.5 else 0elif neighbors == 0:# 如果邻居中没有害虫，当前格子可能会被防治new_grid[i, j] = 0 if np.random.rand() > kill_prob else 1 return new_grid# 绘制元胞自动机状态def plot_grid(grid, title):plt.imshow(grid, cmap='Greys', interpolation='nearest')plt.title(title)plt.show()# 主程序size = 100 # 元胞自动机大小threshold = 5 # 害虫繁殖阈值kill_prob = 0.8 # 防治成功率grid = initialize_grid(size)for i in range(100): # 模拟100个时间步plot_grid(grid, f'Step {i}')grid = update_grid(grid, size, threshold, kill_prob)在这个代码中，我们首先初始化了一个大小为size的元胞自动机，并在中心位置放入一个害虫。
元胞自动机在复杂系统建模中的应用
元胞自动机在复杂系统建模中的应用元胞自动机（Cellular Automata，简称CA）是一种用来描述复杂系统行为的数学模型。
它由一组简单的单元（cell）组成，在一个由相同大小的正方形格子（grid）构成的网格上进行演化。
每个单元可以处于不同的状态，并通过更新规则与其邻居进行交互。
尽管元胞自动机的规则非常简单，但它已被广泛应用于生物、物理、社会科学等领域的复杂系统建模中。
本文将介绍元胞自动机在复杂系统建模中的应用，并探讨其优势和局限性。
元胞自动机最早由美国数学家John von Neumann和Stanislaw Ulam 于20世纪40年代提出。
它广泛应用于不同领域，例如生物学中的细胞生长模拟、物理学中的颗粒传输模拟、社会科学中的城市规划模拟等。
元胞自动机的简单规则和复杂行为之间的关系使其成为复杂系统建模中的强大工具。
首先，元胞自动机在生物学中的应用非常广泛。
生物系统中的许多现象可以通过元胞自动机来模拟和解释。
例如，在细胞生长过程中，细胞与周围细胞进行相互作用，从而形成特定的模式和结构。
通过模拟和研究这些交互作用，科学家可以更好地理解生物系统的发展和演化规律。
元胞自动机还可用于模拟病原体传播、生态系统动力学、遗传算法等生物学问题，为生物学研究提供了新的视角和方法。
其次，元胞自动机在物理学中的应用也非常突出。
在物质传输和分布的模拟中，元胞自动机可以精确地描述粒子之间的相互作用和运动规律。
通过定义单元的状态和更新规则，元胞自动机可以模拟物质在介质中的传输、扩散、聚集等复杂过程。
这种建模方法在材料科学、地球科学、天体物理学等领域得到了广泛应用，为研究人员提供了一种高效而有效的模拟工具。
此外，元胞自动机在社会科学中也有重要的应用。
社会系统是一种充满复杂性和非线性特征的系统，元胞自动机能够较好地刻画其内部的各种相互作用和演化规律。
例如，在城市规划模拟中，通过设定不同的元胞状态和邻居交互规则，可以模拟城市人口密度、交通流动、资源分配等问题，为城市规划者提供决策支持和优化方案。
元胞自动机简单例子
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)=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);% %算法见:孙博文,《分形算法与程序设计：用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);% 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;while 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;endSs(2:end,:)=S;if x<0.5|x>N+0.5;D=1; % out of the rangeelseSs(y,x)=0; % to show the moving particleendset(II,'CData',Ss); % to showpause(0.1);endend模拟卫星云图function CA_sim_cloud;% 使用元胞自动机模拟地球卫星的云图%% reference:% Piazza, E.; Cuccoli, F.;% Cellular Automata Simulation of Clouds in Satellite 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(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)+round(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) & 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)=15;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)% 简单元胞自动机—奇偶规则(模式3)同或运算% N is the size of calculational matrix% Examples:% figure% edwards(200)warning offM=ones(N);M(fix(29*N/59):fix(30*N/59),fix(29*N/59):fix(30*N/59))=0; close allimshow(M,[])for t=1:187;[M,Nu]=jisuan(M);pause(0.1)imshow(M)HH(t)=Nu;endfigure;plot(HH)function [Y,Nu]=jisuan(M);[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,Ymax-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; % 概率pf=6e-5; % 概率faxes;rand('state',0);set(gcf,'DoubleBuffer','on');% S=round((rand(300)/2+0.5)*2);S=round(rand(300)*2);% \copyright: zjliu% Author's email: zjliu2001@Sk=zeros(302);Sk(2:301,2:301)=S;% 红色表示正在燃烧(S中等于2的位置)% 绿色表示绿树(S中等于1的位置)% 黑色表示空格位(S中等于0的位置)C=zeros(302,302,3);R=zeros(300);G=zeros(300);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)+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>7.5)=1;d=find(Ss==1 & Sk==1);for k=1:length(d);r=rand;St(d(k))=round(2*(r<=f)+(r>f));end % for rule (4)Sk=St;R=zeros(302);G=zeros(302);R(Sk==2)=1;G(Sk==1)=1;C(:,:,1)=R;C(:,:,2)=G;set(Ci,'CData',C);set(tp,'string',['T = ',num2str(ti)]) pause(0.2);end。
CA元胞自动机优化模型原代码
CA优化模型原代码:M=load(‘d:\ca\jlwm’)N=load(‘d:\ca\jlwn.asc’)lindishy=load(‘d:\ca\ldfj3.asc’)caodishy=load(‘d:\ca\cdfj3.asc’)gengdishy=load(‘d:\ca\htfj3.asc’)[m,n]=size(M);Xr=[1 1 -1 1 1 1 -1 -1 1 1;1 1 1 1 -1 -1 1 1 1 -1;-1 1 1 1 -1 -1 -1 1 -1 -1;1 1 1 1 1 1 -1 1 1 I; l -1 -1 1 1 -1 -1 -1 1 1;1 -1 -1 1 -1 1 -1 1 -1 -1;-1 1 -1 -1 -1 -1 1 -1 -1 -1;-1 1 1 1 -1 1 -1 1 -1 -1;1 1 -1 1 1 -1 -1 -1 1 1;1 -1 -1 1 1 -1 -1 -1 1 1];caodi=0;lindi=0;gengdi=0;for i=1:mforj=l:nif M(i,j)==4caodi=caodi+1;elseif M(i,j)==3lindi=lindi+1;elseif M(i,j)==2gengdi=gengdi+1;endendendfor i=1:mfor j=1:nif M(i,j)==4if lindishy(i,j)>gengdishy(i,j)if lindishy(i,j)>caodishy(i,j)z=0;for P=max(1,i-1):min(i+1,m)for q=max(j-1,1):min(j+1,n)if (M(p,q)~=0)&&xr(M(p,q),3)==-1z=1;endendendif z== 0caodi=eaodi-1;M(i,j)=3;lindi=lindi+1;endelseif lindishy(i,j)==caodishy(i,j)caoditemp=0;linditemp=0;gengditemp=0;for p=max(1,i-1):min(i+1,m)for q=max(j-1,1):min(i+1,n)if N(p,q)==4caoditemp=caoditemp+1;elseif N(p,q)==3linditemp=linditemp+1;elseif N(p,q)=2gengditemp=gengditemp+1;endendendif linditemp>=max(caoditemp,gengditemp) z=0;for p=max(1,j-1):min(i+1,m)for q=max(j-1,1):min(j+1,n)if(M(p,q)~=0)&&xr(M(p,q),3)==-1;z=1;endendendif Z==0caodi=caodi-1;M(i,j)=3;lindi=lindi+1;endendendelseif lindishy(i,j)==gengdishy(i,j)if lindishy(i,j)>caodishy(i,j)caoditemp=0:linditemp=0;gengditemp=0:for p=max(1,i-1):min(i+1,m)for q=max(j-1,1):min(j+1,n)if N(p,q)==4caoditemp=caoditemp+1;elseif N(p,q)==3linditemp=linditemp+1;elseif N(p,q)==2gengditemp=gengdltemp+1;endendendif linditemp>=gengditempfor p=max(1,j-1):min(i+1,m)for q=max(j-1,1):min(j+1,n)if (M(p,q)~=0)&&xr(M(p,q),3)==-1z=1;endendendif z==0caodi=caodi-1;M(i,j)=3;lindi=lindi+1;endendelseif lindishy(i,j)==caodishy(i,j) caoditemp=0;linditemp=0;gengditemp=0;for p=max(i,i-1):min(i+1,m)for q=max(j-1,1):min(j+1,n)if N(p,q)==4caoditemp=caoditemp+1;elseif N(p,q)==3linditemp=linditemp+1;elseif N(p,q)==2gengditemp=gengditemp+1;endendendif linditemp>=max(caoditemp,gengditemp) z=0;for p=max(1,i-1):min(i+1,m)for q=max(j-1,1):min(j+1,n)if(M(p,q)~=0)&&xr(M(p,q),3)==-1z=1;endendendif z==0caodi=caodi-1;M(i,j)=3;lindi=lindi+1;endendendelseif M(i,j)==2if lindishy(i,j)>gengdishy(i,j)if lindishy(i,j)>caodishy(i,j)z=0;for p=max(1,i-1):min(i+1,m)for q=max(j-1,1):min(j+1,n)if(M(p,q)~=0)&&xr(M(p,q),3)==-1z=1;endendendif z==0M(i,j)=3;Lindi=lindi+1;gengdi=gengdi-1;endelseif lindishy(i,j)==caodishy(i,j) caoditemp=0;linditemp=0;gengditemp=0;for p=max(1,i-1):min(i+1,m)for q=max(j-1,1):min(j+1,n)if N(p,q)==4caoditemp=caoditemp+1;elseif N(p,q)==3linditemp=linditemp+1;elseif N(p,q)==2gengditemp=gengditemp+1;endendendif linditemp>=max(caoditemp,gengditemp) z=0;for p=max(1,i-1):min(i+1,m)for q=max(j-1,1):min(j+1,n)if(M(p,q)~=0)&&xr(M(p,q),3)==-1z=1;endendendif z==0M(i,j)=3;lindi=lindi+1;gengdi=gengdi-1endelseif caoditemp>= gengditemp biaoji=0;for p=max(1,i-1):min(i+1,m)for q=max(j-1,1):min(j+1,n)if(M(p,q)~=0)&&xr(M(p,q),4)==-1 z=1;endendendif z==0M(i,j)=4;caodi=caodi+1;gengdi=gengdi-1;endendelseif lindishy(i,j)<caodishy(i,j)z=0;for p=max(1,i-1):min(i+1,m)for q=max(j-1,1):min(j+1,n)if(M(p,q)~=0)&&xr(M(p,q),4)==-1 z=1;endendendif z==0M(i,j)=4;caodi=caodi+1;gengdi=gengdi-1;endendelseif lindishy(i,j)==gengdishy(i,j) if lindishy(i,j)<caodishy(i,j)z=0;for p=max(1,i-1):min(i+1,m)for q=max(j-1,1):min(j+1,n)if(M(p,q)~=0)&&xr(M(p,q),4)==-1 z=1;endendendif biaoji==0M(i,j)=4;caodi=caodi+1;gengdi=gengdi-1;endelseif lindishy(i,j)>caodishy(i,j) caoditemp=0;linditemp=0;gengditemp=0;for p=max(1,i-1):min(i+1,m)for q=max(j-1,1):min(j+1,n)if N(p,q)==4caoditemp=caoditemp+1;elseif N(p,q)==3linditemp=linditemp+1;elseif N(p,q)==2gengditemp=gengditemp+1; endendendif linditemp>= gengditempz=0;for p=max(1,i-1):min(i+1,m)for q=max(j-1,1):min(j+1,n)if(M(p,q)~=0)&&xr(M(p,q),3)==-1 z=1;endendendif z==0M(i,j)=3;lindi=lindi+1;gengdi=gengdi-1endendelsecaoditemp=0;linditemp=0;gengditemp=0;for p=max(1,i-1):min(i+1,m)for q=max(j-1,1):min(j+1,n)if N(p,q)==4caoditemp=caoditemp+1;elseif N(p,q)==3linditemp=linditemp+1;elseif N(p,q)==2gengditemp=gengditemp+1;endendendif linditemp>=max(caoditemp,gengditemp) z=0;for p=max(1,i-1):min(i+1,m)for q=max(j-1,1):min(j+1,n)if(M(p,q)~=0)&&xr(M(p,q),3)==-1z=1;endendendif z==0M(i,j)=4;lindi=lindi+1;gengdi=gengdi-1endelseif caoditemp>= gengditempz=0;for p=max(1,i-1):min(i+1,m)for q=max(j-1,1):min(j+1,n)if(M(p,q)~=0)&&xr(M(p,q),4)==-1z=1;endendendif z==0M(i,j)=4;caodi=caodi+1;gengdi=gengdi-1;endendendelseifgengdishy(i,j)<caodishy(i,j)z=0;for p=max(1,i-1):min(i+1,m)for q=max(j-1,1):min(j+1,n)if(M(p,q)~=0)&&xr(M(p,q),4)==-1z=1;endendendif z==0M(i,j)=4;caodi=caodi+1;gengdi=gengdi-1;endelseif gengdishy(i,j)==caodishy(i,j) elseif lindishy(i,j)==caodishy(i,j) caoditemp=0;linditemp=0;gengditemp=0;for p=max(1,i-1):min(i+1,m)for q=max(j-1,1):min(i+1,n)if N(p,q)==4caoditemp=caoditemp+1;elseif N(p,q)==3linditemp=linditemp+1;elseif N(p,q)=2gengditemp=gengditemp+1; endendendif caoditemp>= caoditempz=0;for p=max(1,j-1):min(i+1,m)for q=max(j-1,1):min(j+1,n)if(N(p,q)~=0)&&xr(M(p,q),4)==-1; z=1;endendendif Z==0M(i,j)=4;caodi=caodi+1;gengdi=gengdi-1;endendendendendendendfid=fopen(‘d:\ca\lucc’,’at+’)for i=1:mfor j=1:nif M(i,j)>4.5M(i,j)=N(i,j);endfprintf(fid,’%d’, M(i,j)); endfprintf(fid,,’\n’);endfclose(fid);。
- 1、下载文档前请自行甄别文档内容的完整性，平台不提供额外的编辑、内容补充、找答案等附加服务。
- 2、"仅部分预览"的文档,不可在线预览部分如存在完整性等问题,可反馈申请退款(可完整预览的文档不适用该条件!)。
- 3、如文档侵犯您的权益，请联系客服反馈,我们会尽快为您处理(人工客服工作时间：9:00-18:30)。
元胞自动机(CA)代码及应用引言元胞自动机（CA）是一种用来仿真局部规则和局部联系的方法。
典型的元胞自动机是定义在网格上的，每一个点上的网格代表一个元胞与一种有限的状态。
变化规则适用于每一个元胞并且同时进行。
典型的变化规则，决定于元胞的状态，以及其（4或8 ）邻居的状态。
元胞自动机已被应用于物理模拟，生物模拟等领域。
本文就一些有趣的规则，考虑如何编写有效的MATLAB的程序来实现这些元胞自动机。
MATLAB的编程考虑元胞自动机需要考虑到下列因素，下面分别说明如何用MATLAB实现这些部分。
并以Conway的生命游戏机的程序为例，说明怎样实现一个元胞自动机。
●矩阵和图像可以相互转化，所以矩阵的显示是可以真接实现的。
如果矩阵cells的所有元素只包含两种状态且矩阵Z含有零，那么用image函数来显示cat命令建的RGB图像，并且能够返回句柄。
imh = image(cat(3,cells,z,z));set(imh, 'erasemode', 'none')axis equalaxis tight●矩阵和图像可以相互转化，所以初始条件可以是矩阵，也可以是图形。
以下代码生成一个零矩阵，初始化元胞状态为零，然后使得中心十字形的元胞状态= 1。
z = zeros(n,n);cells = z;cells(n/2,.25*n:.75*n) = 1;cells(.25*n:.75*n,n/2) = 1;●Matlab的代码应尽量简洁以减小运算量。
以下程序计算了最近邻居总和，并按照CA规则进行了计算。
本段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);●加入一个简单的图形用户界面是很容易的。
在下面这个例子中，应用了三个按钮和一个文本框。
三个按钮，作用分别是运行，停止，程序退出按钮。
文框是用来显示的仿真运算的次数。
%build the GUI%define the plot buttonplotbutton=uicontrol('style','pushbutton',...'string','Run', ...'fontsize',12, ...'position',[100,400,50,20], ...'callback', 'run=1;');%define the stop buttonerasebutton=uicontrol('style','pushbutton',...'string','Stop', ...'fontsize',12, ...'position',[200,400,50,20], ...'callback','freeze=1;');%define the Quit buttonquitbutton=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]);经过对控件（和CA）初始化，程序进入一个循环，该循环测试由回调函数的每个按钮控制的变量。
刚开始运行时，只在嵌套的while循环和if语句中运行。
直到退出按钮按下时，循环停止。
另外两个按钮按下时执行相应的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 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);% The CA rulecells = (sum==3) | (sum==2 & cells);%draw the new imageset(imh, 'cdata', cat(3,cells,z,z) )%update the step number diaplaystepnumber = 1 + str2num(get(number,'string'));set(number,'string',num2str(stepnumber))endif (freeze==1)run = 0;freeze = 0;enddrawnow %need this in the loop for controls to workend例子1 .Conway的生命游戏机。
规则是：➢对周围的8个近邻的元胞状态求和➢如果总和为2的话，则下一时刻的状态不改变➢如果总和为3 ，则下一时刻的状态为1➢否则状态= 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) + cells(x-1,y+1) + ...cells(3:n,y-1) + cells(x+1,y+1);% The CA rulecells = (sum==3) | (sum==2 & cells);2 .表面张力规则是：➢对周围的8近邻的元胞以及自身的状态求和➢如果总和< 4或= 5 ，下一时刻的状态= 0➢否则状态= 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) + 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 .渗流集群规则：➢对周围相邻的8邻居求和（元胞只有两种状态，0或1 ）。
元胞也有一个单独的状态参量（所谓'记录' ）记录它们之前是否有非零状态的邻居。
➢在0与1之间产生一个随机数r 。
➢如果总和> 0 （至少一个邻居）并且r >阈值，或者元胞从未有过一个邻居，则元胞= 1 。
➢如果总和> 0则设置"记录"的标志，记录这些元胞有一个非零的邻居。
核心代码：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 是图像的尺寸。
最初的图形是由图形操作决定的。
以下程序设定坐标系为一个固定的尺寸，在坐标系里写入文本，然后获得并返回坐标内的内容，并用getframe 函数把它们写入一个矩阵ax = axes('units','pixels','position',[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', [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;经过几十个时间间隔（从MCM Cellular Automata 这个图像开始） ，我们可以得到以下的图像。