One-dimensional short-range nearest-neighbor interaction and its nonlinear diffusion limit

TL;DR

This paper models one-dimensional particle interactions with finite-range repulsion, rigorously derives their macroscopic nonlinear diffusion limits, and validates findings through numerical simulations, connecting microscopic behavior to well-known PDEs like the porous medium equation.

Contribution

It provides a rigorous derivation of nonlinear diffusion equations from microscopic particle models with finite-range interactions in one dimension.

Findings

Derivation of macroscopic nonlinear diffusion equations from microscopic models.

Identification of the porous medium equation as a limit case.

Numerical simulations confirming analytical results.

Abstract

Repulsion between individuals within a finite radius is encountered in numerous applications, including cell exclusion, i.e. avoidance of overlapping cells, bird flocks, or microscopic pedestrian models. We define such individual based particle dynamics in one spatial dimension with minimal assumptions of the repulsion force f as well as their external velocity v and prove their characteristic properties. Moreover, we are able to perform a rigorous limit from the microscopic to the macroscopic scale, where we could recover the finite interaction radius as a density threshold. Specific choices for the repulsion force f lead to well known nonlinear diffusion equations on the macroscopic scale, as e.g. the porous medium equation. At both scaling levels numerical simulations are presented and compared to underline the analytical results.

Equations

$$x_0(t)\le x_1(t)\le \cdots \le x_N(t).$$

$$f:[0,\infty )\to [0,1]\text{ is Lipschitz and nonincreasing,}\mathrm{supp}(f)=[0,1],$$

$$v\in W^{2,\infty }(\mathbb{R}).$$

$$\begin{array}{rl}\mu {\dot{x}}_0& =-F\left(x_1-x_0\right)-F_{ext}(x_0),\\ \mu {\dot{x}}_i& =F\left(x_i-x_{i-1}\right)-F\left(x_{i+1}-x_i\right)-F_{ext}(x_i),1\le i\le N-1,\\ \mu {\dot{x}}_N& =F\left(x_N-x_{N-1}\right)-F_{ext}(x_N).\end{array}$$

$$x\to NRx,t\to N^2\frac{\mu R}{F_0}t.$$

$$\begin{array}{rl}{\dot{x}}_0& =-Nf\left(N(x_1-x_0)\right)+v(x_0),\\ {\dot{x}}_i& =N\left(f\left(N(x_i-x_{i-1})\right)-f\left(N(x_{i+1}-x_i)\right)\right)+v(x_i),1\le i\le N-1,\\ {\dot{x}}_N& =Nf\left(N(x_N-x_{N-1})\right)+v(x_N).\end{array}$$

$${\omega }_i:=N(x_i-x_{i-1}),1\le i\le N,$$

$$\begin{array}{rl}{\dot{\omega }}_1& =N^2\left[2f\left({\omega }_1\right)-f\left({\omega }_2\right)\right]+N\left(v(x_1)-v(x_0)\right),\\ {\dot{\omega }}_i& =N^2\left[2f\left({\omega }_i\right)-f\left({\omega }_{i-1}\right)-f\left({\omega }_{i+1}\right)\right]+N\left(v(x_i)-v(x_{i-1})\right),2\le i\le N-1,\\ {\dot{\omega }}_N& =N^2\left[2f\left({\omega }_N\right)-f\left({\omega }_{N-1}\right)\right]+N\left(v(x_N)-v(x_{N-1})\right),\end{array}$$

$${\dot{x}}_0=-Nf({\omega }_1)+v(x_0),x_i=x_0+\frac{1}{N}\sum _{j=1}^i{\omega }_j.$$

$${\omega }_i(0)={\omega }_{i,0}\ge 0,1\le i\le N,x_0(0)=x_{0,0}.$$

$${\omega }_{\min }{e}^{-\gamma t}\le {\omega }_i(t)\le {\omega }_{\max }{e}^{\gamma t},t\ge 0,1\le i\le N,$$

$${\omega }_{\min }:=\underset{1\le j\le N}{\min }{\omega }_{j,0},{\omega }_{\max }:=\max \left\{1,\underset{1\le j\le N}{\max }{\omega }_{j,0}\right\},\gamma :=\Vert v^{\prime }{\Vert }_{\infty }.$$

$${\dot{Y}}_i(t)=$$

$$+$$

$$\le$$

$${\dot{y}}_i(t)\ge N^2\left[2f\left(e^{-\gamma t}{y}_i(t)\right)-f\left(e^{-\gamma t}{y}_{i-1}(t)\right)-f\left(e^{-\gamma t}{y}_{i+1}(t)\right)\right],2\le i\le N-1.$$

$$\bar{x}:=\frac{1}{N+1}\sum _{i=0}^N{x}_i$$

$$\frac{d}{dt}\bar{x}(t)=\frac{1}{N+1}\sum _{i=0}^{N}v(x_i(t)),$$

$${\mathcal{V}}_x(t):=\frac{1}{N}\sum _{i=0}^N{\left(x_i-\bar{x}\right)}^2$$

