Single to Double Mill Small Noise Transition via Semi-Lagrangian Finite Volume Methods

Single to Double Mill Small Noise Transition via Semi-Lagrangian Finite Volume Methods


We show that double mills are more stable than single mills under stochastic perturbations in swarming dynamic models with basic attraction-repulsion mechanisms. In order to analyse accurately this fact, we will present a numerical technique for solving kinetic mean field equations for swarming dynamics. Numerical solutions of these equations for different sets of parameters will be presented and compared to microscopic and macroscopic results. As a consequence, we numerically observe a phase transition diagram in term of the stochastic noise going from single to double mill for small stochasticity fading gradually to disordered states when the noise strength gets larger. This bifurcation diagram at the inhomogeneous kinetic level is shown by carefully computing the distribution function in velocity space.

1 Introduction

In the last decade, the theoretical and mathematical biology communities have paid a great deal of attention to explain large scale structures in animal groups. Coherent structures appearing from seemingly direct interactions between individuals have been reported in many different species of fish, birds, and insects [40, 53, 50, 7, 23, 10] and many others, see the reviews [17, 12, 44] for more literature in the subject. There are two type of patterns appearing regularly for different species called flocking and milling behavior. In the flocking behavior, individuals agree in moving in certain direction with some spatial shape that changes in time due to the effect of wind, hydrodynamic advantage, own desire, predator presence, or roosting behavior. In the milling case, individuals organise themselves rotating around certain location, that may move in time, forming a thick group in an annular/ring region.

Many individual based (particle) models have been proposed to explain or reproduce these patterns. In general, swarming dynamics have been described by systems of interacting ordinary or stochastic differential equations: see [40, 33, 23] for the biological aspects of the modeling of swarm dynamics and [47, 28, 3, 4, 45, 37, 2, 39, 48, 49, 24] for including effects like the interaction between individuals, self propulsion, roosting, and stochastic forces in these equations. All of these models start by including three basic mechanisms: attraction, repulsion, and reorientation. The form in which each of these effects is including in the modelling depends on the species and particular biological parameters.

Dealing with large particle systems is cumbersome and even prohibitive numerically if the system is very large. Therefore, many authors have proposed to use a mean field approach in which one can derive effective PDEs of Fokker-Planck or Vlasov type [26, 36, 34, 13, 16, 11, 18, 13] from these microscopic systems of equations. As usual in kinetic theory, scaling arguments on these kinetic PDEs lead to macroscopic approximations of the mean field equations as in [13, 36, 18, 51, 25].

A very interesting question that arises naturally is the stability of these patterns under perturbations. This question has been recently analysed in [5, 1, 14] both for flocks and rings in terms of initial data for the model introduced in [28]. However, an stochastic perturbation also leads to very interesting phenomena called phase transition. This phase transition has been reported for the first time in swarming dynamics in the so-called Vicsek model [62], and it has been studied in very detailed at the macroscopic level in [25]. Essentially, in these works they analyse how the system undergoes a transition from an ordered state (pattern) to a disordered state (chaos) by increasing the noise strength in the system.

The aim of our work is twofold. On the one hand, we focus on a very accurate numerical solution technique for the kinetic mean field equations in swarming dynamics. We will describe a splitting scheme with a Semi-Lagrangian solver in space and a semi-implicit finite-volume scheme in velocity space. On the other hand, we make use of this accurate numerical tool to show for the first time, up to our knowledge, a phase transition for an inhomogeneous mean field equations in swarming dynamics. By varying the amplitude of the stochastic forces, we show the influence on milling patterns in the model without noise introduced in [28]. The results of the microscopic, kinetic mean-field, and macroscopic equations are compared and discussed thoroughly. The numerical scheme is an improvement of the schemes used in [56, 43] and of the FVM solver described in [63, 64]. As a conclusion of our study, we have discovered that double mills are quite robust and stable for small stochasticity. Mill solutions immediately leads to double mills solutions under small noise, then fading gradually toward disordered states for large noise strength. This single to double mill noise induced phase transition is the main theoretical biology implication of this work.

The paper is organized as follows. In Section 2, we quickly summarize the microscopic, kinetic mean field, and macroscopic equations derived for the model introduced in [28]. It contains the discussion of different regimes and associated macroscopic approximations ranging from a situation with zero stochastic force, to an intermediate case, where stochastic and propulsion force balance each other, leading finally to the dominating force case. In Section 3, the numerical scheme is described in details and an investigation on the order of convergence of the scheme is presented both for the homogeneous and inhomogeneous cases. Finally, Section 4 show the phase transition of single to double mills for small noise strength and the phase transtion toward disordered state for larger stochastic force. We first conduct a detailed comparison of microscopic and kinetic simulations for both single and double mills with and without noise for showing the accurate comparison of the models.

2 Model Hierarchy: Microscopic, Mean-Field, and Macroscopic Models

We consider classical models for swarming dynamics, which include forces resulting from social interaction between individuals, self-propulsion, and friction as in [28]. The combined effect of self-propulsion and friction is to impose an asymptotic speed for individuals. This means that individuals travel at a typical cruise speed asymptotically, similar to other classical models [62, 25]. Furthermore, a stochastic term given by white noise is included as in [13, 17] to account for random and small error interactions. Moreover, a term describing roosting behaviour is added as discussed and introduced in [18]. The microscopic equations are given by


