當前位置：首頁 > 编程资源 > 编程问答 >内容正文

编程问答

【元胞自动机】基于matlab元胞自动机模拟SEIR传播模型【含Matlab源码 2156期】

發布時間：2023/12/20 编程问答 30 豆豆

生活随笔收集整理的這篇文章主要介紹了【元胞自动机】基于matlab元胞自动机模拟SEIR传播模型【含Matlab源码 2156期】小編覺得挺不錯的,現在分享給大家,幫大家做個參考.

一、元胞自動機簡介

1 元胞自動機發展歷程
最初的元胞自動機是由馮 · 諾依曼在 1950 年代為模擬生物細胞的自我復制而提出的. 但是并未受到學術界重視.1970 年, 劍橋大學的約翰 · 何頓 · 康威設計了一個電腦游戲 “生命游戲” 后, 元胞自動機才吸引了科學家們的注意.1983 年 S.Wolfram 發表了一系列論文. 對初等元胞機 256 種規則所產生的模型進行了深入研究, 并用熵來描述其演化行為, 將細胞自動機分為平穩型, 周期型, 混沌型和復雜型.

2 對元胞自動機的初步認識
元胞自動機（CA）是一種用來仿真局部規則和局部聯系的方法。典型的元胞自動機是定義在網格上的，每一個點上的網格代表一個元胞與一種有限的狀態。變化規則適用于每一個元胞并且同時進行。典型的變化規則，決定于元胞的狀態，以及其（ 4 或 8 ）鄰居的狀態。

3 元胞的變化規則&元胞狀態
典型的變化規則，決定于元胞的狀態，以及其（ 4 或 8 ）鄰居的狀態。

4 元胞自動機的應用
元胞自動機已被應用于物理模擬，生物模擬等領域。

5 元胞自動機的matlab編程
結合以上，我們可以理解元胞自動機仿真需要理解三點。一是元胞，在matlab中可以理解為矩陣中的一點或多點組成的方形塊，一般我們用矩陣中的一點代表一個元胞。二是變化規則，元胞的變化規則決定元胞下一刻的狀態。三是元胞的狀態，元胞的狀態是自定義的，通常是對立的狀態，比如生物的存活狀態或死亡狀態，紅燈或綠燈，該點有障礙物或者沒有障礙物等等。

6 一維元胞自動機——交通規則
定義：
6.1 元胞分布于一維線性網格上.
6.2 元胞僅具有車和空兩種狀態.

7 二維元胞自動機——生命游戲
定義：
7.1 元胞分布于二維方型網格上.
7.2 元胞僅具有生和死兩種狀態.

元胞狀態由周圍八鄰居決定.
規則：

骷髏：死亡；笑臉：生存
周圍有三個笑臉，則中間變為笑臉
少于兩個笑臉或者多于三個，中間則變死亡。

8 什么是元胞自動機
離散的系統: 元胞是定義在有限的時間和空間上的, 并且元胞的狀態是有限.
動力學系統: 元胞自動機的舉止行為具有動力學特征.
簡單與復雜: 元胞自動機用簡單規則控制相互作用的元胞模擬復雜世界.

9 構成要素

（1）元胞 (Cell)

元胞是元胞自動機基本單元：
狀態: 每一個元胞都有記憶貯存狀態的功能.
離散: 簡單情況下, 元胞只有兩種可能狀態; 較復雜情況下, 元胞具有多種狀態.
更新: 元胞的狀態都安照動力規則不斷更新.
（2）網格 (Lattice)
不同維網格

常用二維網格

（3）鄰居 (Neighborhood)

（4）邊界 (Boundary)

反射型：以自己作為邊界的狀態
吸收型：不管邊界（車開到邊界就消失）

（5）規則（狀態轉移函數）
定義：根據元胞當前狀態及其鄰居狀況確定下一時刻該元胞狀態的動力學函數, 簡單講, 就是一個狀態轉移函數.
分類：
總和型: 某元胞下時刻的狀態取決于且僅取決于它所有鄰居的當前狀態以及自身的當前狀態.
合法型: 總和型規則屬于合法型規則. 但如果把元胞自動機的規則限制為總和型, 會使元胞自動機具有局限性.
（6）森林火災

綠色：樹木；紅色：火；黑色：空地。
三種狀態循環轉化：
樹：周圍有火或者被閃電擊中就變成火。
空地：以概率p變為樹木
理性分析：紅為火；灰為空地；綠是樹

