3 Numerical Studies

Models of particle systems contain complex dynamics and static plots often do not best illustrate their behaviour. For this reason, animations of the systems have also been produced to accompany this report. They are available at
https://tom271.github.io/InteractingParticleSystems/

This is a precursor to a webpage that will accompany a future work in which we hope to produce interactive simulations of the system through a Dash app1. Where animations are available, they will be clearly signposted in figure captions. For more details on the implementation, see Appendix D. We begin this section by confirming known results numerically, before starting to answer some of the questions raised in the previous section.

To demonstrate some possible behaviours of the locally scaled particle system (2.1)-(2.2), we simulate it using a smooth herding function(\(G_1\)) with \(\varphi\equiv1\) and \(\sigma=1\) for 1000s. The aim is to corroborate the analytic results for this system. Figure 3.1 shows a histogram of the particles’ positions at \(t=0\) and \(t=1000\). We see that they appear to converge to either \(\mu_+\) or \(\mu_-\), the known stationary measures of the (2.5), depending on the sign of the average velocity of the initial data. This is in agreement with the results of the previous section (and Appendix A). However, if the simulation is ran for longer, we see that these are only metastable states and instead the average velocity of the particles switches between \(+1\) and \(-1\). A similar effect can be produced by increasing the noise or reducing the number of particles. Assuming uniform interaction, we have proved that the stationary measure of the particle system has mean zero. This is numerical corroboration of that fact. This switching of stability is a well-known phenomenon and indeed is behind the strength of many optimisation algorithms such as SGLD [10]. Figure better illustrates this switching by reducing the number of simulated particles. By removing the noise present in (2.1)-(2.2), we can remove this effect. If we assume space homogeneity, set \(\varphi \equiv 1\) and remove the noise the system (\(\sigma=0\)) reduces to \[\begin{align} \mathrm{d} {v}_t^{i,N} & = - v_t^{i,N} \mathrm{d} t + G\left(M(t)\right) \mathrm{d} t , \tag{3.1} \end{align}\] where \[ M(t) = \frac{1}{N}\sum_{j=1}^N v^{j,N}_t \] Summing over \(i\) and dividing by \(N\) both sides of (3.1) one obtains a closed equation for the mean velocity, that is \[\mathrm{d}M(t) = \left(G(M(t)) - M(t)\right)\mathrm{d}t. \] This is in accord with the result for the space homogeneous PDE (2.7) when \(\varphi\equiv1\) and has stationary points where \(M(t) = G(M(t))\). For our examples, this corresponds to \(M(t) =-1, 0 ,+1\).

Regardless of the initial data, the system always converges to $\mu_+$ or $\mu_-$ when $\varphi\equiv1$. Here, initial configurations with positive initial average velocity are coloured green and those with negative initial velocity are orange. Here we simulate $1000$ particles under a uniform interaction with a smooth herding function and $\sigma=1$. After $1000s$ a histogram is used as an approximation for the empirical measure and plotted. The initial data are randomly selected Gaussians with means varying between $-5$ and $5$ while the standard deviations lie between $0.5$ and $3$. Here a smooth herding function $G_1$ was used, and a timestep of $\Delta t =0.01$ in an Euler-Maruyama scheme (see Appendix \@ref(app-implementation))

Figure 3.1: Regardless of the initial data, the system always converges to \(\mu_+\) or \(\mu_-\) when \(\varphi\equiv1\). Here, initial configurations with positive initial average velocity are coloured green and those with negative initial velocity are orange. Here we simulate \(1000\) particles under a uniform interaction with a smooth herding function and \(\sigma=1\). After \(1000s\) a histogram is used as an approximation for the empirical measure and plotted. The initial data are randomly selected Gaussians with means varying between \(-5\) and \(5\) while the standard deviations lie between \(0.5\) and \(3\). Here a smooth herding function \(G_1\) was used, and a timestep of \(\Delta t =0.01\) in an Euler-Maruyama scheme (see Appendix D)

Evolving the dynamics for a longer time suggests an invariant measure with mean $0$ due to the possibility of jumping between states. Here $N=50 , \sigma=1$ and a smooth herding function is used. The grey dashed lines show Gaussian distributions centred at $\pm1$ for reference.

