Monte Carlo stochastic Galerkin methods for non-Maxwellian kinetic models of multiagent systems with uncertainties

In this paper, we focus on the construction of a hybrid scheme for the approximation of non-Maxwellian kinetic models with uncertainties. In the context of multiagent systems, the introduction of a kernel at the kinetic level is useful to avoid unphysical interactions. The methods here proposed, combine a direct simulation Monte Carlo (DSMC) in the phase space together with stochastic Galerkin (sG) methods in the random space. The developed schemes preserve the main physical properties of the solution together with accuracy in the random space. The consistency of the methods is tested with respect to surrogate Fokker–Planck models that can be obtained in the quasi-invariant regime of parameters. Several applications of the schemes to non-Maxwellian models of multiagent systems are reported.


Introduction
Kinetic equations are often studied to describe aggregate trends of large systems of interacting particles and have shown a remarkable effectivity in different research fields, ranging from classical rarefied gas dynamics to socio-economic and traffic flow dynamics.Without reviewing the huge literature on these topics, we mention [6,7,9,12,14,37,42,43] and the reference therein for an introduction to the subject.The contributions have to be further distinguished depending on the type of kernel characterizing the interaction frequency between particles or agents.It is worth mentioning that the introduction of a state-dependent kernel represents an essential tool in kinetic theory to enforce physical properties of rarefied gases [13], whereas it is currently underexplored in less classical applications to multi-agent systems.In this direction, we mention the following recent contributions [18,20,21,23].
The deterministic description of multi-agent phenomena has often to face the lack of essential information on microscopic dynamics, initial states, or boundary conditions.Hence, it is of paramount importance to quantify and control possible deviations from expected trends and patterns due to unavoidable uncertainties in the model parameters and initial distributions.An established idea relies on considering these quantities as random variables influencing the evolution of the kinetic distribution, increasing, therefore, the dimensionality of the problem.In recent years, we experienced a growing interest in the construction of numerical methods for kinetic equations with uncertainties, see the collection [31].Among the most popular techniques for uncertainty quantification, stochastic Galerkin (sG) methods are based on the construction of deterministic solvers and are capable to guarantee spectral convergence in the random field under suitable regularity assumptions [28,34,53].Anyway, their computational cost is generally high due to the curse of dimensionality of kinetic equations and they are highly intrusive with respect to the original formulation of the model.Furthermore, the main physical properties of the solution, like its positivity, entropy dissipation, and hyperbolicity, are lost.Besides sG methods, we find nonintrusive approaches to UQ that do not require significant modifications to the numerical scheme of the deterministic problem and are based on collocation strategies.Therefore, these latter methods are easy to parallelize and do not require any knowledge of the class of probability distributions of random parameters.In this direction, multi-fidelity approaches have been recently developed using control variate techniques, see [4,15,16,24,29,38].
In this work, we follow a different path that is inspired by the novel approach proposed in the seminal work for mean-field equations [8] and further extended to the homogeneous Boltzmann equation in [40].The proposed approach is capable to combine the efficiency of Direct Simulation Monte Carlo (DSMC) methods for nonlinear kinetic equations in the phase space [1,2,35,36] with the accuracy of sG methods in the parameter space.The DSMC-sG method preserves the main physical properties of the kinetic solution along with spectral accuracy in the random space provided minimal regularity assumptions.Anyway, in a non-Maxwellian framework, the numerical formulation of the DSMC-sG method requires the introduction of step functions at the particles' level.As shown in [40] for the variable hard-spheres (VHS) model, this fact can break spectral convergence of sG.For this reason, it has been shown that a mollification of the step function coupled with a thermalization of particles is capable to restore the physical validity of the model together with spectral accuracy in the random space.
In models for collective phenomena, the equilibrium distribution of Boltzmann-type models is unknown and typically only mass is conserved.For this reason, we introduce a surrogate Fokker-Planck model that can be formally derived from the original model in the quasi-invariant limit [45].In the case of non-Maxwellian interactions, we will obtain a nonlocal nonlinear Fokker-Planck class of equations whose equilibrium distribution can be approximated numerically through suitable deterministic methods [39].In particular, we will investigate the effects of a mollification of step functions introduced at the Monte Carlo level, coupled with a correction of nonconserved quantities computed by the approximation of the corresponding surrogate Fokker-Planck model.Numerical tests for kinetic models of wealth distribution and traffic flow have been performed.
The rest of the paper is organized as follows.In Section 2 we introduce non-Maxwellian models for multi-agent systems with random inputs and we formally derive their corresponding Fokker-Planck models.Regularity of the solutions in the random space of the surrogate models has been investigated in Section 2.2.In Section 3 we briefly review the basic features of some existing kinetic models for pure gambling, wealth distribution and traffic flow dynamics with non-Maxwellian kernels.In Section 4 then we construct the DSMC-sG methods and we provide results on the consistency of the method.Finally, in Section 5 several numerical results are presented which show the efficiency and accuracy of the introduced method.

