SetDirectory[NotebookDirectory[]];
Get["DiscreteToPDE.m"];

Bidirectional Pedestrian Flow

(* This is the running example from our paper. *)
pde=DiscreteToPDE[{
T[r,{1,0}]==(1-ρ[1,0])*(1+α*r[2,0]),
T[r,{0,-1}]==(1-ρ[0,-1])*((Subscript[γ, 0])+(Subscript[γ, 1])*b[1,0]),
T[r,{0,1}]==(1-ρ[0,1])*((Subscript[γ, 0])+(Subscript[γ, 2])*b[1,0]),
T[b,{-1,0}]==(1-ρ[-1,0])*(1+α*b[-2,0]),
T[b,{0,1}]==(1-ρ[0,1])*((Subscript[γ, 0])+(Subscript[γ, 1])*r[-1,0]),
T[b,{0,-1}]==(1-ρ[0,-1])*((Subscript[γ, 0])+(Subscript[γ, 2])*r[-1,0])
}/.ρ[i_,j_]->r[i,j]+b[i,j],{x,y}]

(* This is just some syntactic sugar. *)
Format[(Subscript[γ, 12])]:=(Subscript[γ, 1])+(Subscript[γ, 2]);
Format[(Subscript[γ, 21])]:=(Subscript[γ, 1])-(Subscript[γ, 2]);
pde1=pde/.{(Subscript[γ, 1])->((Subscript[γ, 12])+(Subscript[γ, 21]))/2,(Subscript[γ, 2])->((Subscript[γ, 12])-(Subscript[γ, 21]))/2};

(* The mean-field PDEs in conservative form *)
Partial[t][b]==Collect[pde1[[1]],{h,(Subscript[γ, 12]),(Subscript[γ, 21])},PartialIntegrate[#,{r,b},{y,x}]&]

Partial[t][r]==Collect[pde1[[2]],{h,(Subscript[γ, 12]),(Subscript[γ, 21])},PartialIntegrate[#,{b,r},{y,x}]&]

Fernando, Landman, Simpson: “Nonlinear diffusion and exclusion processes with contact interactions” (Physical Review E 81, 011903, 2010)

Simplest case

(* This are the transition probabilities from the simple case (4). *)
(* The result shows that D(C)=1, as claimed in the paper. *)
DiscreteToPDE[{
T[c,{0,1}]==p/4*(1-c[0,1]),
T[c,{0,-1}]==p/4*(1-c[0,-1]),
T[c,{1,0}]==p/4*(1-c[1,0]),
T[c,{-1,0}]==p/4*(1-c[-1,0])},{x,y}]

Contact-maintaining interactions

(* Transition probabilities as given in (6) and (5). *)
(* We use Moore interacting neighbors (defined in neigh). *)
neigh=DeleteCases[Flatten[Table[{i,j},{i,-1,1},{j,-1,1}],1],{0,0}];
maint[v_List]:=Intersection[neigh,v+#&/@neigh];
disc=Function[v,T[c,v]==p/4*(1-c@@v)*With[{P=1-Times@@(1-c@@#&/@maint[v])},u*P+w*(1-P)]]/@{{0,1},{0,-1},{1,0},{-1,0}}

{pde}=DiscreteToPDE[disc,{x,y}]

PartialIntegrate[pde,{c},{x,y},Depth->1]

(* This agrees with the result given in (7). *)
%/.a_Plus/;Length[a]>2:>u+Factor[a-u]

(* This is a special case of the general formula implemented in DiscreteModel: *)
Simplify[(Last/@DiscreteModel[{{0,1},{0,-1},{1,0},{-1,0}},neigh]/.{s->1,p->1,q->1,r->1})-(Last/@disc/.p->1)]

Combined contact-forming or contact-breaking and contact-maintaining interactions

(* This are the transition probabilities as given in (13), *)
(* for 2D square lattice and Moore interacting neighbors. *)
disc=DiscreteModel[{{0,1},{0,-1},{1,0},{-1,0}},neigh]

Timing[{pde}=DiscreteToPDE[disc,{x,y}]]

PartialIntegrate[pde,{c},Depth->1]

(* Verify that our result agrees with Equation (14) in the paper. *)
Expand[(-6 c p+15 (c^2) p-8 (c^3) p-6 (c^4) p+6 (c^5) p-(c^6) p+12 c q-39 (c^2) q+46 (c^3) q-24 (c^4) q+6 (c^5) q-(c^6) q+9 (c^2) r-18 (c^3) r+15 (c^4) r-6 (c^5) r+(c^6) r+s-6 c s+15 (c^2) s-20 (c^3) s+15 (c^4) s-6 (c^5) s+(c^6) s) (4 c u-6 (c^2) u+4 (c^3) u-(c^4) u+w-4 c w+6 (c^2) w-4 (c^3) w+(c^4) w)-(u+(w-u)*(1-c)^4)*(r+(s-p-q+r)*(1-c)^6+(p*(1-9c)+q*(1+9c)-2r)*(1-c)^3)]