在我大四下学期时，毕业设计让我提前接触了美硕研究方向的内容，其中模拟流体的格子-玻尔兹曼方法，利用了网格来模拟流体微粒的方法，让我想到了当初学习数模之时遇到的元胞自动机。
当下很多数模的教材，当讲到一些科学方法时，往往都会直接给出理论计算部分，这样往往容易产生劝退之意。喵姐才疏学浅，但也想以自己的表达方式，以前人的代码基础来谈谈元胞自动机的应用方法。
谈及元胞自动机（Cellular Automata），同学们第一印象或许就是许许多多的小格子，然后这些小格子们以某种函数关系相互相爱相杀，将图像推演下去。
我这么说也不知道有没有说中大家，但这么说其实也差不多，下面我会通过一些例子，简单介绍下这些格子们（cells）是如何“相爱相杀”的。
长文预警
多年以前，世间曾有一片森林，后来我对它用了二向箔，这片森林就成了二维的存在。
在MATLAB里这片森林就是个二维矩阵，1000×1000的位置，每个位置就是一个cell。
并不是说，3D的就不能被元胞自动机模拟了，3D的理论上也是可以的。
不过CA一般多用于二维层面上的模拟，并被认为是离散的动力学模型。并且，网格也不一定非得是四四方方的方块形式，也可以是三角形，多边形等。
有天闲暇，我用MATLAB打开这片森林，想看看cell的数值，里面全是0和2。那么0代表空地，2代表一棵树。
那么假如某棵树着火了，那么在cell里的数值，就应当是1，代表这棵树着火了，并且在下一个时刻，这个cell的数值将变成0，代表树烧没了。
那么凭什么让这棵树附近的树儿们引火上身呢？这里我先介绍一下什么是VonNeumann Neighborhood和Moore Neighborhood。如果是应用VonNeumann邻居，那么则只考虑这棵树前后左右的状态。
若是Moore邻居，则考虑周边八个格子的影响。
你仍然可以考虑周边多个格子的影响，整一些什么Sister Miao Neighborhood什么的，都可以。
在森林火灾模型中，规则就是只要4个邻居cells有一个是1，并且自己是2，那么下个时间节点自己将变成1，再下一个时间节点自己变成0。
森林火灾模型元胞的规则，做到这里，基本上这就完事了，怎么样，是不是很简单？
那么回到我们的绿绿森林，网上现有的源代码里面还设定了0.000005的雷击与0.01的自然生长，即2有0.000005的概率变成1，0有0.01的概率变成2。并且还巧妙地设置了边界连接，那我们就暂且考虑一下这些因素吧，让喵姐在下文详细拆解代码的部分。
Line：1-11
clf
clear all
n=100; %注意line6，生成n×n的零矩阵，就是所谓的lattice，可供我们恣意纵横的棋盘。
Plightning = .0000005; %产生雷击的概率
Pgrowth = .01; %空地长出树木的概率
z=zeros(n,n); %在line3处设置了n，如果你想自定义“棋盘”的话，大可以删除line3，然后直接改这里，这样的话底下的cell规则也要修正一下。
veg=z; %复制一层棋盘，留z作底板
sum=z; %再复制一层z，用作计算该cell周边着火树木的数量
imh = image(cat(3,z,veg,z)); %cat(3,z,veg,z)会生成n×n×3的矩阵，由image生成图像。
axis equal
axis tight
Image（C）会让数组C产生图像，可以用下面这些代码直观地了解一下。
C = [0 2 4 6; 8 10 12 14; 16 18 20 22];
image(C)
colorbar
Line：13
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)
这行代码计算了每个cell的上下左右四个格子中1的数目。
要知道，sum也是一个n×n的矩阵，veg(1:n,[n 1:n-1])中，1:n表示1-100。
[n 1:n-1]可能萌新同学不太理解，且容喵姐在此絮叨一下，[n]很简单，就是1×1的矩阵，[1 2 3]呢？是1×3的矩阵，[1:4]则为[1 2 3 4]，是个1×4的矩阵，那么[1 2 3 1:4]是什么呢？
答案是[1 2 3 1 2 3 4]。
当你看到这里，想必就明白了[n 1:n-1]是一个1×n的矩阵。
而veg(1:n,[n 1:n-1])则是n×n的矩阵。
这行代码萌新理解起来可能比较难，但理解后应当就和多年前大二时期，那个黄昏见晚的下午，坐在大学寝室椅子上的喵姐，之前还一脸懵逼，结果一下子恍然大悟吧。
以这样的形式来计算sum，是考虑到了边界的连接，即火烧到棋盘右侧边界，便会从左侧边界继续烧过来。上下边界同理。
那么同学们准备好恍然大悟了吗？
喵姐真想看看你们恍然大悟的样子，就觉得心有灵犀，大家都好聪明，一点就透。
首先，看矩阵[1:n]，这个一目了然，对吧。
[1:n]=[1:n-1 n]，这个等式对不对？对！
那假若我这样呢？[n 1:n-1]，有什么区别？是不是把n提到左边来了，这是什么意思？
就是veg(1:n,[n 1:n-1])把每行最右边边界那个cell的值提到了最左边第一个，1：n的意思就是第1至n行。
那第二个cell的值呢？答案是来源于左边一个格子。
以此类推。
那么veg(1:n,[n 1:n-1])的意思就很明朗了，该cell值变成左侧一格cell值，注意在这里，是变成左侧一格值，不是加上这个值。
(veg(1:n,[n 1:n-1])==1)则变成了每个cell，若左侧一格cell值为1，则该cell变成1，不是则为0。
好的，理解了这一点，再看(veg(1:n,[2:n 1])==1)，这是把最左侧cell值提到了最右侧去，然后如果是1就保留， 即每个cell右侧一格cell值提到自己身上，如果是1就保留，其他值就归零。
(veg([n 1:n-1], 1:n)==1)，当1：n在括号后侧时，则表示1-n列，此时表示底部边界最后一个cell值提到最顶上第一个，其他每个cell得到头顶上一个cell的值，并且如果这些值为1则保留，为其他值则变成0。
直观一点，可以通过这些简短的代码测试一下：
a=[1 2;
1 2;
1 2;
1 2;
1 2;
1 2;
3 3;
1 2;
1 2;
1 2];
veg=a;
b=veg([10 1:9], 1:2)
(veg([2:n 1],1:n)==1) 表达得到下面一格cell的值，为1则保留，其他值则变成0。
那么
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) ; 则表示4个n×n矩阵相加，这样就有了上下左右四个cell中着火的数目数量。
如果你删除了line3，在line6处设置了自己想要的尺寸，如x（行，即上下方向长度）×y(列，即左右方向宽度)，此处也要改成
sum = (veg(1:x,[y 1:y-1])==1) + (veg(1:x,[2:y 1])==1) + (veg([x 1:x-1], 1:y)==1) + (veg([2:x 1],1:y)==1) ;
Line:14
veg = 2*(veg==2)-((veg==2)&(sum>0|(rand(n,n)<Plightning)))+2*((veg==0) & rand(n,n)<Pgrowth)
先看2*((veg==0) & rand(n,n)<Pgrowth)。
此处，rand可以生成随机数，rand(n,n)则为表示n×n里的每个cell生成随机数，rand(n,n)<Pgrowth，判断该cell生成的随机数是不是比Pgrowth大，假若Pgrowth为0.01，即1%的概率这儿会凭空长树，那么rand(n,n)<0.01，cell中的随机数若是小于0.01的话，那会怎么办呢？该cell会变成1，若是大于0.01，该cell变成0。形象一点的话，如下图。
那么肯定有同学要问了，不是说好了要凭空长出来树吗？要是原来有树了咋办？
所以说，除了rand(n,n)
((veg==2)&(sum>0|(rand(n,n)<Plightning)))
且看veg==2，很好理解，产生个n×n的矩阵，里面树木地区cells，就是里面数值为2的cells变成1，其他cells数值为0。
再看(sum>0|(rand(n,n)0表示该cell周边有火，(rand(n,n)
要是雷劈空地那还好说，倘若是雷劈到树了呢？
没错，((veg==2)&(sum>0|(rand(n,n)<Plightning)))就表示此处cell里面有树，而且他这个时间点它得火。
好了，既然知道了哪些树宝宝要火，下个时间点就得让她们火（即2变成1），如下代码即可实现。
2*(veg==2)-((veg==2)&(sum>0|(rand(n,n)<Plightning)))
多了2*(veg==2)，这会得到整片森林的树，则为n×n的0和2矩阵，即有树处为2，其他处为0。在上面我们已经得到了要火的树宝宝们位置，且在矩阵里以1的形式被标记好了，那么减去这个部分，整片森林中要火的树宝宝瞬间燃烧，即2-1=1。怎么样？元胞自动机就是如同2-1=1一样简单。
考虑到有烧有生，将烧毁部分与生长部分加起来，则得到每一轮演变的完整过程。
veg = 2*(veg==2)-((veg==2)&(sum>0|(rand(n,n)<Plightning)))+2*((veg==0) & rand(n,n)<Pgrowth) ;
我们是不是忘了cell值为1的情况？不，没有，在2*(veg==2)时，上一轮cell为1的树宝宝们直接被0了。
Line:12-17
for i=1:30000
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) ;
set(imh, 'cdata', cat(3,(veg==1),(veg==2),z) ) %这里cdata的意思是生成图像的同时，不替换原有图像。
drawnow
End
到此，本森林火灾模型已经建立完毕。
运行结果：
树木生长时，它看起来就像本公众号那可爱的二维码，迫不及待地等着新同学们来扫码关注一样，不是吗？
改变“棋盘”尺寸的模拟结果：
回到开头，那个被二向箔打击的森林，俯瞰如下图。
新建一个m文件，把bmp图像放到同一文件夹中。
line: 1-4
forest=imread('WoYuanLiangNi.bmp');%导入图像数据至MATLAB
thresh = graythresh(forest);%自动确定二值化阈值；
forgiveyou = im2bw(forest,thresh);%对图像自动二值化
imshow(forgiveyou)%看一下二值化后长啥样
点开forgiveyou数据，发现空白处为1，黑色处为0，所以我要做的就是让空白处全部变成0，黑色处全部变成2。
line1-8:
clc
clear
forest=imread('WoYuanLiangNi.bmp');
thresh = graythresh(forest);
forgiveyou = im2bw(forest,thresh);
%imshow(forgiveyou)
Iforgiveyou = zeros(400,500)+ 2*(forgiveyou==0);
imshow(Iforgiveyou)
事后我点进去Iforgiveyou里面看了下，的确是变成了0和2的矩阵，400×500，是因为我bmp图片的像素就是400×500。
这里我们去掉雷击与自然生长，仅仅在图像中间设定一个燃烧的树。
line:1-21
clc
clear
forest=imread('WoYuanLiangNi.bmp');
thresh = graythresh(forest);
forgiveyou = im2bw(forest,thresh);
%imshow(forgiveyou)
Iforgiveyou = zeros(400,500)+ 2*(forgiveyou==0);
imshow(Iforgiveyou)
z=zeros(400,500);
x=400;
y=500;
sum=z;
veg=Iforgiveyou;
veg(225,445)=1; %这儿~设定一个燃烧的树
imh = image(cat(3,z,veg,z));
for i=1:1000
sum = (veg(1:x,[y 1:y-1])==1) + (veg(1:x,[2:y 1])==1) + (veg([x 1:x-1], 1:y)==1) + (veg([2:x 1],1:y)==1) ;
veg = 2*(veg==2)-((veg==2)&(sum>0)) ;
set(imh, 'cdata', cat(3,(veg==1),(veg==2),z) )
drawnow
end
模拟结果：（并没有理想地烧完）
到这里，森林火灾模型的元胞自动机入门已经结束了，希望各位同学能够有所启发。在设置自己想要的元胞自动机时，需要考虑到各个格子之间的相互影响，写好局部的规则，这往往是很烧头脑的。
如今我们在大学学习或者步入社会，不再像初高中时代填鸭式获得知识。在埃及军区没网的夜晚，我时常抬头看向深邃的夜空，星星清晰明亮，总会让我陷入沉思，感叹前几个世纪开普勒，卡文迪许，牛顿等人，同样是看向这片夜空，他们天才般地得出了行星运动，万有引力等计算公式，被记载于物理课本中，而我却只能看到繁星点点，四处是人间理想，当我们惊叹于前人的惊人的创造时，会发现人类科学的历史就是被这些如同闪耀群星般的伟人推动的。科学历史星河上牛顿的万有引力，爱因斯坦的E=mc的平方耀眼无比，亦或是点点流星如元胞自动机，和其他FEM，FVM，FBM等科学方法，偶然间自己也会莫名有所冲动，想发现一个进步的科学方法，在人类科学的历史星河上，泛出一点属于自己的光芒。
涵盖大学的各个方面，数模入门到轻松拿奖
萌新如何成长成为学霸，升学是保研还是出国
英语该怎么学，其他竞赛怎么办，资料该怎么找...
网友评论