Non-Maxwellian models with uncertain parameters
To introduce the modelling setting we consider a binary interaction model with uncertain mixing [17,38,46].If two sampled particles that are characterized by the pre-interaction states v, w ∈ V ⊆ R interact, then their post-interaction states v , w ∈ V are obtained following the scheme with > 0 a given constant, I 1 , I 2 suitable interaction functions depending on the pre-interaction states and on the random quantity z ∈ R dz , d z ≥ 1.Furthermore, η is a random variable with zero mean and variance σ 2 , and the functions D 1 and D 2 define the local relevance of the diffusion.Under suitable assumptions on the strength of the diffusion it is possible to show that the postinteraction states v , w remain in V.
We adopt a classical kinetic theory approach based on the one-dimensional space-homogeneous Boltzmann equation, that describes the time evolution of the one-body distribution function f = f (t, v, z).The function f (t, v, z) identifies the state of the system, such that f (t, v, z)dv is the fraction of agents characterized by a state comprised between v and v + dv at time t > 0 and parametrised by uncertainties defined in the random vector z with joint distribution p(z).The evolution of f is given by the non-Maxwellian Boltzmann-type model where ϕ : V → R is a test function, while the symmetric function B(v, w, z) denotes the collision kernel that characterizes the collision frequency of agents with states v and w.The notation • expresses the expectation with respect to the random variable η .The model introduced in (2) can be complemented with uncertain initial condition It is worth noting that in the VHS framework of classical kinetic theory for rarefied gas dynamics the collisional kernel is assumed to be function of the relative velocity |v − w|, see [9].The model (2) greatly simplifies in the Maxwellian case corresponding to B(v, w, z) ≡ 1.In the description of collective models in the Maxwellian simplification of considering a constant interaction kernel has been largely considered [22,37,47,49].Within the Maxwellian simplification it is possible to argue on the existence and uniqueness of large time behavior of the resulting model.In some cases, explicit solutions can be obtained, as in the famous one-dimensional Kac model [32].In models for collective phenomena, a precise analytical description of the kinetic emerging equilibrium distribution is very difficult to obtain.A possible way to overcome this difficulty relies on the possibility to study surrogate models, that are approximations of the kinetic model (2) in some limit and whose large time behavior is easily available.