where . is a pairwise radial interaction potential. A classical example is given by the Morse potential , with attraction/repulsion strengths / and radius of interaction /. Other examples for interaction potentials are considered in [19, 15] and the references therein. and are the self propulsion parameter and is a potential defining the roosting force. The parameter describes the amplitude of the stochastic force. The main purposes of the paper are to propose an accurate numerical deterministic scheme for the mean field approximation of this microscopic system and to study numerically the influence of the noise parameter on particular patterns present in this model. In fact, the microscopic model without noise has been shown to be paradigmatic for self-organization, since it exhibits rich dynamical patterns such as flocks, single and double mills, see [28, 22, 13, 18]. The main issue in the present paper is to analyse the qualitative change of the pattern as noise increases. One may expect some of these collective motions to survive for small noise up to some critical value for which the introduced stochasticity destroys the self-organized pattern.

A classical derivation procedure as described, for example, in [13, 9, 52, 58, 32, 27, 57, 6] yields the mean field equations. For the one-particle distribution function and the density

with the normalization condition , one obtains






Different regimes and approximate solutions for macroscopic quantities can be obtained from the above mean field equation using different scaling assumptions ranging from very small to very large stochastic force .

2.1 Vanishing stochastic force

The deterministic case has been treated in detail in many publications. A special feature is the appearance of single and double mill solutions, we refer to [13]. Single mills of the mean field equation

where is given by (4), can be found as follows. We look for a mono-kinetic solution, compare [13, 22]

and are found by integrating against and and closing the equations with the above Ansatz. One obtains

and the momentum equation

on the support of the density . Stationary distributions can be found in the following way. Assuming , we obtain from the hydrodynamic equations above that

on the support of the density . Assuming a rotatory solution given by

compare [13], and looking for radial densities , we obtain that the continuity equation is trivially satisfied and since


Assuming we end up with an integral equation for :


in the support of the density . A numerical investigation of these stationary states of the hydrodynamic equations, called single mills, is performed in [18].

A so called double mill solution can also be derived from an appropriate Ansatz with two opposing kinetic velocities, see [13]. In short, given the spatial profile of a single mill, there is always a double mill solution with the same density but changing the sign of the velocity to part of it. Let us finally remark that in particle simulations of the microscopic system (1)-(2), double mills are observed very often by peforming large perturbations of single mills when . Moreover, these double mills are typically interlaced instead of overlapping on the same spatial region. This means that the velocity of one mill is not totally opposite of the velocity of the other mill at the same spatial point.

2.2 Balance of stochastic and propulsion forces

Assume that and are large for long time scales. The corresponding scaling yields the following equations



To zeroth order in an expansion in the small parameter , we have . The kernel of is given by the solution of

A simple computation gives that the local equilibrium solution reads as

with a normalizing constant and a given spatial density .

Remark \thetheorem

For going to zero goes to

For going to infinity converges to a standard Maxwellian .

To compute one has to proceed to first order. One obtains

This equation is solvable since Then, it is not difficult to check that

with the matrices given by

and the vector fields

Integrating the scaled equation and using the above gives


2.3 Dominating stochastic force

We refer to [13, 18, 8, 31, 38] and consider a scaling with large and long time scales. This yields

A straightforward asymptotic expansion gives to zeroth order where is the standard Maxwellian and

The solution of the stationary problem is given by


with prescribing the total mass. The associated momentum is

To summarize this overview, we expect the stationary solutions of the mean field equations to change as follows when goes from to infinity. For and well prepared inital conditions, one obtains a single mill solution, i.e. a solution in velocity and a radial solution in of equation (6). Double mills appear depending on the initial conditions and typically for large perturbations of the mill solution. As becomes very large, one approaches a Maxwellian solution in velocity and a solution of equation (8) in . These two behaviors will be confirmed by the numerical investigations in Section 4.

On the other hand, the transition from the mill ordered state to the disordered state by increasing the noise is not clear at this point. The macroscopic equations in the middle regime seem to indicate a density profile near a non-isotropic maxwellian due to the matrices in (7). However, the change in the spatial and velocity profile are not clear for small noise when approaching the monokinetic solution. We wil investigate this issue in detail in Section 4 showing that double mills play an important role in this transition.

3 Numerical Scheme

3.1 Discretisation

In the following, we present a second order Semi-Lagrange Finite Volume method for the Fokker-Planck equations. We describe the method in 2-D, see [41] for similar situations in fibre dynamics. An extension to 3-D is straightforward, but computationally expensive, see [42] for similar problems. We split (3) using a Strang splitting [21], that means we consider


where and are given by (4) and (5) respectively. The system is discretised on the grid , where for the spatial domain and are points in the velocity domain. Furthermore, we denote by the discretisation points in time with constant step size . For (10) we use a second order Semi-Lagrange method (see [35, 43, 55]), where we apply an interpolation procedure using cubic Bezier curves. Proceeding in the usual way, we look back in time along the characteristic curves starting in the respective grid point and obtain a grid value for the current time step by interpolating from the old grid values at the characteristic endpoint. On each cell of the grid in space, the solution at a given time is given by the polynomial reconstruction