元胞三種狀態的密度和為1

火轉化為空地的密度等于空地轉換為樹的密度（新長出來的樹等于燒沒的樹）

f是閃電的概率:遠遠小于樹生成的概率；T s m a x T_{smax}T smax
?是一大群樹被火燒的時間尺度
程序實現
周期性邊界條件
購進啊

其中的數字為編號
構建鄰居矩陣

上面矩陣中的數字編號，對應原矩陣相同位置編號的上鄰居編號，一一對應
同樣道理：

（7）交通概念
車距和密度

流量方程

守恒方程

時空軌跡（橫軸是空間縱軸為時間）

紅線橫線與藍色交點表示每個時間車的位置。
如果是豎線則表示車子在該位置對應的時間

宏觀連續模型：

最常用的規則：

紅色條表示速度是滿的。

1 加速規則：不能超過v m a x （ 2 格 / s ） v_{max}（2格/s）v
max（2格/s）
2 防止碰撞：不能超過車距

理論分析：

結果分析: 密度與流量

第一個圖：橫坐標是歸一化后的密度，縱坐標是車流量。第二個圖：理論值與CA的結果

結果分析: 時空軌跡

中間的深色區域是交通堵塞的區域。

二、部分源代碼

% cellular automata
clc;clear;close all;
%%%%%% initialization of parameters
N = 100; % a half length of region
Dn = 0; % the number of death
rand(‘state’,121); % state generating random data
S = zeros(2N+1); % state matrix
r = rand(2N+1);
S(r>0.4) = 4; % uniform distribution of normal people before simulation

vc(1,:) = [1,1,0]; % yellow color state for death or unpeopled
vc(2,:) = [0,1,1]; % black color state for curative
vc(3,:) = [0,0,1]; % blue color state for delitescent
vc(4,:) = [1,0,0]; % red color state for catch an illness
vc(5,:) = [0,1,0]; % green color state for normal people
vc(6,:) = [1,0,1]; % black color state for normal people