Fokker-Planck approximation
The computation of the emerging equilibrium density of the Boltzmann model introduced in (2) is very challenging.An established way to overcome this difficulty relies on the introduction of the quasi-invariant limit [11,45,46] under which it is possible to derive a surrogate Fokker-Planck model for the interaction dynamics.The introduced scaling has connections with the grazing collision limit in classical kinetic theory.The main idea is to introduce a new time scale τ = t and to define Hence, scaling the variance of the random variable η as σ 2 = σ 2 we have that for 1 the interaction dynamics in (1) are quasi-invariant, since v − v 1 and w − w 1.Assuming then ϕ smooth enough and at least ϕ ∈ C 3 0 (V), we can perform the following Taylor expansions Plugging the above expansion in (3) we obtain where R ϕ (f , f ) is a reminder term of the following form Thanks to the smoothness assumptions on ϕ and the boundedness of the third order moment of η , in the limit → 0 we have that and we may assume that f converges to a distribution f (τ, v, z) at least formally, we point the interested reader to [7,45] for related approaches in the Maxwellian and Boltzmann-Povzner frameworks.With a slight abuse of notation, we indicate with f (τ, v, z) the limit distribution as → 0 + , hence f (τ, v, z) is weak solution to the nonlinear nonlocal Fokker-Planck equation complemented with the following boundary conditions for all

Regularity of solutions in the random space
We recall that p(z) : I z → R + is the probability density of the random vector z.We define the weighted norm in L 2 p (V × I z ) as follows In the following we will provide sufficient conditions to guarantee regularity of the solution of the general Fokker-Planck model ( 4).Let us rewrite first the Fokker-Planck models (4) as follows with We have Theorem 1.Given B(v, w, z) > 0, let f (t, v, z) be the solution of the Fokker-Planck model (6).
for all t ≥ 0, provided Proof.We multiply by 2f (t, v, z)p(z) the nonlinear nonlocal Fokker-Planck equation ( 6) and we integrate it over v and z: For the integral I we have For the integral II we have = 0 and D[f ] ≥ 0. Hence, we have Thanks to Gronwall's Lemma we obtain Next, we apply the v derivative to both members of (6) We multiply by 2p(z)∂ v f (t, v, z) both members of the latter equation and we integrate over We have and Therefore, we have and from the Gronwall inequality we get corresponding to (8), from which we conclude the proof.
Theorem 2. Given B(v, w, z) > 0, let f (t, v, z) the solution of the Fokker-Planck model (6) and let us consider the constants where Proof.Let us consider the z derivative of the Fokker-Planck model ( 6) with σ 2 = 0 We multiply by 2p(z)∂ z f and we integrate over Hence, we observe that and thanks to the Young's inequality.Hence, we obtained Thanks to the uniform Gronwall inequality, see [44] p.88, we have Taking into account Theorem 1 we obtain from which we conclude.
< +∞, = 0, 1} exploiting the regularity of f , ∂ v f .Anyway, the estimates are not sharp as the ones obtained for linear equations, see e.g.[30,33].Future research efforts will be dedicated to obtain sharper estimates for nonlinear Fokker-Planck equations.

Examples in non-Maxwellian models for collective phenomena
We briefly present three non-Maxwellian kinetic models for collective phenomena, namely a model for pure gamble [3], a model for wealth distribution [23] where the binary scheme is based on the Cordier-Pareschi-Toscani model [11] and, finally, a variation of the traffic model presented in [47] that includes a speed-dependent interaction kernel.

Pure gambling
In the kinetic models for pure gambling the state space is V = R + .Preliminary Maxwellian models have been introduced in [3] in which the nonlinear Boltzmann-type model (2) with B ≡ 1 has been considered.In the pure gambling processes [19], the entire wealth of two agents is at stake at each interaction and randomly shared between agents.Therefore, assuming that the game is fair, the binary interactions are of the type (1) with where ω is a random variable symmetric with respect to 1/2 and we considered = 1.Furthermore, we consider vanishing diffusion functions In an economic framework, an agent with zero wealth cannot gamble.To mimic this situation, in [23] it has been proposed to modify the classical kinetic gamble model of [3] through an interaction kernel of the following form where the exponent of the kernel is an uncertain quantity, i.e. δ = δ(z).For the introduced gambling rules and in presence of the interaction kernel (9), the wealth density f (t, v, z) satisfies a bilinear non-Maxwellian Boltzmann-type equation that in weak form reads It is worth to observe that, for any δ, equation ( 10) conserved mass and momentum.The mass conservation can be easily observed by taking ϕ(v) = 1, whereas for momentum conservation we consider ϕ(v) = v and we get since −ω and ω − 1 are identically distributed.Assuming now ω ∼ U([0, 1]), it is possible to show that for any 0 < δ < 1 the large time distribution of the model ( 10) incorporated the kernel uncertainties and is a Gamma density of the form The uncertain parameter δ characterizing the interaction kernel has a great influence on the large time behavior of the system.Indeed, it is worth to remark that the variance of and inequalities of the money distribution increase with δ and blow up in the limit δ → 1.

