A multi-scale particle method for mean field equations: the general case

TL;DR

This paper introduces a multi-scale meshfree particle method capable of efficiently approximating mean field equations across different interaction scales, including hyperbolic limits, with applications in swarming, traffic, and granular flow simulations.

Contribution

It extends previous methods to a more general case, including hyperbolic limits, and demonstrates a numerical transition between microscopic and macroscopic modeling approaches.

Findings

High computational efficiency near the macroscopic limit.

Effective handling of both large and small interaction radii.

Versatile application to various mean field problems.

Abstract

A multi-scale meshfree particle method for macroscopic mean field approximations of generalized interacting particle models is developed and investigated. The method is working in a uniform way for large and small interaction radii. The well resolved case for large interaction radius is treated, as well as underresolved situations with small values of the interaction radius. In the present work we extend the approach from [39] for porous media type limit equations to a more general case, including in particular hyperbolic limits. The method can be viewed as a numerical transition between a DEM-type method for microscopic interacting particle systems and a meshfree particle method for macroscopic equations. We discuss in detail the numerical performance of the scheme for various examples and the potential gain in computation time. The latter is shown to be particularly high for situations near the macroscopic limit. There are various applications of the method to problems involving mean field approximations in swarming, tra?c, pedestrian or granular fow simulation.

Figures

Tables

Table 1. Table 1: ℒ 2 superscript ℒ 2 \mathcal{L}^{2} -errors and CPU time for Example 3 with R = 0.02 𝑅 0.02 R=0.02 .

$\#$ Particles	Fokker-Plank	naive	multi-scale	CPU time
	error	error	error	seconds
$200$	0.29	$5.99$	$0.51$	$26$
$400$	0.21	$1.14$	$0.33$	$51$
$800$	0.15	$0.77$	$0.21$	$104$
$1600$	0.07	$0.55$	$0.10$	$216$

Table 2. Table 2: Convergence study for nonlocal Lighthill-Whitham equations with downwind interaction potential and R = 0.2 𝑅 0.2 R=0.2 .

$\#$ particles	naive	multi-scale	CPU time
	error	error	in seconds
$200$	$0.40$	$0.05$	8
$400$	$0.26$	$0.02$	17
$800$	$0.15$	$0.01$	36
$1600$	$0.08$	$0.002$	77

Table 3. Table 3: Convergence study for nonlocal Lighthill-Whitham equations with symmetric interaction potential and R = 0.2 𝑅 0.2 R=0.2 and δ = 0.02 𝛿 0.02 \delta=0.02 .

$\#$ Particles	naive	multi-scale	CPU time
	error	error	in seconds
$200$	$0.45$	$0.13$	$9$
$400$	$0.31$	$0.05$	$19$
$800$	$0.21$	$0.02$	$50$
$1600$	$0.11$	$0.01$	$142$

Table 4. Table 4: Comparison of CPU times between microscopic and multiscale simulations of the hydrodynamic 2D equations. The error analysis is performed at time 18.75 sec.

initial	$\#$ particles	naive	multi-scale	CPU time
spacing		error	error
$0.6$	$4150$	$0.38$	$0.07$	15 min
$0.36$	$11600$	$0.37$	$0.19$	$40$ min
$0.3$	$16600$	$0.29$	$0.16$	$71$ min
$0.24$	$26000$	$0.20$	$0.11$	$90$ min
$0.21$	$34000$	$0.14$	$0.09$	$118$ min
$0.18$	$46000$	$0.11$	$0.07$	$154$ min
$0.1$	$136000$	$0.02$	$0.01$	$603$ min

Table 5. Table 5: Comparison of CPU times between microscopic and multiscale simulations of the hydrodynamic 2D equations for Example 5. The error analysis is performed at time 12.5 sec.

initial	$\#$ particles	naive	multi-scale	CPU time
spacing		error	error
$1$	$1400$	$0.54$	$0.14$	$8$ min
$0.5$	$5700$	$0.36$	$0.18$	$23$ min
$0.35$	$11500$	$0.48$	$0.22$	$52$ min
$0.2$	$35200$	$0.16$	$0.14$	$223$ min

Equations

$${\rho }_R^N(x)=\frac{1}{N}\sum _{j=1}^N{U}_R(x-x_j)$$

$$d{x}_i$$

$$d{v}_i$$

$${\int }_{B_R(0)}{U}_R(x)dx=1.$$

$$U_R=C_R{\left(R-\Vert x\Vert \right)}^2,$$

$$C_R=\frac{(d+1)(d+2)}{2{\tau }_d{R}^{d+2}}$$

$${\int }_{B_R(0)}{U}_R(x)dx=1.$$

$$d{v}_i$$

$$U_R=2C_R{\left(R-|x|\right)}^2,$$

$$U(x)=G(|x|)-L\cdot FG\left(\frac{|x|}{L}\right),$$

$$D=\int U(x)dx=2{\int }_0^{\infty }G(x)dx(1-F{L}^2)>0.$$

$$D=\int U(x)dx=d{\tau }_d{\int }_0^{\infty }G(r)r^{d-1}dx(1-F{L}^{d+1}).$$

