超简单有趣的模拟算法：元胞自动机（CA）原理简介与 matlab 代码实现

元胞自动机的整理与学习笔记：介绍了CA的概念——元胞、状态、元胞空间、邻居、边界、传播规则，并解释了原理。配合matlab代码和详细注释完成实例：奇偶规则、生存游戏、森林火灾模拟，并在森林火灾模拟中考虑了风速等条件

铖铖的花嫁

7721人浏览 · 2023-04-27 11:43:15

铖铖的花嫁 · 2023-04-27 11:43:15 发布

很久之前就就听说了元胞自动机（cellular automata，CA），但是一直没有尝试。得知2023年美赛春季赛只有两道赛题的时候，怕考到这个，所以出一篇元胞自动机的博客，权且当一篇学习笔记。

尝试过后才发现元胞自动机其实并没有那么难，非常好理解，算是万金油的模拟算法，可以作为保底的plan B。

1 原理介绍

1.1 概念

在这里插入图片描述

元胞：模拟的时候就是一个有记忆存储状态功能的格子，是元胞自动机的最基本的组成部分。因此元胞具有以下3个特点：

元胞自动机最基本的单元。
元胞有记忆贮存状态的功能。
所有元胞状态都按照元胞规则不断更新。

状态：每个元胞都有自己的状态，随着时间的推移，状态也会改变。

元胞空间：元胞分布在空间中的格点，可以是1维、2维或n维。

邻居：能影响元胞下一时刻状态的周围元胞。一般来说有三种邻居模式（图是2维情况下的）：

在这里插入图片描述

当然还有其他的邻居模式，不赘述，主要看题目要求。

边界：整个元胞空间不是无穷大的，是需要一个边界范围的。边界也有好几种类型：固定型、周期型、绝热型、映射型……，个人感觉最后还是得看需要解决的问题，具体问题具体分析。

规则：根据元胞自身及邻居元胞的状态，决定下一时刻元胞状态的决定，可以是动力学函数或状态转移方程。

1.2 原理

这个没啥原理，这部分简单讲讲步骤：

Step1：建立假设、简化条件。

Step2：设置初始状态和边界。

Step3：推导状态转方程和条件。

Step4：根据时间迭代。

基于大量的实验，人们发现元胞自动机根据动力学行为能分为 4 类（Wolfram. S.，1986)：

平稳型：自任何初始状态开始，经过一定时间运行后，元胞空间趋于一个空间平稳的构形，这里空间平稳即指每一个元胞处于固定状态。不随时间变化而变化。
周期型：经过一定时间运行后，元胞空间趋于一系列简单的固定结构（Stable Patterns）或周期结构（Perlodical Patterns)。由于这些结构可看作是一种滤波器（Filter），故可应用到图像处理的研究中。
混沌型：自任何初始状态开始，经过一定时间运行后，元胞自动机表现出混沌的非周期行为，所生成的结构的统计特征不再变止，通常表现为分形分维特征。
复杂型：出现复杂的局部结构，或者说是局部的混沌，其中有些会不断地传播。

我们不对每种类型的元胞做实验验证，只提供几种经典和常见的例子。

2 实例

2.1 奇偶游戏

条件：

邻居和为奇数时，中间的数字变成1；
邻居和为偶数时，中间的数字变成2。

%%奇偶规则游戏
clear;clc;
n = 1000;%指定边界长度
Se = zeros(n);
z = zeros(n);
Se(n/2-2:n/2+2,n/2-2:n/2+2)=1;  %初始化中间的点
Ch = imagesc(Se);
axis square;
Sd = zeros(n+2);  %边界
while(true) %死循环
    % 保存现有状态的右2：边界1，右边1
    Sd(2:n+1,2:n+1) = Se;
    % 上下左右加一次
    sumValue = Sd(1:n,2:n+1)+Sd(3:n+2,2:n+1)+Sd(2:n+1,1:n)+Sd(2:n+1,3:n+2);
    % 摸一下
    Se = mod(sumValue,2);
    set(Ch,'cdata',Se);
    pause(0.0001)
end

9倍速效果如下

在这里插入图片描述

还挺好看的

不是个实际例子，但是可以看出元胞自动机的代码模板，反正就是这么一个流程。

2.2 生灭游戏

条件：

如果活细胞周围活细胞数小于2个或者大于3个，则转为死；
如果活细胞周围活细胞数为2个或者3个，则保持活；
如果死细胞周围有3个活细胞，则转为活。

对于活细胞，我们用1表示；对于死细胞，我们用0表示。

n = 200;
p = 0.4;
z = zeros(n);
Se = rand(n)<p;
Sd = zeros(n+2);
Ph = imagesc(Se);
while(true)
    Sd(2:n+1,2:n+1)=Se;
    sumValue = Sd(1:n,1:n)+Sd(1:n,2:n+1)+Sd(1:n,3:n+2)+Sd(2:n+1,1:n)+Sd(2:n+1,3:n+2)+Sd(3:n+2,1:n)+Sd(3:n+2,2:n+1)+Sd(3:n+2,3:n+2);
    for i=1:n
        for j=1:n
            if(sumValue(i,j)==3||(sumValue(i,j)==2&&Se(i,j)==1))
                Se(i,j) = 1;
            else
                Se(i,j) = 0;
            end
        end
    end
    set(Ph,'cdata',Se);
    pause(0.05);