with the Bernstein polynomials and the 16 control points . We used the notations here. Due to the nature of the Bernstein polynomials, the interpolant never leaves the convex hull of the control points. They are chosen appropriately so that we have a third order local approximation for smooth data and a reduced order interpolant without oscillations for non-smooth data.

In order to obtain a cubic Bezier polynomial like (11) with appropriate control values , one computes an interpolating Newton polynomial and then performs a basis transformation from the Newton basis to the Bernstein polynomials. We will demonstrate the procedure for a 1D problem. Assume we have the stencil . Then, for and , we can write the Newton polynomial as follows

where we obtain the coefficients

Figure 1: Convergence of the Semi-Lagrange method with Bezier interpolation for linear advection as in (10). One obtains the predicted third order method.

In the same way, we can write the Standard polynomial basis, the Bernstein polynomial basis, the Bezier coefficients and the Bezier polynomial:

Now a change of basis gives

In order to achieve the right interpolation order, we want the Bezier polynomial to be equal to the Newton polynomial, i.e. . That can be achieved if we set

Since the Bezier polynomial with these control values is the same as the Newton polynomial, we know that we have a local approximation error of fourth order for smooth data, resulting in a third order consistency error for linear advection in the Semi-Lagrange scheme, see figure 1. However, at discontinuities, oscillations will occur. Avoiding them is straight forward, due to the convex hull property. If one of the leaves the interval of neighbouring grid values, we put it back to either of the boundaries or . We note that the mass conservation is destroyed by such a procedure. Thus, additionally a mass conservation procedure as in [43, 29] has been implemented. In literature, similar examples of Semi Lagrange schemes and splitting methods for the Vlasov equation are known, see for example [20, 30, 55, 35, 61]. The contributions mostly differ in the interpolation scheme and the handling of oscillatory solutions. For the advection problem (10), our procedure can be easily extended to 2D by extending the above described vectors in an appropriate way.

Figure 2: A patch of the rectangular grid for FVM.

For the discretization in velocity, we apply a finite volume scheme [46] of second order with a Lax Wendroff type approximation [54] for the advection, which is an enhancement of the methods used in [56, 64]: Denote the midpoints of rectangular cells for (3) by . The boundary between cells and is called with the normal and the boundary midpoint . We further denote the distance between cell midpoints and by , which decomposes into , where are the distances to the boundary midpoint . Since we have an equidistant, rectangular grid in and -direction, is either or and are respectively, see Figure 2. is the set of the indices of neighbour cells to cell . For the FVM nethod, one computes the cell averages over the cells by integrating (9). This gives

where the notation is used. Rewriting

and applying the midpoint rule for the time integration and the divergence theorem for the integration w.r.t , one gets


Now, we have to approximate the quantities and . If we use

we end up with a second order Lax-Wendroff flux function. As any higher order flux function, it will cause oscillations for non-smooth solutions and one has to apply a limiter scheme to remedy this. In the present example, we will mix the Lax-Wendroff flux with a first order upwind advection for non-smooth solutions:


and measures the smoothness of the solution. is the van-Leer limiter [60, 59]

Using this limiting procedure, oscillations are avoided. In order to get a value for , one can expand around in edge normal direction and obtains

for the equidistant grids used here. Together with the fluxes from the other , the second order terms vanish and one obtains an estimate of oder . In total, multiplying by and dividing by , we end up with a second order approximation on the velocity grid. We have to note however, that the discretization is not second order in time, since we did not account for the evaluation of at . This can be remedied by applying the trapezoidal rule in (12) instead of the midpoint rule for the time integral of the diffusion terms, which would lead to a Crank-Nicholson scheme for the diffusive part of the problem. The convergence rate in time does not influence the numerical results that much, since the overall error is still very near to a second order scheme: For a numerical confirmation of the convergence rate of the velocity discretization in the case without interaction and roosting, see Figure 4.

Figure 3: On the left, the diffusion coefficient is set to , on the right, it is . From top to bottom we have: Analytical equilibrium and numerical equilibrium . For , will approach .

In matrix form (12) can be rewritten as


with matrices for the transport coefficients and for the diffusion coefficients. For large values of the diffusion parameter () one will need a very small time step, if the above method (13) is to converge. So instead of (13), we choose the semi-implicit ansatz



can be proven to be strictly diagonally dominant and has positive diagonal entries. We use the conjugate gradient method for solving system (14). For the implementation, it is unavoidable to use some sort of sparse format for the matrices, which decreases computing time by a large factor, since we have to solve this system at every point in every time step.

Figure 4: Logarithmic error plot for different grid sizes for the finite volume method in the space homogeneous case. The coefficient was chosen as . The convergence is second order.

3.2 Numerical investigation of the space homogeneous case

First we test our finite volume scheme for the space homogeneous equation without interaction and roosting force. The equation reads


For the computations we use the numerical values and . The analytical equilibrium solution is given by

with and . For the solution of the time dependent Fokker-Planck equation (15) converges to this equilibrium state, compare figure 3. The convergence to equilibrium is investigated looking at the distance between numerical and analytical equilibrium in . The values of this functional for different grid sizes are displayed in Figure 4.

3.3 Numerical convergence analysis for the full problem