Wealth distribution
In recent years, several kinetic models for wealth distribution have been proposed.Also in this case, the state space is V = R + .In the following, we concentrate on the modelling setting proposed in [11] in the case of interaction with a background distribution.In particular, in [23] it is assumed that elementary wealth changes of an agent are determined by interactions (1) with Therefore, the interaction scheme reads The quantity λ ∈ (0, 1] determines the saving propensity and it is assumed η ≥ −1 + λ.Furthermore, in an economic framework, the probability of transactions in which one player has no wealth to exchange is very rare.To this end, in [23] the authors proposed the kernel with δ ∈ (0, 1] and κ > 0. In the following, we will concentrate on the case δ = δ(z).The resulting non-Maxwellian kinetic model in weak form reads and do not conserve mean and energy.In particular, the following estimates hold Information on the large time behavior can be obtained by relying to a Fokker-Planck model approximating the kinetic model ( 12) where adding boundary conditions of the type (5).The equilibrium distribution is now given by where We point the interested reader to [23] for additional details.

Traffic flow
In kinetic traffic modelling, non-constant interaction kernels have been frequently considered, see e.g.[10, 25-27, 41, 43] and the references therein.In the following, we study the influence of a cross section on a traffic model recently proposed in the Maxwellian framework [47] in which The time evolution of the distribution f (t, v, z) is determined by microscopic binary interactions responsible for speed changes.Given normalized pre-interaction speeds (v, w) ∈ [0, 1] × [0, 1], the post-interaction speeds (v , w ) are determined by (1) with being P (ρ, z) = (1 − ρ) µ(z) ∈ [0, 1], µ(z) > 0, the probability to accelerate with a traffic density ρ ∈ [0, 1].The presence of uncertain quantities in P (ρ, z) is associated with different responses of vehicles in heterogeneous traffic conditions, see [48].Furthermore, we consider D 1 (v) = D(v, ρ) and D 2 (v, z) = 0. Hence, the speed changes are determined by we point the interested reader to [47] for further details on the modeling setting.The choice of the function D has to ensure that v , w ∈ [0, 1] for any v, w ∈ [0, 1].In [47] the authors proposed Hence, the evolution of the density f follows the non-Maxwellian Boltzmann-type equation in weak form Since a priori information on the frequency of interaction is missing, it seem reasonable to introduce an additional uncertain exponent z of the cross section as an uncertain quantity.
From equation (17) we can compute the evolution of mean speed V (t, z) of the system.We fix ϕ(v) = v and the evolution of whose solution depends on the specific interaction kernel considered.For Maxwellian particles, corresponding to the choice α(z) ≡ 0 we recall the results of [47,48] from which we are able to find a close equation for the time evolution of the mean speed V (t, z).
For any other α(z) > 0, the mean speed V (t, z) depends on higher momenta and its evolution cannot be expressed in closed form.To investigate the influence of the uncertain interaction kernel on the evolution of V (t, z), we fix α(z and Let us consider now α(z) = 2.We define the following constants C 1 (z) = P (ρ, z) ≥ 0 and with V (0, z) = V 0 ∈ [0, 1].We point the interested reader to Appendix A for additional details.The impact of the interaction kernel (18) on the introduced traffic model is shown in Figure 1.
As in Section 2.1, the Boltzmann model ( 17) can be approximated through a surrogate Fokker-Planck model in the quasi-invariant limit.In particular, for the introduced traffic model with interaction kernel we get where I has been defined in (15) and B is the interaction kernel as in (18).For the introduced Fokker-Planck equation, at variance with the cases in Sections 3.1-3.2,we cannot compute analytically the equilibrium distribution of the Fokker-Planck model unless α(z) = 0, corresponding to the Maxwellian scenario.

DSMC stochastic Galerkin methods
In this section, we revise the construction of a stochastic Galerkin version of the classical DSMC Algorithm for non-Maxwellian particles, see e.g.[36,37].In more detail, we extend the Direct Simulation Monte Carlo stochastic Galerkin (DSMC-sG) methods, introduced in the gas dynamic framework [40], to the models with uncertain parameters proposed in Section 3. Next, we provide consistency results of the DSMC-sG algorithm with respect to relevant observables and in the reconstruction of the kinetic density.

DSMC-sG for non-Maxwellian models with uncertainties
We first rewrite (2) in strong form to highlight the gain and loss part of the Boltzmann-type equation: where J is the absolute value of the Jacobian of the considered transformation.We denote by Q Σ the operator obtained replacing the kernel B(v, w, z) with B Σ (v, w, z) given by where Σ = max{B(v, w, z)} > 0 is an upper bound for the interaction kernel.By decomposing Q Σ (f, f ) in its gain and loss part we can rewrite the interaction step as Let us now consider a time interval [0, T ] and let us discretize it in n TOT interval of size ∆t > 0. We denote by f n (v, z) the approximation of f (t n , v, z) and we consider the forward Euler scheme where f n+1 is a probability density provided µΣ∆t ≤ 1.
Then, we consider a sample of N particles v i (z, t), i = 1, . . ., N , from the kinetic solution of the Boltzmann model at time t, and we approximate v i (z, t) by its generalized polynomial chaos expansion where {Φ h (z)} M h=0 , M ∈ N, is a set of orthogonal polynomials of degree less or equal to M , orthonormal with respect to the probability density function p(z) where I z is the sample space and δ hk is the Kronecker delta.The choice for the orthogonal polynomials obviously depends on the distribution of the parameters and follows the so-called Wiener-Askey scheme [50,51].In (26), vi,h (t) is the projection of the velocity in the subspace generated by the polynomial of degree h = 0, . . ., M vi,h (t) = To perform collision with a non-Maxwellian kernel, we may rewrite the general binary interaction scheme (1) for two particles v i = v i (z, t), w j = w j (z, t), highlighting the acceptance-rejection process introduced by the classical Nanbu-Babovski method [36] where χ(•) is the indicator function and ξ a uniform random number in (0, 1).Then, we substitute the velocities v i , w i with their gPC approximation v M i , w M j and we project against Φ h (z)p(z)dz on I z for every h = 0, . . ., M .We obtain where V h i,j , Ŵ h i,j are the so-called collisional matrices We stress the fact that the new binary interaction for the projections (30) does not depend on the uncertain parameter z.The DSMC-sG method is summarised in Algorithm 1.

Sround =
x + 1 with probability x − x x with probability 1 − x + x , where x denotes the integer part of x.

Consistency estimates
We want to evaluate the error produced by the DSMC-sG algorithm in the reconstructed distribution function and its moments.In the following, we denote by f (t, v, z) the solution of the Boltzmann equation ( 3) with binary updates (1) and by f (t, v, z) the corresponding mean field approximation, weak solution of the Fokker-Planck equation ( 4).We introduce then the empirical density functions where δ(•) is the Dirac delta function.Being ϕ any test function, we denote by its expectation with respect to the distribution function f (t, v, z), so that we have From the central limit theorem we have the following result [5] Lemma 1.If we denote by E V [•] the expectation with respect to f in the velocity space, for each z the root mean square error satisfies from the polynomial approximation theory [50], we have the following spectral estimate Lemma 2. For any v(z) ∈ H r p (I z ), r ≥ 0, there exists a constant . Next, for any random variable W (z, t) taking values in L 2 p (I z ), we define p (Iz) , and equivalently , see [16,38,40].Then, we have the following result Theorem 3. Let f (t, v, z) be a probability density function, solution of the weak Fokker-Planck equation (4), and f M ,N (t, v, z) the empirical measure obtained from the N -particles sG approximation {v M i (z, t)} i , solution of the time scaled Boltzmann equation (3).If v i (z, t) ∈ H r p (I z ) for every i = 1, . . ., N , and in the quasi-invariant limit → 0, we have the following estimate where ϕ is a test function, C > 0 is a constant independent on M and ξ i = (1−θ)v i +θv M i , θ ∈ (0, 1).
Proof.Thanks to the triangular inequality we have In the quasi-invariant regime → 0, up to the extraction of a subsequence, we have as a consequence we have lim and the first term vanishes in the quasi-invariant limit.The second term can be evaluated exploiting the result of Lemma 1.Therefore, we have Finally, we have and from the mean value theorem ϕ Next, we introduce a uniform grid in the domain V, with ∆v > 0 width of the cell, and we denote by S ∆v ≥ 0 a smoothing function such that ∆v We consider the approximations of the density function obtained by observing that the standard histogram reconstruction corresponds to the choice S ∆v (v) = χ(|v| ≤ ∆v/2)/∆v.Defining we have the following result Theorem 4. The error introduced by the reconstruction of the distribution in the DSMC-sG method, in the grazing limit → 0, satisfies , Proof.Thanks to the triangular inequality we have IV .
In the limit → 0 we have shown that f → f , so the first term vanishes.The second term represents the error introduced by the density reconstruction and is bounded by , where q > 0 depends on the accuracy of the reconstruction.For the last two terms, we observe that with ϕ(•) = S ∆v (v − •).Hence, we can apply the result of Theorem 3 with the just mentioned choice for ϕ.

Numerical results
In this section, we present several numerical tests for the DSMC-sG scheme for the non-Maxwellian models with uncertainties described in Section 3. In all the subsequent tests we will consider N = 10 5 agents and the densities are reconstructed through standard histograms.
In more detail, we first check the consistency of the DSMC-sG approximation of the Boltzmann equation with the exact equilibrium distribution for the kinetic model for gambling.Then we test the convergence to the equilibrium of the Fokker-Planck equation for the discussed wealth distribution and the traffic models.
In the binary interactions (29) we introduced the indicator function χ(•).This term may deteriorate the overall convergence of the DSMC-sG scheme.Coherently with the approach proposed in [40], we introduce the following regularisation where K(β(•)) is a sigmoid function dependent on the parameter β > 0. In particular, we consider With this choice, we note that β 1 is associated with a sharp sigmoid function, on the contrary, a smaller β is linked to a smoother sigmoid.We will return to the influence of the parameter β in the following.
This regularisation induces a different evolution of the relevant observables.Consequently, to keep the exact time evolution of the first two moments together with the spectral convergence, we couple the DSMC-sG scheme with a scaling process of the form where V FP (t, z) and E FP (t, z) are, respectively, the mean velocity and the energy computed from the corresponding surrogate Fokker-Planck model.Similarly, we indicated with V K (t, z) and E K (t, z) mean velocity and energy of the Boltzmann-type model with sigmoid function K(•) (34) in equation (33).The computation of V FP (t, z) and E FP (t, z) follows from the model ( 4) that is solved through standard sG method and for which we can guarantee sufficient regularity under the assumptions of Theorem 2. We highlight how the additional scaling process (35) will be consistent with the original Boltzmann-type model in the regime 1.For clarity purposes, in the rest of the section we will indicate the numerical solution of the Fokker-Planck model f FP (t, v, z) and the numerical solution of the Boltzmann-type model as f (t, v, z), whereas the solution of the Boltzmann-type model with additional scaling process (35) will be denoted by f (t, v, z).For the approximation of the Fokker-Planck model we will implement a standard sG collocation method based on a semi-implicit scheme presented in [39] and further studied in [17,48,52].

Test 1: gambling
We consider the kinetic model for gambling with 1D uncertainty in the collisional kernel.We choose δ(z) = 1 2 z and z ∼ U([0, 1]) in ( 9) fixing κ = 1.Since the random parameter is uniformly distributed, we use the Legendre polynomials in the gPC expansion.In all the simulations, we consider M = 5 Galerkin projections and the time frame [0, T ] with T = 10 discretised with time step ∆t = 0.1.The kinetic density is reconstructed in the interval [0, 10] with ∆v = 0.05.We consider the deterministic initial distribution In this test, we highlight that the equilibrium solution of the Boltzmann-type model can be computed exactly and = 1 and reads see Section 3.1.In Figure 2, we report expected value and variance with respect to the random parameter z of the DSMC-sG approximation.We may observe the good agreement of the considered quantities of interest with the analytical ones.

Test 2: wealth distribution
We consider now the kinetic model for the wealth distribution described in Section 3.2.In particular, we consider the case where the interaction kernel (11) is characterized by δ(z) = z ∼ U([0, 1]) and κ = 1.Therefore, we adopt the Legendre polynomials in the gPC expansion of the velocities.In all the results of this test, we consider a background uniformly distributed as E ∼ U([0.9, 1.1]).Furthermore, we consider the following deterministic initial distribution In Figure 3 we show expectation and variance of f (t, v, z) computed through DSMC-sG method with respect to the analytical equilibrium distribution of the Fokker-Planck model ( 14), for various = 5 × 10 −2 , 10 −1 , 5 × 10 −1 .In the last picture we report also the behavior of the expected mean wealth E z [V (t, z)] for the introduced values of .We consider M = 5 Galerkin projections, ∆t = /10, λ = σ 2 = 0.5 and time frame [0, T ] with T = 10.The kinetic density is reconstructed through standard histogram over the interval [0, 10] with ∆v = 0.05.
In order to show spectral convergence property of the scheme, we consider a reference DSMC-sG evolution of E (t, z) obtained with = 0.1, N = 10 5 , ∆t = 0.1 and sG scheme up to order M = 50.We store the collisional tree generating the reference solution and we check the L 2 convergence of the scheme.In Figure (4) we present the decay of the L 2 error for increasing M obtained from the initial distribution (38).If we consider the original binary dynamics (29), even if the expectation is well described, it can be observed that the spectral accuracy of the method is lost due to discontinuity of the indicator function χ(•).The same test performed for the binary dynamics (33) recovers spectral accuracy.For increasing β 0 the convergence of the scheme is deteriorated, since we approximate a step function.
Coupling now (33) with the process (35), we recover qualitatively consistent approximation of the evolution of relevant quantities of interest together with spectral convergence for moderate values of β > 0, see Figure 5.In this test we solve the Fokker-Planck model (13).14) with the DSMC-sG approximation of f (t, v, z) (regularization without rescaling) and of f (t, v, z) (regularization with rescaling), in terms of expectation in z.Center and right: comparison of V FP , E FP with the DSMC-sG approximation of V (t, z), E (t, z) (regularization without rescaling) and of Ṽ (t, z), Ẽ (t, z) (regularization with rescaling).We consider N = 10 5 , M = 5, ∆t = 0.01 and = 0.1.We fix κ = 1, λ = σ 2 = 0.5

