-
Notifications
You must be signed in to change notification settings - Fork 0
/
RIOTS_contlim_spatiallydependent.m
89 lines (70 loc) · 2.17 KB
/
RIOTS_contlim_spatiallydependent.m
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
function RIOTS_contlim_spatiallydependent(kmr,kmb,kr,kd,tend)
%% Spatially-dependent continuum limit considered in Section 4.2.2
% Use this code and RIOTS_ABM_spatiallydependent.m to reproduce Figures 12, 13
% kmr = rioter motility rate
% kmb = bystander motility rate
% kr = recruitment rate
% kd = defection rate
% tend = total run time
%% Growth rate Equations
% Writes the growth rates for the continuum limit system
% u(1) is (mean) rioter occupancy
% u(2) is (mean) bystander occupancy
Dr=kmr/4; % Effective rioter diffusivity
Db=kmb/4; % Effective bystander diffusivity
m=0;
function [c,f,s]=pdefun(x,t,u,dudx)
c=[1;1];
f=[Dr*(1-u(2))*dudx(1)+Dr*u(1)*dudx(2);Db*(1-u(1))*dudx(2)+Db*u(2)*dudx(1)];
s=[kr*u(2)*(1-(1-u(1))^4)-kd*u(1)*(1-(1-u(2))^4);-kr*u(2)*(1-(1-u(1))^4)+kd*u(1)*(1-(1-u(2))^4)];
end
%% Initial Conditions (comment out conditions not needed)
% Initial conditions with a rioter strip and a bystander strip
function u0=icfun(x)
if x>=80 && x<=90
u0=[1;0];
elseif x>=110 && x<=120
u0=[0;1];
else
u0=[0;0];
end
end
% Initial conditions with a central rioter strip and two bystander strips
% function u0=icfun(x)
% if x>=60 && x<=80||x>=120 && x<=140
% u0=[0;1];
% elseif x>=90 && x<=110
% u0=[1;0];
% else
% u0=[0;0];
% end
% end
%% Boundary Conditions
% Neumann (zero flux) boundary conditions
function [pl,ql,pr,qr]=bcfun(xl,ul,xr,ur,t)
pl=[0;0];
ql=[1;1];
pr=[0;0];
qr=[1;1];
end
%% Spatial and time domains
% Defines 1d x-space 0<=x<=200 and time to run until 'tend'
xval=linspace(0,200,1000);
tval=linspace(0,tend,1000);
%% Solve the PDE
% Solves the continuum limit equations
sol=pdepe(m,@pdefun,@icfun,@bcfun,xval,tval);
r=sol(:,:,1); % Captures rioter distribution solution
b=sol(:,:,2); % Captures bystander distribution solution
%% Plot solutions in x-space
figure(601)
plot(xval,r(end,:),'m--','linewidth',4)
hold on
plot(xval,b(end,:),'c--','linewidth',4)
xlim([0,200])
ylim([0,1])
xlabel('x')
ylabel('r(x,t), b(x,t)')
legend({'r(x,t)' 'b(x,t)'})
hold off
end