$$\frac{d}{dt}{\mathcal{V}}_x(t)=\frac{2}{N}\left[\sum _{i=1}^{N}f({\omega }_i){\omega }_i-\sum _{0\le i\le N}(x_i-\bar{x})\sum _{1\le j\le N,i\ne j}v(x_j)\right],$$

$$x_N(t)-x_0(t)\le e^{\gamma t}{\omega }_{\max },t\ge 0.$$

$$\begin{array}{rl}\sum _{i=0}^N{\dot{x}}_i{\phi }_i& =\sum _{i=1}^{N}f({\omega }_i)N\left({\phi }_i-{\phi }_{i-1}\right)+\sum _{i=0}^N{\phi }_{i}v(x_i).\end{array}$$

$$x_N(t)-x_0(t)=\frac{1}{N}\sum _{i=1}^N{\omega }_i(t),$$

$$\dot{\bar{x}}(t)=0,{\dot{\mathcal{V}}}_x(t)=\frac{2}{N}\sum _{i=1}^{N}f({\omega }_i){\omega }_i\ge 0,\text{for all }t\ge 0.$$

$$\Delta s:=\frac{1}{N},s_i:=i\Delta s,0\le i\le N,$$

$${\partial }_t\omega =-{\partial }_s^{2}f(\omega )+{\partial }_{s}v(x),0<s<1,$$

$$f(\omega (0,t))=f(\omega (1,t))=0,$$

$$\omega \ge 1,s=0,1.$$

$${\partial }_{t}x=-{\partial }_{s}f({\partial }_{s}x)+v(x),0<s<1,$$

$${\partial }_{s}x\ge 1,s=0,1.$$

$$x=X_0(t)+{\int }_0^s\omega (\sigma ,t)d\sigma ,$$

$${\partial }_s\to \omega {\partial }_x,{\partial }_t\to {\partial }_t-\omega {\partial }_{x}f(\omega ){\partial }_x+v{\partial }_x,$$

Full text

One-dimensional short-range nearest-neighbor interaction and its nonlinear diffusion limit

M. Fischer, L. Kanzler, C. Schmeiser

M. Fischer: Universität Wien, Fakultät für Mathematik, Oskar-Morgenstern-Platz 1, 1090 Wien, Austria, & Universität Graz - Institut für Mathematik und Wissenschaftliches Rechnen, Harrachgasse 21, 8010 Wien, Austria, michi.fischer@univie.ac.at

L. Kanzler: CEREMADE, Université Paris Dauphine, Place du Maréchal de Lattre de Tassigny, F-75775 Paris Cedex 16, France, laura.kanzler@dauphine.psl.eu

C. Schmeiser: Universität Wien, Fakultät für Mathematik, Oskar-Morgenstern-Platz 1, 1090 Wien, Austria, christian.schmeiser@univie.ac.at

Abstract

Repulsion between individuals within a finite radius is encountered in numerous applications, including cell exclusion, i.e. avoidance of overlapping cells, bird flocks, or microscopic pedestrian models. We define such individual based particle dynamics in one spatial dimension with minimal assumptions of the repulsion force $f$ as well as their external velocity $v$ and prove their characteristic properties. Moreover, we are able to perform a rigorous limit from the microscopic to the macroscopic scale, where we could recover the finite interaction radius as a density threshold. Specific choices for the repulsion force $f$ lead to well known nonlinear diffusion equations on the macroscopic scale, as e.g. the porous medium equation. At both scaling levels numerical simulations are presented and compared to underline the analytical results.

Keywords: agent-based models, repulsive force, cell-exclusion, nonlinear diffusion limit, porous medium equation

Mathematics subject classification: 82C22, 35R37, 35K55, 92C15

1. Introduction

We consider a chain of particles keeping their order along a straight line and interacting with their neighbors by distance dependent repulsive forces, which vanish above an equilibrium distance. Their movement is further influenced by external forces in the form of a position dependent velocity. For a finite number of such particles the distance between the first and the last particle will remain finite for all time. Our goal is to derive a macroscopic continuum model sharing this property, i.e. for an initial particle density with bounded support there should be a finite upper bound for the length of the support at later times. Such a model, in the form of a nonlinear diffusion equation with drift, will be derived from the particle model by a continuum limit. The macroscopic model is diffusive since we choose a friction dominated (overdamped) microscopic model, being motivated by the dynamics of bacterial colonies living in viscous environments.

Repulsive effects with a finite radius are often used in microscopic particle systems and their corresponding kinetic and macroscopic models. In flocking models, such as the Cucker-Smale- and the Vicsek-model, see [12, 38], they appear as part of an interaction between attraction and repulsion. Examples occur in the modelling of collective behaviour within sheep-herds [30], fish schools [10] and bird flocks [9]. In our setting only neighbouring particles interact if their distance does not exceed a given threshold, which causes the model to combine dynamics of metric interactions, i.e. interactions depending on the distance defined on the state-space, and topological interactions, i.e. interactions depending on the relative separation, which is given by the number of intermediate individuals. For various interacting agent systems, e.g. bird flock dynamics, topological interactions serve as a very realistic model approach, hence its recent interest [3], or see also [20] in a kinetic context.