$${\partial }_{t}f+v\cdot {\nabla }_{x}f=Sf+Lf$$

$$\begin{array}{rl}Sf& ={\nabla }_v\cdot \left(\beta (U_R\star \rho )f\right)+{\nabla }_v\cdot \left(\alpha (U_R\star \rho ){\nabla }_x{U}_R\star \rho f\right).\end{array}$$

$$\begin{array}{rl}Lf& =\gamma {\nabla }_v\cdot \left(vf+\frac{A^2}{2\gamma }{\nabla }_{v}f\right).\end{array}$$

$$F_R\star \rho (x)=\int F_R(x-y)\rho (y)dy$$

$$\rho (x,t):=\int f(x,v,t)dv.$$

$$\int \rho (x,t)dx=1.$$

$${\partial }_t\rho +{\nabla }_x\cdot (\rho u)$$

$${\partial }_{t}u+(u\cdot {\nabla }_x)u+\frac{\delta }{\rho }{\nabla }_x\rho$$

$$\rho u(x,t):=\int vf(x,v,t)dv.$$

$$u=-\frac{1}{\gamma }\left(\alpha (U_R\star \rho ){\nabla }_x{U}_R\star \rho +\beta (U_R\star \rho )+\frac{\delta }{\rho }{\nabla }_x\rho \right)$$

$$\gamma {\partial }_t\rho ={\nabla }_x\cdot \left(\rho \left(\beta (U_R\star \rho )+\alpha (U_R\star \rho ){\nabla }_x{U}_R\star \rho \right)+\delta {\nabla }_x\rho \right).$$

$$\int {\nabla }_x{U}_R(x-y)\rho (y)dy=-\int {\nabla }_y{U}_R(x-y)\rho (y)dy=\int U_R(x-y){\nabla }_y\rho dy\sim {\nabla }_x\rho .$$

$${\partial }_t\rho +{\nabla }_x\cdot (\rho u)$$

$${\partial }_{t}u+(u\cdot {\nabla }_x)u$$

$$\gamma {\partial }_t\rho ={\nabla }_x\cdot (\beta (\rho )\rho )+{\nabla }_x\cdot ((\delta +\rho \alpha (\rho )){\nabla }_x\rho ).$$

$$d{x}_i$$

$$d{v}_i$$

$${\partial }_t\rho ={\partial }_x\cdot (ϵ\rho {\partial }_x{U}_R\star \rho +\rho {\partial }_{x}V+(1-ϵ)\rho U_R\star \rho ).$$

$${\partial }_t\rho ={\partial }_x\cdot (ϵ\rho {\partial }_x\rho +\rho {\partial }_{x}V+(1-ϵ){\rho }^2).$$

Full text

Keywords. Meshfree methods, particle methods, asymptotic preserving methods, interacting particle systems, mean field and hydrodynamic approximations, nonlinear Fokker-Planck equations.

AMS Classification. 82C21, 82C22, 65N06

1 Introduction

Interacting particle models are used in many applications ranging from engineering sciences to social sciences and biology. We refer to [12, 16, 17, 29, 30, 31] for recent work. In certain limits, for large numbers of particles, the solution of the microscopic model can be approximated by the solution of so-called mean field or macroscopic equations. For a mathematical investigation of these limits we refer to [27, 3, 8]. Microscopic interacting particle models require numerically the solution of large systems of ordinary differential equations. Many different numerical approaches can be employed for the numerical solution of the macroscopic approximations. In particular, meshless or particle methods are a popular way to solve these problems. We refer to [4, 26] for the analysis of particle methods, to [24] for a classical particle method and to [20, 51, 52, 55] for different developments of the original idea from [24].

In the present work we extend a particle method developed in [39] which is especially adapted to the solution of macroscopic mean field equations derived from interacting particle systems. The method can be viewed as a numerical transition between a DEM-type method for microscopic interacting particle systems and a meshfree particle method for macroscopic equations. For related approaches we refer to [54, 55]. In macroscopic mean-field models an interaction term appears, which is derived from the microscopic interaction term. Usually, this term has the form of a convolution integral. The particle method approximates these convolution integrals in an appropriate way not using a microscopically large number of macroscopic particles and thus allowing the use of an underresolved meshfree method. In the limit for small $R$ a method which is consistent with the asociated limit equation is obtained. For intermediate values of $R$ one obtains with this procedure an easily calculated correction term. In the present work more general particle systems and more general limit equations are considered compared to [39] where only porous media type limits have been investigated. The present investigations include for example terms leading to hyperbolic limit equations and more general applications.

The paper is organized as follows. Section 2 contains a description of the models under investigation. In particular, the macroscopic hydrodynamic models and their scalar approximations are discussed. Examples range from non-local 1D Lightill Whitham type models to non-local 2D pedestrian flow models. Section 3.1 describes the numerical method and in detail the approximation of the convolution integral. Finally, Section 4 contains the numerical results and a comparison of the methods and their computation times for the above mentioned examples.

2 The models

