没有合适的资源?快使用搜索试试~ 我知道了~
温馨提示
元胞自动机已被应用于物理模拟,生物模拟等领域。 元胞自动机(CA)是一种用来仿真局部规则和局部联系的方法。典型的元胞自动机是定义在网格上的,每一个点上的网格代表一个元胞与一种有限的状态。变化规则适用于每一个元胞并且同时进行。典型的变化规则,决定于元胞的状态,以及其( 4 或 8 )邻居的状态。元胞自动机已被应用于物理模拟,生物模拟等领域。本文就一些有趣的规则,考虑如何编写有效的 MATLAB 的程序来实现这些元胞自动机。
资源推荐
资源详情
资源评论












元胞自动机与
MATLAB
引言
元胞自动机(CA)是一种用来仿真局部规则和局部联系的方法。典型的元
胞自动机是定义在网格上的,每一个点上的网格代表一个元胞与一种有限的状
态。变化规则适用于每一个元胞并且同时进行。典型的变化规则,决定于元胞的
状态,以及其( 4 或 8 )邻居的状态。元胞自动机已被应用于物理模拟,生物
模拟等领域。本文就一些有趣的规则,考虑如何编写有效的 MATLAB 的程序来
实现这些元胞自动机。
MATLAB 的编程考虑
元胞自动机需要考虑到下列因素,下面分别说明如何用 MATLAB 实现这些部分。
并以 Conway 的生命游戏机的程序为例,说明怎样实现一个元胞自动机。
z 矩阵和图像可以相互转化,所以矩阵的显示是可以真接实现的。如果矩阵
cells 的所有元素只包含两种状态且矩阵 Z 含有零,那么用 image 函数来显示
cat 命令建的 RGB 图像,并且能够返回句柄。
imh = image(cat(3,cells,z,z));
set(imh, 'erasemode', 'none')
axis equal
axis tight
z 矩阵和图像可以相互转化,所以初始条件可以是矩阵,也可以是图形。以下
代码生成一个零矩阵,初始化元胞状态为零,然后使得中心十字形的元胞状
态= 1。
z = zeros(n,n);
cells = z;
cells(n/2,.25*n:.75*n) = 1;
cells(.25*n:.75*n,n/2) = 1;
z 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);
z
加入一个简单的图形用户界面是很容易的。在下面这个例子中,应用了三个
按钮和一个文本框。三个按钮,作用分别是运行,停止,程序退出按钮。文
框是用来显示的仿真运算的次数。

%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', ...
'fontsize',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]);
经过对控件(和 CA)初始化,程序进入一个循环,该循环测试由回调函数的每
个按钮控制的变量。刚开始运行时,只在嵌套的 while 循环和 if 语句中运行。
直到退出按钮按下时,循环停止。另外两个按钮按下时执行相应的 if 语句。
stop= 0; %wait for a quit button push
run = 0; %wait for a draw
freeze = 0; %wait for a freeze
while (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) + ...
cells(x-1,y-1) + cells(x-1,y+1) + ...
cells(3:n,y-1) + cells(x+1,y+1);
% The CA rule
cells = (sum==3) | (sum==2 & cells);
%draw the new image
set(imh, 'cdata', cat(3,cells,z,z) )
%update the step number diaplay
stepnumber = 1 + str2num(get(number,'string'));
set(number,'string',num2str(stepnumber))
end
if (freeze==1)
run = 0;
freeze = 0;
end
drawnow %need this in the loop for controls to work
end

例子
1 .Conway 的生命游戏机。
规则是:
¾ 对周围的 8 个近邻的元胞状态求和
¾ 如果总和为 2的话,则下一时刻的状态不改变
¾ 如果总和为 3 ,则下一时刻的状态为 1
¾ 否则状态= 0
核心代码:
x = 2:n-1;
y = 2:n-1;
%nearest neighbor sum
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(3:n,y-1) + cells(x+1,y+1);
% The CA rule
cells = (sum==3) | (sum==2 & cells);
2 .表面张力
规则是:
¾ 对周围的 8 近邻的元胞以及自身的状态求和
¾ 如果总和< 4 或= 5 ,下一时刻的状态= 0
¾ 否则状态= 1
核心代码:
x = 2:n-1;
y = 2:n-1;
%nearest neighbor sum
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(3:n,y-1) + cells(x+1,y+1)+...
cells(x,y);
% The CA rule
cells = ~((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 这个图像开始) ,我们可以得
到以下的图像。
4 .激发介质( BZ reaction or heart)
规则:
¾ 元胞有 10 个不同的状态。状态 0 是体眠。1-5 为活跃状态,、6-9 为是极活跃
状态。
¾ 计算每一个处于活跃的状态的元胞近邻的 8 个元胞。
¾ 如果和大于或等于 3 (至少有三个活跃的邻居) ,则下一时刻该元胞= 1 。
¾ 不需要其它输入,1 至 9 种状态依次出现。如果该时刻状态= 1 那么下一时刻
状态= 2 。如果该时刻状态= 2 ,然后下一时刻状态= 3 ,对于其它的状态
依次类推,直到第 9 种状态。如果状态= 9 ,然后下一状态= 0 并且元胞回
到休息状态。
核心代码:

x = [2:n-1];
y = [2:n-1];
sum(x,y) = ((cells(x,y-1)> 0)&(cells(x,y-1)< t)) + ((cells(x,y+1)> 0)&(cells(x,y+1)< t)) + ...
((cells(x-1, y)> 0)&(cells(x-1, y)< t)) + ((cells(x+1,y)> 0)&(cells(x+1,y)< t)) + ...
((cells(x-1,y-1)> 0)&(cells(x-1,y-1)< t)) + ((cells(x-1,y+1)> 0)&(cells(x-1,y+1)< t)) + ...
((cells(x+1,y-1)> 0)&(cells(x+1,y-1)< t)) + ((cells(x+1,y+1)> 0)&(cells(x+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);
一个 CA 初始图形经过螺旋的变化,得到下图。
5 .森林火灾
规则:
¾ 元胞有 3 个不同的状态。状态为 0是空位,状态= 1 是燃烧着的树木,状态
= 2 是树木。
¾ 如果 4 个邻居中有一个或一个以上的是燃烧着的并且自身是树木(状态为
2 ) ,那么该元胞下一时刻的状态是燃烧(状态为 1 ) 。
¾ 森林元胞(状态为 2 )以一个低概率(例如 0.000005 )开始烧(因为闪电)。
¾ 一个燃烧着的元胞(状态为 1 )在下一时时刻变成空位的(状态为 0 ) 。
¾ 空元胞以一个低概率(例如 0.01 )变为森林以模拟生长。
¾ 出于矩阵边界连接的考虑,如果左边界开始着火,火势将向右蔓延,右边界
同理。同样适用于顶部和底部。
核心代码:
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) ;
注意环形连接是由序标实现的。
剩余22页未读,继续阅读
资源评论

- lijian9912023-05-13非常有用的资源,可以直接使用,对我很有用,果断支持!
- oppoaply2022-10-10资源值得借鉴的内容很多,那就浅学一下吧,值得下载!

passionSnail
- 粉丝: 235
- 资源: 4534
上传资源 快速赚钱
我的内容管理 展开
我的资源 快来上传第一个资源
我的收益
登录查看自己的收益我的积分 登录查看自己的积分
我的C币 登录后查看C币余额
我的收藏
我的下载
下载帮助


安全验证
文档复制为VIP权益,开通VIP直接复制
