Monte Carlo gPC methods for diffusive kinetic flocking models with uncertainties

TL;DR

This paper develops Monte Carlo gPC numerical schemes for kinetic flocking models with uncertainties, effectively capturing phase transitions and noise effects while preserving positivity and achieving high accuracy.

Contribution

It introduces a novel combination of Monte Carlo and stochastic Galerkin methods for uncertain kinetic flocking models, ensuring positivity and high accuracy.

Findings

Validated methods on kinetic alignment models with noise.

Demonstrated accurate capture of phase transition phenomena.

Showed influence of uncertainties on smoothing and confidence bands.

Abstract

In this paper we introduce and discuss numerical schemes for the approximation of kinetic equations for flocking behavior with phase transitions that incorporate uncertain quantities. This class of schemes here considered make use of a Monte Carlo approach in the phase space coupled with a stochastic Galerkin expansion in the random space. The proposed methods naturally preserve the positivity of the statistical moments of the solution and are capable to achieve high accuracy in the random space. Several tests on a kinetic alignment model with self propulsion validate the proposed methods both in the homogeneous and inhomogeneous setting, shading light on the influence of uncertainties in phase transition phenomena driven by noise such as their smoothing and confidence band.

Figures

Equations

$$\begin{array}{rl}& {\partial }_{t}f+v\cdot {\nabla }_{x}f={\nabla }_v\cdot \left[\mathcal{B}[f]f+{\nabla }_v(Df)\right],\end{array}$$

$$\mathcal{B}[f](\theta ,x,w,t)=S(\theta ,v)+{\int }_{{\mathbb{R}}^{d_x}}{\int }_{{\mathbb{R}}^{d_v}}P(\theta ,x,x_{\ast })(v-v_{\ast })f(\theta ,x_{\ast },v_{\ast },t)d{v}_{\ast }d{x}_{\ast },$$

$$P(\theta ,x,x_{\ast })=\delta (x-x_{\ast }),S(\theta ,v)=\alpha (\theta )(|v{|}^2-1)v,D=D(\theta ),$$

$${\partial }_{t}f+v\cdot {\nabla }_{x}f={\nabla }_v\cdot \left[\alpha (\theta )(|v{|}^2-1)vf+{\rho }_f(v-u_f)f+D(\theta ){\nabla }_{v}f\right],$$

$${\rho }_f(\theta ,x,t)={\int }_{{\mathbb{R}}^{d_v}}f(\theta ,x,v,t)dv,{\rho }_f(\theta ,x,t)u_f(\theta ,x,t)={\int }_{{\mathbb{R}}^{d_v}}vf(\theta ,x,v,t)dv.$$

$$\xi (\theta ,v,t)=\Phi (\theta ,v)+(U\ast f)(\theta ,v,t)+D(\theta )\log f(\theta ,v,t),$$

$$\Phi (\theta ,v)=\alpha (\theta )\left(\frac{|v{|}^4}{4}-\frac{|v{|}^2}{2}\right),$$

$${\partial }_{t}f(\theta ,v,t)={\nabla }_v\cdot \left(f(\theta ,v,t){\nabla }_v\xi (\theta ,v,t)\right).$$

$$\begin{array}{rl}\mathcal{E}(\theta ,t)=& {\int }_{{\mathbb{R}}^{d_v}}\left(\alpha (\theta )\frac{|v{|}^4}{4}+(1-\alpha (\theta ))\frac{|v{|}^2}{2}\right)f(\theta ,v,t)dv-\frac{1}{2}|u_f{|}^2\\ & +D(\theta ){\int }_{{\mathbb{R}}^{d_v}}f(\theta ,v,t)\log f(\theta ,v,t)dv,\end{array}$$

$$u_f(\theta ,t)={\int }_{{\mathbb{R}}^{d_v}}vf(\theta ,v,t)dv.$$

$$f^{\infty }(\theta ,v)=C\exp \left\{-\frac{1}{D(\theta )}\left[\alpha (\theta )\frac{|v{|}^4}{4}+(1-\alpha (\theta ))\frac{|v{|}^2}{2}-u_{f^{\infty }}\cdot v\right]\right\},$$

$${\partial }_{t}f={\nabla }_v\cdot \left[\alpha (\theta )(|v{|}^2-1)vf+(v-u_f(\theta ,t))f+D(\theta ){\nabla }_{v}f\right],$$

$$P(\theta ,x,x_{\ast })=\frac{H(\theta )}{{\left(1+|x-x_{\ast }{|}^2\right)}^{\gamma }},H\ge 0,$$

$$\theta :(\Omega ,F)\to (I_{\theta },{\mathcal{B}}_{I_{\theta }}),$$

$$\mathbb{E}\left[{\Phi }_h(\theta ){\Phi }_k(\theta )\right]={\int }_{I_{\theta }}{\Phi }_h(\theta ){\Phi }_k(\theta )\Psi (\theta )d\theta ={\delta }_{hk},$$

$$f(\theta ,x,v,t)\approx f^M(\theta ,x,v,t)=\sum _{h=0}^M{\hat{f}}_h(x,v,t){\Phi }_h(\theta ),$$

$${\partial }_t{f}^M+v\cdot {\nabla }_v{f}^M={\nabla }_v\cdot \left[\mathcal{B}[f^M]f^M+{\nabla }_v\left(D{f}^M\right)\right].$$

$$\begin{array}{rl}& {\partial }_t{\hat{f}}_h(x,v,t)+v\cdot {\nabla }_x{\hat{f}}_h(x,v,t)=\\ & {\nabla }_v\cdot \left[\sum _{k=0}^M{S}_{hk}{\hat{f}}_k(x,v,t)+\sum _{k=0}^M{P}_{hk}[f^M]{\hat{f}}_k(x,v,t)+\sum _{k=0}^M{D}_{hk}{\nabla }_v{\hat{f}}_k(x,v,t)\right],\end{array}$$

$$S_{hk}(v)={\int }_{I_{\theta }}S(\theta ,v){\Phi }_h(\theta ){\Phi }_k(\theta )\Psi (\theta )d\theta ,$$

$$P_{hk}[\hat{f}]={\int }_{I_{\theta }}{\int }_{{\mathbb{R}}^{d_x\times d_v}}P(\theta ,x,x_{\ast })(v-v_{\ast })f^M(\theta ,x_{\ast },v_{\ast },t){\Phi }_h(\theta ){\Phi }_k(\theta )\Psi (\theta )d{v}_{\ast }d{x}_{\ast }d\theta$$