Figure 3.2: Evolving the dynamics for a longer time suggests an invariant measure with mean \(0\) due to the possibility of jumping between states. Here \(N=50 , \sigma=1\) and a smooth herding function is used. The grey dashed lines show Gaussian distributions centred at \(\pm1\) for reference.

Green hexagonal particles are those moving with positive velocity. Orange circles shows negative velocity. The left column of histograms shows the approximated particle distribution at time \(t\), while the right column shows the approximate distribution up to time \(t\), that is calculated on \([0,t]\). We can see the shifts in the stability from the majority of the particles being green (positive velocity) to orange (negative velocity). The orange lines on the histograms show the invariant measure predicted by the analysis of the space homogeneous PDE (2.7) when \(\varphi\equiv1\).

In the absence of noise, we can time the convergence to equilibrium for different initial data as we know that no switching is possible. Once a baseline has been set for uniform interaction, we can compare this with non-uniform interactions, \(\varphi\not\equiv1\). To do this we use the following initial setup. We start the particles in two diametrically opposite clusters of equal width: one containing one third of the total particles with positive velocity, and the other containing two thirds of the particles with negative velocity. This allows for an interesting setup—most particles begin movement with negative velocity yet the average velocity can be positive. We give each particle in the small cluster velocity \(1.8\), while in the larger cluster all particles have velocity \(-0.2\). The initial average velocity of the system is thus \(\frac{7}{15}\approx 0.467\). \[\begin{align*} x^{i,3N}_0 &\sim U\left[\pi-\frac{\pi}{10},\pi+\frac{\pi}{10}\right], &&v^{i,3N}_0 = -0.2 \qquad \text{ for } 1\leq i\leq 2N,\\ x^{i,3N}_0 &\sim U\left[-\frac{\pi}{10},\frac{\pi}{10}\right], &&v^{i,3N}_0 = 1.8\qquad \text{ for } 2N+1\leq i\leq 3N. \end{align*}\] With this initial distribution, the system is very far from the stationary distribution of the PDE. To control the influence of the interaction function, we use the hard cutoff indicator function described in Figure . By varying the value of \(\gamma\), the range of interaction can be controlled. Note that this interaction is not uniformly bounded above zero, and as such does not satisfy the assumptions (2.3). For the herding function we use the step function common in the modelling literature. \[\begin{align*} \varphi(x) &= \mathbb{I}_{[-\pi\gamma,\pi\gamma]}(x)\,,\gamma\in \bigg\lbrack-\frac{1}{2},\frac{1}{2}\bigg\rbrack,\\ G(u) &= \frac{u+\beta \mathrm{sgn}(u)}{1+\beta}\, , \beta>0 \end{align*}\] As a guide, doubling the radius of interaction \(\gamma\), gives the fraction of the torus which can be seen by each particle, and \(\gamma= \frac{1}{2}\) is equivalent to \(\varphi\equiv1\).

Figure shows how this affects the speed of convergence. Here the number of particles \(N\) has no discernible effect on the dynamics. Surprisingly, we see here that the average velocity can settle into a periodic orbit not centred at \(-1,0\) or \(1\), as we might na"ively expect. The system does not converge to \(\pm1\) either in law or in the sense of ergodic averages—for an explanation of the difference, see Appendix C. This cannot occur in the noisy particle system and indeed is unexpected. However to see how this can occur, consider a system containing just three particles. There are three possible behaviours From the two particles case, we know the particles will either not interact and come to rest or they can coalesce and move in unison. Both of these are possible in a three particle system however a third option appears. Two of the particles can coalesce and quickly converge to velocity \(\pm1\). When this occurs, decay to rest is no longer possible. Upon meeting the third particle, it can either be absorbed into the cluster or the interaction can be insufficient and it gets left behind. The cluster containing only two particles will then accelerate back to its preferred velocity while the third particle returns to rest—until the cluster comes back around. This is what produces these at first counterintuitive periodic patterns. The same is happening here for very weak interactions. See online for the three particle case and the one present in Figure 3.3.