end

二倍速效果如下：

在这里插入图片描述

丑丑的

最后趋于平衡了

2.3 森林火灾

条件：

对于树木：它可能会被雷劈，导致树木着火变成火树；它可能会被附近的火树点着，变成火树；也可能保持现状
对于火树：它可能会烧成灰，变成空地；也可能会继续烧，保持现状
对于空地：它可能会长出树木，变成树；也可能保持现状
对于火树传播：火树附近1格的树有概率变成火树

各种概率简单起见我自己定了，如果是比赛，可能需要假设和计算概率方程。

2.3.1 普通森林大火

此时我们不考虑其他因素，简单地将火烧成空地的概率定位0.5，树被雷劈的概率定义为0.000001，空地长树的概率定为 0.01，树被旁边火烧到的概率定为0.7。

clear;clc;
%火灾
n = 300;                           % 定义表示森林的矩阵大小
k = 30000;                         % 迭代次数
Pground = 0.2;                     % 从着火变成空地的概率
Plight = 1e-6; Pgrowth = 2e-3;      % 定义闪电和生长的概率  
P2=0.7; %旁边有火，树着火的概率
% UL = [n,1:n-1]; DR = [2:n,1];      % 定义上左，下右邻居
veg=zeros(n,n)+2;                    % 初始化表示森林的矩阵
imh = image(cat(3,veg,veg,veg));   % 可视化表示森林的矩阵
Sd = zeros(n+2);                 %边界
% veg = 空地为0 着火为1 树木为2
for i=1:k
    Sd(2:n+1,2:n+1) = veg;
    sumValue = (Sd(1:n,2:n+1)==1)+(Sd(2:n+1,1:n)==1)+(Sd(2:n+1,3:n+2)==1)+(Sd(3:n+2,2:n+1)==1);
    for p=1:n
        for q=1:n
            if(veg(p,q)==2 && ((sumValue(p,q)>0 && rand()<P2)||rand()<Plight))
                %首先要是树，而且需要邻居有火，就会一定概率着火;或者被雷劈了,就会直接着火
                veg(p,q)=1;
            elseif(veg(p,q)==1&&rand()<Pground)
                %如果是火且满足概率，则变为空地
                veg(p,q) = 0;
            elseif(veg(p,q)==0&&sumValue(p,q)==0&&rand()<Pgrowth)
                %如果是空地，且周围没有火，那么以一定概率长成树
                veg(p,q) = 2;
            end
        end
    end
    set(imh, 'cdata', cat(3,(veg==1),(veg==2),zeros(n)))
    drawnow  
end

7倍速效果如下图

在这里插入图片描述

还是能看出来森林大火和树生长维持动态平衡的。

2.3.2 东风森林大火

在这里插入图片描述
如上图所示，如果是北风，不同颜色的火树（圆圈）对风方向2格的树木的影响是这样的（X代表是否会被火树影响）。可以看出传播方向越边缘的树，越难被火树波及，在火树林正前方的树很容易被点着。

其实没难多少，就是加了个判断和概率，计算的时候再算一下吹风的时候会按照风的方向影响2格的树木，这个已经写在代码注释中了~~

clear;clc;
%火灾
n = 300;                           % 定义表示森林的矩阵大小
k = 30000;                         % 迭代次数
Pground = 0.2;                     % 从着火变成空地的概率
Plight = 1e-6; Pgrowth = 2e-3;      % 定义闪电和生长的概率  
Windposb = 2e-1;                    % 定义风的效果
P2=0.7;                             %旁边有火，树着火的概率
veg=zeros(n,n)+2;                    % 初始化表示森林的矩阵
imh = image(cat(3,veg,veg,veg));   % 可视化表示森林的矩阵
Sd = zeros(n+4);                 %边界
% veg = 空地为0 着火为1 树木为2
for i=1:k
    Sd(3:n+2,3:n+2) = veg;
    % 周围有火
    sumValue = (Sd(2:n+1,3:n+2)==1)+(Sd(3:n+2,2:n+1)==1)+(Sd(3:n+2,4:n+3)==1)+(Sd(4:n+3,3:n+2)==1);
    sumWindValue = (Sd(1:n,5:n+4)==1)+(Sd(2:n+1,5:n+4)==1)+(Sd(3:n+2,5:n+4)==1)+(Sd(4:n+3,5:n+4)==1)+(Sd(5:n+4,5:n+4)==1);
    for p=1:n
        for q=1:n
            if(veg(p,q)==2 && ((sumValue(p,q)>0 && rand()<P2)||rand()<Plight||(rand()<Windposb && sumWindValue(p,q)>0)))
                %首先要是树，而且需要邻居有火，就会一定概率着火;或者被雷劈了,就会直接着火;或者是东风让东火烧到了
                veg(p,q)=1;
            elseif(veg(p,q)==1&&rand()<Pground)
                %如果是火且满足概率，则变为空地
                veg(p,q) = 0;
            elseif(veg(p,q)==0&&sumValue(p,q)==0&&rand()<Pgrowth)
                %如果是空地，且周围没有火，那么以一定概率长成树
                veg(p,q) = 2;
            end
        end
    end
    set(imh, 'cdata', cat(3,(veg==1),(veg==2),zeros(n)))
    drawnow  
end