We consider a hierarchy of models ranging from microscopic interacting particle systems and their mean field approximations to hydrodynamic and scalar macroscopic approximations.

2.1 Microscopic and mean field equations for interacting particle system

The starting point for our model is a general microscopic model for $N$ interacting particles [37]. We define the empirical density as

[TABLE]

and consider systems of equations of the form

[TABLE]

with $i=1,\ldots,N$, $(x_{i},v_{i})\in\mathrm{IR}^{d}\times\mathrm{IR}^{d}$, $d\in\mathbb{N}$. $\alpha:\mathrm{IR}^{+}\rightarrow\mathrm{IR}^{+},\beta:\mathrm{IR}^{+}\rightarrow\mathrm{IR}^{d}$, $\gamma,A\geq 0$ and $W_{t}$ a d-dimensional Brownian motion. $U_{R}$ is a sufficiently smooth, repulsive potential with support in $B_{R}(0)$ and positive integral which is normalized to $1$, that means

[TABLE]

We assume that $U_{R}$ approximates the delta distribution as $R$ goes to [math], $U_{R}(x)=\frac{1}{R^{d}}U_{1}(\frac{x}{R})$. Possible extensions are given in the remarks below.

Remark 1**.**

We note that physically speaking $U_{R}$ is used here on the one hand to define the empirical density, which is used in the definition of the coefficients $\alpha$ and $\beta$. On the other hand it is used as interaction potential in the term $\nabla_{x}\rho^{N}_{R}$. For the following calculations we could as well use different functions $U_{R}$ at the different places in the equation.

Example 1**.**

As an example for a repulsive interaction potential we choose for $\|x\|\leq R$

[TABLE]

where

[TABLE]

and $\tau_{d}$ is the volume of the unit sphere in $\mathrm{IR}^{d}$. The coefficients are chosen such that

[TABLE]

We note that non-symmetric potentials could be chosen as well.

Remark 2**.**

In 1D the present formulation includes one sided interaction potentials $V_{R}$ approximating the Heaviside function, i.e. $\partial_{x}V_{R}=U_{R}$. An equation of the form

[TABLE]

turns into

[TABLE]

with $\beta(\rho)=\rho\alpha(\rho)$. Compare the systems considered in [49, 6].

Example 2**.**

For the potential $V_{R}$ in the remark above with $\partial_{x}V_{R}=U_{R}$ one might choose potentials,such that $U_{R}$ is given as in example 1 or non-anticipating examples like

[TABLE]

if $x<0$ and [math] otherwise.

Remark 3**.**

For the following considerations, one might also consider more general attractive-repulsive potentials. For example in 1-D the following class of potentials can be considered, compare [43, 6]:

[TABLE]

where $G:\mathrm{IR}_{+}\rightarrow\mathrm{IR}_{+}$ is a monotone decaying, integrable function, decaying to 0 as $x$ goes to infinity. This includes, for example, the Morse potential ($G(x)=exp(-x),x>0$) used among other applications for modelling the swarming behaviour of birds, see [10]. However, for the following, the potential should lead to a spreading behaviour of the solutions goverened for long times by a diffusive equation, in our case a porous media type equation, see [45, 6, 43]. For the above class of potentials a purely repulsive potential is characterized by $F,L<1$, see [43]. However, this is not a necessary condition for the spreading behaviour of the solution. The analysis in [43] shows that one still observes a spreading behaviour, if $F<1$ and $1-F\cdot L^{2}>0$, i.e. also for $L>1$ with $L<1/\sqrt{F}$. We note that in this case

[TABLE]

In $d$ dimensions the corresponding formula for a rotationally symmetric potential $U$ defined by (5) is

[TABLE]

This is larger than [math], if $1-FL^{d+1}>0$.

Remark 4**.**

*Velocity dependent interactions $U_{R}=U_{R}(x_{j}-x_{i},v_{j}-v_{i})$ might be included, as well as further explicit dependence of $\alpha$ and $\beta$ on $x_{i},v_{i}$. For example, one could include an exterior potential $V(x)$ via $\beta(x)=\nabla_{x}V(x)$. *

For $N$ going to infinity, one can derive in the limit of a large number of particles the associated mean field equation [8, 10, 3, 44]. One obtains for the distribution function $f=f(x,v,t)$ of the particles the mean field equation

[TABLE]

with force term

[TABLE]

and diffusion term

[TABLE]

Here, the convolution is defined as

[TABLE]

and the density as

[TABLE]

We normalize

[TABLE]

2.2 Hydrodynamic and scalar macroscopic models

Multiplying the mean field equation with $1$ and $v$ and closing the equations by approximating the distribution function with a function having mean $u$ and variance $\delta\rho$, one obtains the continuity and momentum equations

[TABLE]

with the momentum

[TABLE]

Neglecting time derivatives and inertia terms in the $u$-equation we approximate the velocity $u$ as

[TABLE]

and obtain the following scalar equation for the density

[TABLE]

2.3 Localized models

Since $U_{R}$ approximates for small values of $R$ a $\delta$ distribution, one obtains formally

