matlab——元胞自动机

（主要参考zhou lvwen教学视频）
元胞自动机是离散的动力学系统，即元胞定义在有限的时间空间，且个数有限
数学表示A=(L,d,s,N,f)
L: 元胞网格空间
d:空间维数
s:有限离散状态集合（例如0，1）
N：某邻域内所有元胞出现的可能性
f：局部映射或局部规则
元胞是其基本单元，每一个都具有记忆储存状态的功能
几种不同邻居关系

边界条件：周期型，定值型，吸收型，反射型
规则：
总和型：看邻居总的状态
合法型：
实例：森林火灾

代码

% simulate forest fire with cellular automata
n = 300;%定义表示森林矩阵大小
Plight = 5e-6; %闪电概率
Pgrowth = 1e-2;%树生长概率
UL = [n 1:n-1]; %上左邻居处理
DR = [2:n 1];
veg=zeros(n,n);%初始化矩阵
imh = image(cat(3,veg,veg,veg));%可视化（R，G，B表示三层矩阵）
% veg = empty=0 burning=1 green=2三种不同状态
for i=1:3000%循环开始
    %nearby fires?
    sum =            (veg(UL,:)==1) + ...
        (veg(:,UL)==1)     +      (veg(:,DR)==1) + ...
                     (veg(DR,:)==1);%求上下左右四个邻居
                 %按照规则更新森林矩阵：树=总-着火+新生

    veg = 2*(veg==2) - ...
          ( (veg==2) & (sum>0 | (rand(n,n)<Plight)) ) + ...
          2*((veg==0) & rand(n,n)<Pgrowth) ;
     
    set(imh, 'cdata', cat(3,(veg==1),(veg==2),zeros(n)) )%可视化
    drawnow
end

实例：单轨交通问题
代码

function [rho, flux, vmean] = ns(rho, p, L, tmax, animation, spacetime)
%
% NS: This script implements the Nagel Schreckenberg cellular automata based
% traffic model. Car move forward governed by NS algorithm:
% USAGE: flux = ns(rho, p, L, tmax, isdraw)
%        rho       = density of the traffic密度即p(rou)
%        p         = probability of random braking随机减速概率
%        L         = length of the load
%        tmax      = number of the iterations时间步
%        animation = if show the animation of the traffic
%        spacetime = if plot the space-time after the simuation ended.
%        flux      = flux of the traffic
%


if nargin == 0; 
    rho = 0.25; p = 0.25; L = 100; tmax = 100; 
    animation = true; spacetime = true;
end

vmax = 5;                  % maximun speed
% place a distribution with density
ncar = round(L*rho);%ncar=L*rho
rho = ncar/L;

x = sort(randsample(1:L, ncar));%1:ncar 随机个不重复的数
v = vmax * ones(1,ncar);   % start everyone initially at vmax

if animation; h = plotcirc(L,x,0.1); end

flux = 0;                  % number of cars that pass through the end
vmean = 0;
road = zeros(tmax, L);

for t = 1:tmax
    % acceleration加速规则
    v = min(v+1, vmax);
    
    %collision prevention防止碰撞
    gaps = gaplength(x,L); % determine the space vehicles have to move车距
    v = min(v, gaps-1);
    
    % random speed drops随机减速
    vdrops = ( rand(1,ncar)<p );
    v = max(v-vdrops,0);
    
    % update the position
    x = x + v;
    passed = x>L;          % cars passed at time r
    x(passed) = x(passed) - L;% periodic boundary conditions周期边界条件
    
    if t>tmax/2
        flux = flux + sum(v/L); %flux = flux + sum(passed); 
        vmean = vmean + mean(v);
    end
    road(t,x) = 1;
    
    if animation; h = plotcirc(L,x,0.1,h); end
end
flux = flux/(tmax/2);
vmean = vmean/(tmax/2);

if spacetime; figure;imagesc(road);colormap([1,1,1;0,0,0]);axis image; end

% -------------------------------------------------------------------------

function gaps = gaplength(x,L)
% 
% GAPLENGTH: determine the gaps between vehicles
%
ncar = length(x);
gaps=zeros(1, ncar);
if ncar>0
    gaps = x([2:end 1]) -x;
    gaps(gaps<=0) = gaps(gaps<=0)+L;
end

% -------------------------------------------------------------------------

function h = plotcirc(L,x,dt,h)
W = 0.05;  R = 1;
ncar = length(x);

theta = [0 : 2*pi/L : 2*pi];
xc = cos(theta);     yc = sin(theta);
xinner = (R-W/2)*xc; yinner = (R-W/2)*yc;
xouter = (R+W/2)*xc; youter = (R+W/2)*yc;

xi = [xinner(x);  xinner(x+1); xouter(x+1); xouter(x)];
yi = [yinner(x);  yinner(x+1); youter(x+1); youter(x)];
if nargin == 3
    color = randperm(ncar);
    h = fill(xi,yi, color); hold on
    plot(xinner,yinner, 'k', xouter,youter, 'k','linewidth',1.5)
    plot([xinner; xouter], [yinner; youter],'k','linewidth',1.5)
    axis image; axis((R+2*W)*[-1 1 -1 1]); axis off
else
    for i=1:ncar;  set(h(i),'xdata',xi(:,i),'ydata',yi(:,i));  end
end
pause(dt)

总结可以看出：元胞自动机具有如下特点
1.离散的空间与时间
2.离散有限状态（树的着火，空地，生长；有车，无车）
3.元胞同质（没有特殊）
4.局部作用与计算
使用：交通问题，传染病问题等
注意：
元胞自动机属于网格模型，注意其余粒子模型区别