Mean-field control variate methods for kinetic equations with uncertainties and applications to socio-economic sciences
MMean-field control variate methods for kinetic equations withuncertainties and applications to socio-economic sciences
Lorenzo Pareschi ∗ , Torsten Trimborn † , Mattia Zanella ‡ February 5, 2021
Abstract
In this paper, we extend a recently introduced multi-fidelity control variate for theuncertainty quantification of the Boltzmann equation to the case of kinetic models arisingin the study of multiagent systems. For these phenomena, where the effect of uncertaintiesis particularly evident, several models have been developed whose equilibrium states aretypically unknown. In particular, we aim to develop efficient numerical methods basedon solving the kinetic equations in the phase space by Direct Simulation Monte Carlo(DSMC) coupled to a Monte Carlo sampling in the random space. To this end, exploitingthe knowledge of the corresponding mean-field approximation we develop novel mean-field Control Variate (MFCV) methods that are able to strongly reduce the varianceof the standard Monte Carlo sampling method in the random space. We verify theseobservations with several numerical examples based on classical models , including wealthexchanges and opinion formation model for collective phenomena.
Keywords: uncertainty quantification, kinetic equations, mean field approximations,control variate methods, Monte Carlo methods, stochastic sampling, multi-fidelity meth-ods
Contents ∗ Department of Mathematics and Computer Science, University of Ferrara, Italy( [email protected] ) † Institut f¨ur Geometrie und Praktische Mathematik, RWTH Aachen University, Germany( [email protected] ) ‡ Department of Mathematics ”F. Casorati”, University of Pavia, Italy ( [email protected] ) a r X i v : . [ m a t h . NA ] F e b Mean Field Control Variate DSMC methods 13
In recent years kinetic theory emerged as a sound theoretical framework to describe a widerange of collective phenomena in socio-economy and life sciences for systems composed by asufficiently large number of agents [9, 12, 13, 18, 19, 21, 46, 48]. For an introduction to thesetopics we refer to the recent surveys and monographs [4, 5, 32, 37].Nevertheless, the design of realistic models for the description of human behavior has toface the lack of first principles and the dynamics are often inferred from empirical observationsbased on experimental data [3, 7, 26]. Especially in the field of socio-economic applicationsthe precise form of the microscopic interactions is largely unknown. Thus, one typically con-structs microscopic social forces which are able to qualitatively fit the macroscopic behavior ofthe system. In the context of kinetic modelling, this issue can be translated in structural un-certainties present both in initial observations and interaction rules, which can be consideredin the form of uncertain parameters of the model depending on random quantities. Hence,the quantification of uncertainties in such social microscopic interactions on the observablebehavior is of major importance.The introduced uncertainties inevitably increase the dimensionality of the problems. There-fore, we need to develop new numerical methods to efficiently quantify the impact of unknownquantities on the overall dynamics. Among the various methods for uncertainty quantifica-tion we find intrusive stochastic Galerkin (sG) methods which provide spectral convergencetowards the solution of the problem under suitable regularity assumptions [27, 43, 53]. BesidesG methods we find non-intrusive approaches for UQ which does not require strong mod-ification of the numerical scheme for the deterministic problem, like stochastic collocationmethods. Such methods are non-intrusive, easy to parallelize [52] and, in principle, do notrequire any knowledge of the probability distribution of the uncertain parameter. We referthe interested reader to [17, 22, 27–29, 34, 41, 42, 51] for an additional overview on numericalmethods for uncertainty quantification of hyperbolic and kinetic equations.In this work we concentrate on a non-intrusive approach based on a stochastic Monte Carlo(MC) sampling for Boltzmann-type equations. In in comparison to sG methods, techniquesbased on MC sampling have a lower impact on the curse of dimensionality [8, 24, 30, 31, 37].In details, following the methodology recently introduced in [15,16] for rarefied gas dynamics,we develop a variance reduction method based on a control variate approach which exploits amicro-macro-type decomposition of the solution. In the simplest setting proposed in [15] as asurrogate model to reduce the variance of the MC estimator, the corresponding Maxwelliansteady state solution has been used. In many socio-economic applications, however, theequilibrium states of the Boltzmann models are unknown and therefore it is an open problemthe determination of a suitable surrogate model that can be used as control variate.2n the present paper, we propose to use as surrogate model the corresponding mean fieldmodel of Fokker-Planck type obtained as an approximation of the original Boltzmann modelin a grazing-type limit [37, 49]. In fact there are several reasons for adopting this controlvariate. First, if the mean-field model’s steady state is known, we can use it directly as acontrol variate. Second, more in general, we can use the whole solution of the Fokker-Planckmodel as a control variate. This is possible thanks to the recently introduced structurepreserving schemes for Fokker-Planck equations that preserve exactly the steady state of theequation [39]. The latter choice can be generally applied whereas for the former the steadystate of the mean field model needs to be known. We call this novel variance reductionmethod, mean field control variate (MFCV) method.Furthermore, our mean field control variate approach is not based on a deterministic solverof the Boltzmann model but on a Direct Simulation Monte Carlo (DSMC) solver. This isbeneficial since DSMC solver are widely used in order to solve Boltzmann type equations[11, 35, 37] and can be easily applied to generalized Boltzmann model simply by exploitingthe rules defining the microscopic dynamics. This is an additional difficulty since our methodneeds to couple the deterministic solution of the mean field model with the DSMC solver of theBoltzmann equation. In this respect, our MFCV method can be regarded as a hybrid methodin the phase space. The crucial assumptions of the MFCV method are the following: First, themean field model has to be less expensive to solve than the original kinetic model. Secondly,the collision regime of the space homogeneous Boltzmann model needs to be sufficiently closeto the grazing collision regime.The rest of the manuscript is organized as follows: In the next section we introduce ageneral Boltzmann model well suited for socio-economic applications. Furthermore, we recallthe derivation of the corresponding mean field model and discuss some examples. In Section3 we introduce the basics of DSMC methods and shortly discuss the MC sampling method foruncertainty quantification. Then we present in Section 4 the novel MFCV method and em-phasize its potential advantages in reducing the variance of the estimator used for uncertaintyquantification. Several examples in the subsequent section demonstrate the advantages of theMFCV method in comparison to the standard MC technique. We finish this study with ashort discussion of the present method and comment on further developments.
Let us consider a general binary interaction model with uncertain mixing [17, 47]. The pairof interacting agents is characterized by the pre-interaction states v, w ∈ V ⊆ R and thepost-interaction states v (cid:48) , w (cid:48) ∈ V are obtained as follows v (cid:48) = v + ε [( p ( z ) − v + q ( z ) w ] + D ( v, z ) η ε ,w (cid:48) = w + ε [ p ( z ) v + ( q ( z ) − w ] + D ( w, z ) η ε , (1)where ε > p i , q i , i = 1 , , are suitable interaction functions dependingon a random variable z ∈ Ω ⊆ R d z , d z ≥
1. Furthermore, η ε is a random variable with zeromean and variance σ ε and the function D ( · , z ) is the local relevance of the diffusion. We willassume that the post-interaction states w (cid:48) , v (cid:48) remain in the set V up to the introduction ofsuitable conditions on η ε , see [46].To describe the evolution of a large system of agents undergoing binary interaction weadopt a kinetic approach. Hence, we introduce the distribution function f = f ( t, w, z ), such3hat f ( t, w, z ) dw is the fraction of agents, at time t ≥
0, characterized by a state comprisedbetween w and w + dw and parametrized by the uncertainty z . The evolution of f is givenin terms of the Boltzmann-type model for Maxwellian interactions ddt (cid:90) V ϕ ( w ) f ( t, w, z ) dw = 1 ε (cid:28)(cid:90) (cid:90) V ( ϕ ( w (cid:48) ) − ϕ ( w )) f ( t, w, z ) f ( t, v, z ) dwdv (cid:29) , (2)being ϕ : V → R any observable quantity which may be expressed as a function of themicroscopic state w of the agents. The symbol (cid:104)·(cid:105) denotes the expectation with respect to η ε . The model (2) can be also complemented with uncertainties on the initial condition f (0 , w, z ) = f ( w, z ).We can recast model (2) in symmetric form as follows ddt (cid:90) V f ( t, w, z ) ϕ ( w ) dw = 12 ε (cid:28)(cid:90) (cid:90) V ( ϕ ( w (cid:48) ) + ϕ ( v (cid:48) ) − ϕ ( w ) − ϕ ( v )) f ( t, w, z ) f ( t, v, z ) dwdv (cid:29) (3)Taking ϕ ( w ) = 1 is easily seen that the number of agents is conserved in time. Theevolution of the mean is obtained for ϕ ( w ) = w which gives ddt m f ( t, z ) = 12 (cid:90) (cid:90) V [( p ( z ) + p ( z ) − v + ( q ( z ) + q ( z ) − w ] f ( t, w, z ) f ( t, v, z ) dw dv, and the mean is conserved for p ( z ) + p ( z ) = 1, q ( z ) + q ( z ) = 1, indeed at the microscopiclevel the binary mean is conserved provided (cid:104) w (cid:48) + v (cid:48) (cid:105) = v + w + ε [( p ( z ) + p ( z ) − v + ( q ( z ) + q ( z ) − w ] = v + w. Notice that any nonconserved moment of the distribution function f explicitly depends on z . An extensive qualitative study of the previously introduced Boltzmann model is very chal-lenging. In particular, the asymptotic behavior of socio-economic Boltzmann models like (2)are unknown. Therefore we employ the quasi-invariant interaction limit [13, 36, 46] to derivea mean-field approximation of the Boltzmann model.In kinetic theory a mean field model can be derived by the grazing limit [20, 44, 49] of theBoltzmann equation. The main idea is to introduce a scaling parameter such that we havea small change in the kinetic variable after each interaction, while keeping the macroscopicproperties of the model unchanged.Let us introduce a time scaling parameter ε > τ = εt, f ε ( τ, w, z ) = f ( τ /ε, w, z ) . (4)Then, the distribution f ε is solution to ddτ (cid:90) V ϕ ( w ) f ε ( τ, w, z ) dw = 12 ε (cid:28)(cid:90) (cid:90) V ( ϕ ( w (cid:48) ) + ϕ ( v (cid:48) ) − ϕ ( v ) − ϕ ( w )) f ε ( τ, w, z ) f ε ( τ, v, z ) dwdv (cid:29) . (5)4ence, scaling the variance of the introduced random variables as σ ε = εσ we can observethat for ε (cid:28) w (cid:48) − w and v (cid:48) − v are small. Assuming now ϕ sufficiently smooth and at least ϕ ∈ C ( V ) we can perform thefollowing Taylor expansions ϕ ( w (cid:48) ) − ϕ ( w ) = ( w (cid:48) − w ) ∂ w ϕ ( w ) + 12 ( w (cid:48) − w ) ∂ w ϕ ( w ) + 16 ( w (cid:48) − w ) ∂ w ϕ ( ¯ w ) ,ϕ ( v (cid:48) ) − ϕ ( v ) = ( v (cid:48) − v ) ∂ v ϕ ( v ) + 12 ( v (cid:48) − v ) ∂ v ϕ ( v ) + 16 ( v (cid:48) − v ) ∂ v ϕ (¯ v ) , where ¯ w ∈ (min { w, w (cid:48) } , max { w, w (cid:48) } ), ¯ v ∈ (min { v, v (cid:48) } , max { v, v (cid:48) } ). Plugging the above ex-pansion in (5) we have ddτ (cid:90) V ϕ ( w ) f ε ( τ, w, z ) dw = 12 ε (cid:104) (cid:28)(cid:90) (cid:90) V [( w (cid:48) − w ) ∂ w ϕ ( w ) + ( v (cid:48) − v ) ∂ v ϕ ( v )] f ε ( τ, w, z ) f ε ( τ, v, z ) dv dw (cid:29) + 12 (cid:28)(cid:90) (cid:90) V [( w (cid:48) − w ) ∂ w ϕ ( w ) + ( v (cid:48) − v ) ∂ v ϕ ( v )] f ε ( τ, w, z ) f ε ( τ, v, z ) dv dw (cid:29) + R εϕ ( f ε , f ε ) , where R εϕ ( f ε , f ε ) is a reminder term with the following form R εϕ ( f ε , f ε ) = 16 ε (cid:28)(cid:90) (cid:90) V [( w (cid:48) − w ) ∂ w ϕ ( w ) + ( v (cid:48) − v ) ∂ v ϕ ( v )] f ε ( τ, w, z ) f ε ( τ, v, z ) dv dw (cid:29) . Thanks to the assumed smoothness we can argue that ϕ and its derivatives are bounded in V . Furthermore, since η ε has bounded moment of order three (cid:10) | η ε | (cid:11) < + ∞ , we can observethat in the limit ε → + we have | R εϕ ( f ε , f ε ) | → . Therefore, in the limit ε → + , it can be shown that f ε converges, up to subsequences, toa distribution function ˜ f = ˜ f ( τ, w, z ) which is weak solution to the following Fokker-Planckequation ∂ τ ˜ f ( τ, w, z ) + ∂ w (cid:20)(cid:18)(cid:90) V P ( v, w, z ) ˜ f ( τ, v, z ) dv (cid:19) ˜ f ( τ, w, z ) (cid:21) = σ ∂ w ( D ( w, z ) ˜ f ( τ, w, z )) , (6)where P ( v, w, z ) = 12 [( p ( z ) + q ( z ) − w + ( p ( z ) + q ( z )) v ]provided the following boundary conditions are satisfied for all z ∈ Ω − (cid:18)(cid:90) V P ( v, w, z ) ˜ f ( τ, v, z ) dv (cid:19) ˜ f ( τ, w, z ) + σ ∂ w ( D ( w, z ) ˜ f ( τ, w, z ) (cid:12)(cid:12)(cid:12)(cid:12)(cid:12) w ∈ ∂V = 0 D ( w, z ) ˜ f ( τ, w, z ) (cid:12)(cid:12)(cid:12)(cid:12)(cid:12) w ∈ ∂V = 0 . (7)5he asymptotic analysis of equation (6) is considerable simpler compared to the originalBoltzmann-type model and convergence towards a unique equilibrium distribution can beobtained under suitable hypotheses, see [10, 21, 45]. Therefore, we have obtained a surrogatemodel with reduced complexity whose large time behavior can be more easily studied. Inparticular, the steady state distribution ˜ f ∞ ( w, z ) of (6) is determined by imposing σ ∂ w ( D ( w, z ) ˜ f ∞ ( w, z )) = ∂ w (cid:20)(cid:18)(cid:90) V P ( v, w, z ) ˜ f ∞ ( v, z ) dv (cid:19) ˜ f ∞ ( w, z ) (cid:21) . (8)In fact, for many models the solution of the differential equation (8) is known analytically. Inthe next paragraph we present some relevant examples of socio-economic Boltzmann modelswhich fit in the previously introduced general model. We shortly present two socio-economic models namely a rather general opinion formationmodel [46] and the Cordier-Pareschi-Toscani (CPT) model [12]. As we will see, the modelsare characterized by different steady states in their mean-field approximation.We consider first a model for opinion formation where V = [ − ,
1] and the binary interac-tion rules can be framed in (1) in the symmetric case, i.e. putting in evidence the additional de-pendence of interaction by the opinion variable p = q = q ( | v − w | , z ), p = q = p ( | v − w | , z ),with additionally q ( | v − w | , z ) = 1 − p ( | v − w | , z ). Hence, the binary interaction scheme sim-plifies to v (cid:48) = v + εp ( | w − v | , z )( w − v ) + D ( v, z ) η ε ,w (cid:48) = w + εp ( | v − w | , z )( v − w ) + D ( w, z ) η ε , (9)where the function 0 ≤ p ( | v − w | , z ) ≤ | v − w | . In the present case, in order to produce post-interaction opinionsin the interval [ − , η ε should be such that for all z ∈ Ω and v, w ∈ V we have − − v − εp ( | v − w | , z )( w − v ) ≤ D ( v, z ) η ε ≤ − v − εp ( | v − w | , z )( w − v ) . In particular for D ( w, z ) = D ( w ) = 1 − w and p ( | v − w | , z ) = p ( z ) we obtain the followingbound | η | (1 + | w | ) ≤ (1 − ε max { p ( z ) } ). In the quasi-invariant opinion limit described aboveand provided the boundary conditions (7) are satisfied we reduce to study the following model ∂ t ˜ f ( t, w, z ) + ∂ w (cid:20)(cid:18)(cid:90) V p ( | v − w | , z )( v − w ) ˜ f ( t, v, z ) dv (cid:19) ˜ f ( t, w, z ) (cid:21) = σ ∂ w ( D ( w, z ) ˜ f ( t, w, z )) . In this case, in the interaction scheme (9) we easily see that the mean opinion m ( z ) isconserved in time and, if we consider also uncertainties in the initial distribution, the steadystate distribution of the Fokker-Planck model reads˜ f ∞ ( w, z ) = C ( z ) (1 + w ) − p ( z ) m ( z )2 σ (1 − w ) − − p ( z ) m ( z )2 σ exp (cid:26) − p ( z )(1 − m ( z ) w ) σ (1 − w ) (cid:27) , (10)6here C ( z ) is a normalization factor such that (cid:82) V ˜ f ∞ ( w, z ) dw = 1. Other form of equilibriumdistribution can be determined, for example the choice D ( w ) = √ − w produce a Beta-typesteady state of the form˜ f ∞ ( w, z ) = 2 − σ (cid:16) m ( z ) σ , − m ( z ) σ (cid:17) (1 + w ) m ( z ) σ − (1 − w ) − m ( z ) σ − , (11)where B( · , · ) is the Beta function. We refer to [38, 46] for a detailed discussion.The second example that we mention has been presented in [13] and models wealth ex-changes between agents composing a simple economy. The wealth variable is here allowed totake values on the positive half line therefore V = R + . Again, the considered binary inter-actions can be framed in (1) in the symmetric case, i.e. p = q = q ( z ) and p = q = λ ( z )with q ( z ) = 1 − λ ( z ). Furthermore, we consider here D ( w, z ) = w . The uncertain parameter λ ( z ) ∈ [0 ,
1] determines the proportion of wealth that a single agents wants to invest, thequantity 1 − λ ( z ) is the so-called saving propensity. Thus, the binary scheme reads for wealthexchanges reads v (cid:48) = (1 − ελ ( z )) v + ελ ( z ) w + vη ε w (cid:48) = (1 − ελ ( z )) w + ελ ( z ) v + wη ε , (12)with | η ε | ≤ (1 − ε ). Additionally, we assume that there is no uncertainty in the initialconditions of our Boltzmann-type model. Hence, the corresponding mean field model providedreads ∂ t ˜ f ( t, w, z ) + λ ( z ) ∂ w (cid:104) ( m ˜ f ( z ) − w ) ˜ f ( t, w, z ) (cid:105) = σ ∂ w ( w ˜ f ( t, w, z )) . We observe that the uncertain mean wealth m ˜ f ( z ) is conserved in time. Hence, the equilibriumstate can be computed and reads˜ f ∞ ( w, z ) = ( µ ( z ) − µ ( z ) Γ( w ) w µ ( z ) exp (cid:32) − ( µ ( z ) − m ˜ f ( z ) w (cid:33) , (13)where Γ( · ) denotes the gamma function and µ ( z ) := 1 + λ ( z ) σ . Notice that in this case thesteady state exhibits tails with polynomial decay determined by the uncertain quantity µ ( z ). In this section we introduce a MC-DSMC method to solve the uncertain Boltzmann equationwhere both the physical variables as well as the uncertain parameters are solved by MonteCarlo approximations. The realization of the method represents the starting point for theconstruction of our mean-field control variate strategy.
The efficient computation of the highly non-linear Boltzmann model is a major task and hasbeen tackled by several studies [14, 35, 40]. In deterministic methods the multi-dimensionalintegral of the collision kernel needs to be approximated by a quadrature formula which suf-fers the so-called curse of dimensionality. Additionally, preservation of the main physical7uantities is a true challenge at the discrete level that makes the scheme design model de-pendent. On the other hand, Monte Carlo methods for kinetic equations naturally employthe microscopic dynamics to satisfy the physical constraints and are much less sensitive tothe curse of dimensionality [8]. The most popular examples of Monte Carlo methods for theBoltzmann equation are the classical Direct Simulation Monte Carlo (DSMC) methods byBird and Nanbu [6, 33, 35]. The convergence of the methods to the solution of the Boltzmannequation has been rigorously proven in [2, 50]. For a detailed introduction to DSMC solvers,especially for socio-economic Boltzmann type equations of the form (2), we refer to [37].In order to summarize the classical DSMC algorithms for simplicity we focus on the casewithout uncertainty. We are interested in the evolution of the density f = f ( w, t ) solutionof (5) with initial condition f (0 , w ) = f ( w ). We recall the symmetrized version of thesimulation algorithm originally proposed by Nambu [33], in a similar way one can considerNanbu’s algorithm [6]. Let us consider a time interval [0 , T ], and let us discretize it in n t intervals of size ∆ t . Algorithm 1 (DSMC method) .
1. Compute the initial sample particles { w i , i = 1 , . . . , N } ,by sampling them from the initial density f ( w ) for n = 0 to n t − given { w ni , i = 1 , . . . , N }◦ set N c = Sround( N ∆ t/ ◦ select N c pairs ( i, j ) uniformly among all possible pairs,- perform the collision between i and j , and compute w (cid:48) i and w (cid:48) j according to the collision law- set w n +1 i = w (cid:48) i , w n +1 j = w (cid:48) j ◦ set w n +1 i = w ni for all the particles that have not been selected end for Here, by Sround( x ) we denote the stochastic rounding of a positive real number x Sround( x ) = (cid:26) (cid:98) x (cid:99) + 1 with probability x − (cid:98) x (cid:99)(cid:98) x (cid:99) with probability 1 − x + (cid:98) x (cid:99) where (cid:98) x (cid:99) denotes the integer part of x .The kinetic distribution, as well as its moments, is then recovered from the empiricaldensity function f N ( t, w ) = 1 N N (cid:88) i =1 δ ( w − w i ( t )) , (14)where δ ( · ) is the the Dirac delta and { w i ( t ) , i = 1 , . . . , N } are the samples at time t ≥
0. Forany test function ϕ , if we denote by( ϕ, f )( t ) = (cid:90) V ϕ ( w ) f ( t, w ) dw, we have ( ϕ, f N )( t ) = 1 N N (cid:88) i =1 ϕ ( w i ( t )) . (cid:82) V f ( t, w ) dw = 1 we have that ( ϕ, f ) = E V [ ϕ ], where E V [ · ] is theexpectation of the observable quantity ϕ with respect to the density f . In the sequel, wewill also make use of the notation E [ · ] to denote the expectation in the random space ofuncertainties. We will implicitly assume that for multidimensional variables, as in the case ofstatistical samples, the expected values are done with respect to each variable.Thanks to the central limit theorem the following result holds [8]: Lemma 1.
The root mean square error is such that for each t ≥ E V (cid:104) (( ϕ, f ) − ( ϕ, f N )) (cid:105) / = σ ϕ N / , where σ ϕ = V ar V [ ϕ ] with V ar V [ ϕ ]( t ) = (cid:90) V ( ϕ ( w ) − ( ϕ, f )( t )) f ( t, w ) dw. If we are interested in the evaluation of the error produced in the approximation of struc-tured quantities like the reconstruction of the distribution function we can operate as follows:we introduce a uniform grid in V ⊆ R where each cell has width ∆ w > S ∆ w ≥ w (cid:90) V S ∆ w ( w ) dw = 1 . Then, we consider the approximation of the empirical density (14) obtained by f N, ∆ w ( t, w ) = 1 N N (cid:88) i =1 S ∆ w ( w − w i ( t )) . (15)In the simplest case, we have S ∆ w ( w ) = χ ( | w | ≤ ∆ w/ / ∆ w , where χ ( · ) is the indicatorfunction, that corresponds to the standard histogram reconstruction.The numerical error of the reconstructed DSMC solution (15), can be estimated from (cid:107) f ( t, · ) − f N, ∆ w ( t, · ) (cid:107) L p ( V,L ( V )) ≤ (cid:107) f ( t, · ) − f ∆ w ( t, · ) (cid:107) L p ( V ) + (cid:107) f ∆ w ( t, · ) − f N, ∆ w ( t, · ) (cid:107) L p ( V,L ( V )) where f ∆ w ( t, w ) = (cid:90) V S ∆ w ( w − w ∗ ) f ( t, w ∗ ) dw ∗ and we defined (cid:107) g (cid:107) L p ( V,L ( V )) = (cid:107) E V (cid:2) g (cid:3) / (cid:107) L p ( V ) . (16)Now, the second term can be evaluated from Lemma 1 taking ϕ ( · ) = S ∆ w ( w − · ), w ∈ V , andgives (cid:107) f ∆ w ( t, · ) − f N, ∆ w ( t, · ) (cid:107) L p ( V,L ( V )) = (cid:107) σ S (cid:107) L p ( V ) N / where σ S ( t, w ) = V ar V [ S ∆ w ( w − · )]( t ) . (17)9inally, the first term is bounded by (cid:107) f ( t, · ) − f ∆ w ( t, · ) (cid:107) L p ( V ) ≤ C f (∆ w ) q , accordingly to the accuracy used in the reconstruction. For example, q = 1 in the case ofsimple histogram reconstruction, since (cid:90) V (cid:12)(cid:12)(cid:12)(cid:12)(cid:12) f ( t, w ) − w (cid:90) w +∆ w/ w − ∆ w/ f ( t, w ∗ ) dw ∗ (cid:12)(cid:12)(cid:12)(cid:12)(cid:12) p dw ≤ C ∆ w p , where C is a bound for the derivative of f ( t, w ). Therefore, in the general case, we have thefollowing result. Theorem 1.
The error introduced by the reconstruction function (15) satisfies (cid:107) f ( t, · ) − f N, ∆ w ( t, · ) (cid:107) L p ( V,L ( V )) ≤ (cid:107) σ S (cid:107) L p ( V ) N / + C f (∆ w ) q , where C f depends on the q derivative in w of f ( t, w ) and σ S is defined in (17) . Remark 1.
It is worth to remark that in the above consistency estimates we neglected theerrors due to time discretization. By adapting the arguments in [1] to the present case it ispossibile to show that Nanbu’s scheme converges in law almost surely to the solution of thetime discretized space homogeneous Boltzmann equation (2) in the deterministic case. Weleave a detailed analysis of the convergence properties in connection with the specific form ofthe uncertain binary interaction (1) to further researches.
In order to analyze the uncertainty of our model (3) one may not be only interested in theexpected value and variance of the distribution function f but rather in any quantity ofinterest (QoI) which can be computed from the distribution function. For that purpose weintroduce the operator q [ f ]( t, w, z ) which is, in the simplest setting the identity, i.e. q [ f ] = f .More generally, the QoI defined by q [ f ] is a functional of f , like, for example, the momentsof the kinetic solution.Let Ψ( z ) be the probability distribution function of the uncertainty z ∈ Ω, hence we definethe expected value of the operator q [ f ]( t, w, z ) by E [ q [ f ]]( t, w ) = (cid:90) Ω q [ f ]( t, w, z )Ψ( z ) dz, whose variance is defined as V ar [ q [ f ]]( t, w ) = (cid:90) Ω ( q [ f ]( t, w, z ) − E [ q [ f ]]( t, w )) Ψ( z ) dz. In socio-economic applications, together with the moments of the distribution f thatare linked with observable quantities, we can define other operators characterizing the mainfeatures of the emerging distribution. The first example is given by the tail distribution1 − F ( t, w, z ) := (cid:90) ∞ w f ( t, v, z ) dv. L ( F ( w, z )) = (cid:82) w vf ∞ ( v, z ) dv (cid:82) ∞ vf ∞ ( v, z ) dv , (18)where F ( w, z ) = (cid:90) w f ∞ ( v, z ) dv, as follows G ( z ) = 1 − (cid:90) L ( x, z ) dx. In the present setting the QoI defined by q [ f ] needs to be approximated by the DSMC solver ofthe Boltzmann equation. Therefore, given uncertain random samples { w ( t, z ) , w ( t, z ) , ..., w N ( t, z ) } ,the corresponding empirical density f N is defined by (14).In particular, the particle approximation of the Lorenz curve at time t reads L ( F ( t, w, z )) ≈ L N ( F N ( t, w, z )) := (cid:80) w i ( t,z ) ≤ w w i ( t, z ) (cid:80) Nk =1 w k ( t, z ) , with F N ( t, w, z ) := card { w i ( t, z ) ≤ w, i = 1 , ..., N } . We shortly recall the standard MC sampling method for uncertainty quantification of a kineticequation of type (2). Assume f ( t, w, z ), w ∈ V , solution of a PDE with uncertainties only inthe initial distribution f ( w, z ), z ∈ Ω ⊆ R d z . The MC sampling method can be regarded asthe simplest method for UQ and can be formulated as follows Algorithm 2 (MC-DSMC algorithm) . The algorithm can be compactly summarized in thefollowing steps.1.
Sampling:
Sample M independent identically distributed (i.i.d.) initial distribution f k, = f ( t, w, z k ) , k = 1 , ..., M from the random initial data f .2. Solving:
For each realization of f k, the underlying kinetic equation (2) is solved nu-merically by a DSMC solver for the kinetic equation. We denote the solution at time t n by f k,nN , k = 1 , ..., M where N is the sample size of the DSMC solver for the kineticequation.3. Reconstruction:
Estimate a statistical moment of the quantity of interest q [ f n ] ofthe random solution field, e.g. for the expected value with the sample mean of theapproximated solution q [ f nN ] E [ q [ f n ]] ≈ E M [ q [ f nN ]] := 1 M M (cid:88) k =1 q [ f k,nN ] . f N ( t, w, z ) = 1 N N (cid:88) i =1 δ ( w − w i ( t, z )) , being { w i ( t, z ) , i = 1 , . . . , N } the samples of the particles at time t ≥ w i ∈ L (Ω).For example, if q [ f ] = ( ϕ, f ) we have q [ f N ]( z, t ) = ( ϕ, f N )( z, t ) = 1 N N (cid:88) i =1 ϕ ( w i ( z, t )) , and the following result holds Lemma 2.
The root mean square error of the MC-DSMC method satisfies E (cid:104) E V [( E [( ϕ, f )] − E M [( ϕ, f N )]) ] (cid:105) / ≤ ν ( ϕ,f ) M / + σ ϕ,M N / , where ν ϕ,f ) = V ar [( ϕ, f )] and σ ϕ,M = E M [ σ ϕ ] with σ ϕ = V ar V [ ϕ ] .Proof. The above estimate follows from E (cid:104) E V [( E [( ϕ, f )] − E M [( ϕ, f N )]) ] (cid:105) / ≤ E (cid:104) ( E [( ϕ, f )] − E M [( ϕ, f )]) (cid:105) / + E V (cid:104) ( E M [( ϕ, f )] − E M [( ϕ, f N )]) (cid:105) / ≤ ν ( ϕ,f ) M / + σ ϕ,M N / , where we used the fact that by Cauchy-Schwartz and Lemma 1 we have E V (cid:104) ( E M [( ϕ, f )] − E M [( ϕ, f N )]) (cid:105) ≤ M M (cid:88) k =1 E V [(( ϕ, f )( z k ) − ( ϕ, f N )( z k )) ]= 1 N (cid:32) M M (cid:88) k =1 σ ϕ ( z k ) (cid:33) . Let us consider the reconstructed distribution with uncertainty f N, ∆ w ( t, w, z ) = 1 N N (cid:88) i =1 S ∆ w ( w − w i ( t, z )) , (19)and let us focus on the accuracy of the expectation of the solution E [ f ]. We can estimate (cid:107) E [ f ]( t, · ) − E M [ f N, ∆ w ]( t, · ) (cid:107) L p ( V,L (Ω ,L ( V ))) ≤ (cid:107) E [ f ]( t, · ) − E [ f ∆ w ]( t, · ) (cid:107) L p ( V ) + (cid:107) E [ f ∆ w ]( t, · ) − E M [ f N, ∆ w ]( t, · ) (cid:107) L p ( V,L (Ω)) , (cid:107) g (cid:107) L p ( V,L (Ω ,L ( V ))) = (cid:107) E [ E V [ g ]] / (cid:107) L p ( V ) . (20)We can bound the first term as in the case without uncertainty (cid:107) E [ f ]( t, · ) − E [ f ∆ w ]( t, · ) (cid:107) L p ( V ) ≤ C E [ f ] (∆ w ) q whereas the second term is bounded using Lemma 2 with ϕ ( · ) = S ∆ w ( w − · ) to get (cid:107) E [ f ∆ w ]( t, · ) − E M [ f N, ∆ w ]( t, · ) (cid:107) L p ( V,L (Ω ,L ( V )) ≤ (cid:107) ν ( S,f ) (cid:107) L p ( V ) M / + (cid:107) σ S,M (cid:107) L p ( V ) N / , where ν S,f ) = V ar [( S ∆ w ( w − · ) , f )] . (21)As a consequence we have shown that Theorem 2.
The error introduced by the reconstruction function (19) in the MC-DSMCmethod satisfies (cid:107) E [ f ]( t, · ) − E M [ f N, ∆ w ]( t, · ) (cid:107) L p ( V,L (Ω ,L ( V ))) ≤ (cid:107) ν ( S,f ) (cid:107) L p ( V ) M / + (cid:107) σ S,M (cid:107) L p ( V ) N / + C E [ f ] (∆ w ) q where ν S,f ) is defined in (21) and σ S,M = E M [ σ S ] with σ S defined in (17) . Remark 2.
In the above consistency estimates we used the error norm defined in (20) . Afrequently used norm in uncertainty quantification based on Monte Carlo strategies is theexpectation of the error defined as (cid:107) g (cid:107) L (Ω ,L p ( V,L ( V ))) = E (cid:104) (cid:107) E V [ g ] / (cid:107) L p ( V ) (cid:105) / . (22) The two norm for p (cid:54) = 2 differs and are related by Jensen’s inequality (cid:107) g (cid:107) L p ( V,L (Ω ,L ( V ))) ≤ (cid:107) g (cid:107) L (Ω ,L p ( V,L ( V ))) . We refer to [15] for a more detailed discussion.
In order to improve the accuracy of standard MC sampling methods, we introduce a novelclass of mean field control variate methods [15]. The key idea is to take advantage of thereduced cost of the mean field model which approximates the asymptotic behavior of theoriginal Boltzmann model. This enables us to reduce the variance of the MC estimate.More precisely we consider two control variates approaches obtained by the the mean fieldapproximation: the direct numerical solution of the mean field model and its correspondingsteady state.The consistency of the mean field approximation (6) with the integro-differential equation(5) has been previously discussed in Section 2. We denote the solution of the time scaledBoltzmann equation (5) by f ε = f ε ( t, w, z ) and the solution of the corresponding mean field13pproximation (6) by ˜ f = ˜ f ( t, w, z ). The grazing collision regime of the time scaled kineticequation (5) corresponds to a small scaling parameter ε . Thus, up to an extraction of asubsequence, we indicate lim ε → f ε ( t, w, z ) = ˜ f ( t, w, z ) , being f ε solution of (5). Therefore, also the equilibrium distributions are compatible, i.e.lim t →∞ lim ε → f ε ( t, w, z ) = ˜ f ∞ ( w, z ) . As we observed, the steady state ˜ f ∞ ( w, z ) of the Fokker-Planck model is analytically com-putable in several cases, whereas the steady state of the Boltzmann model (5) is unknown.The parameter dependent control variate method with λ ∈ R can be formulated by intro-ducing the quantity q λ [ f ε ] = q [ f ε ] − λ ( q [ ˜ f ] − E [ q [ ˜ f ]]) . (23)It is straightforward to observe that the expected value satisfies E [ q λ [ f ε ]] = E [ q [ ˜ f ]] . We can state the following
Theorem 3.
The optimal value λ ∗ which minimizes the variance of (23) is given by λ ∗ := C ov [ q [ f ε ] , q [ ˜ f ]] V ar [ q [ ˜ f ]] , (24) where C ov [ · , · ] denotes the covariance. The corresponding variance of q λ ∗ [ f ε ] is then V ar [ q λ ∗ [ f ε ]] = (cid:16) − ρ q [ f ε ] ,q [ ˜ f ] (cid:17) V ar [ q [ f ε ]] , (25) where ρ q [ f ε ] ,q [ ˜ f ] := C ov [ q [ f ε ] , q [ ˜ f ]] (cid:113) V ar [ q [ f ε ]] V ar [ q [ ˜ f ]] ∈ ( − , , is the Pearson’s correlation coefficient between q [ f ε ] and q [ ˜ f ] . In particular, we have lim ε → C ov [ q [ f ε ] , q [ ˜ f ]] V ar [ q [ ˜ f ]] = 1 , lim ε → V ar [ q λ ∗ [ f ε ]] = 0 . Proof.
We aim to choose the coefficient λ such that the variance V ar [ q λ [ f ε ]] of the controlvariate formulation is minimized. We have V ar [ q λ [ f ε ]] = V ar [ q [ f ε ]] − λ C ov [ q [ f ε ] , q [ ˜ f ]] + λ V ar [ q [ ˜ f ]] . (26)Then we can differentiate (26) with respect to λ to obtain that (24) minimizes the variance.The second part of the theorem follows immediately since lim ε → f ε = ˜ f holds in the quasi-invariant regime for f ε solution to (5). 14quation (25) reveals that the mean field control variate approach may lead to a strongvariance reduction provided that the correlation coefficient is close to one.In the simplest setting the control variate formulation in (23) can be modified using thesteady state ˜ f ∞ ( w, z ) of the mean-field model (6) q λ [ f ε ] = q [ f ε ] − λ ( q [ ˜ f ∞ ] − E [ q [ ˜ f ∞ ]]) . (27)For the mean field control variate steady state (27) similar results holds in the large timelimit lim t →∞ lim ε → C ov [ q [ f ] , q [ ˜ f ∞ ]] V ar [ q [ ˜ f ∞ ]] = 1 , lim t →∞ lim ε → V ar [ q λ ∗ [ f ]] = 0 , where now λ ∗ := C ov [ q [ f ε ] , q [ ˜ f ∞ ]] V ar [ q [ ˜ f ∞ ]] . In practice it is only possible to compute the optimal λ ∗ numerically. Furthermore, it is ofparamount importance to be able to compute E [ q [ ˜ f ]] or E [ q [ ˜ f ∞ ]] exactly or with very smallerror in order to keep advantage of the control variate approach. In a MC setting, we use M realizations of our random variable z to define the Mean FieldControl Variate (MFCV) estimator E [ q λ ∗ [ f ε ]] ≈ E M [ q λ ∗ [ f ε ]] = E M [ q [ f ε ]] − λ ∗ N ( E M [ q [ ˜ f ]] − E [ q [ ˜ f ]]) , (28)where λ ∗ N = Cov M [ q [ f ε ] , q [ ˜ f ]] V ar M [ q [ ˜ f ]] , (29)and E [ q [ ˜ f ]] denotes the exact value of the expectation of the considered QoI or its numericalapproximation with negligible error. Furthermore, we introduced the following notations E M [ q [ f ε ]] := 1 M M (cid:88) k =1 q [ f kε ] , E M [ q [ ˜ f ]] := 1 M M (cid:88) k =1 q [ ˜ f k ] V ar M [ q [ ˜ f ]] := 1 M − M (cid:88) k =1 ( q [ ˜ f k ] − E [ q [ ˜ f ]]) ,Cov M [ q [ f ε ] , q [ ˜ f ]] := 1 M − M (cid:88) k =1 ( q [ f kε ] − E M [ q [ f ε ]]) ( q [ ˜ f k ] − E [ q [ ˜ f ]]) , with f kε and ˜ f k solutions of the Boltzmann-type model and of the mean-field model, respec-tively, relative to k th realization of the random variable z .The MFCV algorithm based on a DSMC algorithm reads: Algorithm 3 (MFCV-DSMC algorithm) . The main steps of the MFCV-DSMC method foruncertainty in the initial data can be summarized as follows:1.
Sampling:
Sample M independent identically distributed (i.i.d.) samples of the initialdistribution f k, = f ( w, z k ) , k = 1 , ..., M from the random initial data f ( w, z ) . . Solving:
For each realization f ,k , k = 1 , ..., M • Compute the control variate ˜ f k,n , k = 1 , ..., M at time t n solving with a suitabledeterministic method the mean field model (6) (or using the steady state ˜ f ∞ ) andcompute E [ q [ ˜ f n ]] (or E [ q [ ˜ f ∞ ]] ) with negligible error. • Solve the kinetic equation (5) by a DSMC solver with sample size N . We denotethe solution at time t n by f k,nε,N , k = 1 , ..., M .3. Estimating: • Estimate the optimal value of λ ∗ at time t n by (29) and denote it as λ ∗ ,nM . • Compute the the expectation of any quantity of interest q [ f ε,N ] of the random so-lution field with the mean-field control estimator E λ ∗ M [ g [ f ε,N ]] = E M [ q [ f nε,N ]] − λ ∗ ,nM ( E M [ q [ ˜ f n ]] − E [ q [ ˜ f n ]]) . The error bound of the MFCV-DSMC method may improve significantly in comparisonto the standard MC method thanks to the relation (25) once ρ q [ f ε ] ,q [ ˜ f ] ≈ λ ∗ , wehave the following Lemma 3.
The root mean square error of the MFCV-DSMC method satisfies E (cid:20) E V (cid:20)(cid:16) E [( ϕ, f ε )] − E λ ∗ M [( ϕ, f ε,N )] (cid:17) (cid:21)(cid:21) / ≤ (cid:16) − ρ ϕ,f ε ) , ( ϕ, ˜ f ) (cid:17) / ν ( ϕ,f ε ) M / + σ ϕ,M N / , where ν ϕ,f ε ) = V ar [( ϕ, f ε )] and σ ϕ,M = E M [ σ ϕ ] with σ ϕ = V ar V [ ϕ ] .Proof. Let us observe that E [( ϕ, f ε )] = E [( ϕ, f λ ∗ ε )] , E λ ∗ M [( ϕ, f ε, N ) = E M [( ϕ, f λ ∗ ε,N )]and then, using Theorem 3, we get ν ϕ,f λ ∗ ε ) = V ar [( ϕ, f λ ∗ ε )] = (cid:16) − ρ ϕ,f ε ) , ( ϕ, ˜ f ) (cid:17) V ar [( ϕ, f ε )] . The conclusion follows from Lemma 2 together with the identity( ϕ, f λ ∗ ε ) − ( ϕ, f λ ∗ ε,N ) = ( ϕ, f ε ) − ( ϕ, f ε,N ) . In the case of the reconstruction function (19), we have (cid:13)(cid:13)(cid:13) E [ f ε ]( t, · ) − E λ ∗ M [ f ε,N, ∆ w ]( t, · ) (cid:13)(cid:13)(cid:13) L p ( V,L (Ω ,L ( V ))) ≤ (cid:107) E [ f ε ]( t, · ) − E [ f ε, ∆ w ]( t, · ) (cid:107) L p ( V ) + (cid:13)(cid:13)(cid:13) E [ f ε, ∆ w ]( t, · ) − E λ ∗ M [ f ε,N, ∆ w ]( t, · ) (cid:13)(cid:13)(cid:13) L p ( V,L (Ω ,L ( V ))) , where the first term is bounded as in Theorem 2 and the second term can be bounded usingLemma 3 with ϕ ( · ) = S ∆ w ( w − · ). Thus we proved the following result.16 heorem 4. The error introduced by the reconstruction function (19) in the MFCV-DSMCmethod satisfies (cid:13)(cid:13)(cid:13) E [ f ε ]( t, · ) − E λ ∗ M [ f ε,N, ∆ w ]( t, · ) (cid:13)(cid:13)(cid:13) L p ( V,L (Ω ,L ( V ))) ≤ (cid:13)(cid:13)(cid:13)(cid:13)(cid:16) − ρ S,f ε ) , ( S, ˜ f ) (cid:17) / ν ( S,f ε ) (cid:13)(cid:13)(cid:13)(cid:13) L p ( V ) M / + (cid:107) σ S,M (cid:107) L p ( V ) N / + C E [ f ] (∆ w ) q (30) where and ν S,f ε ) is defined as in (21) and σ S,M is defined in Theorem 2.
As a consequence when the solution of the full model is close to the solution of the controlvariate the statistical error due to the uncertainty vanishes. This justifies the use of a largernumber of samples in the state space in agreement with the reconstruction used in order tobalance the last two error terms in (30).
In this section, we present several numerical examples of the mean field control variate methodwith application to different kinetic models in socio-economic sciences. The presented testcases consider uncertain initial data and uncertain interaction coefficients as well. We startpresenting relevant tests on the steady state control variate approach that, we remark, hasno impact on the simulation costs, since the expected value of steady states can be computedoffline with arbitrary accuracy. Subsequently, we consider the general mean field control vari-ate approach where the mean-field model is solved using a second order structure preservingfinite difference method [39]. In this case, we cannot perform the computation of E [ q [ · ]] of-fline. Hence, the computational costs play a relevant role and we will discuss the connectionbetween simulation costs and performance. We consider in this section the case where the steady state ˜ f ∞ ( w, z ) of the mean-field modelis analytically known and is used as control variate. We will refer to this method, where thecontrol variate term is evaluated off line, as Mean Field Steady State Control Variate (MFCV-S). If not otherwise indicated the quantity of interest is here the distribution itself q [ f ε ] = f ε .The reconstruction is performed using standard first order histogram approximation with N Z = 100 grid points. In order to compute with negligible error E [ ˜ f ∞ ] we adopt a stochasticcollocation approach with 20 collocation nodes. Test 1: Opinion model with uncertainty
We consider the kinetic model for opinion formation resulting from binary interactions (9)with uncertainties present on the initial distribution or on the interaction strength. Let usfirst assume that the uncertainty is present in the initial distribution of the problem. Moreprecisely, we assume that z ∼ U ([0 , f ( w, z ) is given by f ( w, z ) = w ∈ (cid:20)
14 ( z − ,
14 ( z + 2) (cid:21) . (31)17 -4 -3 -2 Figure 1:
Test 1 . Convergence of the expected value of the solution E [ f ε ], in the case ofthe Boltzmann model for opinion formation (A) with uncertain initial data, for the classicalMC sampling method at time t = 2. The solutions are averaged over 50 runs to reducestatistical fluctuations. The DSMC solver has been implemented with N = 2 × samplesand ε = 10 − .Hence, the first moment m ( z ) = z/ p ( | v − w | , z ) = 1 , D ( v ) = (cid:112) − v , (A)and thus we obtain a steady state of the Fokker-Planck equation of the form ˜ f ∞ ( w, z ) givenby (11). The Boltzmann equation is solved with N = 2 × particles up to t = 5 where thesolution approaches the steady state.In Figure 1, we report the L error of E [ f ε ] at t = 5 for different number of samples M from which we can clearly deduce the MC convergence rate of . The convergence towards E [ ˜ f ∞ ] is given in the top row of Figure 2 for different ε = 10 − , − and different number ofsamples M = 20 (left plot) and M = 1280 (right plot). The bottom row of Figure 2 showsthat using the MFCV-S approach, in comparison to the classical MC method, for M = 20 weobtain already a good fit to the mean field steady state for ε = 10 − .In Figure 3 (left) we report the L error of the expected density obtained by the MC andMFCV-S methods for increasing number of samples at time t = 5. We obtain an improvementin accuracy for the MFCV-S between one and two orders of magnitude using the same numberof samples.As a further test for opinion dynamics we consider the kinetic model resulting from binaryinteractions (9) with p ( | w − v | , z ) = 34 + z , z ∼ U ([ − , , D ( w ) = 1 − w , (B)so that the resulting steady state of the Fokker-Planck model is the Maxwellian-like distribu-tion (10). The initial data in this case is (31) in the deterministic setting z = 0 and the finaltime is t = 20.In Figure 3 (right) we report the L error of expected probability distribution functioncomputed by the MC and MFCV-S method at the final time for different number of samples.18 Figure 2:
Test 1 . Approximation of E [ f ε ], in the case of the Boltzmann model for opinionformation (A) with uncertain initial data. The DSMC solver has been implemented with N = 2 × . We considered different scalings ε = 10 − , − and two different number ofsamples M = 20 (left) and M = 1280 (right). The cross markers represent E [ ˜ f ∞ ]. (Top)standard MC sampling; (Bottom) MFCV-S method. -4 -3 -2 -1 -6 -5 -4 -3 -2 Figure 3:
Test 1 . Error of the MFCV-S estimate and classical MC method in the case of theBoltzmann model for opinion formation for increasing samples M . We considered N = 2 × in the DSMC solver. (Left) Solution at t = 5 of case (A) with uncertain initial data; (Right)Solution at t = 20 of case (B) with uncertain interaction parameters.19 Figure 4:
Test 2 . Expected distribution E [ f ε ] and expected Lorentz curve for the MC,MFCV-S and the mean field model for problem (A). The DSMC solver of the Boltzmannmodel used N = 2 × particles. Left: Expected distribution function E [ f ] computed with M = 10 samples for the MC and MFCV-S method. Right: Expected Lorenz curve computedwith M = 20 samples for the MC and MFCV-S method.As in case (A) we obtain an improvement between one and two orders of accuracy for theMFCV-S method compared to the classical MC method. Test 2: Wealth model with uncertainty
We study two test cases of the CPT model defined by the binary interaction scheme (12).First, we consider uncertainty in the initial condition and secondly in the saving propensity.We will consider as computational domain the interval [0 , z ∈ U ([0 , f ( w, z ) defined by f ( w, z ) = w ∈ (cid:104) z , z (cid:105) . Furthermore, we consider λ ( z ) = 1 , D ( w ) = w, (A)so that the large time behavior of the Fokker-Planck model ˜ f ∞ ( w, z ) is given by (13) with m ( z ) = 1 + z E [ f ε ] and of the expectedvalue of the Lorenz curve, see Section 3.2, at the final time t = 30 for ε = 10 − . We canobserve a better agreement of the solution of the MFCV-S method with the expected steadystate solution of the mean field model than the solution obtained by the MC solver.In Figure 5 we compare the expected QoI computed by the MC and MVCV-S solver fordifferent number of samples M plotted against time. We can observe that the MC methodneeds at least 8 times more samples in order to reach the same accuracy than the MFCV-Smethod with 10 samples. 20 -8 -7 -6 -5 -4 -3 -2 -6 -5 -4 -3 Figure 5:
Test 2 . The L error for the expected QoI computed for problem (A) by the MCand MFCV-S method against time. On the left we show the expected probability densityfunction and Gini coefficient, whereas on the right the expected tail index and Lorenz curve.The DSMC solver has been implemented with N = 2 × particles.Next, we consider uncertainty in the interaction λ ( z ) = 12 + z , z ∼ U ([ − , , D ( w ) = w. (B)The initial condition is uniformly distributed on [0 , m ˜ f ≡ L error for E [ f ] at final time t = 30 computed by the MC andMFCV-S method for M = 10. We deduce that similar to the previous test cases the MFCV-Smethod improves the accuracy in comparison to the MC method between one and two ordersof magnitude. The MFCV method relies on the numerical solution of the Fokker-Planck equation (6) whichis able to capture the transient behavior much better than the MFCV-S method. In fact theFokker-Planck equation (6) is computationally less expensive than the original kinetic modelsince the expensive integral operator is replaced by the simpler differential Fokker-Planckoperator. In order to solve the mean field model we adopted the second order structurepreserving method for nonlocal Fokker-Planck equations developed in [39]. Indeed, for ourcontrol variate method it is of immanent importance that the scheme is able to capture theasymptotic behavior of the model with arbitrary accuracy.We denote the numerical approximation of the mean field model at time t n by ˜ f n ∆ w andthus the control variate estimate for the QoI q [ f nε ] = q [ f ε ( t n )] becomes E λ ∗ M [ q N [ f nε ]] = E M [ q N [ f nε ]] − λ ∗ ,nM ( E M [ q [ ˜ f n ∆ w ]] − E [ q [ ˜ f n ∆ w ]]) . (32)For the structure preserving method we select a semi implicit time integrator and a Gaussquadrature rule to evaluate the fluxes as defined in [39]. The boundary conditions are set tozero flux boundary conditions. As pointed out previously, the computational cost to evaluate(32) is not negligible due to the time dependence of the control variate. In fact, in order21 -5 -4 -3 -2 Figure 6:
Test 2 . L error for the expected distribution E [ f ε ] for problem (B) computed bythe MC and MFCV-S method for different number of samples M at t = 30. We considered N = 2 × in the DSMC solver.to have an accurate estimate of E [ q [ ˜ f n ∆ w ]] the Fokker-Planck equation needs to be solvedfor a sufficiently large number of random samples M MF at each time step and estimate E [ q [ ˜ f n ∆ w ]] ≈ E MF [ q [ ˜ f n ∆ w ]]. In the sequel, we discuss the issue of computational costs in moredetails. Estimating the computational cost
The costs of the solution of the Boltzmann equation at one time step can be estimated bythe number of samples in physical space N times the number of samples in random space M .We have Cost( E M [ q N [ f ε ]]) = N M.
For the mean field model with N MF grid points and M MF samples in random space we get:Cost( E M MF [ q N MF [ ˜ f ]]) = N MF M MF . Additionally, the total computation costs depend on the selected time step. The time stepof the Boltzmann model is allowed to take value ∆ t ∈ (0 , ε ] in the DSMC solver. As in allour simulations we choose the maximum allowed value ∆ t = ε . Consequently, since we use asemi-implicit method for the mean-field solver with time step ∆ t MF we have ε = ∆ t = ∆ t MF k , k ≥ . We aim to control the cost of the mean field model by the total cost of the Boltzmann model,namely Cost( E M [ q N [ f ]]) ≥ Cost( E M MF [ q N MF [ ˜ f ]] . Thus, we have to ensure that
N M ≥ N MF M MF k -6 -5 -4 -5 -4 -3 Figure 7:
Test 3 . The L error of E [ f ] computed by the MFCV and MFCV-S method fordifferent number of samples M at t = 0 .
1. We consider the same opinion formation model ofTest 1. The left hand side corresponds to the setting of case (A) and the right to case (B).is satisfied. Provided that the number of points and particles in physical space are fixed, weobtain the following upper bound on M MF M MF ≤ k N MN MF . (33)In practice, for a given ε , N and M , we define the number of grid points N MF in thedeterministic mean-field solver which provides the largest value of the time step ∆ t MF toensure stability. The maximum allowed number of samples in the mean-field model is thenchosen accordingly to (33). Test 3: Opinion model with uncertainty
In order to show the performance of the MFCV method we select the same setting as definedin Section 5.1. The number of grid points of our mean field model is set to N MF = 20.Furthermore, the number of particles of the Boltzmann model is fixed to N = 2 × and wefix k = 1. Hence, for M MF and arbitrary M we get accordingly to (33) the estimate M MF ≤ × M
20 = 10 M. (34)For our simulations we have chosen M MF = 10 which is accordingly to (34) exactly theupper bound for M = 10. Thus, the cost of the MFCV method are comparable to the costof the corresponding MC method.In Figure 7 we compare the MFCV method with the MFCV-S for different numbers ofsamples M . The QoI is E [ f ε ]. As expected we gain accuracy and the error of the MFCVmethod is below the error of the MFCV-S method in transient times, we considered here t = 0 .
1. 23 -4 -3 -6 -5 -4 -3 -2 -1 Figure 8:
Test 4 . The L error of E [ f ] computed by the MFCV and MFCV-S method fordifferent number of samples M at t = 1. We consider the same wealth model of Test 2. Theleft hand side corresponds to the setting of case (A) and the right to case (B). Test 4: Wealth model with uncertainty
For the CPT model we consider N MF = 100 and we choose N = 5 × , M MF = 5 × and k = 1. Thus the upper bound (33) is satisfied for all samples M ≥ L error of E [ f ε ] computed by the MFCV or MFCV-S methodat fixed times but for different number of samples M . We deduce that the MFCV method isable to be more accurate than the MFCV-S method also in short times. Test 5: Bounded confidence opinion model with uncertainty
In this last test case we consider the opinion model introduced in section 2.2 with boundedconfidence interaction rule [25,38]. The uncertainty is present in the function P which weightsthe compromise tendency. We consider P ( | w − v | ) = χ ( | w − v | < z ) , z ∼ U ([1 , , D ( w ) = 1 − w . (35)In comparison to the previous test cases the steady sate distribution of the mean field modelis unknown. Therefore, we cannot apply the MFCV-S method but only the MFCV method.The number of grid points of our mean field model is set to N MF = 40, the number of particlesof the Boltzmann solver is fixed to N = 2 × and the number of samples for the mean fieldmodel is given by M = 5 × . The simulations have been conducted with ε = 2 . × − and the initial distribution is given by f (0 , w ) = C exp (cid:40) − (cid:18) w + 12 (cid:19) (cid:41) + C exp (cid:40) − (cid:18) w − (cid:19) (cid:17)(cid:41) , with normalization constant C > L error of the E [ f ] plotted for increasing number of samples M (left plot) and for fixed number of samples in the time interval [0 , -6 -5 -4 -3 -2 -1 -6 -5 -4 -3 -2 -1 Figure 9:
Test 5. L error for E [ f ] using the bounded confidence model of opinion, fordifferent number of samples (left) and evolution in time (right). The simulation has beenconducted with N = 2 × particles for the Boltzmann solver and the scaling ε = 2 . × − . We introduced a novel variance reduction technique for uncertainty quantification in Boltz-mann type equations of interest in socio-economic applications. The method relies on acontrol variate strategy in order to speed up the convergence of the standard Monte Carlomethod sampling method in the random space. In contrast to the classical Boltzmann equa-tion of rarefied gas dynamics we cannot use as in [15] the knowledge of the local equilibriumstate to design control variate methods. For this purpose, we considered the mean field ap-proximation of the Boltzmann model as surrogate model for variance reduction. Therefore,the mean field control variate (MFCV) method makes use of the less expensive solution ofthe mean field approximation to accelerate the Monte Carlo convergence. In the physicalspace, the mean field model has been computed with a deterministic structure preservingmethod whereas the kinetic Boltzmann model has been solved by a standard DSMC method.Thus, the novel MFCV method can be regarded as a hybrid method in the physical space.Whenever the steady state of the mean field model is known it is possible to consider a sim-plified control variate strategy that uses the steady state as surrogate model. This latterapproach, even if less accurate in transient regimes, leads to a strong computational costreduction since the control variate can be evaluated off line. The numerical results confirmthat the MFCV methods outperform the standard MC method in all applications considered.Although we have focused in this work on space-homogeneous Boltzmann type equations forsocio-economic applications it is possible to apply the new MFCV method to a larger classof kinetic equations. Besides applications to other Boltzmann equations, several extensionsare actually under study. Among these, the case of space non-homogeneous Boltzmann typeequations and the introduction of multiple control variates as in [16].
Acknowledgments
This research is funded by the Deutsche Forschungsgemeinschaft (DFG, German ResearchFoundation) under Germany’s Excellence Strategy – EXC-2023 Internet of Production –390621612 and supported also by DFG HE5386/15. T. T. acknowledges the support by25he ERS Prep Fund - Simulation and Data Science. The work was partially funded by theExcellence Initiative of the German federal and state governments. L. P. would like to thankthe Italian Ministry of Instruction, University and Research (MIUR) to support this researchwith PRIN Project 2017, No. 2017KKJP4X - Innovative numerical methods for evolutionarypartial differential equations and applications. This paper was written within the activities ofthe GNFM and GNCS of INDAM. The research was partially supported by the Italian Min-istry of University and Research (MUR): Dipartimenti di Eccellenza Program (2018–2022) -Dept. of Mathematics ”F. Casorati”, University of Pavia.
References [1] H. Babovsky. A convergence proof for Nanbu’s Boltzmann simulation scheme.
Eur. J.Mech. B Fluids , 8(1):550–591, 1989.[2] H. Babovsky and R. Illner. A convergence proof for Nanbu’s simulation method for thefull Boltzmann equation.
SIAM J. Numer. Anal. , 26(1):45–65, 1989.[3] M. Ballerini, N. Cabibbo, R. Candelier, A. Cavagna, E. Cisbani, I. Giardina, A. Orlandi,G. Parisi, A. Procaccini, M. Viale, and Zdravkovic. Empirical investigation of starlingflocks: a benchmark study in collective animal behaviour.
Anim. Behav. , 76:201–215,2008.[4] N. Bellomo, P. Degond, and E. Tadmor, editors.
Active Particles, Volume 1 . Modelingand Simulation in Science, Engineering and Technology. Birkh¨auser Basel, 2017.[5] N. Bellomo, P. Degond, and E. Tadmor, editors.
Active Particles, Volume 2 . Modelingand Simulation in Science, Engineering and Technology. Birkh¨auser Basel, 2019.[6] G. Bird. Direct simulation and the Boltzmann equation.
Phys. Fluids , 13(11):2676–2681,1970.[7] M. Bongini, M. Fornasier, M. Hansen, and Maggioni. Inferring interaction rules fromobservations of evolutive systems I: The variational approach.
Math. Mod. Meth. Appl.Sci. , 27(5):909–951, 2017.[8] R. E. Caflisch. Monte Carlo and quasi-Monte Carlo methods.
Acta Numer. , 7:1–49,1998.[9] J. A. Carrillo, M. Fornasier, G. Toscani, and F. Vecil. Particle, kinetic, and hydrodynamicmodels of swarming. In
Mathematical Modeling of Collective Behavior in Socio-Economicand Life Sciences , pages 297–336. Springer, 2010.[10] J. A. Carrillo and G. Toscani. Exponential convergence toward equilibrium for homoge-neous Fokker–Planck-type equations.
Math. Meth. Appl. Sci. , 21(13):1269–1286, 1998.[11] C. Cercignani.
Rarefied Gas Dynamics: From Basic Concepts to Actual Calculations ,volume 21. Cambridge University Press, 2000.[12] S. Cordier, L. Pareschi, and C. Piatecki. Mesoscopic modelling of financial markets.
J.Stat. Phys. , 134(1):161–184, 2009. 2613] S. Cordier, L. Pareschi, and G. Toscani. On a kinetic model for a simple market economy.
J. Stat. Phys. , 120(1-2):253–277, 2005.[14] G. Dimarco and L. Pareschi. Numerical methods for kinetic equations.
Acta Numer. ,23:369–520, 2014.[15] G. Dimarco and L. Pareschi. Multi-scale control variate methods for uncertainty quan-tification in kinetic equations.
J. Comput. Phys. , 388:63–89, 2019.[16] G. Dimarco and L. Pareschi. Multiscale variance reduction methods based on multi-ple control variates for kinetic equations with uncertainties.
Multiscale Model. Simul. ,18(1):351–382, 2020.[17] G. Dimarco, L. Pareschi, and M. Zanella. Uncertainty quantification for kinetic modelsin socio–economic and life sciences. In
Uncertainty Quantification for Hyperbolic andKinetic Equations , pages 151–191. Springer, 2017.[18] B. D¨uring, L. Pareschi, and G. Toscani. Kinetic models for optimal control of wealthinequalities.
Eur. Phys. J. B , 91(10):265, 2018.[19] B. D¨uring and G. Toscani. Hydrodynamics from kinetic models of conservativeeconomies.
Physica A , 384(2):493–506, 2007.[20] G. Furioli, A. Pulvirenti, E. Terraneo, and G. Toscani. The grazing collision limit of theinelastic Kac model around a L´evy-type equilibrium.
SIAM J. Math. Anal. , 44(2):827–850, 2012.[21] G. Furioli, A. Pulvirenti, E. Terraneo, and G. Toscani. Fokker-planck equations in themodeling of socio-economic phenomena.
Math. Mod. Meth. Appl. Sci. , 27(01):115–158,2017.[22] I. Gamba, S. Jin, and L. Liu. Error estimate of a bi-fidelity method for kinetic equationswith random parameters and multiple scales.
Int. J. Uncertain. Quan. , page to appear,2020.[23] J. L. Gastwirth. The estimation of the Lorenz curve and Gini index.
Rev. Econ. Stat. ,pages 306–316, 1972.[24] M. B. Giles. Multilevel Monte Carlo methods.
Acta Numer. , 24:259–328, 2015.[25] R. Hegselmann and U. Krause. Opinion dynamics and bounded confidence: models,analysis and simulation.
J. Artif. Soc. Soc. Simul. , 5(3), 2002.[26] M. Herty, A. Tosin, G. Visconti, and M. Zanella. Reconstruction of traffic speed dis-tributions from kinetic models with uncertainties. In A. Tosin and G. Puppo, editors,
Mathematical Descriptions of Traffic Flow: Micro, Macro and Kinetic Models , SEMA-SIMAI Springer Series. Springer, 2021.[27] J. Hu and S. Jin. A stochastic Galerkin method for the Boltzmann equation with uncer-tainty.
J. Comput. Phys. , 315:150–168, 2016.2728] J. Hu, L. Pareschi, and W. Yubo. Uncertainty quantification for the BGK model of theBoltzmann equation using multilevel variance reduced Monte Carlo methods. PreprintarXiv:2004.07638, 2020.[29] S. Jin and L. Pareschi, editors.
Uncertainty Quantification for Hyperbolic and KineticEquations , volume 14 of
SEMA SIMAI Springer Series . Springer International Publish-ing, 2017.[30] L. Liu and X. Zhu. A bi-fidelity method for the multiscale Boltzmann equation withrandom parameters.
J. Comput. Phys. , 402:108914, 2020.[31] S. Mishra, C. Schwab, and J. ˇSukys. Multi-level Monte Carlo finite volume meth-ods for nonlinear systems of conservation laws in multi-dimensions.
J. Comput. Phys. ,231(8):3365–3388, 2012.[32] G. Naldi, L. Pareschi, and G. Toscani, editors.
Mathematical modeling of collective behav-ior in socio-economic and life sciences . Modeling and Simulation in Science, Engineeringand Technology. Birkh¨auser Basel, 2010.[33] K. Nanbu. Direct simulation scheme derived from the Boltzmann equation. I. Monocom-ponent gases.
J. Phys. Soc. Jpn. , 49(5):2042–2049, 1980.[34] L. Pareschi. An introduction to uncertainty quantification for kinetic equations andrelated problems. In G. Albi, S. Merino-Aceituno, A. Nota, and M. Zanella, editors,
Trails in Kinetic Theory: Foundational Aspects and Numerical Methods , SEMA-SIMAISpringer Series. Springer, 2021.[35] L. Pareschi and G. Russo. An introduction to Monte Carlo method for the Boltzmannequation. In
ESAIM: Proceedings , volume 10, pages 35–75. EDP Sciences, 2001.[36] L. Pareschi and G. Toscani. Self-similarity and power-like tails in nonconservative kineticmodels.
J. Stat. Phys. , 124(2):747–779, 2006.[37] L. Pareschi and G. Toscani.
Interacting Multiagent Systems: Kinetic Equations andMonte Carlo Methods . Oxford University Press, 2013.[38] L. Pareschi, G. Toscani, A. Tosin, and M. Zanella. Hydrodynamic models of preferenceformation in multi-agent societies.
J. Nonlin. Sci. , 29(6):2761–2796, 2019.[39] L. Pareschi and M. Zanella. Structure preserving schemes for nonlinear Fokker–Planckequations and applications.
J. Sci. Comput. , 74(3):1575–1600, 2018.[40] L. Pareschi and M. Zanella. Monte Carlo stochastic Galerkin methods for the Boltzmannequation with uncertainties: Space-homogeneous case.
J. Comput. Phys. , 423:109822,2020.[41] B. Peherstorfer, K. Willcox, and M. Gunzburger. Survery of multifidelity methods inuncertainty propagation, inference, and optimization.
SIAM Rev. , 60(3):550–591, 2018.[42] M. P. Pettersson, G. Iaccarino, and J. Nordstrom.
Polynomial Chaos Methods for Hy-perbolic Partial Differential Equations , volume 10. Springer, 2015.2843] G. Po¨ette, B. Despr´es, and D. Lucor. Uncertainty quantification for systems of conser-vation laws.
J. Comput. Phys. , 228(7):2443–2467, 2009.[44] G. Pomraning. The Fokker-Planck operator as an asymptotic limit.
Math. Mod. Meth.Appl. Sci. , 2(01):21–36, 1992.[45] G. Toscani. Entropy production and the rate of convergence to equilibrium for theFokker-Planck equation.
Quarterly of Applied Mathematics , 57(3):521–541, 1999.[46] G. Toscani. Kinetic models of opinion formation.
Commun. Math. Sci. , 4(3):481–496,2006.[47] A. Tosin and M. Zanella. Boltzmann-type models with uncertain binary interactions.
Commun. Math. Sci. , 16(4):963–985, 2018.[48] T. Trimborn, L. Pareschi, and M. Frank. Portfolio optimization and model predictivecontrol: A kinetic approach.
Discrete & Continuous Dynamical Systems-B , 24(11):6209–6238, 2019.[49] C. Villani. On a new class of weak solutions to the spatially homogeneous Boltzmannand Laundau equations.
Arch. Ration. Mech. Anal. , 143(3):273–307, 1998.[50] W. Wagner. A convergence proof for Bird’s direct simulation Monte Carlo method forthe Boltzmann equation.
J. Stat. Phys. , 66(3-4):1011–1044, 1992.[51] D. Xiu.
Numerical Methods for Stochastic Computations . Princeton University Press,2010.[52] D. Xiu and J. S. Hesthaven. High-order collocation methods for differential equationswith random inputs.
SIAM J. Sci. Comput. , 27(3):1118–1139, 2005.[53] X. Zhu, E. M. Linebarger, and D. Xiu. Multi-fidelity stochastic collocation method forcomputation of statistical moments.