Repulsion forces between individuals are furhter highly relevant in modelling size exclusion effects. On the microscopic level, this has been studied for pattern formation in bacterial colonies [14, 15, 39]. In general, such modelling of size exclusion is a macroscopic alternative to models based on e.g. cellular automata, compare [32], or microscopic asymmetric exclusion process, see [6], in which size exclusion has no influence on the diffusive term, typically a linear diffusive-term is derived. The modelled repulsion can also be seen as a cut-off potential. On a microscopic level this is often used for better computational speed and e.g. steadier movement of pedestrians as in the optimal-step-model, see [33]. A further investigation of size-exclusion effects between interacting particles, even grouped in species, can be found in [4, 5]. Performing the many particle limit under usage of the matched asymptotic expansion technique leads to a system of nonlinear cross-diffusion equations for each species. Remarkable is the effect the size-exclusion has on the diffusivities on the macroscopic scale: while excluded-volume interactions between particles of the same species favour enhancement of the diffusion rates, its diminishment of is caused by particles of the other species.

This nonlinear diffusion term can be microscopically interpreted as particles trying to reach a desired density. Density thresholds similar to velocity thresholds are an important topic in pedestrian dynamics, see [1, 19], gaining new focus in the context of social distancing [26], and cannot be calibrated in other macroscopic models, see [19]. Moreover, nonlinear diffusion equations are widely studied in the form of the well-known porous medium equation. Similar to this work, in [29] various versions of the porous medium equation have been derived in the many particle limit, starting from systems of ordinary differential equations modelling their interaction. While in [29] the macroscopic limit was achieved in a mean-field sense with varying interaction intensities, our approach relies on the definition of the discrete density, measuring the inverse of the difference of two neighboring particle positions. For an overview over the theory of porous medium equations we refer the reader to [36, 37]. While usual porous medium models consider diffusion without threshold, our approach can be seen as a generalisation of the nonlinear diffusion of a gas expanding in porous media, where the repulsion of its particles is limited by its force only having finite range (see Section 3).

Related limits from microscopic particle dynamics on the real line to macroscopic equations can be found in  [13, 17]. In [17] the authors derive a nonlinear diffusive equation coupled with an aggregation term from the corresponding dynamics of the particle positions. On the one hand, it is very similar to our work in methods and assumptions on the nonlinearity. On the other hand, our assumptions of a purely external forcing coupled with the nonlinear interactions show a different aspect of this problem, together with a stronger convergence result.

On a macroscopic level nonlinear diffusion equations with degenerate diffusivity lead to moving boundary problems. It is closely related to the Stefan-problem, see classical books [27, 31]. Analytically exact solutions are still matter of current research, see [11, 18] and also numerically challenges occur e.g. also in cancer research, see [35, 21]. On an agent-based level, cancer growth was investigated in [28]. More recently, in  [24], an individual-based mechanical model describing the dynamics of two contiguous cell populations with different characteristic behaviours was formulated and its formal continuum model was presented. The latter leading to the aforementioned free boundary problem with a nonlinear diffusive term, which exhibits travelling wave solutions.

This work is structured as follows. In Section 2 the microscopic model in the form of an ODE system is formulated and its characteristic properties are derived. Section 3 contains a formal derivation of the macroscopic model and a discussion of its qualitative properties. This includes the derivation of a Eulerian formulation of the model, which is originally written in terms of Lagrangian coordinates. In Section 4 the macroscopic limit is carried out rigorously, providing also a existence of solutions for the continuum model. Some of the formal results are illustrated by numerical simulations in Section 5. Finally, we conclude this article with Section 6.

2. The microscopic model – individual based dynamics

Consider a chain of $N+1$ point particles with time dependent positions $x_{i}(t)\in\mathbb{R}$, $0\leq i\leq N$, such that

[TABLE]

Neighboring particles $i$ and $i+1$ interact by a distance dependent repulsive force

$F(x_{i+1}-x_{i})$, written as $F(r)=F_{0}f(r/R)$ in terms of the dimensionless function $f$, which satisfies

[TABLE]

i.e. there is no interaction between neighbors further apart than the equilibrium distance $R>0$. Additionally, the particles are under the action of a position dependent external force, written as $F_{ext}(x)=F_{0}v(\frac{x}{NR})$, where $v$ fulfills

[TABLE]

Balancing these forces with friction against a nonmoving environment (with friction coefficient $\mu>0$) leads to the ODE system

[TABLE]

We introduce a nondimensionalization by

[TABLE]

This is a diffusive macroscopic rescaling (by the factors $N$ and, respectively, $N^{2}$) of the natural microscopic scaling. The scaled system reads

[TABLE]