$$D_{hk}={\int }_{I_{\theta }}D(\theta ){\Phi }_h(\theta ){\Phi }_k(\theta )\Psi (\theta )d\theta .$$

$$\mathbb{E}[f(\theta ,x,v,t)]\approx {\hat{f}}_0(x,v,t),Var[f(\theta ,x,v,t)]\approx \sum _{h=0}^M{\hat{f}}_h(x,v,t)\mathbb{E}[{\Phi }_h^2(\theta )]-{\hat{f}}_0^2(x,v,t).$$

$${\partial }_t\hat{\text{f}}+v\cdot {\nabla }_x\hat{\text{f}}={\nabla }_v\cdot \left[(\text{S}+\text{P}[\hat{\text{f}}])\hat{\text{f}}+\text{D}{\nabla }_v\hat{\text{f}}\right].$$

$$\Vert \hat{\text{f}}{\Vert }_{L^2}:={\left[{\int }_{{\mathbb{R}}^{d_x}\times {\mathbb{R}}^{d_v}}\left(\sum _{h=0}^M{\hat{f}}_h^2(x,v,t)\right)dvdx\right]}^{1/2}.$$

$$\begin{array}{r}\Vert f^M{\Vert }_{L^2(\Omega )}={\int }_{I_{\theta }}{\int }_{{\mathbb{R}}^{d_x}\times {\mathbb{R}}^{d_v}}{\left(\sum _{h=0}^M{\hat{f}}_h{\Phi }_h(\theta )\right)}^{2}dvdx\Psi (\theta )d\theta .\end{array}$$

$$\Vert f^M{\Vert }_{L^2(\Omega )}=\Vert \hat{\text{f}}{\Vert }_{L^2}.$$

$$\Vert \hat{\text{f}}{\Vert }_{L^2}^2\le e^{t(C_{\mathcal{B}}+2C_D)}\Vert \hat{\text{f}}(0){\Vert }_{L^2}^2.$$

$$\begin{array}{rl}& {\int }_{{\mathbb{R}}^2}\left[{\partial }_t\left(\frac{1}{2}{\hat{f}}_h^2\right)+v{\partial }_x\left(\frac{1}{2}{\hat{f}}_h^2\right)\right]dvdx\\ & ={\int }_{{\mathbb{R}}^2}{\partial }_v\left[\sum _{k=0}^M(S_{hk}+P_{hk}){\hat{f}}_k+{\partial }_v\sum _{k=0}^M{D}_{hk}{\hat{f}}_k\right]{\hat{f}}_{h}dvdx.\end{array}$$

$$\begin{array}{rl}& \sum _{k=0}^M{\int }_{{\mathbb{R}}^2}{\hat{f}}_h{\partial }_v\left((S_{hk}+P_{hk}){\hat{f}}_k\right)dvdx\\ & =\sum _{k=0}^M{\int }_{{\mathbb{R}}^2}\left({\hat{f}}_k{\hat{f}}_h{\partial }_v(S_{hk}+P_{hk})+(S_{hk}+P_{hk}){\hat{f}}_h{\partial }_v{\hat{f}}_k\right)dvdx\\ & =-\sum _{k=0}^M{\int }_{{\mathbb{R}}^2}(S_{hk}+P_{hk}){\partial }_v\left({\hat{f}}_h{\hat{f}}_k\right)dvdx-\sum _{k=0}^M{\int }_{{\mathbb{R}}^2}{\hat{f}}_k{\partial }_v\left((S_{hk}+P_{hk}){\hat{f}}_h\right)dvdx,\end{array}$$

$$\begin{array}{rl}2\sum _{h,k=0}^M{\int }_{{\mathbb{R}}^2}{\hat{f}}_h{\partial }_v\left((S_{hk}+P_{hk}){\hat{f}}_k\right)dvdx& =-\sum _{h,k=0}^M{\int }_{{\mathbb{R}}^2}(S_{hk}+P_{hk}){\partial }_v({\hat{f}}_k{\hat{f}}_h)dvdx\\ & =\sum _{h,k=0}^M{\int }_{{\mathbb{R}}^2}{\hat{f}}_h{\hat{f}}_k{\partial }_v(S_{hk}+P_{hk})dvdx.\end{array}$$

$$\sum _{h,k=0}^M{\int }_{{\mathbb{R}}^2}{\hat{f}}_h{\partial }_v\left((S_{hk}+P_{hk}){\hat{f}}_k\right)dvdx\le \frac{C_{\mathcal{B}}}{2}\sum _{h,k=0}^M{\int }_{{\mathbb{R}}^2}|{\hat{f}}_h{\hat{f}}_k|dvdx,$$

Full text

**Keywords: uncertainty quantification, stochastic Galerkin, collective behaviour.

**

MSC: 35Q83, 65C05, 65M70.

1 Introduction

Uncertainty quantification (UQ) for partial differential equations describing real world phenomena gained in recent years lot of momentum in various communities. Without intending to review the very huge literature on this topic we mention [3, 18, 22, 34, 40, 44, 45, 46, 48, 49] and the references therein. One of the main advantages of this approach relies in its capability to provide a sound mathematical framework to reproduce realistic experiments through the introduction of the stochastic quantities, reflecting our incomplete information on some features on the systems’ modelling. This argument is especially true for the description of emergent social structures in interacting agents’ systems in socio-economic and life sciences. Common examples are the emergence of consensus phenomena in opinion dynamics, flocking and milling patterns in swarming of animals or humans and the formation of stable wealth distributions in economic systems, see [39]. It is worth observing how for these models we can have at most statistical information on initial conditions and on the modeling parameters, which are in practice substituted by empirical social forces embedding a huge variability, see for example [5].

From a mathematical viewpoint, the kinetic equations we are interested in are nonlinear Vlasov-Fokker-Planck equations depending on random inputs taking into account uncertainties in the interaction terms or in the boundary conditions. These equations arise in the description of collective phenomena modeling the evolution of a distribution function $f=f(\theta,x,v,t)$, $t\geq 0$, $x\in\mathbb{R}^{d_{x}}$, $v\in\mathbb{R}^{d_{v}}$, $d_{x},d_{v}\geq 1$, and $\theta\in I_{\theta}\subseteq\mathbb{R}$ a random field, according to

[TABLE]

where $\mathcal{B}[\cdot]$ is a non–local operator of the form

[TABLE]

being $P\geq 0$ for all $v,v_{*}\in\mathbb{R}^{d_{v}}$ and $D=D(\theta)\geq 0$ is an uncertain constant diffusion that depends on the introduced uncertainty.

