小组成员:张逸凡,邓国章,李衡逸
设计思想: 如同之前在方案设计当中提及的,本实验设计集中在于还原原本课堂中进行的关于瘟疫传播模拟的同时,模拟现实当中社区间固定的人员流通会导致何种不同的瘟疫传播,同时探究在这样的情况前提之下实行各种阻隔瘟疫传播的措施对瘟疫传播影响有何种变化。 为了满足以上的需求,我们首先需要对原本处于一个大区域内的人员进行分隔,在此实验当中我们将区域分隔为了九个小型社区,在这之后为了模拟实现日常社区生活当中人们通常总是遵循着一定的行动规律,我们设定社区间存在着固定的人员流通模式,例如1号社区的人员只会在1号与3,7号社区之间行动。 我们同样参考了著名的SIR模型,在此模型内将人群分为 可感染,感染,已康复三种状态,再来拟合不同状态间的转移概率。 SIR模型的建立基于以下三个假设: ⑴不考虑人口的出生、死亡、流动等种群动力因素。人口始终保持一个常数,即N(t)≡K。 ⑵一个病人一旦与易感者接触就必然具有一定的传染力。假设 t 时刻单位时间内,一个病人能传染的易感者数目与此环境内易感者总数s(t)成正比,比例系数为β,从而在t时刻单位时间内被所有病人传染的人数为βs(t)i(t)。 ⑶ t 时刻,单位时间内从染病者中移出的人数与病人数量成正比,比例系数为γ,单位时间内移出者的数量为γi(t)。 但是应注意到,模型对人群的分类不够细致,没有明确考虑隔离的因素。而现实中对疑似病人的隔离是控制疫情传播的有效手段。模型没有引入反馈机制,在预测过程中,单纯依据已有数据预测未来较长一段时间的数据,必然会使准确度降低。此外,微分方程组求解较为困难,且对初值比较敏感,这对模型的稳健性是一个很大的影响。 遵循以上三点虽然会带来一定的局限性,但在进行改进,;如导入反馈机制,加入我们预定讨论的固定行为模式等条件之后能更好的控制实验的变量,得到更为稳定的结果以供观察分析。 在这之后我们通过调整各项参数来实现对之前设定的各个探究项目的实验,如封城,减少人员流动,降低传染率等等,并观察结果。
数学思想: 运用抽象和等效替代的思想,用一定区域内自由运动的,类似气体分子的球体的运动抽象替代一定区域内人群的运动和接触。 数学抽象是在探究数学问题时常用的一种数学方法,基本上可划分为四种类型: 1.弱抽象.即从原型中选取某一特征(侧面)加以抽象,使原型内涵减少,结构变弱,外延扩张,获得比原结构更广的结构;2.强抽象.即通过在原型中引人新特征,使原型内涵增加,结构变强,外延收缩,获得比原结构内容更丰富的结构;3.构象化抽象.即根据数学发展的逻辑上的需要,构想出不能由现实原型直接抽取的、完全理想化的数学对象,作为一种新元素添加到某种数学结构系统中去,使之具有完备性,即运算在此结构系统中畅行无阻.如在实变函数论中,勒贝格可积函数与平方可积函数概念的引人; 4.公理化抽象.即根据数学发展的需要,构想出完全理想化的新的公理(或基本法则),以排除数学悖论,使整个数学理论体系恢复和谐统一非欧几何学平行公理、非阿基米德公理等都是公理化抽象的产物.这种抽象思维的法则可称为:‘.公理更新和谐化法则”;此处我们主要使用的是弱抽象原理,将人员与社区从具体的生物与地理概念抽象为单位与简单的平面区域。 等效替代法实际上主要是用于物理范围内的科学方法,是在保证某种效果(特性和关系)相同的前提下,将实际的、陌生的、复杂的物理问题和物理过程用等效的、简单的、易于研究的物理问题和物理过程代替来研究和处理的方法。但我们认为在这里这是对其数学概念意义的一种创新体现,单位与人员,区域与社区,以及传染可能性与传染概率,日常人员行动模式与区域间的人员固定交换都体现着等效替代法的主要思想。 与上述思想类似的理想模型也是我们应用了的方法之一,通过建立模型来揭示原型的形态、特征和本质的方法称为理想模型法。把复杂的问题简单化,摒弃次要因素,抓住主要因素,对实际问题进行理想化处理。有时为了更加形象地描述所要的现象、问题,还需要引入一些模型。 而根据大数原理,即随机现象的大量重复中往往出现几乎必然的规律。此法则的意义是:风险单位数量愈多,实际损失的结果会愈接近从无限单位数量得出的预期损失可能的结果。因此我们的模拟应运行多次,消除偶然性的影响,使得结果更具有代表性(演示视频是最显著的一次运行结果)。
技术方案(源代码): Mover[] movers = new Mover[450]; float InfectChance=0.2; int DteleCool=10; //任意人移动到其他区域的间隔。越小交通频率越大 int teleCool=DteleCool; float QAwidth=720,QAheight=720; // int I=1,Su=449,R=0; int[] Icount=new int[3600]; int[] Sucount=new int[3600]; int[] Rcount=new int[3600]; //用这三个数组存储按时间顺序的三种人的数量 PVector S=new PVector(0,0),G=new PVector(0,0); //用于显示移动轨迹 void settings() { size(1080, 720); smooth(); }
void setup() { hint(ENABLE_STROKE_PURE); background(0); for (int i = 0; i < movers.length; i++) { movers[i] = new Mover(); } movers[0].state=1; }
void draw() { noStroke(); fill(0); rect(0, 0, QAwidth, QAheight); stroke(255); strokeWeight(30); line(0,QAheight/3,QAwidth,QAheight/3); line(0,2QAheight/3,QAwidth,2QAheight/3); line(QAwidth/3,0,QAwidth/3,QAheight); line(2QAwidth/3,0,2QAwidth/3,QAheight); //分割网格 noStroke(); fill(0); rect(720, 0, 1080-QAwidth, QAheight); for (int i = 0; i< movers.length; i++) { for (int j = 0; j< movers.length; j++) {
if(movers[i].pos.dist(movers[j].pos)<20) { PVector a = movers[i].SocialDis(movers[j]); movers[i].applyForce(a.mult(2)); } else { PVector a = movers[i].SocialDis(movers[j]); movers[i].applyForce(a); } float c=random(0,1); if(c<InfectChance) if(movers[i].pos.dist(movers[j].pos)<15) { if(movers[i].state==1&&movers[j].state==0) { movers[j].state=1; Su--; I++; } } } if(movers[i].state==1)movers[i].lifetime--; if(movers[i].lifetime<0){ if(movers[i].state==1) { I--; R++; } movers[i].state=2; } movers[i].applyForce(movers[i].dir); movers[i].update(); movers[i].show(); movers[i].checkEdges();} fetch(); fetch(); fetch(); fetch(); fetch(); textSize(32); fill(0, 102, 153); text("Susceptible: "+Su, 730, 30); text("Infectious: "+I, 730, 60); text("Removed: "+R, 730, 90); float t=second(); if(t){ } DrawGraph(730,120,t); if(I0) text(“Pandemic is Over”, 730, 490); ;
}
void fetch()//移动时调用 { int t =millis(); if(teleCool0) { teleCool=DteleCool; int a1x,a1y; int a2x,a2y; a1x=movers[0].areax; a1y=movers[0].areay; a2x=int(random(0,3)); a2y=int(random(0,3)); while(a2xa1x&&a1y==a2y) { a2x=int(random(0,3)); a2y=int(random(0,3)); }; int j = int(random(movers.length)); //float detectP=random(1); //if(detectP<0.9) //这一段来决定是否限制感染者出行和限制概率 j=int(random(movers.length)); S=movers[j].pos; movers[j].teleport(a1x,a1y,a2x,a2y); strokeWeight(10); stroke(244,255,30,100*teleCool/DteleCool); line(a1x,a1y,a2x,a2y); noStroke(); G.x=movers[j].pos.x; G.y=movers[j].pos.y;
} else { teleCool–; print(teleCool+"\n"); strokeWeight(10); stroke(244,255,30,100*teleCool/DteleCool); line(S.x,S.y,G.x,G.y); noStroke(); } }
void DrawGraph(float x,float y,float t)//绘制统计图表 { push(); Sucount[second()]=Su; Icount[second()]=I; Rcount[second()]=R; translate(x,y); for(int i=0;i<t+1;i++){ fill(127,207,255); rect(i300/t,0,600/t,200Sucount[i]/movers.length); fill(125,127,127); rect(i300/t,200Sucount[i]/movers.length,600/t,200Rcount[i]/movers.length); fill(255,27,27); rect(i300/t,(200Sucount[i]/movers.length)+(200Rcount[i]/movers.length),600/t,200*Icount[i]/movers.length); } pop(); }
class Mover {
PVector pos; PVector vel; PVector acc; float maxSpeed; float awareness; int state=0; float lifetime=200; float r; PVector dir = PVector.random2D(); int areax=int(random(0,3)),areay=int(random(0,3)); // int areax=1,areay=1;
Mover() { pos = new PVector(random(areax*QAwidth/3+20,(areax+1)QAwidth/3-20), random(areayQAheight/3+20,(areay+1)*QAheight/3-20)); vel = new PVector(0, 0); acc = new PVector (0, 0); maxSpeed = 10; r = 10; }
void update() { pos.add(vel); vel.add(acc); vel.limit(1.5); acc.mult(0); }
void show() { noStroke(); if(state0) fill(127,207,255); if(state1) fill(255,27,27); if(state==2) fill(125,127,127); ellipse(pos.x, pos.y, r, r); }
void checkEdges() {
if (pos.x > (areax+1)*QAwidth/3-21) { vel.x=-vel.x; dir.x=-dir.x; } else if (pos.x < areax*QAwidth/3+19) { vel.x=-vel.x; dir.x=-dir.x; } if (pos.y > (areay+1)*QAheight/3-21) { vel.y=-vel.y; dir.y=-dir.y; } else if (pos.y < areay*QAheight/3+19) { vel.y=-vel.y; dir.y=-dir.y; }}
void applyForce(PVector force) { //m*a = f j假设每个小球的质量为单位1 acc.add(force.mult(0.01));}
PVector SocialDis(Mover m) { PVector force = PVector.sub(m.pos, pos); float distance = force.mag(); distance = constrain(distance, 0.001,QAwidth); force.normalize(); force.mult(-150/((distance)*(distance)));//假设每个小球的质量为单位1,按照万有引力的公式计算小球之间的斥力 return force; } void teleport(int a1x,int a1y,int a2x,int a2y) { PVector workpos = new PVector(random(a2x*QAwidth/3+20,(a2x+1)*QAwidth/3-20), random(a2y*QAheight/3+20,(a2y+1)*QAheight/3-20)); PVector middle=PVector.lerp(pos,workpos,0.2); areax=a2x; areay=a2y; pos=workpos; fill(244,255,30); ellipse(pos.x, pos.y, 2*r, 2*r); ellipse(workpos.x, workpos.y, 2*r, 2*r); }}