Test 3: traffic flow
In this last test, we consider the traffic model of Section 3.3, affected by an uncorrelated 2D random parameter z = (z 1 , z 2 ) with p(z) = p 1 (z 1 )p 2 (z 2 ).In particular, we consider z 1 affecting µ(z 1 ) in the interaction function I(v, w, z 1 ) defined in (15) and z 2 affecting α(z 2 ) in the kernel B(|v − w|, z 2 ) defined in (18).Under these assumptions, the gPC expansion of the velocities being {Φ h (z 1 )} M1 h=0 and {Φ k (z 2 )} M2 k=0 the polynomials orthogonal with respect to the distributions p 1 (z 1 ) and p 2 (z 2 ), respectively.Substituting v M1,M2 i (z 1 , z 2 , t) into the binary interaction (16) and and proceeding as in Section 4.1, we obtain with the following collision matrix We consider the following deterministic initial distribution To assess the impact of the single uncertain parameters on the dynamics we first consider the case µ(z 1 ) = 1 + 2z 1 and α(z 2 ) = 2z 2 with uncorrelated uncertainties z 1 , z 2 ∼ U([0, 1]).In Figure 6 we show the DSMC-sG approximation of the solution of Fokker-Planck model for traffic (23) in terms of expected value and variance in z = (z 1 , z 2 ) of the distribution function and of the macroscopic quantities.We considered two different densities ρ = 0.4 and ρ = 0.6 and the Fokker-Planck is solved on a grid of N v = 51 points such that ∆v = 0.02 and ∆t = ∆v/2.As before, the  ]), with the DSMC-sG reconstruction of traffic distributions (first and second column) and of mean velocity (third column) for various = 0.05, 0.1, 0.5.We considered ρ = 0.4 (top row) and ρ = 0.6 (bottom row).We set for the DSMC-sG method N = 10 5 , M 1 = M 2 = 5, ∆t = , and the deterministic solver for Fokker-Planck is such that ∆v = 0.02 and ∆t = ∆v/2.The time frame is t ∈ [0, T ], T = 300.
DSMC-sG provides a good approximation in the limit 1 of the solution of the Fokker-Planck model.
We show the L 2 convergence of the DSMC-sG scheme in Figure 7.In details, we considered the case with binary interactions (29) in the left plot, whereas the case with regularization of the step function as in (33), with β = 0.01, is presented in the right plot.The error has been computed with respect to a reference DSMC-sG evolution of E (t, z) with ∆t = = 0.1, N = 10 5 , and M 1 = M 2 = 50.As before, in this test we store the collisional tree of the reference solution and we check L 2 convergence for increasing M 1 , M 2 .The error is presented here in log 10 and we can clearly observe spectral accuracy in the case with regularization.
Finally, coupling (33) with the process (35), we recover a qualitatively consistent approximation of the evolution of relevant quantities of interest in the case of traffic flow model, see Figure 8.