As a follow-up question to the progress of the analytical understanding of real phenomena, the existence of phase transitions driven by noise represents a deeply fascinating issue, see [17, 19, 23, 25, 28, 29] and the references therein. The notion of phase transition has been fruitfully borrowed from thermodynamics to highlight the phase change of the system under specific stimuli. In particular, it is of interest the emergence of patterns in the collective dynamics for critical strengths of noise as in the classical Kuramoto model [8, 14, 37, 42] or in collective behavior [6, 28]. In the present paper we concentrate on a flocking model for interacting agents with self-propulsion and diffusion where the parameters are assumed to be stochastic. This model has been recently investigated in [1, 6, 7, 10] in absence of uncertainties. The introduction of uncertain quantities points in the direction of a more realistic description of the underlying processes and helps us to compute possible deviations from the prescribed deterministic behavior.

Suitable numerical methods that preserve the positivity of the distribution function are developed and are based on the so-called of Monte Carlo generalized polynomial chaos (MCgPC) methods. The introduced class of schemes is based on a Monte Carlo approach in the phase space that is coupled with stochastic Galerkin decomposition in the random space [34, 35, 39, 46]. This method has been recently proposed in [18] in the case of zero diffusion and it will be extended for models incorporating noise in the present setting. Furthermore, beside the natural positivity preservation, the resulting methods are spectrally accurate in the random space for sufficiently regular uncertainties. Furthermore, the adoption of a spectrally accurate methods in the random space is particularly efficient if compared with other non-intrusive approaches.

The rest of this paper has the following structure: in Section 2 we introduce in detail the class of kinetic flocking models of interest, the issue of phase transition will be treated with a particular focus on the relation between self-propulsion strength and noise intensity in the case they both depend on uncertain quantities. Section 3 focuses on the construction of Monte Carlo generalized polynomial chaos methods, a reduction in the computational complexity is here achieved through a Monte Carlo mean-field algorithm discussed in previous works. Finally, Section 4 is devoted to numerical test for the validation of the proposed schemes. Here continuous and Monte Carlo schemes will be compared to show the effectiveness of the MCgPC approach.

2 Phase transitions in kinetic flocking models with uncertainties

Kinetic models for aggregation-diffusion dynamics encountered extensive investigation in recent years [6, 7, 11, 13, 15, 16, 24, 26]. This class of models describes the aggregate behaviour of large systems of self-propelled particles for which stable patterns emerge asymptotically. Nevertheless, at the present times, the study of the influence of ineradicable uncertainties in the modeling parameters is quite an unexplored area. We mention in this direction [3, 18, 30, 31, 44].

Let us consider first a kinetic flocking model for self-propelled particles with both state-independent interactions and uncertain diffusion given by (LABEL:eq:MF_general) with the following choice

[TABLE]

being $\delta(x-x_{*})$ the Dirac delta distribution centered in $x_{*}\in\mathbb{R}^{d_{x}}$. The model of interest is now characterized by a localised Vlasov-Fokker-Planck equation with random inputs for the dynamics of $f=f(\theta,x,v,t)$ and having the following form

[TABLE]

where $\alpha(\theta)\geq 0$, $D(\theta)\geq 0$ are respectively self–propulsion strength and intensity of the diffusion operator, and where $u_{f}$ is the momentum of the system which is not conserved in time due to the presence of the self-propulsion term

[TABLE]

It seems worth stressing that, in the present setting, the role of the additive diffusion operator $\Delta_{v}f$ in (LABEL:eq:MF_general), representing the impact of unpredictable events, is strongly different from that of the uncertainties introduced in the definition of the model parameters. The velocity diffusion term is, in fact, representative of all the possible modifications of the dynamics not modeled in a structural way, whereas $\theta\in I_{\theta}$ summarizes all source of modifications in the prescribed dynamics. We point the reader to [44] for a more focussed discussion in a related setting.

In the following we will treat separately the case space-homogeneous case and then we will provide some insights on the non-localised inhomogeneous setting.

2.1 Space homogeneous problem

In [6] the authors investigated the space homogeneous version of (3) in the deterministic setting. It has been proven that a phase transition between the so-called polarized and unpolarized motion takes place as the noise intensity $D$ increases and for a specific range of the values of the self-propulsion strength $\alpha$. At a difference with the cited results, the present formulation of the model embeds from the very beginning the presence of uncertain quantities.

Let us consider for simplicity the case of unitary mass, so that $f$ is a probability density for all times $t\geq 0$. First, we observe that (3) can be rewritten as a gradient flow. In fact, if we define

[TABLE]

with $U(v)$ a interaction potential, given by $U(v)=\tfrac{|v|^{2}}{2}$, and $\Phi(v)$ a confining potential of the form

[TABLE]

the equation reads

[TABLE]

A free energy functional which dissipates along solution is defined for all $\theta\in I_{\theta}$ as follows

[TABLE]

with

[TABLE]

Stationary solutions of the space–homogeneous problem satisfy the identity $\nabla_{v}\xi(\theta,t)=0$ for all $\theta\in I_{\theta}$ and have the form

[TABLE]

being $C>0$ a normalization factor. It is possible to prove the following result, see [6, Theorem 2.1].

Theorem 1**.**

The space–homogeneous nonlinear Fokker–Planck equation

[TABLE]

exhibits a phase transition in the following sense:

• •

For small enough diffusion $D(\theta)\rightarrow 0$ for all $\theta\in I_{\theta}$ there is a function $u=u(D(\theta))$ with $\lim_{D\rightarrow 0}u(D(\theta))=1$, such that $f^{\infty}(\theta,v)$ with $u=(u(D(\theta)),0,\dots,0)$ is a stationary solution of the original problem.

• •

For large enough diffusion coefficient $D$ the only stationary solution is the symmetric distribution given by (4) with $u_{f}\equiv 0$.

This behavior is reminiscent of the one observed for the Vicsek model [6, 23], where agents move with a constant speed and interact with their neighbors via a local alignment force and are subject to noise action. In particular, a critical noise intensity has been discovered whose value determines a phase transition between the ordered states and a chaotic state characterized by a null asymptotic average velocity of the system of agents. The critical noise value is theoretically demonstrated in [23] for the Vicsek model while for the present model there is a strong numerical evidence for its existence [6].

2.2 Space inhomogeneous case

In the non-localized inhomogeneous setting sharp results on analogous phase transitions are sill non present in the literature due to the additional difficulties of the models. In the zero diffusion limit and without self-propulsion forces, small variations of the arguments in [15] prove that for a Cucker-Smale type interaction