Up to now, we only tested the numerical scheme for spatial and velocity domain separately. In this section a numerical convergence analysis is performed for the two regimes of large diffusion and no diffusion , for the full problem in space and velocity. Starting from a Single Mill initial condition in all cases, we will compare the stationary distributions for different grid sizes. Let be the density of the numerically computed stationary distribution function and be the reference density. To measure the errors we compute the distance


The reference solutions are given on a finer grid than the numerical solutions. We will interpolate the numerical data with a 3rd order accurate procedure to the reference grid and compute the distance there. The occurring integrals will then be evaluated with the midpoint rule on each grid cell.


In this case, we do not have an analytical equilibrium. However, a very good approximation of the stationary density can be obtained by taking the solution of the fixed-point equation

Figure 5: The top row contains the original data of the numerical solution for , below there is the interpolated solution, which is then compared to . The colour scaling is the same as for the reference density in figure 7.

This equation has been obtained in the previous section in the limit . For the stationary density is well approximated by the limit density. The computational domain is given by . Table 1 shows the number of grid points in each of the 4 dimensions, the grid size in the spatial domain, and the error.

Figure 6: The blue line shows the second order accuracy of the numerical scheme in for . The green and red lines are references for convergence rates of first or second order, respectively.
Figure 7: On the left: radial plots of the numerical densities for the different grid sizes and the reference density for . On the right: the reference solution on its native grid.
15 22 30 45
8.5714 5.7143 4.1379 2.7273
0.0049827 0.0044338 0.0021669 0.0008489
Table 1: Grid size and error for .

The step sizes in the velocity domain are decreased in the same way as the spatial grid for the convergence plots. In fact, we use as many grid points for the spatial as for the velocity grid. The respective densities can be found in figure 5. We get the logarithmic error plots in figure 6. Radial plots for the densities and the reference density can be found in figure 7. Note, that the colour scaling is the same for reference density in figure 7 and the densities in figure 5.


In this case without diffusion, we obtain a single mill solution as equilibrium solution. Since there is no analytical solution to compare, we take as a reference solution the numerical solution with a fine resolution. In this investigation, we use the solution with 60 grid points in each of the 4 directions as the reference solution. The computational domain will be given by . We denote the density of the fine resolved numerical solution by and compare them to the densities obtained from the lower resolutions, in the norm stated in (16). Table 2 shows the number of grid points in each of the 4 dimensions, the grid size in the spatial domain, and the error.

Figure 8: The top row contains the original data of the numerical solution for , below there is the interpolated solution, which is then compared to . The color scaling is the same as for the reference density in figure 10.
Figure 9: The blue line shows the errors of the numerical scheme in for . The green and red lines show the errors corresponding to a convergence rate of first or second order, respectively. Since the error line is parallel to the red line, we conclude that the numerical scheme is of second order for .
15 22 30 45
7.1429 4.7619 3.4483 2.2727
0.0075987 0.0045355 0.0025808 0.0009124
Table 2: Grid size and error for .

The respective densities can be found in figure 8. From this, we get the logarithmic error plot in figure 9. Radial plots for the densities and the reference solution can be found in figure 10.

Figure 10: On the left: radial plots of the numerical densities for the different grid sizes and the reference density for . On the right: the reference solution on its native grid.

4 Numerical investigation of milling States

In this section we look at numerical solutions for the full problem (3) and compare them to the microsopic solutions of equation (1) and the macroscopic approximations derived in the last section. For we use the same numerical values as in the last section. Moreover, the constants of the interaction potential are chosen as in [18] as . Depending on the initial conditions special stationary solutions, so called single or double mills, appear for diffusion coefficient . They are obtained as stationary states of the interacting particle system as well as the kinetic equation. We discuss these stationary states numerically and concentrate on the investigation of the behaviour of these solutions for increasing noise parameter .

4.1 Milling solutions

The single mill solution in the microscopic context is a structure, where all the particles with positions go around the origin with velocities

In the kinetic case we display the spatial distribution of the density . The velocities are visualized by looking at the velocity distributions at fixed spatial points . If a single mill is obtained, will become concentrated at

For comparison, we generate a distribution function from the microscopic results, by counting particles in a 4D-histogram, which is then normalized to mass 1 as for the kinetic solution. We note, that a large amount of particles is needed to generate a reasonably smooth histogram. In figure 11, we plotted 400 microscopic particle positions with arrows indicating their velocities. The actual computation was carried out with 13000 interacting particles, where we have chosen 400 particles randomly for plotting the results. The histogram is based on the position and velocity data of all 13000 particles. The spatial density obtained from the kinetic equation and pictures of the velocity distributions at corresponding fixed points obtained from microscopic and kinetic solution are also shown in figure 11. The results coincide well, although the peaks in the kinetic velocity distribution are broader than the microscopic ones, due to the diffusivity of the kinetic solver.

Figure 11: Comparison of single mills from microscopic equations (top) and from kinetic equations (bottom). To the left, we have the particle positions and kinetic particle distribution, to the right, there are the velocity distributions at fixed points on the grid. Note, that the grid in velocity space is different for the microscopic and kinetic case, so the positions do not coincide exactly.

The single mill is obtained for for a specific type of initial conditions. One has to prepare the system in a state which is not very far from the single mill solution. In the kinetic case, this can be realised by using the initial condition

is chosen in such a way, that . One can achieve a corresponding initial condition in the microscopic case by placing particles randomly at positions with and assigning random velocities with the corresponding deviation from the direction .