[TABLE]

Thus, one obtains from the hydrodynamic equations the damped isentropic Euler equations with exterior potential

[TABLE]

From the scalar equations one obtains with the same procedure as before a nonlinear Fokker-Planck equation of the form

[TABLE]

A simple example with a solution which converges to a stationary distribution is given by the following

Example 3**.**

We consider for $d=1$ the following coefficients $\alpha(\rho,x)=\epsilon,\beta(\rho,x)=\partial_{x}V(x)+(1-\epsilon)\rho,\epsilon\in(0,1]$, $A=0,\delta=0$. This leads to the microscopic problem

[TABLE]

and the scalar equation

[TABLE]

The scalar localized macroscopic approximation is

[TABLE]

with the stationary solution

[TABLE]

where $C_{\epsilon}\in\mathrm{IR}$ is determined by the normalization $\int\rho_{\infty}(x)dx=1$. $\epsilon=1$ gives the well known solution of the Fokker-Planck porous media equation $\rho_{\infty}=(C-V(x))_{+}$, see [1]. For $\epsilon$ tending to [math] one approaches a stationary distribution function, which is for a convex potential $V$ given by $\rho_{\infty}=(-V^{\prime}(x))_{+}$ for $C<x$ and [math] otherwise. $C$ is in both cases chosen such that the normalization condition is fulfilled.

We consider further examples from traffic and pedestrian flow in one- and two space dimensions.

Example 4**.**

A simple microscopic traffic model is given as follows. We model the acceleration by choosing $\gamma>0$ and $\beta=-\gamma$. Breaking interactions are modelled by an interaction potential $V_{R}$ given by a smooth version of the Heaviside function such that $\partial_{x}V_{R}=U_{R}$, where $U_{R}$ is a smooth version of the $\delta$-function. $U_{R}$ is concentrated in the negative half plane to model the fact that the interaction is essentially restricted to an interaction with the predecessors. We call $V_{R}$ a downwind potential in this case. Finally, we choose $\alpha=\gamma$ and $A=0$. We obtain with $x_{i},v_{i}\in\mathrm{IR},i=1\cdots N$

[TABLE]

*The corresponding mean field hydrodynamic equations are *

[TABLE]

and the associated scalar nonlocal viscous Lighthill-Whitham model is

[TABLE]

We remark that nonlocal Lighthill-Whitam model has been investigated in [25, 2]. There, for one-sided downwind monotone convolution kernels concentrated in the negative half plane existence and uniqueness of solutions have been shown as well as a maximum principle. In particular, the density is not exceeding the maximal density, which is in this case eqal to $1$ for suitable initial conditions. The local limits are

[TABLE]

and the Lighthill Whitham equation

[TABLE]

We note that for the one-sided downwind kernel we obtain an approximation to order $\mathcal{O}(R^{2})$ in the following way. Define

[TABLE]

using $U_{R}\star\rho=\rho+D_{R}\partial_{x}\rho$ we obtain

[TABLE]

This is a stable equation for $\delta\geq 0$. We note that using the potential defined in Example 2 we obtain $D_{1}=\frac{1}{4}$. If instead of using a down-wind interaction potential we consider symmetric potentials we obtain $D_{1}=0$ and the approximation to order $\mathcal{O}(R^{3})$ is

[TABLE]

[TABLE]

In this case we need $\delta>0$ in order to obtain a stable convergence of the solutions of equation (17) as $R$ goes to [math], compare [46]. A numerical investigation of non-local Lighthill-Whitham type equations with symmetric or upwind potentials can be found in [2]. In this case there is no maximum principle and the density might exceed the ’maximal density’.

Example 5**.**

For the 2D case we consider a model for pedestrian flow. For $x_{i},v_{i}\in\mathrm{IR}^{2},i=1,\ldots,N$ and $W(\rho)=u_{max}(1-\frac{\rho}{\rho_{max}})$ we consider the microscopic model

[TABLE]

*with symmetric interaction kernel $U_{R}$ and with *

[TABLE]

where

[TABLE]

The hydrodynamic mean field limit is

[TABLE]

and the scalar limit is

[TABLE]

where

[TABLE]

The local limits are

[TABLE]

and

[TABLE]

together with

[TABLE]

which is a viscous form of the Hughes model. For a detailed modelling of the interactions between pedestrians we refer for example to [17]. For rigorous results on the Hughes model and the approximation of the model via particle systems we refer to [18, 35, 19].

In the following, our goal will be to develop a meshfree particle method for equations (7), (8) for different ranges of parameters and, in particular, for the limit equations (9), (10).

3 Numerical method

For a general description of the particle method used here and further references on the subject, we refer to [39]. Here, we concentrate on the approximation of the interaction term. We follow the approach in [39]. However, additionally to [39] we have to treat here the terms leading to hyperbolic limit equations in a suitable way.

3.1 Approximation of the interaction term using particle methods

The key point of our method is to approximate the integrals

[TABLE]