[TABLE]

unconditional alignment, i.e. convergence to a profile travelling with the same mean velocity, emerges in the case $\gamma\leq 1/2$ without conditions on the initial configuration of the system or on its dimensionality both at the kinetic and particle level [15, 20]. In this case the asymptotic distribution is a Dirac delta in the velocity space, meaning that for large times the agents share the same velocity and that they form a group with frozen mutual distances. At the particle level, the case $\gamma>1/2$ leads to a conditional flocking, i.e. alignment emerges for provable initial configurations of the system [20]. Initial density conditions to achieve asymptotic alignment for $\gamma>1/2$ are still not proven, in general we may expect flocking for high mean-density distributions [16].

In the following we will present a numerical insight for the case of a VFP model with Cucker-Smale type interaction forces (5) in dimension $d_{x}=1$ and $d_{v}=1$, self-propulsion forces and uncertain diffusion coefficient. In particular, numerical results obtained with accurate numerical schemes highlight that a phase transition can occur also in this regime, shedding light on the deep interplay between alignment forces and noise strength in more general settings with not localized alignment [20, 38].

3 Stochastic Galerkin methods for kinetic equations

We introduce Stochastic Galerkin (SG) numerical methods with applications to the nonlinear Vlasov-Fokker-Planck (VFP) equation (LABEL:eq:MF_general). We discuss the class of stochastic Galerkin (SG) methods and, in particular, we concentrate on the generalized Polynomial Chaos (gPC) decomposition [22, 46, 47, 49]. These methods gained increased popularity in recent years since they guarantees spectral accuracy in the random space under suitable regularity conditions.

In our schemes, we exploit Monte Carlo methods for the approximation of the numerical solution of the (VFP) in the phase space taking advantage of particle based reformulation of the problem that converges in distribution to the solution for an increasing number of agents. The core idea then is to apply SG-gPC techniques for the efficient approximation of the random field in the resulting MC approximation. We highlight that the combination of the two approaches leads to a positive approximation of statistical moments of the solution even in the nonlinear case.

In the following we will derive a SG-gPC scheme for the continuous problem and we will consider specific formulation for the so-called Monte Carlo gPC methods (MCgPC), see [18].

3.1 Preliminaries on SG-gPC expansion

In this section we derive a SG-gPC approximation for the uncertain nonlinear VFP equation (LABEL:eq:MF_general) with nonlocal drift $\mathcal{B}[\cdot]$ having the structure (2). Let $(\Omega,F,P)$ be a probability space where as usual $\Omega$ is the sample space, $F$ is a $\sigma-$algebra and $P$ a probability measure, and let us define the random variable

[TABLE]

with $I_{\theta}\subseteq\mathbb{R}$ and $\mathcal{B}_{I_{\theta}}$ the Borel set, and with known probability density function $\Psi(\theta):I_{\theta}\rightarrow\mathbb{R}^{+}$. Hence, we consider the linear space $\mathbb{P}_{M}$ of polynomials of degree up to $M$, which is generated by a family of orthonormal polynomials $\{\Phi_{h}(\theta)\}_{h=0}^{M}$ such that

[TABLE]

being $\delta_{hk}$ the Kronecker delta function. Assuming that $\Psi(\theta)$ has finite second order moment, we can approximate the distribution $f\in L^{2}(\Omega,\mathcal{F},P)$ in terms of the following chaos expansion

[TABLE]

being $\hat{f}_{h}(x,v,t)=\mathbb{E}\left[f(\theta,x,v,t)\Phi_{h}(\theta)\right]$ for all $h=0,\dots,M$. The SG-gPC formulation for the approximated distribution $f^{M}=f^{M}(\theta,x,v,t)$ solution of (LABEL:eq:MF_general) is then

[TABLE]

Thanks to the orthogonality of the polynomial basis of $\mathbb{P}_{M}$ we obtain a coupled system of $M+1$ purely deterministic PDEs for the evolution of each projection $\hat{f}_{h}(x,v,t)$, $h=0,\dots,M$ that reads

[TABLE]

having defined the following matrices for $h,k=0,\dots,M$

[TABLE]

The statistical quantities of interest may be defined in terms of the coefficients of the expansion. In particular we have

[TABLE]

We introduce the vector $\hat{\textbf{f}}=\left(\hat{f}_{0},\dots,\hat{f}_{M}\right)$ and the $(M+1)\times(M+1)$ matrices $\textbf{S}=\{S_{hk}\}_{h,k=0}^{M}$, $\textbf{P}[\hat{\textbf{f}}]=\{P_{hk}\}_{h,k=0}^{M}$, $\textbf{D}=\{D_{hk}\}_{h,k=0}^{M}$, whose components are given by (7)-(9), we can reformulate (6) in more compact form as follows

[TABLE]

In the following we indicate with $\|\hat{\textbf{f}}\|_{L^{2}}$ the standard $L^{2}$ norm of the vector $\hat{\textbf{f}}(x,v,t)$

[TABLE]

We can easily observe that since the norm of $f^{M}(\theta,x,v,t)$ in $L^{2}(\Omega)$ is

[TABLE]

thanks to the orthonormality of the basis $\{\Phi_{k}\}_{k=0}^{M}$ in $L^{2}(\Omega)$ it follows that

[TABLE]

We can prove the following stability result. For notation simplicity we consider the case $d_{x}=d_{v}=1$.

Theorem 2**.**

Assume that there exists $C_{\mathcal{B}},C_{D}>0$ such that $\|\partial_{v}(S_{hk}+P_{hk})\|_{L^{\infty}}\leq C_{\mathcal{B}}$ and $D_{hk}\leq C_{D}$ for all $h,k=0,\dots,M$, then

[TABLE]

Proof.

We multiply (6) by $\hat{f}_{h}$ and integrate over $\mathbb{R}^{2}$ to obtain

[TABLE]

After integration by parts the transport term vanishes and the right hand term reads

[TABLE]

and thanks to the symmetry of S and P we have

[TABLE]

Since $\|\partial_{v}(S_{hk}+P_{hk})\|_{L^{\infty}}\leq C_{\mathcal{B}}$ we have

[TABLE]

and from Cauchy-Schwarz we have

[TABLE]

Furthermore, since $D_{hk}\leq C_{D}$ we have

[TABLE]

and we obtain

[TABLE]

and thanks to Gronwall’s Lemma, we conclude. ∎