Conclusion
In this work, we studied an extension of a recently introduced DSMC-sG hybrid approach [8,40] for uncertainty quantification of kinetic equations to non-Maxwellian Boltzmann-type models for multi-agent systems.The proposed method combines a DSMC solver in the physical space with a stochastic Galerkin method in the random space and is based on a generalized Polynomial Chaos expansion of statistical samples of a DSMC solver.The DSMC-sG solution of non-Maxwellian models with uncertainties requires a suitable reformulation of classical DSMC solvers.The class of kinetic models of interest can be formally approximated by surrogate Fokker-Planck-type models in the quasi-invariant regime.For these models, the regularity in the random space has been  Similarly, for the second term of (44) we have  Hence, we get V − (t, z) ≤ V (t, z) ≤ V + (t, z) where where Hence, we consider the case α(z) = 2.We define the following constants C 1 (z) = P (ρ, z) ≥ 0 and C 2 (z) = 2 (P 2 (ρ, z) + 1 − P (ρ, z)) > 0 and from (19) we get where M 3 (t, z) = Arguing as before, we get Therefore we get V − (t, z) ≤ V (t, z) ≤ V + (t, z) where + e (C1(z)+C2(z))t 1 with V (0, z) = V 0 ∈ [0, 1].

Figure 1 :
Figure 1: Upper and lower bounds for the time evolution of the mean velocity V (t, z) for α(z) = 1 (red line) and α(z) = 2 (blue line), for different values of the traffic density ρ and the parameter µ(z).The initial velocity is V 0 = 1 − ρ.