We shall mostly work with a reformulation in terms of the new unknowns

[TABLE]

satisfying

[TABLE]

which is coupled to

[TABLE]

This system will be considered subject to initial conditions

[TABLE]

We start with stating global existence and boundedness of the solution.

Theorem 1**.**

Let $f$ and $v$ satisfy (1) and (2) respectively and let $0\leq\omega_{i,0}<\infty$ for all $i\in 1,\dots,N$. Then there exists a unique global solution $(\omega_{1},\ldots,\omega_{N})\in C^{1,1}\left([0,\infty)\right)^{N}$ of (6), (8), which satisfies

[TABLE]

where we defined

[TABLE]

Proof.

Global existence and uniqueness follow from the Lipschitz continuity of $f$ and $v$. In order to show the upper bound of $\omega_{i}(t)$, we define $Y_{i}(t):=e^{-\gamma t}\omega_{i}(t)$ for $2\leq i\leq N-1$ and estimate

[TABLE]

Assume that for a time $t_{0}$ we have $Y_{i}(t_{0})=\omega_{max}$. We remember $f\geq 0$ and ${\rm supp}(f)=[0,1]$, stated in (1). Hence, we have $f(\omega)=0$ for $\omega\geq\omega_{max}$ and we can conclude $\dot{Y}_{i}(t_{0})\leq 0$, which gives the upper bound $\omega_{i}(t)\leq e^{\gamma t}\omega_{max}$. With analogous arguments we this upper bound for the boundary values $\omega_{1}$ and $\omega_{N}$ can be derived.

In order to show the bound from below we define in a similar manner $y_{i}(t):=e^{\gamma t}\omega_{i}(t)$, from which we obtain the estimate

[TABLE]

We assume that for a time $t_{0}$ we have $y_{i}(t_{0})=\omega_{min}$ and $y_{j}(t_{0})\geq\omega_{min}$ for $1\leq j\leq N$. Due to the monotonicity of $f$ we see that $\dot{y}_{i}(t)\geq 0$ has to hold, which gives the lower bound $\omega_{i}(t)\geq\omega_{\min}e^{-\gamma t}$. In an analogous way the lower bound for the boundary points $\omega_{1}$ and $\omega_{N}$ can be shown, which completes the proof. ∎

Theorem 2**.**

Let $(x_{0},\ldots,x_{N})$ be a solution of (4), satisfying (5), (8). Then

the center of mass

[TABLE]

moves with respect to the average of all particle velocities:

[TABLE] 2. 2)

the change in time of the variance

[TABLE]

is given by

[TABLE] 3. 3)

the distance between the leftmost and the rightmost particle remains bounded in bounded time-intervals:

[TABLE]

Proof.

The scalar product of (4) with $(\varphi_{0}(t),\ldots,\varphi_{N}(t))$ gives, after summation by parts,

[TABLE]

This immediately implies 1) with $\varphi_{0}=\cdots=\varphi_{N}=1$. The choice $\varphi_{i}(t)=\frac{2}{N}(x_{i}(t)-\bar{x})$, $i=0,\ldots,N$, gives the time evolution of $\mathcal{V}_{x}$ proving 2). Statement 3) is a consequence of

[TABLE]

and of the upper bound in Theorem 1. ∎

Remark 3*.*

We would like to point out that for the case where the particles are not exposed to an external force, i.e. $v\equiv 0$, the center of mass $\bar{x}$ is conserved and the variance is non-decreasing in time:

[TABLE]

Hence, the particle positions expand around their center of mass, but not too much.

3. The macroscopic model

The continuum model in Lagrangian coordinates:

Interpreting the particle index as a discrete Lagrangian variable, the connection to the continuum is made by the definitions

[TABLE]

a discretization of the Lagrangian coordinate $s\in[0,1]$. Assuming the existence of a function $\omega(s,t)$, such that $\omega_{i}(t)\approx\omega(s_{i},t)$ as $N\to\infty$, the formal limit of the second equation in (6) gives

[TABLE]

a nonlinear diffusion equation with the diffusivity $-f^{\prime}(\omega)\geq 0$, which is bounded by the Lipschitz continuity of $f$. The limits of the first and third equation in (6) lead to

[TABLE]

equivalent to

[TABLE]

This looks like incomplete information on the boundary. However, it is sufficient in view of the degenerate diffusivity. If, on the one hand, $\omega>1$ next to the boundary, then the diffusivity vanishes there, the solution does not change with $t$, and the boundary condition is satisfied. If, on the other hand $\omega(0+,t)\leq 1$ (or, respectively, $\omega(1-,t)\leq 1$) then the solution needs to take the boundary value 1.

Similarly, the continuum limit of the particle positions $x(s,t)$ in (4), satisfying $\partial_{s}x=\omega$, solves the Neumann type problem

[TABLE]

Eulerian coordinates – the particle density:

Our next goal is to write the equation (11) in terms of the Eulerian coordinate $x$ instead of the Lagrangian coordinate $s$. This produces a moving boundary problem posed on the $x$-interval $[X_{0}(t),X_{1}(t)]:=[x(0,t),x(1,t)]$. The coordinate transformation can be written as

[TABLE]

implying

[TABLE]

where (11) and (14) have been used. Consequently the Eulerian version of (11) reads

[TABLE]

This is equivalent to the conservation law

[TABLE]

for the macroscopic particle density $\rho:=1/\omega$, which is complemented by the boundary conditions

[TABLE]

and by the dynamics of the moving boundaries, determined from (14):

[TABLE]

Jump discontinuities:

As a consequence of the degeneracy of the diffusivity $D(\rho):=-\rho^{-2}f^{\prime}(1/\rho)$ for $\rho\leq 1$, the nonlinear diffusion equation (16) supports jumps between values $\rho(x_{*}-,t)=1$ and $\rho(x_{*}+,t)<1$. The velocity of the jump location $x_{*}(t)$ is given by the Rankine-Hugoniot condition

[TABLE]

With the obvious changes also the case $\rho(x_{*}-,t)<1$ and $\rho(x_{*}+,t)=1$ can be considered. Such jumps typically separate regions where $\rho\geq 1$ from regions where $\rho<1$. In the case $\rho(X_{0,1}(t),t)=1$, the moving boundary equation (17) can be seen as a special case of (18), with $\rho$ continued by zero outside of $[X_{0},X_{1}]$. Otherwise, when $\rho(X_{0,1}(t),t)<1$, the boundary does not move. Moreover, note that if $v\equiv 0$ and with the additional natural assumption $\partial_{x}\rho(x_{*}(t)-,t)\leq 0$, the formula above implies $\dot{x}_{*}\geq 0$, i.e. the jump moves towards the region of lower density.

Special cases for $v\equiv 0$

: We below mentioned special cases will be considered in the case of no external force on the particles, i.e. $v\equiv 0$.

First we note that the statements above do not cover the situation of initial data with a smooth transition between $\rho>1$ and $\rho<1$. Consider $v\equiv 0$ and an initial datum $\rho_{0}(x)=1-cx$, $c>0$, and the choice $f(\omega)=(1-\omega)_{+}$. We expect that a discontinuity develops with location $x_{*}(t)$ starting at $x_{*}(0)=0$. Approximating the denominator on the right hand side of (18) by its value at $t=0$, we obtain

[TABLE]

and therefore $x_{\star}(t)\approx\sqrt{2t}$ for small $t$. This behavior with infinite initial velocity also occurs in the Stefan-problem, see e.g. [2, Chapter 1, Example 1]. It can also be seen on a microscopic level in Figure 1(a) for $x_{0}$ and $x_{N}$.

Another special case, again for $v\equiv 0$, is initial data with $\rho_{0}>1$ in a bounded interval, $\rho_{0}=1$ outside of it, and $f(\omega)=(1/\omega-1)_{+}^{m}$, $m\in\mathbb{R}$. For $m\geq 1$ the shifted density $\rho-1$ solves the porous medium equation with initial data with bounded support, a problem very well studied (see e.g. [37]). In particular, in that case supp$(\rho-1)$ will grow, and the long-time behavior of the solution is given by an explicitly computable self-similar Barenblatt profile. For $m\in[0,1)$ we are in the case where $\rho-1$ solves the fast diffusion equation, which also is well investigated, see e.g. [8]. For $m<0$ this phenomenon is known as super-fast diffussion. For a survey of results regarding these types of nonlinear diffusion we refer the reader to [36]. If the initial data are such that the flux $\partial_{x}f(1/\rho)$ initially vanishes at the boundary of the support, then a waiting time phenomenon occurs, where the edges of the support start moving at a positive time. Indeed, for initial data where the right-most point is given by $x_{*}(0)$, such that $\rho(x_{*}(0)-)>1$ and $\rho(x_{*}(0)+)=1$, we can calculate under consideration $f\left(1/\rho\right)=(\rho-1)_{+}^{m}$ similar to (18) the velocity of the jump location

[TABLE]

One can see clearly that the smaller the exponent $m$ and the closer the left-sided limit $\rho_{0}(x_{*}(0)-)$ is to the value 1, the flatter has to be the initial density $\rho_{0}$ left of the boundary point $x_{*}(0)$ in order to observe the aforementioned waiting time phenomenon.

Less clear is the situation with initial data of the form $\rho_{0}(x)=1-cx$, $c>0$ and a general nonlinearity $f$. We conjecture that, on the one hand, a discontinuity develops with infinite initial speed as above, whenever $x=o(\partial_{x}f(1/\rho))$ as $x\to 0-$, and on the other hand, the discontinuity only appears after a waiting time for $\partial_{x}f(1/\rho)=o(x)$ as $x\to 0-$. However, we are not aware of any rigorous results on these questions.

Decay to equilibrium:

The dynamics of (16) dissipates the $L^{2}$-norm. After continuation of $\rho$ by zero outside of $[X_{0},X_{1}]$, we obtain

[TABLE]

which is non-positiv, if $v$ is non-decreasing. Moreover, if we again consider the special case $v\equiv 0$, we observe that the dissipation vanishes for $\rho<1$ or $\rho$ independent from $x$. During the evolution we expect to see intervals $I_{+}(t)$, where $\rho\geq 1$, separated from intervals $I_{-}(t)$, where $\rho<1$, by moving jump discontinuities, where $\rho=1$ are the boundary conditions for the intervals $I_{+}$. Therefore we expect that equilibria have intervals $I_{+,\infty}$, where $\rho\equiv 1$, separated by intervals $I_{-,\infty}$, where $\rho<1$ and otherwise arbitrary. The number of $I_{+,\infty}$-intervals might be smaller than that of $I_{+}(0)$-intervals, since these intervals might merge by collisions of the moving jump discontinuities. As soon as this coarsening process is over, the limit $I_{+,\infty}=[a,b]$ of each interval $I_{+}(t)$ can be predicted. Let $[c,d]$ be big enough to contain $[a,b]$ with $\rho(c,t)=\rho_{0}(c),\rho(d,t)=\rho_{0}(d)<1$. Then $a$ and $b$ can be computed from the conservation of mass and of the center of mass:

[TABLE]

4. The rigorous macroscopic limit

The macroscopic limit will be carried out in the individual based model in terms of the unknowns $\omega_{i}$, as in (6). However, since the velocity of $x_{0}$ may be unbounded as $\Delta s=1/N\to 0$ (see (7)), we represent the variables $x_{i}$ in terms of their average $\bar{x}$, satisfying (9). Therefore we consider the initial value problem

[TABLE]

coupled to

[TABLE]

with

[TABLE]

The last relation is obtained by averaging the expression for $x_{i}$ in (7), which allows to express $x_{0}$ in terms of $\bar{x}$. For the initial data we still assume

[TABLE]

This problem is equivalent to (6), (7), and therefore the results of Theorem 1 remain valid, i.e. the existence and uniqueness of a global solution of (21)–(23), satisfying

[TABLE]

where we recall the notation $\gamma:=\|v^{\prime}\|_{\infty}$. The connection to the continuum is made by the definition of the piecewise constant interpolants

[TABLE]

With these definitions, (22), (23) imply

[TABLE]

For passing to the macroscopic limit in the nonlinearities $f$ and $v$, some regularity of $\omega$ and $x$ will be needed. Since the observations of the preceding section show that jump discontinuities of $\omega$ have to be expected, bounded variation is the best regularity we can hope for.

Lemma 4**.**

Let $(\omega_{1},\ldots,\omega_{N},\bar{x},x_{0},\ldots,x_{N})$ be a solution of (21)–(23), with initial data satisfying (24) and $TV(\omega_{\Delta s}(\cdot,0))$ bounded idependently from $\Delta s$. Then

[TABLE]

where $\omega_{\Delta s}$ is defined in (25).

Proof.

The result for $x_{\Delta s}(\cdot,t)$ is an immediate consequence of its monotonicity and of Theorem 2, 3). The bound on the total variation of $\omega_{\Delta s}$ will be obtained by calculating its time derivative:

[TABLE]

where we used summation by parts for the formulation of the term $\Sigma_{1}$. By separately checking the cases that the signs are equal or different we observe

[TABLE]

For dealing with the second term $\Sigma_{2}$, we note that

[TABLE]

where the remainder term contains mean values $\tilde{x}_{i}\in(x_{i},x_{i+1})$ and $\hat{x}_{i}\in(x_{i-1},x_{i})$. The boundedness of $\omega_{\Delta s}$ and Assumption (2) imply

[TABLE]

for $c(t)>$ independent from $\Delta s$. Combining our results, we arrive at

[TABLE]

which gives a bound on $TV(\omega_{\Delta s}(\cdot,t))$ on every finite time interval by the Gronwall inequality. ∎

These total variation bounds can be used to also get some regularity in time.

Lemma 5**.**

With the assumptions of Lemma 4, $\partial_{t}\omega_{\Delta s},\partial_{t}x_{\Delta s}\in L^{\infty}_{loc}([0,\infty);W^{-1,\infty}((0,1)))$ uniformly as $\Delta s\to 0$.

Proof.

For a test function $\varphi\in W^{1,\infty}_{0}((0,1))$ we define

[TABLE]

and compute

[TABLE]

This leads to the estimate

[TABLE]

where we have used the Lipschitz constant $L_{f}$ of $f$, Lemma 4, and Theorem 1.

For the time derivative of $x_{\Delta s}$ we proceed similarly:

[TABLE]

∎