We highlight how, as for classical spectral methods, the introduced SG decomposition leads to the loss of important physical properties of the solution like the positivity. A positivity preserving scheme is instead obtained if the matrix $\textbf{S}+\textbf{P}=\{S_{hk}+P_{hk}\}_{h,k=0}^{M}$ and the matrix $\textbf{D}=\{D_{hk}\}_{h,k=0}^{M}$ are diagonal. This case corresponds to the situation where no uncertainties are present in self-propelling and diffusion terms and in case of interaction with a background distribution that does not encode any uncertainty. Under the introduced assumptions and with uncertain initial distribution of agents $f(\theta,x,v,0)=f_{0}(\theta,x,v)$ the SG-gPC decomposition of the original problem reads for all $h=0,\dots,M$

[TABLE]

Therefore, we have to solve a set of $M+1$ decoupled VFP equations. In this case we can apply structure preserving type schemes to ensure the positivity of the statistical moments, accurate description of the large time behavior of each projection and entropy dissipation, see [4, 13, 40, 41, 50].

To tackle efficiently the general fully nonlinear case in the following we will introduce a novel scheme that exploits the spectral convergence in the random space of SG-gCP methods and that naturally preserves the positivity of the numerical distribution.

3.2 Particle based SG-gPC formulation of kinetic equations

As for classical spectral methods, the solution of the coupled SG system (6) for $f^{M}$ looses its positivity and then a clear physical meaning. In order to overcome the difficulty recently the so-called Monte Carlo gPC (MCgPC) scheme has been proposed. These methods combine the advantages of a Monte Carlo approach for the approximation of $f$ in the phase space and they conserve spectral accuracy in the random space. We refer to [18] for the study a wide range of mean-field equations describing the emergence of patterns and ordered behavior in large interacting systems in the zero diffusion limit.

The Monte Carlo (MC) method is a probabilistic particle method that describe the evolution of density functions by resorting to the computation of interactions in a finite set of randomly chosen particles of which it is known a priori the dynamics, see [21, 39]. In MC methods the updated distribution function in the phase space is typically reconstructed from particles as post-processing. Several approaches are possible for the reconstruction step, which can estimated through a histogram, weighted integration rule methods or by a convolution of the empirical particle distribution with a suitable mollifier, see [33]. All the mentioned approximations preserve the positivity of the obtained numerical distribution function. Concerning the order of accuracy for MC methods, we have a convergence rate of the order $\mathcal{O}(\frac{1}{\sqrt{N}})$ where $N$ is the number of considered samples, see [12].

Therefore, in order to set-up a MC scheme for the general equation (LABEL:eq:MF_general) we need to find a system of stochastic differential equations which converge in distribution to the correct solution of the problem. This problem has been tackled in a more general setting in [9] and equation (LABEL:eq:MF_general) can be derived from the following second order system of SDEs for $(x_{i}(\theta,t),v_{i}(\theta,t))\in\mathbb{R}^{d_{x}}\times\mathbb{R}^{d_{v}}$, $i=1,\dots,N$, with random inputs

[TABLE]

with initial positions and velocities $(x_{i}(\theta,0),v_{i}(\theta,0))$. In (10) we denoted by $\{W_{i}\}_{i=1}^{N}$ a vector of $N$ independent Wiener processes and independent with the random varlable $\theta$. Notice that the random variable $\theta$ plays the role of a parameter, thus for each sampling from the distribution $\Psi(\theta)$, we can apply the mean-field results developed in [9] to conclude that (3) is the mean field limit of (10).

It is worth to remark that in the case of vanishing self-propulsion, i.e. $\alpha(\theta)\equiv 0$, and for symmetric interactions, i.e. $P(\theta,x_{i},x_{j})=P(\theta,x_{j},x_{i})$, the above system of SDEs preserves the initial mean velocity of the system in the limit $N\rightarrow+\infty$. Indeed we have

[TABLE]

to conclude due to a law of large numbers argument. In general, the mean velocity is not conserved due to the presence of the self-propulsion forces and appropriate methods must be developed to catch the correct transient regimes.

In the following we summarize the assumptions to ensure the mean-field convergence. Let us introduce the empirical measure associated to the particle system (10) weighted by the distribution of the uncertainty $\Psi(\theta)$, i.e.

[TABLE]

The transition to chaos for the particle system (10) follows from a small variation of [9, Theorem 1.1] under non restrictive assumptions on the a kernel $K(x_{i}-x_{j},v_{i}-v_{j})=-P(\theta,|x_{i}-x_{j}|)(v_{i}-v_{j})$ and the drift $S(\theta,v_{i})$. The proof is based on the following argument: the $N$ interacting process $(x_{i}(\theta,t),v_{i}(\theta,t))_{i=1}^{N}$ behaves in the limit $N\rightarrow+\infty$ like the process defined by the set $(\bar{x}_{i}(\theta,t),\bar{v}_{i}(\theta,t))_{i=1}^{N}$ solution of the McKean-Vlasov equations

[TABLE]

where $f$ is the probability law of $(\theta,\bar{x}_{i}(\theta,t),\bar{v}_{i}(\theta,t))$and $\{W_{i}\}_{i=1}^{N}$ are the Wiener processes characterizing (10). The processes are identically distributed and independent from the random variable $\theta$, and by the Itô formula their law $f$ evolves according a VFP equation of type (LABEL:eq:MF_general). We summarize this in the following result representing a direct extension of Theorem 1.1 in [9] that holds for all $\theta\in I_{\theta}$.

Theorem 3** ([9, Theorem 1.1]).**

Let $f_{0}(\theta,x,v)$ be a Borel probability measure and let $(x_{i}(\theta,0),v_{i}(\theta,0))$ for all $i=1,\dots,N$ be $N$ independent variables for all $\theta\in I_{\theta}$. We assume that the drift and the antisymmetric kernel satisfy that there exist constants $A,L,p>0$ such that

[TABLE]

for all $x,y\in\mathbb{R}^{d_{x}}$, $v,w\in\mathbb{R}^{d_{v}}$. If the particle system (10) and the processes (11) have global solutions on the finite time interval $[0,T]$ with initial data $({x}_{i}(\theta,0),{v}_{i}(\theta,0))$ such that

[TABLE]

with $f(\theta,\cdot,\cdot)=\textrm{law}(\bar{x}_{i},\bar{v}_{i})$ for $\theta\in I_{\theta}$, then there exists a constant $C>0$ such that

[TABLE]

If there exists $p^{\prime}>p$ such that

[TABLE]

then for all $0<\epsilon<1$ there exists a constant $C>0$ such that

