1.元胞自动机定义和适用范围
不同于一般的动力学模型,元胞自动机不是由严格定义的物理方程或函数确定,而是用一系列模型构造的规则构成。凡是满足这些规则的模型都可以算作是元胞自动机模型。因此,元胞自动机是一类模型的总称,或者说是一个方法框架。其特点是时间、空间、状态都离散,每个变量只取有限多个状态,且其状态改变的规则在时间和空间上都是局部的。
元胞自动机——应用于森林火灾和传染病场景
最近接触了元胞自动机模型,做了一些资料搜查,并进行学习,推荐这篇文章澳洲变燠洲,考拉成烤拉!澳大利亚山火为什么难以控制?
以下对所学进行记录。
森林火灾元胞自动机原理
在元胞自动机模型中,空间被离散成网格,每一个网格被称为元胞。森林火灾元胞有三种状态:树,火(正在燃烧的树)和空(空地)状态。元胞下一时刻状态的更新规则如下:
树变火:一棵树,其上下左右若有一个状态为火,下一刻就会变成火。或者一棵树遇上闪电,下一刻就会变成火。由于遇上闪电着火的概率Plight很小。
火变空:火在下一时刻会变成空。
空变树:空地下一时刻会以一个很小的概率Pgrowth长出新树。
改进模型会考虑树的对角位置有没有着火。或者会考虑风向(比如吹西风(火从东吹向西),火的西边着火的机率会变大(顺风),火的东边着火的几率变小(逆风)),这里盗个图:
图a是基础元胞自动机,图b是考虑对角的元胞自动机,图c是吹西风的元胞自动机。
本人这里对基础元胞自动机,考虑对角情况的元胞自动机,考虑吹西风的元胞自动机三种模型进行仿真,仿真结果如下:
P=[]; for m=1:10 clear D T fire_time_lightning fire_time_itself aspect tdata Index; %% Orginal %the first 3 dimension is RGB,R is the fire,G is the tree. %Black is the meaning of no tree. global n D T Y fire_time_lightning fire_time_itself fire_time_demend pull aspect count_1 tdata Pull_times n=500; % the length of the forest matrix D=zeros(n); T=zeros(n); Y=zeros(n,n,3); % draw the picture by matrix Y(RGB) fire_time_lightning=0; fire_time_itself=0; fire_time_demend=0; times=1; pull=1/times; % the rate of pull the fire out aspect=ceil(rand(1)*4); % 1 is right,2 is up,3 is left and 4 is down. count_1=0; tdata=[]; % the day by each fire happened. Pull_times=0; f1=1/1000; % f1 is the probability in ceil when it being struck by lightning. f2=1/500; % f2 is the probability in ceil when it being fired itself. Z=Terrain(); [scale_b,S]=Forest(Z); Tem=S; Yi=imshow(Y); set(gcf,'DoubleBuffer','on'); % set up the double cache to prevent the flash in palying animation constantly t=0; tp=title(['T = ',num2str(t)]); %ap=title(['aspect = ',num2str(aspect)]); %while 1 % t=t+1; for t=1:2400 %% Fire in the early time if rem(t,50)==0&&t<1000 Fire_Demand(S); end %% Lightning if rand<f1 %Is there happen sth with lightning? OriginFireLightning(S); set(Yi,'CData',Y); end %% Fire itself if t>10 OriginFireCritical(S,t,f2); set(Yi,'CData',Y); end %% FireRule1 [a,b]=FireRule1(S); %% FireRule2 S=FireRule2(S,a,b); set(Yi,'CData',Y); set(tp,'string',['T = ',num2str(fix(t/(24))),'D = ',num2str(t),'h ']) %set(ap,'string',['aspect = ',num2str(aspect)]) pause(2e-6) end scale_n=sum(sum(Y(:,:,2))); fire_count=fire_time_itself+fire_time_lightning; Lost_area=scale_b-scale_n; Lost_rate_all=Lost_area/scale_b; Lost_rate_average=Lost_rate_all/(fire_count); Tem=Tem-S; Tem=Tem.*Tem; %Lost_Value=sum(sum(Tem)); %N=[pull Lost_rate_all Lost_rate_average fire_count] P=[P;fire_time_demend Pull_times Lost_rate_all fire_time_itself fire_time_lightning] %figure (2) %plot(tdata(:,1),tdata(:,2)); end
版本:2014a
需要完整代码或代写加QQ 1564658423