In the deterministic system ($\sigma=0$), unexpected behaviours can occur. For $\gamma=0.05$, we see the average velocity settle into a periodic orbit, and for $\gamma=0.1$, a convergence to $-1$. In a), we see the time taken for the average velocity to converge. Although skewed by the periodic trajectory, it is possible to see the dependence on $\gamma$. In b), the colour refers to the size of $\gamma$. In all these simulations, we use a step herding function and the opposing clusters initial setup.In the deterministic system ($\sigma=0$), unexpected behaviours can occur. For $\gamma=0.05$, we see the average velocity settle into a periodic orbit, and for $\gamma=0.1$, a convergence to $-1$. In a), we see the time taken for the average velocity to converge. Although skewed by the periodic trajectory, it is possible to see the dependence on $\gamma$. In b), the colour refers to the size of $\gamma$. In all these simulations, we use a step herding function and the opposing clusters initial setup.

Figure 3.3: In the deterministic system (\(\sigma=0\)), unexpected behaviours can occur. For \(\gamma=0.05\), we see the average velocity settle into a periodic orbit, and for \(\gamma=0.1\), a convergence to \(-1\). In a), we see the time taken for the average velocity to converge. Although skewed by the periodic trajectory, it is possible to see the dependence on \(\gamma\). In b), the colour refers to the size of \(\gamma\). In all these simulations, we use a step herding function and the opposing clusters initial setup.

When there are only three particles in the system, the periodic behaviours seen in Figure 3.3 become much more intuitive. Two clusters move together, but the interaction is not sufficient to attract the third. The bottom plot shows the average velocity of all the particles (red). Faint lines show the velocities of each individual particle. Here we use a hard cutoff interaction with \(\gamma=0.05\) and a smooth herding function \(G_1\).

Also surprising here is that the initial average velocity of the ensemble, \(M(0)\), does not govern the convergence as it did for both the PDE and the particle system when \(\varphi \equiv 1\). We see that there is a radius \(\gamma = 0.05\) for which the system converges to velocity \(-1\)—at odds with its initial average velocity, \(M(0)=\frac{7}{15}\). This occurs as the two initial clusters do not interact quick enough. An intuitive reasoning for this is that, seen separately, each cluster is converging to a different value. The cluster at \(0\) has positive velocity and thus accelerates towards average velocity \(1\), while the cluster at \(\pi\) has negative velocity and accelerates towards \(-1\). This can cause the average velocity of the entire ensemble to become negative before sufficient interaction occurs. Thus by the time the interactions begin, the system as whole will converge to \(-1\). We see this happening whenever \(\gamma \leq 0.1\).

Now that the role of the interaction function is better understood, we take a similar approach in our study of the herding function \(G\). In particular, we are interested in how the presence of a discontinuity affects the dynamics. As such, we take the smoothed herding function \[ G_{\alpha}(u) = \frac{\mathrm{atan}(\alpha u)}{\mathrm{atan}(\alpha)}\, , \alpha > 0 . \] The parameter \(\alpha\) allows us to control the gradient of the herding function. As $, G(u) $ tends towards a step function. Figure 3.4 shows the effect of three different values: \(\alpha = 1, 5,10\). It can be seen that as the herding function gets steeper, periodic orbits become more prevalent. The system behaves unlike we expect from any of the analytic results.

As the steepness of the herding function increases, we see more and more erratic behaviour for small $\gamma$. Here the herding function is smooth, $G_{\alpha}$ with $\alpha=1,5,10$. The interaction function has a hard cutoff.As the steepness of the herding function increases, we see more and more erratic behaviour for small $\gamma$. Here the herding function is smooth, $G_{\alpha}$ with $\alpha=1,5,10$. The interaction function has a hard cutoff.As the steepness of the herding function increases, we see more and more erratic behaviour for small $\gamma$. Here the herding function is smooth, $G_{\alpha}$ with $\alpha=1,5,10$. The interaction function has a hard cutoff.