[TABLE]

for all $0\leq t\leq T$ and $N\geq 1$.

We have denoted with $\mathbb{E}_{\theta,W}[\cdot]$ the expectation taken with respect to the distribution of the Wiener processes and the random variable $\theta$. We highlight how the provided bounds hold for all $\theta\in I_{\theta}$ and in particular we take expectation with respect to the introduced stochastic parameter. It is easily verified that for the introduced self-propulsion force and the kernels of interest the assumptions are met. Therefore the mentioned theorem gives the following quantitative result on the convergence of the empirical measure $f^{(N)}$ to the distribution $f$: let $\varphi$ a Lipschitz map on $\mathbb{R}^{d_{x}}\times\mathbb{R}^{d_{v}}$ then

[TABLE]

where $\varepsilon(N)\rightarrow 0$ for $N\rightarrow+\infty$ with a rate $\varepsilon(N)$ defined in (12)-(13) depending on the assumptions on the solutions in Theorem 3 and the law of large numbers to deal with the last term. These arguments are standard in interacting particle systems and analogous to [9] and the references therein.

In the next section we derive a SG-gPC decomposition of (10) so that we will preserve the exponential convergence in the random space with respect to all the uncertain quantities.

3.2.1 Stochastic Galerkin scheme for the particle system

In the following we consider the gPC decomposition of the microscopic dynamics (10). Space and velocity variables of the $i$th agent, for all $i=1,\dots,N$, are approximated by $(x_{i}^{M},v_{i}^{M})$ where

[TABLE]

being as before $\{\Phi_{k}(\theta)\}_{k=0}^{M}$ an orthonormal basis of $L^{2}(\Omega)$.

Next we give explicit representation of the SG-gPC expansion at the particle level. To obtain the SG-gPC decomposition for the particle system we rewrite (10) in terms of $x_{i}^{M}$, $v_{i}^{M}$ for all $i=1,\dots,N$, whose projection in the linear spaces of degree $h=0,\dots,M$ reads

[TABLE]

and again thanks to the orthogonality of the polynomial basis we obtain for all $i=1,\dots,N$ the following coupled system of SDEs describing the evolution of each component of the original variables in the phase space

[TABLE]

where

[TABLE]

and

[TABLE]

It is worth to remark that in the case of vanishing self-propulsion $s_{hk}(v_{i}^{M})=0$ for all $h,k=0,\dots,M$ and for symmetric interactions $p_{hk}^{ij}=p_{hk}^{ji}$ we recover the conservation of the mean velocity of the system also in the SG decomposition.

The convergence of the SG-gPC expansion (6) for sufficiently regular function follows from standard results in polynomial approximation theory, we recall for example [27]. In general, thanks to the property of the introduced polynomial basis we have

Theorem 4**.**

We consider $\theta\in I_{\theta}$ with distribution $\Psi(\theta)$ and the basis $\{\Phi_{h}\}_{h=0}^{M}$ of orthonormal polynomial basis in $L^{2}(\Omega)$ . Then, for any $g(\theta,t)\in L^{2}(\Omega)$

[TABLE]

where $g^{M}=\sum_{k=0}^{M}\left(\int_{I_{\theta}}g(\theta,t)\Phi_{k}(\theta)d\Psi(\theta)\right)\Phi_{k}(\theta)$.

In the present setting, thanks to Theorem 4 it is reasonable to expect the convergence of the empirical measure

[TABLE]

to the empirical measure $f^{(N)}$ as $M\rightarrow\infty$ for all $t\geq 0$ and for $(x_{i}^{M},v_{i}^{M})$ solution to (15). Hence, from the mean-field convergence showed in Theorem 3 we would guarantee that $f^{(N,M)}$ converges to $f(\theta,x,v,t)$ solution of the initial VFP problem by taking first the limit of the Galerkin modes $M\to\infty$ and then the limit of large number of particles $N\to\infty$. We can prove the following result following the footprint of (14) and Theorem 4 applied to $x_{i}(\theta,t)$ and $v_{i}(\theta,t)$.

Theorem 5**.**

Let us define the following empirical measure

[TABLE]

Hence, provided that $x_{i}^{M}(\theta,t)\in L^{2}(\Omega)$, $v_{i}^{M}(\theta,t)\in L^{2}(\Omega)$ we have

[TABLE]

in $P(P(\Omega\times\mathbb{R}^{d_{x}}\times\mathbb{R}^{d_{v}}))$.

Proof.

To prove the convergence in $P(P(\Omega\times\mathbb{R}^{d_{x}}\times\mathbb{R}^{d_{v}}))$ we consider a sufficiently regular test function in all variables $(\theta,x,v)$, that we can assume of the form $\varphi_{1}(\theta)\varphi_{2}(x,v)$ without loss of generality, and we compute

[TABLE]

providing the convergence result for $M\rightarrow+\infty$. ∎

In the last section we will give numerical evidence of this result.

3.2.2 Monte Carlo gPC scheme

We now approximate the limiting stochastic kinetic equation taking advantage of the particle reformulation of the problem. In fact, since the solution of the system of SDEs (10) converges in distribution to the solution of the original problem (LABEL:eq:MF_general) for $N\rightarrow+\infty$, we can approximate the original dynamics by means of a Monte Carlo (MC) method in the phase space. The main drawback of this approach lies in the computational cost $\mathcal{O}(M^{2}N^{2})$, since at each time step and for each gPC projection each agent modifies its velocity in a genuine nonlinear way.

A significant reduction in terms of computational cost can be achieved through a mean field MC evaluation of the interaction dynamics as originally proposed in [2], see [18] for the UQ framework. Thanks to this approach we have an efficient algorithm for transport and interaction in the phase space and we can reconstruct the expected solution from the particle system from positions and velocities at the microscopic level, which is considered in the SG-gPC setting as in Section 3.2.1. This approach has been recently analysed in connection to other problems in [36].

Algorithm 1** (MCgPC for nonlocal nonlinear VFP equations).**

Consider $N$ samples $(x_{i},v_{i})$ with $i=1,\dots,N$ from the initial $f_{0}(x,v)$, and fix $S\leq N$ a positive integer;

• 2.

Perform gPC expansion up to order $M\geq 0$ over the set of $S\leq N$ particles to obtain the projections $(\hat{x}_{i,h},\hat{v}_{i,h})$, $h=0,\dots,M$.

• for $n=0$ to $T-1$