with $F_{R}=U_{R}$ and $F_{R}=\nabla_{x}U_{R}$. We denote the Voronoi cell around particle $i$ given by the particle locations of the other particles by $V_{i}$. $|V_{i}|$ denotes the volume of this cell. Then, a naive or microscopic approximation of the integral terms would be to use

[TABLE]

where $N$ is now the number of macroscopic particles used in the particle method. For small values of $R$ this results in the following problem. We consider a situation where the method is underresolved, that means where one uses a small number of macroscopic particles in contrast to the large number of microscopic particles described by the macroscopic equations. Then, the value evaluated from (23) will be zero due to the large distances between the macroscopic particles and the relatively small value of $R$. However, the actual value of the integral will not be zero even for very small $R$ due to the corresponding ’infinite’ number of microscopic particles described by the macroscopic equations.

We resolve this problem using a higher order approximation of the integral. This yields in the limit for small $R$ a method for the limiting nonlinear Fokker-Planck equations, even if the number of macroscopic particles is still small. We use an approximation of the density given by

[TABLE]

where $V_{j}$ denotes the Voronoi cell associated to particle/mesh point $x_{j}$ and $\chi$ denotes the characteristic function. The approximation $\sigma_{j}$ of the first derivative $\nabla_{x}\rho(x_{j})$ is determined via a least squares approximation using the neighbouring points. Then, we obtain for the integral (22)

[TABLE]

We approximate first the integral around the center point $x_{i}$. To approximate the integral we distinguish between the cases $|V_{i}|>|B_{R}|$ and $|V_{i}|<|B_{R}|$. For $|V_{i}|>|B_{R}|$, we proceed as follows.

[TABLE]

This is equal to 1 for $F_{R}=U_{R}$ due to the normalization of the potential. For $F_{R}=\nabla_{x}U_{R}$ the expression is equal to 0, since

[TABLE]

Moreover,

[TABLE]

This expression is for $F_{R}=U_{R}$ equal to $-\sigma_{i}\cdot\mu$, where $\mu=\int_{B_{R}(0)}yU_{R}(y)dy$ is the mean value of $U_{R}$. For $F_{R}=\nabla_{x}U_{R}$ we obtain

[TABLE]

If $|V_{i}|<|B_{R}|$ then we first compute $R_{V_{i}}$ such that $|B_{R_{V_{i}}}|=|V_{i}|$. Then,

[TABLE]

Thus we obtain $\int_{B_{R_{V_{i}}}(0)}U_{R}(y)dy$ for $F_{R}=U_{R}$. For $F_{R}=\nabla_{x}U_{R}$ we get

[TABLE]

Finally, we compute

[TABLE]

This is equal to $-\sigma_{i}\cdot\int_{B_{R_{V_{i}}}(0)}yF_{R}(y)dy$ for $F_{R}=U_{R}$. For $F_{R}=\nabla_{x}U_{R}$ we get

[TABLE]

The integrals over the Voronoi cells with the points $x_{j},j\neq i$ as centerpoints are approximated in the following way. A simple second order approximation is given by the midpoint rule

[TABLE]

and

[TABLE]

Altogether, one obtains

[TABLE]

and

[TABLE]

with the correction factors

[TABLE]

[TABLE]

for $|V_{i}|<|B_{R}|$ and

[TABLE]

[TABLE]

for $|V_{i}|<|B_{R}|$. We note that in order to obtain a stable approximation we have to guarantee that $\alpha_{i}^{\prime}$ and $\beta_{i}\prime$ are non-negative. Alltogether, we have the following algorithm

Algorithm 1**.**

For each $x_{i}$ check for particles inside the interaction-region $B_{R}(x_{i})$. 2. 2.

Compute for all particles inside this region the size of the Voronoi cells $|V_{j}|,j=1,\ldots,N.$ 3. 3.

Use the neighbouring particles inside the neighbourhood region of radius $h$ (not only the interaction region) to compute an approximation $\sigma_{i}$ of the first derivative $\nabla_{x}\rho(x_{i})$. 4. 4.

Compute $F_{R}\star\rho(x_{i})$ according to the above formulas.

Remark 5**.**

For $|V_{j}|<<|B_{R}|$ the approximation behaves like the microscopic interaction approximation, since $\alpha_{i},\beta_{i}$ go to [math]. For $|V_{j}|>>|B_{R}|$, i.e. the underresolved situation, we have an approximation which behaves like a solution method for the macroscopic equations since $\alpha_{i},\beta_{i}$ go to $1$ and the other terms vanish.

Remark 6**.**

More accurate approximations of the integrals are possible at the expense of a more complicated approximation, compare [39].

Remark 7**.**

In a situation as in the above remark with $\beta\neq 0$ we need an upwind procedure to stabilize the numerical approximation. A first order upwind procedure amounts to adding numerical diffusion proportional to $|V_{i}|.$ In comparison, the above factor $\alpha_{i}^{\prime}$, describing the physical diffusion in case of an unsymmetric potential, is proportional to $R$. In case $|V_{i}|\sim R$, the two diffusion coefficients are of the same order. Thus, the physical diffusion is of the same order as the numerical diffusion. In order to capture the effects of the physical diffusion a higher order upwinding procedure would be necessary to reduce the numerical diffusion.