Figure 12: Comparison of double mills from microscopic equations (top) and from kinetic equations (bottom). To the left, we have the particle positions and kinetic particle distribution, to the right, there are the velocity distributions at fixed points on the grid. Note, that the grid in velocity space is different for the microscopic and kinetic case, so the positions do not coincide exactly.

A double mill is, from the microscopic point of view, a situation, where one part of the particles at the positions is circulating around the origin with speed

while the other part at around the same location is going in the opposite direction . The kinetic spatial density looks like the one for the single mill solution, but the velocity distributions at fixed points show two symmetric peaks in at the sphere . Computing a histogram from the microscopic solution as before, we compare the results in figure 12. Again, the two solutions agree reasonably well. Note, that it is not possible to obtain a double mill with the macroscopic model because of the monokinetic closure, see [13].

The double mill is again obtained for by choosing appropriate initial conditions. For example, one may use the following distribution function for the kinetic equation:

is chosen as before. The corresponding initial values in the microscopic context is obtained by placing particles randomly in the same region in space and assign random velocity vectors , which are limited by . Note, that the initial condition to get a double mill pattern is less restrictive than the one for the single mill. This indicates that a double mill is in a certain sense more stable or robust than a single mill for small stochasticity. This fact will become clearer investigating situations with increasing noise parameter in the next subsection.

Figure 13: Solutions of the microscopic system with 400 particles: At the top on the left, one can see a single mill for , on the right, we have a double mill for . Below, is raised first to 0.123, where we still have a double mill until we end up in an unordered state at .

4.2 Influence of the noise amplitude: Phase Transition

In this section the microscopic and kinetic equations are investigated for different values of starting with a single and a double mill initial condition. Figure 13 shows the behaviour of the microscopic system for three values of the noise coefficient: For , we obtain the expected single mill or double mill depending on the initial conditions, then for the stationary state is still a double mill independent of the initial state being a single or a double mill initial distribution. For , we have an unordered state.

Figure 14: Solutions of the mean field equation: On the left, we have the spatial density and on the right, we have the velocity distribution at some fixed positions in the spatial domain. From top to bottom there are results for different values of the diffusion coefficient . The constants of the interaction potential are chosen according to [18]: .

The crucial point is, that starting from the single mill state for , there is a transition from single to double mills for small values of , before milling structures vanish as increases. This can be observed for the solutions of the kinetic equation as well. Figure 14 shows spatial density and velocity distributions at fixed points in space, for several values of for single mill initial conditions. The evaluation points in space are indicated above the respective velocity distribution. For one observes the milling solutions. In this case a -type velocity distribution is observed for the fixed spatial points and a circular density distribution with zero density in the center, as before. However, yields a double mill type stationary state, as can be seen by comparing the result to the double mill for in figure 12. Further increase of yields the disappearance of milling structures.

Moreover, for we compare the microscopic and the kinetic solutions. We generate again a histogram from the particle data and give the comparison to the kinetic result in figure 15. Again, the solutions match, with a slightly higher diffusivity in the kinetic result due to the numerical scheme.

Figure 15: Comparison of results from kinetic equation (top) and from microscopic equations (bottom). In each case, a single mill initial condition was used and the diffusion/noise parameter was set to .

Additionally, we plot in Figure 16 the velocity distribution functions at one fixed spatial point for a wider range of values of , with separate plots for single mill and double mill initial conditions. One can see, that the single mill state fades into a double mill and then into a normal distribution. The double mill however persists for small values of .

Figure 16: Here, we look at at one of the fixed spatial points. On the left, we have a single mill, which fades away gradually with increasing diffusion coefficient . On the right, we have the same situation starting from a double mill.

In order to examine this behaviour more closely, we computed for each numerical solution of the kinetic equation the distance to a single mill


Figure 17 shows the values of this functional depending on the diffusion parameter for initial conditions yielding either single or double mills. We can observe the presumed behaviour, where the double mill persists for some range of and then fades away as the single mill does.

Figure 17: Distance of the solution from the single mill solution (17) for different values of . The blue line is obtained for single mill initial condition, the green one for a double mill initial condition.
Figure 18: Comparison of spatial densities from microscopic and kinetic equations. We show radial plots of them for different values of : From left to right, we have , and .

In order to further validate our results for the kinetic equation, we have also compared the spatial densities of microscopic, kinetic and macroscopic equations. In figure 18 one can find radial plots of the spatial densities of the mean-field equation and the microscopic system for different values of . For we have the above-mentioned circular distribution with zero density in the center, which fades away for until there is a normal distribution for and higher values of . Furthermore, we compare our microscopic and kinetic solutions to the macroscopic single mill result computed in [18] for in figure 19. We get good agreement of both the spatial density and the tangential moment . The mean-field solution is a bit more diffusive than the corresponding microscopic and macroscopic ones. We note again that the macroscopic solution is not able to reproduce the double mill solution, since a mono-kinetic closure is used in section 2.1. For large the radial density of the mean field solution is shown together with the corresponding macroscopic solution of the diffusion problem from section 2.3 and the microscopic solution in figure 20.

Figure 19: On the left: Radial plot of the densities obtained from microscopic, mean field and macroscopic single mill solution for . On the right: radial plot of the corresponding tangential moment.
Figure 20: Radial plot of the densities obtained from microscopic, mean field and macroscopic solution for .

