流行病传染模拟_matlab

function covid19(~) n = 1000; % 初始人口规模 infect = 10; % 新患病个体之间的间隔 birth = .020; % 出生率 mortality = .050; % 死亡率 virulence = .020; % 传染距离 duration = 50; % 感染时间 speed = 1; % 显示更新之间的间隔 barrier = .70; % 屏障长度 randu = @(n,m) 2*rand(n,m)-1; % Uniform random in [-1 1] % Fixed parameters bw = .02; % 屏障宽度 ms = 24; % markersize epidemic = 5; % durations per epidemic stepsize = .03; tmax = 2000; fontsize = 10; fontname = get(groot,'fixedwidthfontname'); % fixed width % --------------------------------------------------------- A = []; data = []; deaths = []; figs = gcf; n0 = n; if nargin == 1 thumbnail return end [paws,stop,census] = init_gui; restart_callback while get(stop,'value') == 0 % Census count = histc(A.type,1:5); count(6) = deaths; % Update display if mod(t,speed) == 0 for k = 1:n set(A.p(k),'xdata',A.x(k),'ydata',A.y(k)) end xlabel(t) title(n) set(census,'string',num2str(count,'%3d')) drawnow limitrate end % Skip rest of loop if no more sick or pause or runaway simulation if (t > 0 && count(4) == 0) || get(paws,'value') == 1 || t >= tmax continue end if t > 0 data(:,end+1) = count; end % Introduce new virus if t <= epidemic*duration && mod(t,infect) == 0 introduce_sick end % Time step A.x = A.x + stepsize*A.u; A.y = A.y + stepsize*A.v; A.age = A.age + 1/10; t = t + 1; % Bounce off sides k = find(abs(A.x) > 1-stepsize); A.x(k) = sign(A.x(k))*(1-stepsize); A.u(k) = -A.u(k); k = find(abs(A.y) > 1-stepsize); A.y(k) = sign(A.y(k))*(1-stepsize); A.v(k) = -A.v(k); % Bounce off barrier k = find(abs(A.y) >= (1-barrier) & abs(A.x) < 2*bw); A.x(k) = 2*bw*sign(A.x(k)); A.u(k) = -A.u(k); % Possibly infectious encounters z = A.x + A.y*1i; d = triu(abs(z - z.')); [r,q] = find(0 < d & d < virulence); for k = 1:length(r) rk = A.type(r(k)); qk = A.type(q(k)); if rk == 4 && qk <= 3 A.type(q(k)) = 4; set(A.p(q(k)),'color',colors(4)); A.t(q(k)) = t; end if qk == 4 && rk <= 3 A.type(r(k)) = 4; set(A.p(r(k)),'color',colors(4)); A.t(r(k)) = t; end end % births if rand < birth*n/100 narg = findobj('tag','n+'); n_callback(narg); A.age(n) = 0; end % deaths nv = count(4)+count(5); if rand < mortality*nv/10 narg = findobj('tag','n-'); n_callback(narg,nv) end % young -> adult at 20 f = (A.type == 1) & (A.age >= 20); A.type(f) = 2; set(A.p(f),'color',colors(2)) % adult -> mature at 60 f = (A.type == 2) & (A.age >= 60); A.type(f) = 3; set(A.p(f),'color',colors(3)) A.u(f) = 0; A.v(f) = 0; % sick -> immune after duration time steps f = (A.type == 4) & (A.t < t - duration); A.type(f) = 5; set(A.p(f),'marker','o', ... 'markersize',ms/4, ... 'linewidth',1, ... 'color',colors(5)); end close(figs) % ---------------------------------------------------------------------- function restart_callback(~,~) axis([-1 1 -1 1]) axis square noticks t = 0; n0 = n; shg drawnow limitrate cla make_barrier(barrier,bw) t = 0; % time deaths = 0; data = zeros(6,0); narg = findobj('tag','n+'); A.x = randu(n,1); % initial position, [-1,1] A.y = randu(n,1); A.u = randu(n,1); % constant velocity, [-1,1] A.v = randu(n,1); A.type = 2*ones(n,1); % adults A.age = 20 + 22*randn(n,1); % 20 < age < 42 A.t = zeros(n,1); % time stamp A.p = zeros(n,1); % graphics pointer c = colors(2); for k = 1:n A.p(k) = line(A.x(k),A.y(k), ... 'linestyle','none', ... 'marker','.', ... 'color',c, ... 'markersize',ms); end end function [paws,stop,census] = init_gui if isequal(get(gcf,'name'),'covid19-1') clf else set(gcf,'name','covid19-1','numbertitle','off') end uicontrol('style','pushbutton', ... 'units','normal', ... 'position',[.86 .87 .10 .05], ... 'background','w', ... 'string','重新', ... 'fontsize',fontsize, ... 'fontweight','bold', ... 'callback',@restart_callback); paws = uicontrol('style','toggle', ... 'units','normal', ... 'position',[.86 .52 .10 .05], ... 'background','w', ... 'string','停止', ... 'fontweight','bold', ... 'fontsize',fontsize); stop = uicontrol('style','toggle', ... 'units','norm', ... 'position',[.86 .80 .10 .05], ... 'string','关闭', ... 'background','w', ... 'fontweight','bold', ... 'fontsize',fontsize); help_callback = @(~,~) doc('covid19'); uicontrol('style','push', ... 'units','norm', ... 'position',[.86 .73 .10 .05], ... 'string','帮助', ... 'background','w', ... 'fontsize',fontsize, ... 'fontweight','bold', ... 'callback',help_callback); blog_callback = @(~,~) ... web('https://blog.csdn.net/zzx2016zzx/'); uicontrol('style','push', ... 'units','norm', ... 'position',[.86 .66 .10 .05], ... 'string','博客', ... 'background','w', ... 'fontsize',fontsize, ... 'fontweight','bold', ... 'callback',blog_callback); uicontrol('style','push', ... 'units','norm', ... 'position',[.86 .59 .10 .05], ... 'string','画图', ... 'background','w', ... 'fontsize',fontsize, ... 'fontweight','bold', ... 'callback',@plot_callback); types = '年轻人 成年人 中年人 患病 免疫 死亡 '; uicontrol('style','frame', ... 'units','normal', ... 'position',[.835 .20 .155 .25], ... 'background','w') uicontrol('style','text', ... 'units','normal', ... 'position',[.84 .22 .09 .21], ... 'fontname',fontname, ... 'fontsize',fontsize, ... 'fontweight','bold', ... 'background','w', ... 'horiz','left', ... 'string',types) census = uicontrol('style','text', ... 'units','normal', ... 'position',[.938 .22 .05 .21], ... 'fontname',fontname, ... 'fontsize',fontsize, ... 'fon