Figure 3.4: As the steepness of the herding function increases, we see more and more erratic behaviour for small \(\gamma\). Here the herding function is smooth, \(G_{\alpha}\) with \(\alpha=1,5,10\). The interaction function has a hard cutoff.

To aid in the stochastic case, it will be pertinent to examine the existence of spatially inhomogeneous distributions in the deterministic case. Starting with the two clusters of uneven density described above, we evolve the dynamics for \(50s\) seconds before counting the number of clusters that remain. We do this by looking at the histogram at the final time. The ``eyeball metric’’ was chosen after considering other methods such as \(k\)-means, Jenks optimisation or counting the minima of a kernel density estimate for its ease of adaptation to the torus. Table 3.1 below shows the number of clusters present after \(50s\). We see that as \(\gamma\) increases, the number of clusters present decreases. Note that despite the clusters, this system is still converging to a uniform distribution in the sense of ergodic averages.

Table 3.1:
N=24 N=168 N=360 N=456
\(\gamma=0.05\) 3 2 3 2
\(\gamma=0.15\) 2 2 2 1
\(\gamma=0.3\) 1 1 1 1
\(\gamma=0.45\) 1 1 1 1

In the deterministic setting, we have seen that a larger interaction radius speeds up convergence and if the radius is too small it is possible to converge to a different value than that expected. We also see that stable clusters form—as expected for a deterministic system. Surprisingly, we saw periodic average velocities emerge in many cases when the interaction radius was small. The herding function also affects convergence, with the presence of a steep gradient also causing periodic average velocities. Now, we turn our attention back to the stochastic system to see whether these phenomena persist in the presence of noise. We immediately expect to lose these periodic average velocities, as well as any clusters present.

To see the influence of noise, Figure 3.5 shows the effect on average velocity when \(\sigma=0.1\). Again we begin the evolution from two diametrically opposed clusters initially travelling with deterministic speeds \(-0.2\) and \(1.8\). It can be seen that no periodic orbits persist and all systems converge to \(\pm1\) except when \(\gamma=0\). Furthermore, smaller values of \(\gamma\) cause convergence to slow down, similar to the deterministic case. The difference here is clearer as there are no periodic orbits skewing the dataset.

One realisation of the locally scaled particle model with the step herding function; hard cutoff interaction function with $\gamma = 0, 0.05,0.1,\dots,0.5$ and $N=24, 72, 120, \dots, 408$. The initial distribution is two opposing clusters travelling in opposite directions. This is the same setup as Figure \@ref(fig:interaction) but with $\sigma=0.01$. We see that there no periodic orbits persist, unlike the deterministic case. The colour of the lines corresponds to the interaction radius $\gamma$.One realisation of the locally scaled particle model with the step herding function; hard cutoff interaction function with $\gamma = 0, 0.05,0.1,\dots,0.5$ and $N=24, 72, 120, \dots, 408$. The initial distribution is two opposing clusters travelling in opposite directions. This is the same setup as Figure \@ref(fig:interaction) but with $\sigma=0.01$. We see that there no periodic orbits persist, unlike the deterministic case. The colour of the lines corresponds to the interaction radius $\gamma$.

Figure 3.5: One realisation of the locally scaled particle model with the step herding function; hard cutoff interaction function with \(\gamma = 0, 0.05,0.1,\dots,0.5\) and \(N=24, 72, 120, \dots, 408\). The initial distribution is two opposing clusters travelling in opposite directions. This is the same setup as Figure 3.3 but with \(\sigma=0.01\). We see that there no periodic orbits persist, unlike the deterministic case. The colour of the lines corresponds to the interaction radius \(\gamma\).

To see if the volatile behaviour produced by increasing the gradient of \(G\) persists in the presence of noise, we now repeat the earlier experiment but here with \(\sigma=0.01\). Figure 3.6 shows the result: when there is noise in the system, its effect is much reduced. Note however that for \(\gamma = 0.1\), the system still settles into an unexpected state before the noise pushes the average velocity towards \(-1\).