Since, for $1\leq q<\infty$, $BV([0,1])\subset L^{q}((0,1))\subset W^{-1,\infty}((0,1))$, where the first inclusion is compact [25, Corollary 3.49], we conclude from Lemmas 4 and 5 and from [34] that $\{\omega_{\Delta s},\Delta s>0\}$ and $\{x_{\Delta s},\Delta s>0\}$ are relatively compact in $L^{p}_{loc}\left([0,\infty)\times[0,1]\right)$ for every $p<\infty$. This finally leads to the following convergence result:

Theorem 6**.**

Let $\omega_{\Delta s}(\cdot,0)\in BV([0,1])$ and let (24), (1), and (2) hold. Let $\omega_{\Delta s}$ be defined by (25) in terms of the solution of (21)–(23). Then $\lim_{\Delta s\to 0}\omega_{\Delta s}=\omega$ and $\lim_{\Delta s\to 0}x_{\Delta s}=x$ in $L^{p}_{loc}([0,\infty)\times[0,1])$ for any $1\leq p<\infty$, restricting to appropriate subsequences. The limit $\omega\in L^{\infty}((0,\infty)\times(0,1))$ is a weak solution of the initial value problem for (11), satisfying the boundary conditions (12) weakly with respect to time. The function $x$ is determined from

[TABLE]

Proof.

For a test function $\varphi\in C_{0}^{\infty}([0,\infty)\times(0,1))$ we test (21) against $\varphi_{i}(t)$ (defined in (27)), noting that $\varphi_{1}=\varphi_{N}=0$ for $\Delta s$ small enough. After an integration by parts with respect to $t\geq 0$ and summation by parts with respect to $i$ we obtain

[TABLE]

where the last equation follows from

[TABLE]

and from

[TABLE]

Restricting to appropriate subsequences we have $\omega_{\Delta s}\to\omega$ in $L^{p}_{loc}\left([0,\infty)\times[0,1]\right)$ as $\Delta s\to 0$, and we may pass to the limit in (29):

[TABLE]

where $\omega_{0}(s)$ is the limit of the initial data $\omega_{\Delta s}(0)$, which exists due to the assumption $\omega_{\Delta s}(0)\in BV(\mathbb{R})$. This is the weak formulation of the initial value problem for (11).

Since $x_{0}$ and $x_{N}$ are bounded, we can pass to the limit in the weak formulations of

[TABLE]

(see (4)) to obtain the boundary conditions (12). Finally, (28) is obtained by passing to the limit in (26). ∎

Remark 7*.*

For the case $v\equiv 0$ an existence theory for the continuous problem written in terms of $x(s,t)$ (see (14)) can also be carried out by interpreting the problem as gradient flow for the energy functional

[TABLE]

The basic theory (see e.g. [16]), however, only gives $x\in C([0,\infty);L^{2}((0,1)))$ and not much information on $\omega=\partial_{s}x$.

5. Numerical simulations

5.1. Microscopic model:

We illustrate the previous statements with numerical experiments in $x$ and $\omega$. We solve the systems (4) and (6) for the following two choices of nonlinearities

[TABLE]

with an implicit Euler algorithm to conserve the characteristic properties. We discretize as follows,

[TABLE]

including the boundary values $x_{0}$, $x_{N}$, based on classical ideas as in [7]. System (6) can be calculated in every timestep from the results of (31) to avoid coupling.

For $v\equiv 0$ and $f=f_{1}$ we simulate $N=20$ agents and chose a time-stepping of $\Delta t=0.1\Delta s^{2}$, in which a typical parabolic CFL-condition in incorporated. The non-linearity in $f$ is solved with a fixed-point approach over $n=40$ iterations as proposed in [22]. The results can be found in Figure 1.

The Min-Max-principle in Theorem 1, visualized in Figure 1(d), already shows for $N=20$ the relation with a parabolic system which has in general a smoothing effect, as can be seen in Figure 1(b). However, this effect does not occur if the points $x_{i}$ are too far apart, which is visualized in Figure 1(c). Indeed, in that case groups of particles remain unmoved, since their metric distance is above the given interaction threshold.

Figure 2 corresponds to simulations for $f=f_{1}$ and non-trivial $v$, where the two simulated choices of $v$ can be seen in Figure 2(a). In Figure 2(b) the agents are accelerating if their position is around $x=1$, caused by the enhanced value of $v_{1}$ . However, the distance (depending on $\Delta s$ and $N$) to the relative leading agent does not decrease due to the repulsion. For $v_{2}$ (orange in Figure 2(a)) the agents congregate around position $x=1$, before the repulsion acts due to the resulting higher densities, which can be seen in Figure 2(c).

5.2. Macroscopic model

We investigate problem (16) on an open domain $\mathbb{R}$ with $\rho_{0}>0$ and discretize it as $\rho(i\Delta x,j\Delta t)=\rho^{j}_{i}$ explicit in time via

[TABLE]

In the case $f=f_{1}$ we discretized the Laplacian by the usual finite difference approximation and used a first-order upwind scheme depending on the sign of $v$.