3.2 Special cases

The above formulas give for radially symmetric potentials

[TABLE]

and

[TABLE]

where $\tau_{d}$ is the volume of the unit sphere in $\mathrm{IR}^{d}$. $\alpha_{i}^{\prime}$ and $\beta_{i}$ are [math] in the radially symmetric case. Rewriting gives

[TABLE]

with

[TABLE]

for $z\in[0,1]$ and $A(z)=1,z\geq 1$. Moreover,

[TABLE]

with

[TABLE]

for $z\in[0,1]$ and $B^{\prime}(z)=1,z\geq 1$.

Remark 8**.**

We note that $\beta_{i}^{\prime}$ is positive, if the potential fulfills

[TABLE]

or

[TABLE]

This is the case for any repulsive potential. Also potentials like the Morse potential fulfilling the condition in Remark 3 fulfill the condition, since

[TABLE]

This expression is negative for all $z\in[0,\infty]$, if $F<1$ and $1-FL^{d+1}>0.$

For the quadratic potential $U_{R}$ described in the first section in Remark 1 the above formulas give the following expressions for $\alpha_{i},\beta_{i}^{\prime}$:

[TABLE]

with

[TABLE]

for $z\in[0,1]$ and $A(z)=1,z\geq 1$. Moreover,

[TABLE]

with

[TABLE]

for $z\in[0,1]$ and $B^{\prime}(z)=1,z\geq 1$. We note that $B^{\prime}\geq 0$. Moreover, $\alpha_{i}^{\prime}=\beta_{i}=0$.

In 1D this gives

[TABLE]

and the following expressions

[TABLE]

and

[TABLE]

for $z\in[0,1]$. In 2-D we have

[TABLE]

and the following expressions

[TABLE]

and

[TABLE]

for $z\in[0,1]$.

Remark 9**.**

Considering the 1-D case, the potential in Example 2, i.e.

[TABLE]

for $-R<x<0$ and [math] otherwise, leads to the expressions

[TABLE]

with

[TABLE]

and $A(z)=1$ for $z\geq 1$ and

[TABLE]

with

[TABLE]

and $A^{\prime}(z)=1$ for $z\geq 1$. We note that $A^{\prime}\geq 0$.

4 Numerical results

In this section we present a series of numerical experiments for the hydrodynamic (7) and scalar equations (8) and their localized approximations (9) and (10). In particular, we will investigate and compare the different schemes for the non-local equations for situations near the local limit. One observes the following: if a well resolved situation with a large number of particles compared to the interaction radius is considered, then the results of the multi-scale and the naive or microscopic approximation coincide. In this case the multiscale method behaves similiar to a microscopic (DEM-type) simulation. If, however, a strongly underresolved situation is considered, i.e. the number of grid particles is small compared to the interaction radius, then the numerical solution using the microscopic approximation of the integral deviates strongly from the correct solution. On the contrary, using the multi-scale method with the correction factors introduced above, a good approximation of the correct solution is obtained. In this case the multi-scale method is essentially a meshfree numerical method for macroscopic equations. Since all methods require approximately the same amount of computation time per grid-particle, the above observations can be rephrased as follows: for situations with relatively small interaction radius we obtain, using the multiscale method, a reduction in computation time by several orders of magnitude. On the other hand, for relatively large interaction radius, the computation times of the microscopic approximation and of the multiscale method are similar.

4.1 Numerical results in 1-D

4.1.1 Test case (Example 3)

We study 1-D movement of the particles under the influence of a confining potential. The equations are described in Example 3. More precisely, we look at equations (11) and its limit equation (12). The interaction potential is given by (3). We consider the case of small $\epsilon<<1$. The porous media case $\epsilon=1$ has been treated in Ref. [39]. We choose $V(x)=|x|^{2}/2$ and $\gamma=1$. We choose different values for the interaction radius $R$. For this example the particles are assumed to be distributed on a fixed equidistant grid, i.e. we use here an Eulerian approach. Situations with arbitrary particle locations and a Lagrangian approach are considered in the following subsections. The initial density is given by $\rho(x,0)=1$ for $x\in[-0.5,0.5]$ and [math] elsewhere. As time proceeds the density converges to a stationary solution depending on the value of $R$. For the localized limit equation this stationary solution is for small values of $\epsilon$ approximated by $\rho_{\infty}=-x$ for $-\sqrt{2}<x<0$ and $\rho_{\infty}=0$ otherwise, such that $\int\rho_{\infty}(x)dx=1$. In order to obtain the steady state solution we have solved the equation until $t=8$.

In Figure 1 we plot the solutions for a fixed number of particles $N=1600$ and different values of $R$, considering in this way the different algorithms in the well-resolved, in intermediate cases and in the underresolved case. We plot the limiting Fokker-Planck solution and the microscopic and multi-scale solutions. In the well resolved case, i.e. here for larger values of $R$, the microscopic approximation of the integral and the multi-scale solution give similiar results deviating from the solution of the limiting Fokker-Planck equation. For the under-resolved case, for smaller values of $R$, the microscopic approximation does not give the correct results as discussed above, in contrast to the the multi-scale method. In particular, one observes that using the microscopic discretization for an underresolved situation, the numerical evolution terminates before the stationary state is reached, if the compact supports of the interaction kernels of the finite number of particles do not overlap any more.