Generate $N$ Brownian paths $\{\eta_{i}\}_{i=1}^{N}=\{W_{i}^{n+1}-W_{i}^{n}\}_{i=1}^{N}\sim\mathcal{N}(0,1)$, $\{\eta_{i}(0)\}_{i=1}^{N}=0$

• for $i=1$ to $N$
• a)

sample $S$ particles $j_{1},\dots,j_{S}$ uniformly without repetition among all particles;

• b)

compute the position and velocity change

[TABLE]

• end for

Reconstruction $\mathbb{E}[f(\theta,x,v,n\Delta t)]=\int_{I_{\theta}}f(\theta,x,v,n\Delta t)\Psi(\theta)d\theta$ .

• end for

We sketch in Figure 1 the two approaches for the approximation of statistical quantities of the nonlinear nonlocal stochastic VFP problems of interest. In the left branch we find the standard SG approach where first we consider the gPC approximation of the original problem, which generates a coupled system of PDEs that can be solved through deterministic methods to obtain the evolution of expected quantities. On the right branch we describe the introduced MCgPC procedure, based on a particle reformulation of the problem which converges in distribution in the phase space to the solution of the problem. The advantage of considering a gPC scheme for the microscopic system lies in the preservation of the typical spectral convergence in the random space of the method. We highlight how thanks to the adoption of the computational strategy in Algorithm 1 the overall cost becomes $\mathcal{O}(M^{2}SN)$, $S\ll N$.

For the reconstruction of expected quantities, in the present manuscript we consider the histogram of position and velocity of the set of particles in the phase space, we point the reader to [33] for possible alternatives. Thanks to the MC approach the resulting method preserves the positivity of the expected distribution function.

Remark 1**.**

In the case $S=N$ we obtain the typical convergence rate of order $\mathcal{O}(1/\sqrt{N})$ of Monte Carlo in the phase space, where $N$ is the number of particles and spectral convergence in $M$ in the random space. Due to the presence of the fast evaluation of interactions, the case $S<N$ induces an additional error of order $\mathcal{O}(\sqrt{1/S-1/N})$ in the microscopic dynamics. Let $f^{(S,M)}(\theta,x,v,t)$ the empirical distribution density estimated at time $t>0$ from Algorithm 1 and let $J=\{j_{1},\dots,j_{S}\}$, $|J|=S$, be the vector of indexes sampled uniformly without repetition in $\{1,\dots,N\}$. Then for any test function $\varphi$ we have

[TABLE]

whose leading error is given by $I_{1}$ since $S\ll N$ and, for $N\gg 0$, we may observe that the accuracy with which we approximate $f^{(N,M)}$ with $f^{(S,M)}$ is $\mathcal{O}\left(\sqrt{\dfrac{1}{S}-\dfrac{1}{N}}\right)$ thanks to a central limit theorem-type argument.

4 Numerical tests

In this section we present several numerical examples based on (3) both in the homogeneous and inhomogeneous case. We test the effectiveness of the MC-gPC scheme through several tests based on VFP equations. In all test the integration of the system of stochastic differential equations (10) is performed through a standard Euler-Maruyama method whereas the solution of the system of PDEs derived from the SG procedure is solved through a standard central scheme coupled with a fourth order Runge-Kutta integration. In the whole section we will consider a uniform noise, therefore Gauss-Legendre polynomial basis are chosen in the gPC setting. Numerical investigations on the influence of uncertainties in phase transition phenomena are presented through the section. Finally, we explore the non localized inhomogeneous model with Cucker-Smale type interactions.

4.1 Test 1: Space homogeneous case

In this first test we consider the space homogeneous problem in 1D with uncertain diffusion parameter of the form

[TABLE]

where $\theta\sim\mathcal{U}([-1,1])$ and $\lambda\leq 1$ is such that $D(\theta)\geq 0$ for all $\bar{D}\geq 0$. In this first test we consider a constant self propelling strength $\alpha(\theta)=\alpha\geq 0$.

The evolution of the density function $f(\theta,v,t)$, $v\in\mathbb{R}$ is ruled by the following PDE

[TABLE]

whose SG-gPC approximation is given for all $h=0,\dots,M$ by

[TABLE]

being

[TABLE]

At the particle level we obtain from (15) the following coupled system of SDEs for the evolution of the particles’ velocities

[TABLE]

being

[TABLE]

and $\left\{W_{i}\right\}_{i=1}^{N}$ defines a set of $N$ independent Wiener processes. Furthermore, we indicated with $\hat{u}_{h}$ the $h$th projection of the average velocity of the system.

The long time solution is given by (4) as discussed in Section 2.1. We consider as deterministic initial condition the Gaussian distribution

[TABLE]

At the particle level the initial velocities are a sample of $N\gg 0$ velocities from the distribution $f_{0}(v)$. In Figure 2 we compare the numerical expected distribution of the SG-gPC system (17) and the reconstructed expected density of the MCgPC scheme at time $T=50$ for two uncertain diffusion coefficients (16) with $\bar{D}=0.2$, $\bar{D}=0.8$ and $\lambda=0.1$ in both cases. It is easily observed how for an increasing number of particles the statistical quantities are well approximated by the scheme.

In Figure 3 we study the convergence of the temperature of the system

[TABLE]

in $L^{2}(\Omega)$ obtained through the MCgPC algorithm for an increasing order of the gPC expansion. We considered $N=10^{4}$ particles and two uncertain diffusion parameters of the form (16) with $\bar{D}=0.2$ and $\bar{D}=0.8$, $\lambda=0.1$ and $\theta\sim\mathcal{U}([-1,1])$. Time integration has been performed up to time $T=50$ with time step $\Delta t=10^{-2}$.

To complete the overview of the one dimensional setting we computed in Figure 4 the large time behaviour of the expected average velocity $\mathbb{E}[u_{f}(\theta,T)]$ for several $\bar{D}=\frac{\kappa}{10}$, $\kappa=0,\dots,10$ with constant self-propulsion term $\alpha=1$. We considered an uncertain diffusion of the form (16) with $\lambda=10^{-1}$ (top row) and $\lambda=10^{-3}$ (bottom row). In each figure we marked in dashed red the confidence band computed in terms of the approximated standard deviation $\sqrt{\textrm{Var}(u_{f})}$.

We can easily observe the regions of maximal sensitivity with respect to the presence of uncertainties. In particular, for high diffusion values the variability of the expected average velocity vanishes and we may argue that the phase transition predicted in [6] is actually a quite stable pattern in the space homogeneous regime. Nevertheless, the averaging of uncertain quantities acts as a smoothing factor of the phase transition as we can clearly observe in Figure 4-(b) and 4-(d). For a vanishing influence of $\theta\in I_{\theta}$ given by $\lambda\rightarrow 0$ in (16) the transition becomes sharper coherently with the deterministic case, see figure 4-(d). Summarizing the main effect of the uncertainties is the smoothing of the transition point making it less sharper and abrupt than in the deterministic cases.