We choose $\Delta x=0.001$ for a sharp visualisation of the shock and again $\frac{\Delta t}{\Delta x^{2}}=0.1$ in compliance with typical CFL-conditions, following classical literature, see [23]. All initial-data chosen for the following simulations are typical sums of bump-functions of the form

[TABLE]

and constants.

In Figure 3(a) we see the combination of the diffusive effect and a drift to the left of system (16) for the choices $f=f_{1}$ and $v\equiv 1$. The diffusive effect acts stronger than the external forcing $v$.

Figure 3(b) corresponds to simulations with $f=f_{1}$ and $v\equiv 0$ were the initial data were chosen such that they continuously pass through the threshold-value 1. The position of the jump $x_{*}$, which is formed for positive time, was calculated by (18). We used for the discretization of (18) an explicit Euler algorithm for solving the ODE. Knowing $x_{*}$ is monotone increasing and considering the left- and right-handed limits, we used in the numerator a downwind- and for the denominator an upwind-approach. We see that the discretisation (32) creates the correct shock-speed. Additionally continuous initial data lead to satisfying results since as mentioned, the discontinuity only occurs for $t=0$ and $x_{*}$ is Lipschitz for $t>0$.

Figure 3(c) corresponds to the discussion regarding decay to equilibria for the case $v\equiv 0$ and $f=f_{1}$, which can be found at the end of Section 3. We visualized a simulation for showing two colliding plateaus $\rho_{1,2}$ (with $\rho_{0}=\rho_{1}+\rho_{2}$), which are spreading in space, since each of their densities exceeds the threshold value 1, until they merge. The shocks move towards each other until the collision at time $t=0.031$ happens. The dynamics stop, once the final interval-length $[a,\,b]$ is reached, where the boundaries $a$ and $b$ are calculated by (19)-(20).

In Figure 3(d) we visualized the in Section 3 mentioned waiting-time phenomenon, which occurs in porous medium equations. We chose $f=f_{2}$, no external forcing $v\equiv 0$ and initial data

[TABLE]

Therefore, $\partial_{x}f(1/\rho)$ initially vanishes at the boundary of the support of $\rho-1$, which is given by $x=\pm\frac{1}{4}$ and is marked with a dashed line.

5.3. Comparison of the Microscopic and the Macroscopic Model

We conclude the numerical simulations with an experiment, which compares the microscopic and the macroscopic dynamics in order to show consistency between the two scales. We start with a continuous initial density $\rho_{0}+0.5$ as in Equation (33) (with $h=3$, $b=0.3$ and $m=0$). We calculate the corresponding discrete values $x_{i}$ via (15). In Figure 4 we see the solutions of the corresponding microscopic dynamics (6), solved numerically by (32), and the one of the corresponding macroscopic dynamics (16), solved via (31), which are plotted beside each other.

6. Conclusion

In this work we modelled particles interacting in a certain radius and driven by an external force. We introduce the scaled distance $\omega$, which can also be seen as a derivative in Lagrangian coordinates $s$. From this interpretation, it was possible to define a microscopic density $\rho$ by inverting $\omega$.

For the microscopic systems (4) and (6) we could establish an existence and uniqueness result, together with a maximum-principle (Theorem 1). This result as well as observations of the particle’s positions (Theorem 2) are underlined with simulations in Section 5.1. Properties of (6) for $x$ and $\omega$ can be transferred to the density $\rho$.

On a macroscopic level our main focus was on (16), a conservation law for which we could derive jump conditions for discontinuous initial data, which can be related to the moving boundary of a Stefan problem. Moreover, we discussed specific choices for the repulsive force $f$, which lead to well known nonlinear diffusion equations as the porous medium equation or the fast diffusion equation. A rigorous limit from (6) to (LABEL:om-IVP) could be established in Theorem 6. Due to the non-linearity $f$, passing to the limit is non-trivial. In order to conclude by finding a weakly convergent sequence for system (6), we used bounds on the total variation with respect to the spatial variable of the solution and a compact interpolation theorem [34].

While this work has focused solely on diffusion in one spatial dimension, there are several natural extensions of the model that arise from practical considerations: for instance, source/sink terms, as well as aniosotropic effects that could arise in higher dimension. One motivating biological example that incorporate these effects is Bacteria growth:

[TABLE]

where $\alpha\rho(\rho_{*}-\rho)$ models the growth or death of bacteria depending on the carrying capacity $\rho_{*}>0$. In Figure 5 corresponding simulations for $f=f_{1}$, $\alpha=0.5$ and $\rho_{*}=1.5$ can be found. We leave the physical study of this and related models also in higher spatial dimensions as well as their rigorous derivation to future work.

Acknowledgments:

The authors are grateful to three anonymous referees, whose detailed comments led to a significant improvement of this article. Moreover, the authors thank Diane Peurichard for fruitful discussions forming the starting point for these investigations.

This work has been supported by the Austrian Science Fund, grants no. W1245 and F65. L.K. received funding by the European Commission under the Horizon2020 research and innovation programme, Marie Sklodowska-Curie grant agreement No 101034255.