Finally, we investigate the convergence rate of the method numerically. We consider a situation near the localized limit with $R=0.02$, where the limit equation gives the correct solution up to order $10^{-2}$ and is used as the reference solution. We compute the error of the particle methods for the scalar model (11) comparing the numerical results to the numerically determined stationary solution. The scalar problem is solved with microscopic and multi-scale approximation of the interaction term. In Figure 2 and Table 1 the $\mathcal{L}^{2}$-errors versus the number of particles are plotted. One observes the deterioration of the method using the microscopic approximation of the interaction term for smaller numbers of particles. The multi-scale method is able to treat all ranges with a similiar accuracy. We remark that the computation time for the different methods is approximately the same for the same number of particles. Looking at Table 1 we observe that the multiscale method with $N=200$ particles yields for the present example the same error as the microscopic method with $N=1600$ particles. Thus, the computation time for the naive, microscopic method is approximately an order of magnitude larger than the time for the multi-scale method.

4.1.2 Non-local traffic flow (Example 4)

We consider equations (14) and (16) with different initial conditions and different choices of the interaction potential. For this example, we consider a fully Lagrangian approach. First we consider a symmetric potential as in Example 1 and initial conditions leading to a rarefaction wave solution. We choose

[TABLE]

The solution of the limit equation (16) is a rarefaction wave

[TABLE]

We choose and $\delta=0.1$. In Figure 3 the rarefaction solution is compared at time $t=2$ for the scalar model (16) and for (14) with microscopic interaction approximation and multi-scale approximation. We use a fixed number of particles $N=800$ and an interaction radius ranging from $R=0.002$ to $R=0.4$. One observes a good coincidence of microscopic and multi-scale approximation for large $R$ and a stronger deviation, the smaller the value of $R$ is chosen. In this situation the influence of larger values of $R$ on the exact solution is a small increase of the smearing of the solution. We note that in this situation using a one-sided downwind interaction potential as in Example 2 and $\delta=0$ gives similar results.

Second we consider a shock solution of the Riemann problem for equation (14) and (16). The initial values are now chosen as

[TABLE]

The solution of the limit equation (16) is a shock

[TABLE]

We consider first the case of a one-sided downwind potential chosen as in Example 2 and chose $\delta=0$. In Figure 4 we compare again different values of $R$ for a fixed number $N=800$ of particles. We observe, that in the well resolved case for larger values of $R$ the microscopic approximation and the multi-scale solution give a smeared out shock solution deviating from the solution of the limiting equation which is a moving shock. For the under-resolved case with smaller values of $R$, the microscopic approximation does not give the correct results, in particular, the speed of the computed wave is wrong and equal to the one of the advection problem obtained after setting $U_{R}\star\rho$ in (14) equal to [math]. This is in contrast to the multi-scale method which computes the correct wave speed.

For the investigation of the convergence of the method we consider a situation with $R=0.2$. The scalar problem is solved with microscopic and multi-scale approximation of the interaction term. We compute the error of the different methods for different numbers of grid-particles comparing the numerical results to the numerical solution of (14) determined with a fine grid with $N=6400$. For  $6400$ particles the difference between multi-scale and naive method is again of the order $10^{-2}$. In Figure 5 and Table 2 the $\mathcal{L}^{2}$-errors versus the number of particles are plotted. One observes again the deterioration of the method using the microscopic approximation of the interaction term for smaller numbers of particles. The multi-scale method is able to treat all ranges with a reasonable accuracy for the first order method.

Similiar to the example in the last subsection comparable errors are obtained using $N=200$ particles and a computation time of $8s$ for the multi-scale method and $N=1600$ particles and a computation time of $77s$ for the microscopic method. This yields again an order of magnitude gain in computation time.

Finally we investigate the above example with a symmetric interaction potential and a regularization $\delta=0.02$. In Figure 6 we compare again different values of $R$ for a fixed number $N=800$ of particles. We observe, as before that in the well resolved case for larger values of $R$ the microscopic approximation and the multi-scale solution give an oscillating solution deviating from the moving shock solution of the limit equation. For the under-resolved case with smaller values of $R$, the microscopic approximation does not give the correct results as in the case of the downwind potential, whereas the multi-scale method computes the correct solution.

For the investigation of the convergence of the method we consider as before $R=0.2$. We compute the error of the particle method for different numbers of grid-particles comparing the numerical results to the numerical solution of (14) determined with a fine grid with $N=6400$. The scalar problem is solved with microscopic and multi-scale approximation of the interaction term. In Figure 7 and Table 3 the $\mathcal{L}^{2}$-errors versus the number of particles are plotted. One observes a similiar behaviour as in the non-symmetric potential case. We remark that the method developed here is still working in the present case, since we have regularized the equations with $\delta>0$ obtaining a dissipative limit. For $\delta$ much smaller, the limit would be dominated by dispersion and the method is not supposed to work more efficienly than the naive, microscopic approximation.