As in Figure \@ref(fig:steepherding), here the herding function is smooth, $G_{lpha}$ with $lpha=1,5,10$ however here we introduce some noise, $\sigma=0.01$. The interaction function has a hard cutoff across a range of $\gamma$. Here we see that the noise removes the erratic periodic behaviours regardless of gradient. However note that there is an effect for some $\gamma$, particularly when $\alpha=10$.As in Figure \@ref(fig:steepherding), here the herding function is smooth, $G_{lpha}$ with $lpha=1,5,10$ however here we introduce some noise, $\sigma=0.01$. The interaction function has a hard cutoff across a range of $\gamma$. Here we see that the noise removes the erratic periodic behaviours regardless of gradient. However note that there is an effect for some $\gamma$, particularly when $\alpha=10$.As in Figure \@ref(fig:steepherding), here the herding function is smooth, $G_{lpha}$ with $lpha=1,5,10$ however here we introduce some noise, $\sigma=0.01$. The interaction function has a hard cutoff across a range of $\gamma$. Here we see that the noise removes the erratic periodic behaviours regardless of gradient. However note that there is an effect for some $\gamma$, particularly when $\alpha=10$.

Figure 3.6: As in Figure 3.4, here the herding function is smooth, \(G_{lpha}\) with \(lpha=1,5,10\) however here we introduce some noise, \(\sigma=0.01\). The interaction function has a hard cutoff across a range of \(\gamma\). Here we see that the noise removes the erratic periodic behaviours regardless of gradient. However note that there is an effect for some \(\gamma\), particularly when \(\alpha=10\).

When there is noise in the system, we can compare the rate at which clusters break to the deterministic system where it is possible for clusters to persist. To do this, we calculate the discrete \(\ell^1\) distance between the uniform distribution and a histogram of all particle positions at time \(t\). This will give a quantitative measure of how quickly clusters disperse in the stochastic regime. Figure shows the difference in \(\ell^1\) distance for a \(N=87\) and $= 0.02,0.04,0.09 $. It can be seen that regardless of \(\gamma\), the noise in the velocity variable quickly breaks any clusters present. Note that the \(\ell^1\) distance does not converge to 0 as expected as there is not enough particles in the system for the histogram to ever be a good approximation to the uniform distribution at any given time point. These results suggest that spatially inhomogeneous invariant distributions do not exist for the particle system even when \(\varphi\not \equiv 1\).

A comparison of $\ell^1$ distance between a histogram of particle positions and the uniform distribution. The deterministic system (orange) retains clusters and therefore a high distance, while even a small amount of noise (here, $\sigma=0.05$) disperses the particle positions.A comparison of $\ell^1$ distance between a histogram of particle positions and the uniform distribution. The deterministic system (orange) retains clusters and therefore a high distance, while even a small amount of noise (here, $\sigma=0.05$) disperses the particle positions.A comparison of $\ell^1$ distance between a histogram of particle positions and the uniform distribution. The deterministic system (orange) retains clusters and therefore a high distance, while even a small amount of noise (here, $\sigma=0.05$) disperses the particle positions.

Figure 3.7: A comparison of \(\ell^1\) distance between a histogram of particle positions and the uniform distribution. The deterministic system (orange) retains clusters and therefore a high distance, while even a small amount of noise (here, \(\sigma=0.05\)) disperses the particle positions.

3.1 On the Difference Between Local and Global Scaling