% Sc = ones(2N+1,2N+1,3); % color state matrix
Sc = mark_color(S,vc);
%%%%%% initialization of Figure
set(gcf,‘Position’,[200,60,600,700]);
axes(‘Position’,[0.1,0.4,0.6,0.5]); % mapping
I1 = imshow(Sc);
th = title(‘The spread days of the disease: 1’,‘FontSize’,14);
axes(‘Position’,[0.78,0.4,0.2,0.5]); hold on;
ylim([-0.1,1.05]); xlim([-0.1,1]);
% box on;
axis off
rectangle(‘Position’,[0,0.35,0.3,0.10],‘FaceColor’,vc(1,:));
text(0.4,0.4,{‘D’});
rectangle(‘Position’,[0,0.05,0.3,0.10],‘FaceColor’,vc(2,:));
text(0.4,0.1,{‘R’});
rectangle(‘Position’,[0,0.65,0.3,0.10],‘FaceColor’,vc(3,:));
text(0.4,0.70,{‘E’});
rectangle(‘Position’,[0,0.5,0.3,0.10],‘FaceColor’,vc(4,:));
text(0.4,0.55,{‘I’});
rectangle(‘Position’,[0,0.8,0.3,0.10],‘FaceColor’,vc(5,:));
text(0.4,0.85,{‘S’});
rectangle(‘Position’,[0,0.2,0.3,0.10],‘FaceColor’,vc(6,:));
text(0.4,0.25,{‘Q’});
set(gca,‘Xtick’,[],‘Ytick’,[]);
axes(‘Position’,[0.1,0.08,0.6,0.2]); % stats
Bh(1) = bar(1,sum(sum(S1)),‘FaceColor’,vc(2,:)); hold on;
Bh(2) = bar(2,sum(sum(S2)),‘FaceColor’,vc(3,:));
Bh(3) = bar(3,sum(sum(S==3)),‘FaceColor’,vc(4,:));
ylabel(‘number of people’);
xlim([0,4]);set(gca,‘Xtick’,[0:4])
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%%% iterative calculation
S(N+N/2,N+1) = 3; % center cell is catched an illness
S(N+(N-10)/2,N+1) = 3; % center cell is catched an illness
S(N-N/2,N+1) = 3; % center cell is catched an illness
S(N+1,N+1) = 3; % center cell is catched an illness
S(N+1,N+N/2) = 3; % center cell is catched an illness
S(N+1,N-N/2) = 3; % center cell is catched an illness
SL(S3) = 0; % time for the state of catching an illness
SGL = zeros(2*N+1); % 隔離時間
Sc = mark_color(S,vc);
set(I1,‘CData’,Sc);
Yao = 60;
cure=0;
ge=0;
ill=0;
for ktime = 1:200; % unit: day
% 潛伏者 5-10天有 80%的概率變成發病者，也就是狀態2有80%概率變成狀態3 ，
% 有20%的概率變成狀態1
s2 = find(S2);
for k = 1:length(s2);
if SL(s2(k))>4.5 & SL(s2(k))<10.5; % 5-10天
if rand<0.8;
S(s2(k)) = 3; % 2–>3, 80%
ill=ill+1;
SL(s2(k)) = 1; % time updata
else
S(s2(k)) = 1; % 2–>1, 20%
cure=cure+1;%治好人數更新
SL(s2(k)) = 1; % time updata
end
else
SL(s2(k)) = SL(s2(k))+1; % updating data
end
end
s3 = find(S3);
%%隔離死亡
s5 = find(S5);
for k = 1:length(s5);
if SGL(s5(k))>0.5 & SGL(s5(k))<10.5; % 5-10天
if rand<0.4;
S(s5(k)) = 0; % 2–>0, 80%
Dn = Dn +1;
% ill=ill+1;
SGL(s5(k)) = 0; % time updata
else
% S(s5(k)) = 1; % 2–>1, 20%
SGL(s5(k)) = SGL(s5(k))+1; % time updata
end
else
if rand<0.5;
S(s5(k)) = 0; % 2–>3, 80%
Dn = Dn +1;
% ill=ill+1;
SGL(s5(k)) = 0; % time updata
else
% S(s5(k)) = 1; % 2–>1, 20%
SGL(s5(k)) = SGL(s5(k))+1; % time updata
end
% SL(s2(k)) = SL(s2(k))+1; % updating data
end
end
%%%隔離狀態
for k = 1:length(s3);
if ktime>60;
if rand<min(SL(s3(k))*0.05+0.05,0.4);
S(s3(k)) = 5;
ge = ge+1;
if SL(s3(k))<10.5 % 如果在患病 10天內發現
if rand<0.8 & rand<min(1,Yao/ge)
S(s3(k)) = 1; % 5–>1. % 50%被治好, 藥物的期望0.9，也就是有90%可以被治好
cure=cure+1; %治好人數更新
SGL(s3(k)) = 0; % 隔離時間清空
else
SGL(s3(k)) = SGL(s3(k))+1; % 隔離時間更新
end
elseif SL(s3(k))>10.5 % 患病10天到20天
if rand<0.2 & rand<rand<min(1,Yao/ge);
S(s3(k)) = 1; % 5–>1. 有20%被治好，藥物的期望0.9，也就是有90%可以被治好
cure=cure+1;%治好人數更新
SGL(s3(k)) = 0; % 隔離時間清空
else
SGL(s3(k)) = SGL(s3(k))+1; % 隔離時間更新
end
elseif SGL(s3(k))<10.5 % 在隔離5天內
if rand<0.2 & rand<rand<min(1,Yao/ge)
S(s3(k)) = 1; % 5–>1. 有20%被治好，藥物的期望0.9，也就是有90%可以被治好
cure=cure+1;%治好人數更新
SGL(s3(k)) = 0; % 隔離時間清空
else
SGL(s3(k)) = SGL(s3(k))+1; % 隔離時間更新
end
elseif SGL(s3(k))>10.5 % 在隔離5天以上
if rand<0.5 & rand<rand<min(1,Yao/ge)
S(s3(k)) = 1; % 5–>1. 有50%被治好，藥物的期望0.9，也就是有90%可以被治好
cure=cure+1;%治好人數更新
SGL(s3(k)) = 0; % 隔離時間清空
else
% S(s3(k)) = 0; % 5–>1. 有50%未被治好，藥物沒有作用而死亡
% SGL(s3(k)) = 0; % 隔離時間清空
SGL(s3(k)) = SGL(s3(k))+1; % 隔離時間更新
end
end
else
SL(s3(k)) = SL(s3(k))+1; % 未隔離者發表時間更新
end
end