4.2 Comparison of numerical algorithms in 2D

In this section we solve the hydrodynamic system (18) using the multi-scale algorithm and the algorithm with naive, microscopic evaluation of the integral in 2D. We compare the results with each other and with the local hydrodynamic problem (20). We consider a fully Lagrangian approach.

4.2.1 Hydrodynamic system for test case

We consider a two dimensional domain $[0,45]\times[0,50]$. We initially generate particles in $[0,30]\times[0,50]$. We simplify the hydrodynamic system (18) considering a fixed vector $\hat{e}(x)=(1,0)^{T}.$ The maximal velocity is $u_{max}=1$ and the maximal density is chosen as $\rho_{max}=1$. A symmetric interaction potential and $\delta=0.1$ are chosen. Moreover, $\alpha=1,\gamma=1$. We use $N=4150$ particles and the interaction radius $R$ equal to $0.2$. The initial condition is given by

[TABLE]

In Figure 8 a.) the solutions of the nonlocal equations (18) with microscopic and multi-scale interaction aproximation for the interaction radius $R=0.2$ are plotted. Moreover, the time development of the normalized total mass in the domain $x\in[0,100]$ is plotted for microscopic and multi-scale interaction approximation in Figure 8 b.).

In Table 4 we compute the error and the CPU times of naive and multi-scale method for different numbers of grid-particles, see also Figure 9. The reference solution is computed by using a spacing of $\Delta x=0.11$ and approximately $136000$ particles. For this very fine resolved case the difference of naive and multi-scale solution is approximately equal to $10^{-2}$. One observes in this physical situation a smaller but still relevant gain in computation time.

4.2.2 Hydrodynamic systems for pedestrian flow (Example 5)

In our final example we investigate the hydrodynamic systems (18) and (20). Compared to the last subsection we consider a more complicated configuration defined in Ref. [40], compare also ([22]), and a full coupling to the eikonal equation.

Pedestrians are initialized on the left of the domain and evacuated towards the exits on the right as shown in Figure 10. As initial value we choose a constant value of $\rho=1$ in the region $[0,30]\times[0,50]$. In the center of the computational domain an obstacle is located. For the eikonal equation we use $\phi=0$ on the two exits and $\phi=\infty$ on all walls as boundary conditions.

We choose the maximum velocity $u_{max}=2$ and the maximum density $\rho_{max}=10$. In this case, we vary the initial average distance between grid points from $0.15$ to $1$, i.e. the number of grid particles varies between $1440$ and $61800$. Moreover, we choose the following parameter: $\gamma=500,\alpha=1000$.

First we consider an underresolved situation with relatively small value $R=0.2$ and $\Delta x=0.5$, i.e a number of particles of approximately $N=5650$. In Figure 11 we plot the solution at fixed time using the localized equation and the nonlocal equation with multiscale method and microscopic integration. In this case the solution computed via the multi-scale method and the one computed from the localized equation coincide, whereas the microscopic method gives strongly different results. In Figure 12 we show a comparison of solutions obtained from the microscopic scheme with decreasing discretization sizes or increasing number of particles.

Figure 13 considers a well-resolved case with $\Delta x=0.2$ and $R=0.4$. In this case we observe good coincidence of the solutions of microscopic an multi-scale approximations.

Finally, in Table 5 and Figure 15 we compute the error and the CPU times of naive and multi-scale method for different numbers of grid-particles. The errors are determined along a line with $y=37$ and $x\in[25,55]$. The relative $\mathcal{L}^{2}$-errors are given as well as the computation times in minutes. The reference solution is computed by using a spacing of $\Delta x=0.15$ and approximately $62000$ particles. For this fine resolved case the difference of naive and multi-scale solution is approximately equal to $10^{-2}$. Looking at Table 5 and the multiscale error with $1400$ particles and the naive error with $35200$ particles one observes in this more complex situation a gain in computation time by more than an order of magnitude. comparing the naive computation with $35200$ particles with the multi-scale simulation with $5700$ particles there is still a gain of an order of magnitude.

5 Conclusion and Outlook

We have extended a multi-scale meshfree particle method for macroscopic mean field approximations working uniformly for a large range of interaction parameters $R$. The well resolved case for large $R$ is treated as well as underresolved situations for small values of $R$. The method can be considered as a numerical transition from a microscopic system simulation to a macroscopic averaged simulation. Application of the method are shown for pedestrian flow simulations. The potential gain in computation time compared to a microscopic simulation is depending on the situation under consideration wit a potential gain of several orders of magnitude. In the situations considered in this paper we have obtained a gain in computation time of approximately one order of magnitude. Situations with more complex geometries or moving boundaries as well as the investigation of attractive-repulsive potentials will be treated in future work, compare [41].

Acknowledgment

This work is supported by the German research foundation, DFG grant KL 1105/20-1.