If the system is scaled globally, as in (2.12)-(2.13), we see a difference in speed of convergence. Figure 3.8 shows that for smaller values of \(\gamma\), the system has a much larger convergence time. The systems appear to settle into a steady state dependent on the value of \(\gamma\), before moving out of this state and towards \(-1\). Many trajectories settle into a periodic pattern very close to \(-1\)—these are caused by one or two particles remaining out of a large cluster as discussed in the deterministic case.As \(\gamma\) gets larger, the local scaling converges to the global scaling and we see less difference in the trajectories (cf. Figure 3.3. This is thought to be because for this interaction function, the interaction is always less than 1. Hence, the strength of the interaction in the local scaling will always be greater than that in the global scaling, due to the monotonicity of \(G\).

If the model is scaled globally, low values of $\gamma$ have a much larger convergence time. Here the system was evolved for $500s$ from an initial configuration of two clusters of uneven density. For higher interaction radii, behaviour is very similar. Some systems persist in a periodic trajectory for even larger times (tested up to $1000s$). Note the logarithmic scale.If the model is scaled globally, low values of $\gamma$ have a much larger convergence time. Here the system was evolved for $500s$ from an initial configuration of two clusters of uneven density. For higher interaction radii, behaviour is very similar. Some systems persist in a periodic trajectory for even larger times (tested up to $1000s$). Note the logarithmic scale.

Figure 3.8: If the model is scaled globally, low values of \(\gamma\) have a much larger convergence time. Here the system was evolved for \(500s\) from an initial configuration of two clusters of uneven density. For higher interaction radii, behaviour is very similar. Some systems persist in a periodic trajectory for even larger times (tested up to \(1000s\)). Note the logarithmic scale.

3.2 Simulation of Kinetic Models

Simulating the kinetic models is an ongoing work. Many standard numerical techniques seem unsuitable for this model. The PDE (2.5) is degenerate, in the sense that the noise is not elliptic. By only acting on the velocity, we must be very careful not to introduce any artificial diffusion into the space variable that is a common pitfall of many standard techniques such as upwind schemes (Figure 3.9 illustrates this for the one-way wave equation). Doing so would surely break any clusters that may form under the dynamics. As mass is also conserved in this system (See Appendix A), we also must be careful not to use any techniques that cause mass loss or allow negative solutions. To combat this, we have begun adapting 2DChebClass, a pseudospectral method originally developed for the solution of problems arising in dynamical density functional theory (DDFT) [11, 12]. These provide exponential accuracy in the number of grid points, meaning the mesh can be much coarser than in a typical finite volume scheme. The drawback of these methods is that they are limited in the geometry of the domain in which they can solve problems. As here we are only considering a two dimensional model on \(\mathbb{T}\times\mathbb{R}\), this is not an issue.

Even for smooth initial data, a first order upwind scheme exhibits severe artificial diffusion. Both plots show the one way wave equation $u_t-u_x=0$ with an indicator initial profile (a) and a Gaussian profile (b).Even for smooth initial data, a first order upwind scheme exhibits severe artificial diffusion. Both plots show the one way wave equation $u_t-u_x=0$ with an indicator initial profile (a) and a Gaussian profile (b).

Figure 3.9: Even for smooth initial data, a first order upwind scheme exhibits severe artificial diffusion. Both plots show the one way wave equation \(u_t-u_x=0\) with an indicator initial profile (a) and a Gaussian profile (b).

Simulating the kinetic model \@ref(eq:nl) using _2DChebClass_. Here we use a smooth herding function, uniform interaction ($\varphi=1$) and initial data $f_0 = \exp(-(v-1)^2/2)(\exp(\sin(x))+2)$. The equation is solved on $\mathbb{T}\times[-5,5]$, with 10 mesh points in both dimensions. The solution appears to contradict the non-negativity, however this is just an artefact of the interpolating method used for plotting---the solution calculated is always non-negative.

Figure 3.10: Simulating the kinetic model (2.5) using 2DChebClass. Here we use a smooth herding function, uniform interaction (\(\varphi=1\)) and initial data \(f_0 = \exp(-(v-1)^2/2)(\exp(\sin(x))+2)\). The equation is solved on \(\mathbb{T}\times[-5,5]\), with 10 mesh points in both dimensions. The solution appears to contradict the non-negativity, however this is just an artefact of the interpolating method used for plotting—the solution calculated is always non-negative.

The simulation shown above

References

[10] Nagapetyan, T., Duncan, A. B., Hasenclever, L., Vollmer, S. J., Szpruch, L. and Zygalakis, K. (2017). The true cost of stochastic gradient langevin dynamics.

[11] Goddard, B., Nold, A. and Kalliadasis, S. (2017). 2DChebClass.

[12] Nold, A., Goddard, B. D., Yatsyshin, P., Savva, N. and Kalliadasis, S. (2017). Pseudospectral methods for density functional theory in bounded and unbounded domains. Journal of Computational Physics 334 639–64.

  1. See https://dash.plotly.com/↩