4.2 Test 2. Space inhomogeneous case

In this section we focus on inhomogeneous models. First we consider the localised case for which a phase transition is expected to happen similarly to the homogeneous case, see [6], even if this has not yet been proved. To explore other possible alternatives we will also consider the case of Cucker-Smale type interactions for which no analogous theoretical results regarding phase transitions currently exist. In all the tests of the present section we consider a sample of $N=10^{5}$ particles for the MCgPC scheme.

4.2.1 Test 2A. Localized interaction case

We consider a space dependent interaction function of Dirac delta form, i.e. $P(x,x_{*})=\delta(x-x_{*})$. With this choice agents interact only is they share the same position. In this case, the MCgPC scheme translates in a second order system of SDEs of the form

[TABLE]

where $s_{hk}$ has the same definition of the space homogeneous case. As before we consider an uncertain diffusion parameter $D(\theta)$ and a constant self-propulsion $\alpha=1$.

We introduce a discretization of the phase space $\{x_{i}\}_{i=1}^{N_{x}}$, $\{v_{i}\}_{i=1}^{N_{v}}$ such that $x_{i+1}-x_{i}=\Delta x$, $v_{i+1}-v_{i}=\Delta v$ made by $N_{x}=20$ gridpoints for the space variable and by $N_{v}=40$ points for the velocity variable. At time $t=0$ we consider the following multivariate Gaussian distribution

[TABLE]

we considered $\mu_{x}=0$, $\mu_{v}=1$, $\sigma_{v}=\frac{1}{10}$, $\sigma_{x}=\frac{1}{100}$. Let us consider the evolution of the gPC system of PDEs derived in (6) with localized interactions in (8). At the numerical level we perform a second order Strang splitting approach. The transport step is integrated through a third order Runge-Kutta scheme with a fifth order WENO reconstruction [43]. In the following we compare the behavior of the continuous approximation with the one obtained through the MCgPC scheme where, to localize the interaction, we considered as a smoothing factor the indicator function of the numerical cell $[x_{i},x_{i+1}]$.

In Figure 5 and Figure 6 we compare the evolution of the expected density over the time interval $[0,5]$ for the inhomogeneous kinetic equation (LABEL:eq:MF_general) with local interactions. In particular, we present the expected density obtained through numerical integration of the SG-gPC system of PDEs with the one approximated through the MCgPC scheme with fast evaluation of the interaction term through a mean-field MC approach described in Algorithm 1 and interaction subset of $S=10$ particles. The uncertain diffusion has the form (16) with $\bar{D}=0.2$, $\bar{D}=0.8$ and $\lambda=0.1$. Periodic boundary conditions has been considered. The results are in agreement with the space homogeneous case since the interaction are localised, it is easily observed how a phase transition occurs switching from an ordered state to a chaotic isotropic state for increasing values of the diffusion. In Figure 7 the average in position distributions in velocity are compared between the MCgPC and gPC. These results corroborate the appearance of a smoothed phase transition: a shifted distribution with positive average velocity for small values of the average diffusion is observed while for larger diffusion the distribution becomes almost symmetric.

4.2.2 Test 2B. Cucker-Smale type interactions

In this section we consider the evolution of statistical quantities obtained from (LABEL:eq:MF_general) with space dependent interaction kernel of the Cucker-Smale (CS) type

[TABLE]

This choice of interaction has been introduced in [20] and received a lot of attention in the recent literature on kinetic models for emergent behavior. Without intending to review the huge related literature we mention [15, 16, 32] for an introduction on the topic. Recent efforts on CS dynamics in the UQ setting are [3, 34]. As we already mentioned, at the present times there are no analytical results regarding possible phase transitions for this highly nonlinear kernel in the VFP setting. The determination of sharp estimates on the possible phase transitions is an open problem which goes beyond the purpose of the present manuscript. Hence, to complete the overview on the introduced methods we give a numerical insight on the features of the resulting dynamics.

In order to observe the action of uncertainty in the system we consider a stochastic diffusion parameter $D(\theta)=\bar{D}+\lambda\bar{D}\theta$, with $\lambda=10^{-1}$ and $\theta\sim\mathcal{U}([-1,1])$ in the cases $\bar{D}=0.2$ and $\bar{D}=0.8$. In Figure 8-9 we report the expected dynamics of the VFP model (LABEL:eq:MF_general) with CS interactions and stochastic diffusion obtained through the MCgPC scheme. Periodic boundary conditions have been taken into account. In particular, a sample of $N=10^{5}$ particles is considered for the reconstruction of the expected density in the interval $[-2,2]\times[-3,3]$ with $N_{x}=20$, $N_{v}=40$ gridpoints. In the random space we considered $M=3$ modes for the gPC expansion. We can observe how the expected density moves from an ordered state for small values of diffusion to a chaotic isotropic for sufficiently big values of diffusion. We postpone to later works a more detailed insight on these phenomena for general interaction kernels.

Conclusion

The introduction of uncertainties in swarming dynamics is of paramount importance for the large-scale description of realistic phenomena. We concentrated in this manuscript on the case of Vlasov-Fokker-Planck equations for emerging collective behaviour with phase transitions depending on the strength of self-propelling and noise terms. We extended to these models the construction of nonnegative gPC schemes for kinetic-type equations. The proposed approach takes advantage of a Monte Carlo approximation of the kinetic equation in phase space which is coupled with a stochastic Galerkin gPC expansion. Several numerical tests were proposed both in the homogeneous and inhomogeneous settings, proving the effectiveness of the approach and shedding light on the action of uncertainties in terms of the stability of the phase transitions, We generically observe that they lead to a smoothing of the sharp transition values. Several extensions of the present work are possible from both the analytical and numerical viewpoints and in connection with more general kernels in the collision/interaction dynamics.

Acknowledgments

JAC was partially supported by the EPSRC grant number EP/P031587/1. MZ is member of the National Group on Mathematical Physics (GNFM) of INdAM (Istituto Italiano di Alta Matematica), Italy. This work has been written within the activities of the Excellence Project CUP: E11G18000350001 of the Department of Mathematical Sciences ”G. L. Lagrange” of Politecnico di Torino funded by MIUR (Italian Ministry for Education, University and Research).