4.3 Computational remarks

All the computations were carried out on the cluster of the University of Kaiserslautern on an Intel Xeon E5 2670 with eight cores at 2.6 GHz each. The solver for the Fokker-Planck equation is written in Fortran and with OpenMP parallelisation. We note that for the mean field equation, we only need to compute distances between points of the grid to determine the influence of the interactions. Since the grid is fixed, these distances can be precomputed. The methods discussed here can be extended to three dimensions in space and velocity domain, but then memory consumption and computing time will become important issues. For example, the storage of the precomputed pairwise distances between grid points in 3 dimensions takes a considerable amount of memory.


JAC acknowledges support from the MINECO Spanish projects MTM2011-27739-C04-02 (FEDER), 2009-SGR-345 from Agència de Gestió d’Ajuts Universitaris i de Recerca-Generalitat de Catalunya, the Royal Society through a Wolfson Research Merit Award, and the Engineering and Physical Sciences Research Council (UK) grant number EP/K008404/1. The work of the last two authors has been supported by the German research foundation, DFG grant KL 1105/17-1 and KL 1105/18-1.


  1. Department of Mathematics, Imperial College London, South Kensington Campus, London SW7 2AZ, UK (
  2. Technische Universität Kaiserslautern, Department of Mathematics, Erwin-Schrödinger-Straße, 67663 Kaiserslautern, Germany ({klar,roth}
  3. Fraunhofer ITWM, Fraunhoferplatz 1, 67663 Kaiserslautern, Germany
  4. Technische Universität Kaiserslautern, Department of Mathematics, Erwin-Schrödinger-Straße, 67663 Kaiserslautern, Germany ({klar,roth}


  1. Albi, G., Balagué, D., Carrillo, J. A., von Brecht, J.: Stability Analysis of Flock and Mill Rings for Second Order Models in Swarming. SIAM J. Appl. Math., 74 (2014), pp. 794–818.
  2. Ballerini, M., Cabibbo, N., Candelier, R., Cavagna, A., Cisbani, E., Giardina, L., Lecomte, L., Orlandi, A., Parisi, G., Procaccini, A., Viale, M., Zdravkovic, V.: Interaction ruling animal collective behavior depends on topological rather than metric distance: evidence from a field study. Proc Natl Acad Sci USA, 105 (2008) , pp. 1232-1237.
  3. Barbaro, A., Taylor, K., Trethewey, P.F., Youseff, L., Birnir, B.: Discrete and continuous models of the dynamics of pelagic fish: application to the capelin. Mathematics and Computers in Simulation, 79 (2009), pp. 3397–3414.
  4. Barbaro, A., Einarsson, B., Birnir, B., Sigurthsson, S., Valdimarsson, H., Palsson, O.K., Sveinbjornsson, S., Sigurthsson, T.: Modelling and simulations of the migration of pelagic fish. ICES J. Mar. Sci., 66 (2009), pp. 826–838.
  5. Bertozzi, A. L., von Brecht, J., Sun, H., Kolokolnikov, T., Uminsky, D.: Ring Patterns and their Bifurcations in a Nonlocal Model of Biological Swarms, to appear in Comm. Math. Sci., (2014).
  6. Bolley, F., Canizo, J.A., Carrillo, J.A.: Stochastic Mean-Field Limit: Non-Lipschitz Forces & Swarming. Math. Mod. Meth. Appl. Sci., 21 (2011), pp. 2179–2210.
  7. Bonabeau, E., Dorigo, M., Theraulaz, G.: Swarm Intelligence: From Natural to Artificial Systems. Intelligence: From Natural to Artificial Systems (Oxford University Press, New York, 1999);
  8. Bonilla, L., Götz, T., Klar, A., Marheineke, N., Wegener, R.: Hydrodynamic limit for the Fokker-Planck equation of fiber lay–down models. SIAM Appl. Math., 68 (2007), pp. 648-655.
  9. Braun, W., Hepp, K.: The Vlasov Dynamics and Its Fluctuations in the 1/N Limit of Interacting Classical Particles. Commun. Math. Phys., 56 (1977), pp. 101-113.
  10. Camazine, S., Deneubourg, J.-L., Franks, N.R., Sneyd, J., Theraulaz, G., Bonabeau, E.: Self-Organization in Biological Systems. Princeton University Press (2003).
  11. Cañizo, J.A., Carrillo, J.A., Rosado, J.: A well-posedness theory in measures for some kinetic models of collective motion. Math. Mod. Meth. Appl. Sci., 21 (2011), pp. 515-539.
  12. Cañizo, J.A., Carrillo, J.A., Rosado, J.: Collective Behavior of Animals: Swarming and Complex Patterns. Arbor, 186 (2010), pp. 1035–1049.
  13. Carrillo, J.A., D’Orsogna, M.R., Panferov, V.: Double milling in self-propelled swarms from kinetic theory. Kinetic and Related Models, 2 (2009), pp. 363-378.
  14. Carrillo, J. A., Huang, Y., Martin, S.: Nonlinear stability of flock solutions in second-order swarming models. Nonlinear Anal. Real World Appl., 17 (2014), pp. 332–343.
  15. Carrillo, J. A., Huang, Y., Martin, S.: Explicit flock solutions for quasi-morse potentials. to appear in European J. Appl. Math., (2014).
  16. Carrillo, J.A., Fornasier, M., Rosado, J., Toscani, G.: Asymptotic Flocking Dynamics for the kinetic Cucker-Smale model. SIAM J. Math. Anal., 42 (2010), pp. 218–236.
  17. Carrillo, J.A., Fornasier, M., Toscani, G., Vecil, F.: Particle, Kinetic, and Hydrodynamic Models of Swarming, in Naldi, G., Pareschi, L., Toscani, G. (eds.) Mathematical Modeling of Collective Behavior in Socio-Economic and Life Sciences, Series: Modelling and Simulation in Science and Technology, Birkhauser, (2010), pp. 297–336.
  18. Carrillo, J.A., Klar, A., Martin, S., Tiwari, S.: Self-propelled interacting particle systems with roosting force. Math. Mod. Meth. Appl. Sci., 20 (2010), pp. 1533–1552.
  19. Carrillo, J.A., Martin, S., Panferov, V.: A new interaction potential for swarming models. Physica D, 260 (2013), pp. 112–126.
  20. Carrillo, J. A., Vecil, F.: Nonoscillatory interpolation methods applied to Vlasov-based models. SIAM J. Sci. Comput., 29 (2007), pp. 1179–1206.
  21. Cheng, C., Knorr, G.: The integration of the Vlasov equation in configuration space, J.Comp.Phys., 22 (1976), pp. 330–351.
  22. Chuang, Y.L., D’Orsogna, M.R., Marthaler, D., Bertozzi, A.L., Chayes, L.: State transitions and the continuum limit for a 2D interacting, self-propelled particle system. Physica D, 232 (2007), pp. 33-47.
  23. Couzin, I.D., Krause, J., Franks, N.R., Levin, S.A.: Effective leadership and decision making in animal groups on the move. Nature, 433 (2005), pp. 513-516.
  24. Couzin, I.D., Krause, J., James, R., Ruxton, G. and Franks, N.: Collective memory and spatial sorting in animal groups. Journal of Theoretical Biology, 218 (2002), pp. 1-11.
  25. Degond, P., Frouvelle, A., Liu, J.-G.: Macroscopic limits and phase transition in a system of self-propelled particles. J. Nonlinear Sci., 23 (2013), pp. 427–456.
  26. Degond, P., Motsch, S.: Continuum limit of self-driven particles with orientation interaction. Math. Models Methods Appl. Sci., 18 (2008), pp. 1193-1215.
  27. Dobrushin, R.: Vlasov equations. Funct. Anal. Appl., 13 (1979), pp. 115-123.
  28. D’Orsogna, M.R., Chuang, Y.L., Bertozzi, A.L., Chayes, L.: Self-propelled particles with soft-core interactions: patterns, stability, and collapse. Phys. Rev. Lett., 96 (2006), 104302
  29. Douglas, J., Huang, C., Pereira, F.: The Modified Method of Characteristics with Adjusted Advection. Numer. Math., 83 (1999), pp. 353–369.
  30. Filbet, F., Sonnendrücker, E., Bertrand, P.: Conservative numerical schemes for the Vlasov equation. J. Comput. Phys., 172 (2001), pp. 166–187.
  31. Götz, T., Klar, A., Marheineke, N., Wegener, R.: A stochastic model and associated Fokker–Planck equation for the fiber lay-down process in nonwoven production processes. SIAM Appl. Math, 67 (2007), pp. 1704-1717
  32. Golse, F.: The mean field limit for the dynamics of large particle systems. Journées équations aux dérivées partielles, 9 (2003), pp. 1–47.
  33. Grégoire, G., Chaté, H.: Onset of collective and cohesive motion. Phy. Rev. Lett., 92 (2004), 025702
  34. Ha, S.-Y., Liu, J.-G.: A simple proof of the Cucker-Smale flocking dynamics and mean-field limit. Comm. Math. Sci., 7, (2009), pp. 297-325
  35. Sonnendrücker, E., Roche, J., Bertrand, P., Ghizzo, A.: The Semi-Lagrangian Method for the Numerical Resolution of the Vlasov Equation. J. Comp. Phys. 149 (1999) pp. 201–220.
  36. Ha, S.-Y., Tadmor, E.: From particle to kinetic and hydrodynamic descriptions of flocking. Kinetic and Related Models, 1 (2008), pp. 415-435.
  37. Hemelrijk, C.K. and Kunz, H.: Density distribution and size sorting in fish schools: an individual-based model. Behavioral Ecology, 16 (2005), pp. 178-187.
  38. Herty, M., Klar, A., Motsch, S., Olawsky, F.: A smooth model for fiber lay-down processes and its diffusion approximations. Kinetic and Related Models 2 (2009), pp. 480-502.
  39. Hildenbrandt, H, Carere, C., and Hemelrijk, C. K.: Self-organised complex aerial displays of thousands of starlings: a model. Behav. Ecol., 21 (2010), pp. 1349–1359.
  40. Huth, A. and Wissel, C.: The Simulation of the Movement of Fish Schools. Journal of Theoretical Biology, 152 (1992), pp. 365-385
  41. Klar, A., Reuterswärd, P., Seaïd, M.: A semi-Lagrangian method for the Fokker–Planck equation of fiber dynamics. J. Sci. Comput., 38 (2009), pp. 349-367
  42. Klar, A., Maringer, J., Wegener, R.: A 3D model for fiber lay-down in nonwoven production processes. Math. Mod. Meth. in Appl. Sci, 22(9) (2012).
  43. Klar, A., Reuterswärd, P., Seaïd, M.: A semi-Lagrangian method for a Fokker-Planck equation describing fiber dynamics. J.Sci.Comp., 38 (2009), pp. 349–367.
  44. Kolokolnikov, T., Carrillo, J. A., Bertozzi, A. L., Fetecau, R., Lewis, M.: Emergent behaviour in multi-particle systems with non-local interactions. Phys. D, 260 (2013), pp. 1–4.
  45. Kunz, H. and Hemelrijk, C. K.: Artificial fish schools: collective effects of school size, body size, and body form. Artificial Life, 3 (2003), pp. 237-253.
  46. Leveque, R.J.: Finite Volume Methods for Hyperbolic Problems. Cambridge University Press (2002).
  47. Levine, H., Rappel, W.J., Cohen, I.: Self-organization in systems of self-propelled particles. Phys. Rev. E, 63 (2000), 017101.
  48. Li, Y.X., Lukeman, R., Edelstein-Keshet, L.: Minimal mechanisms for school formation in self-propelled particles. Physica D, 237 (2008), pp. 699-720.
  49. Li, Y.X., Lukeman, R., Edelstein-Keshet, L.: A conceptual model for milling formations in biological aggregates. Bull Math Biol., 71 (2008), pp. 352-382.
  50. Mogilner, A., Edelstein-Keshet, L., Bent, L., Spiros, A.: Mutual interactions, potentials, and individual distance in a social aggregation. J. Math. Biol., 47 (2003), pp. 353-389.
  51. Motsch, S., Tadmor, E.: A new model for self-organized dynamics and its flocking behavior. J. Stat. Phys., 144 (2011), pp. 923–947.
  52. Neunzert, H.: The Vlasov equation as a limit of Hamiltonian classical mechanical systems of interacting particles. Trans. Fluid Dynamics, 18 (1977), pp. 663-678.
  53. Parrish, J., Edelstein-Keshet, L.: Complexity, pattern, and evolutionary trade-offs in animal aggregation. Science, 294 (1999), pp. 99-101.
  54. Quarteroni, A., Valli, A.: Numerical Approximation of Partial Differential Equations. Springer, Second Edition (1997).
  55. Qiu, J.-M., Christlieb, A.: A Conservative high order semi-Lagrangian WENO method for the Vlasov Equation. Journal of Computational Physics, 229 (2010), pp. 1130–1149.
  56. Roth, A., Zharovsky, E., Klar, A., Simeon,B.: A Semi-Lagrangian Method for 3-D Fokker Planck Equations for Stochastic Dynamical Systems on the Sphere. to appear in J. Sci. Comput., (2014).
  57. Spohn, H.: Large scale dynamics of interacting particles. Texts and Monographs in Physics, Springer (1991).
  58. Spohn, H.: Kinetic equations from Hamiltonian dynamics: Markovian limits. Rev. Modern Phys., 52 (1980), pp. 569–615.
  59. Toro, E.: Riemann Solvers and Numerical Methods for Fluid Dynamics. Springer, Third Edition (2009).
  60. Van Leer, B.: Towards the ultimate conservative difference scheme II. Monotonicity and conservation combined in a second order scheme. J. Comp. Phys., 14 (1974), pp. 361–370.
  61. Vecil, F., Laffitte, P., Rosado, J.: A numerical study of attraction/repulsion collective behavior models: 3D particle analyses and 1D kinetic simulations. Physica D., 260 (2013), pp. 127–144.
  62. Vicsek, T., Czirók, A., Ben-Jacob, E., Cohen, I., Shochet, O.: Novel type of phase transition in a system of self-driven particles. Phys. Rev. Lett., 75 (1995), pp. 1226–1229.
  63. Zharovsky, E., Simeon, B.: A space-time adaptive approach to orientation dynamics in particle laden flow. Procedia Computer Science, 1 (2010), pp. 791–799.
  64. Zharovski, E., Moosaie, A., LeDuc, A., Manhart, M., Simeon, B.: On the numerical solution of a convection-diffusion equation for particle orientation dynamics on geodesic grid. IMACS Journal Applied Numerical Mathematics, 62 (2012), pp. 1554–1566.
Comments 0
Request Comment
You are adding the first comment!
How to quickly get a good reply:
  • Give credit where it’s due by listing out the positive aspects of a paper before getting into which changes should be made.
  • Be specific in your critique, and provide supporting evidence with appropriate references to substantiate general statements.
  • Your comment should inspire ideas to flow and help the author improves the paper.

The better we are at sharing our knowledge with each other, the faster we move forward.
The feedback must be of minumum 40 characters
Add comment
Loading ...
This is a comment super asjknd jkasnjk adsnkj
The feedback must be of minumum 40 characters
The feedback must be of minumum 40 characters

You are asking your first question!
How to quickly get a good answer:
  • Keep your question short and to the point
  • Check for grammar or spelling errors.
  • Phrase it like a question
Test description