Uncertainty quantification for kinetic models in socio-economic and life sciences
Giacomo Dimarco, Lorenzo Pareschi, Mattia Zanella
TL;DR
This paper reviews the application of uncertainty quantification methods to kinetic models in socio-economic and life sciences, emphasizing the importance of handling parameter uncertainties in models describing social phenomena.
Contribution
It introduces various kinetic models based on nonlinear Vlasov and Fokker-Planck equations and discusses numerical approaches for uncertainty quantification that maintain key solution features.
Findings
Uncertainty quantification is essential in socio-economic kinetic models.
Numerical methods can preserve main features of kinetic solutions under uncertainty.
Models successfully describe social phenomena like opinion formation and wealth distribution.
Abstract
Kinetic equations play a major rule in modeling large systems of interacting particles. Recently the legacy of classical kinetic theory found novel applications in socio-economic and life sciences, where processes characterized by large groups of agents exhibit spontaneous emergence of social structures. Well-known examples are the formation of clusters in opinion dynamics, the appearance of inequalities in wealth distributions, flocking and milling behaviors in swarming models, synchronization phenomena in biological systems and lane formation in pedestrian traffic. The construction of kinetic models describing the above processes, however, has to face the difficulty of the lack of fundamental principles since physical forces are replaced by empirical social forces. These empirical forces are typically constructed with the aim to reproduce qualitatively the observed system behaviors,…
Peer Reviews
No public reviews on file for this paper yet. If you reviewed it on a platform where reviews are public (OpenReview, ICLR, NeurIPS, ICML), you can paste yours below so the community can read it here.
Videos
No videos yet. Explain this paper in a talk, walkthrough, or lecture? Add one.
Taxonomy
TopicsOpinion Dynamics and Social Influence · Complex Network Analysis Techniques · Ecosystem dynamics and resilience
Uncertainty quantification for kinetic models
in socio–economic and life sciences
Giacomo Dimarco Department of Mathematics and Computer Science, University of Ferrara, Via Machiavelli 35, 44121 Ferrara, Italy ([email protected]).
Lorenzo Pareschi Department of Mathematics and Computer Science, University of Ferrara, Via Machiavelli 35, 44121 Ferrara, Italy ([email protected]).
Mattia Zanella Department of Mathematical Sciences, Politecnico di Torino, Corso Duca degli Abruzzi 24, Torino, Italy ([email protected]).
Abstract
Kinetic equations play a major rule in modeling large systems of interacting particles. Recently the legacy of classical kinetic theory found novel applications in socio-economic and life sciences, where processes characterized by large groups of agents exhibit spontaneous emergence of social structures. Well-known examples are the formation of clusters in opinion dynamics, the appearance of inequalities in wealth distributions, flocking and milling behaviors in swarming models, synchronization phenomena in biological systems and lane formation in pedestrian traffic. The construction of kinetic models describing the above processes, however, has to face the difficulty of the lack of fundamental principles since physical forces are replaced by empirical social forces. These empirical forces are typically constructed with the aim to reproduce qualitatively the observed system behaviors, like the emergence of social structures, and are at best known in terms of statistical information of the modeling parameters. For this reason the presence of random inputs characterizing the parameters uncertainty should be considered as an essential feature in the modeling process. In this survey we introduce several examples of such kinetic models, that are mathematically described by nonlinear Vlasov and Fokker–Planck equations, and present different numerical approaches for uncertainty quantification which preserve the main features of the kinetic solution.
Contents
1 Introduction
Kinetic models describing the collective behavior of a large group of interacting agents have attracted a lot of interest in the recent years in view of their potential applications to various fields, like sociology, economy, finance and biology [1, 2, 3, 5, 7, 22, 23, 24, 33, 34, 35, 40, 51, 71, 73]. One of the major difficulties in applying the classical toolbox of kinetic theory to these systems is the lack of fundamental principles which define the microscopic dynamic. In addition, experimental results are typically non reproducible and, as a consequence, the model construction is dictated by its ability to describe qualitatively the system behavior and the formation of emergent social structures. A degree of uncertainty is therefore implicitly embedded in such models, since most modeling parameters can be assigned only as statistical information from experimental results [6, 8, 15, 50, 64].
From a mathematical viewpoint, the kinetic models we will consider in the present survey are characterized by nonlinear Vlasov–Fokker–Planck equations with random inputs taking into account uncertainties in the initial data, in the interaction terms and/or in the boundary conditions. The models describe the evolution of a distribution function , , , , , and a random field, accordingly to
[TABLE]
where is a linear operator describing the agents’ dynamics with respect to the variable,typically , is a non–local operator of the form
[TABLE]
and , for all , is a function describing the local relevance of the diffusion. We refer to [27, 41, 79, 85, 86] for an introduction to the subject in relation with kinetic theory. In the rest of the chapter, to avoid unnecessary difficulties, we will mainly restrict to the case of a one-dimensional random input distributed as . In the homogeneous case , the kinetic models are characterized by nonlinear Fokker-Planck equations.
The classic Fokker-Planck equation with uncertainties
The most classical example is represented by the linear Fokker-Planck model obtained for corresponding to
[TABLE]
where
[TABLE]
are the (conserved) mean velocity and the temperature of the particles. In the above expressions we assumed an uncertain initial data such that for all . The stationary solution in this case is represented by a Maxwellian distribution with uncertain momentum and temperature given by
[TABLE]
Opinion formation with uncertain interaction
A kinetic Fokker-Planck model of opinion formation for , where denote the two extremal opinions, corresponds to the choices [73, 82]
[TABLE]
In the above nonlocal interaction term is a function taking into account uncertainties in the compromise propensity between the agents’ opinions.
In the simple case and deterministic initial data, the model preserves the mean opinion and we can analytically compute the steady state distribution
[TABLE]
with a normalization constant.
Wealth distribution with uncertain diffusion
If we now consider a measure of the agents’ wealth, a Fokker-Planck model describing the wealth evolution of agents is obtained taking [33, 73]
[TABLE]
where the term characterizes the uncertain strength of diffusion. An explicit expression of the steady state distribution is given in the case
[TABLE]
where is the Gamma function and is the so–called Pareto exponent, which is now dependent on the random input.
Swarming models with uncertainties
As a final example we consider a kinetic model for the swarming behavior [12, 22, 23, 24, 31, 38, 52, 47, 58, 61]. In particular we focus on a model with self–propulsion and uncertain diffusion, see [9, 10]. The dynamics for the density of agents in position with velocity is described by the Vlasov-Fokker-Planck equation (LABEL:eq:MF_general) characterized by
[TABLE]
where
[TABLE]
with a localization kernel, a self–propulsion term and the uncertain noise intensity.
In the space–homogeneous case , stationary solutions have the form
[TABLE]
with a normalization constant and
[TABLE]
We stress that, in all the above reported examples, uncertainty may be present in other modeling parameters by further increasing the dimensionality and the complexity of the kinetic model.
The development of numerical methods for kinetic equations presents several difficulties due to the high dimensionality and the intrinsic structural properties of the solution. Non negativity of the distribution function, conservation of invariant quantities, entropy dissipation and steady states are essential in order to compute qualitatively correct solutions. Preservation of these structural properties is even more challenging in presence of uncertainties which contribute to increase the dimensionality of the problem. We refer to [43, 62, 81] for recent surveys on numerical methods for kinetic equations in the deterministic case.
For this reason we will focus on the construction of numerical methods for uncertainty quantification (UQ) which preserves the structural properties of the kinetic equation and, in particular, which are able to capture the correct steady state of the problem with arbitrary accuracy. We will discuss different numerical approaches based on the major techniques used for uncertainty quantification. In the deterministic case, similar approaches for nonlinear Fokker-Planck equations were previously derived in [18, 17, 29, 65, 70, 80]. Related methods for the case of nonlinear degenerate diffusion equations were proposed in [14, 28] and with nonlocal terms in [19, 21]. We refer also to [4] for the development of methods based on stochastic approximations and to [56] for a recent survey on schemes which preserve steady states of balance laws and related problems.
The simplest class of methods for quantifying uncertainty in partial differential equations (PDEs) are the stochastic collocation methods. Stochastic collocation methods are non-intrusive, so they preserve all properties of the deterministic numerical scheme, and easy to parallelize. In Section 3 we describe the structure preserving methods recently developed in [74, 75] together with a collocation approach and show how the resulting schemes preserve non negativity, conservation and entropy dissipation. In addition they capture the steady states with arbitrary accuracy and may achieve high convergence rates (spectral convergence for smooth solutions). Next in Section 4, we consider the closely related class of statistical sampling methods, most notably Monte Carlo (MC) sampling. In order to address the slow convergence of MC methods, we discuss here the development of Monte Carlo methods based on a Micro–Macro decomposition approach introduced in [44]. These methods preserve the structural properties of the kinetic problem, are capable to significantly reduce the statistical fluctuations of standard Monte Carlo and increase their computational efficiency by reducing the number of statistical samples in time. Section 5 is devoted to stochastic Galerkin methods based on generalized Polynomial Chaos (gPC). Although these deterministic methods may achieve high convergence rates for smooth solutions, they suffer from the disadvantage that they are highly intrusive and that increase the computational complexity of the problem. As a consequence, the main physical properties of the solution are typically lost at a numerical level. For this class of methods we show how to construct generalized Polynomial Chaos schemes based on the Micro–Macro formalism which preserve the steady states of the system [45]. Finally, in Section 6 several numerical applications to problem in socio-economy and life sciences are presented.
2 Preliminaries
In this section we recall some analytical properties of the considered kinetic models which will be useful for the development of the different numerical methods. Except in some simple case, a precise analytic description of the global equilibria of equation (LABEL:eq:MF_general) is very difficult [25, 26, 84]. A deeper insight into the large time behavior can be achieved by resorting to the asymptotic behavior of the corresponding space homogeneous models, leading to nonlinear Fokker–Planck type equations [83].
2.1 Fokker-Planck type equations
In the space homogeneous case the distribution function reduces to , , , and is solution of the following problem
[TABLE]
where
[TABLE]
together with an initial datum and suitable boundary conditions on .
We review in the present stochastic setting the classical results for the trend to equilibrium of the problem (2.1) in the simplified case of a one-dimensional problem with a linear drift term, i.e.
[TABLE]
Conservation of mass is imposed on the previous equation by considering suitable boundary conditions [73]. The stochastic stationary solution of equation (2.3) is given by the solution of
[TABLE]
The stochastic Fokker–Planck equation (2.3) may be rewritten in the equivalent forms
[TABLE]
which corresponds to the stochastic Landau form, whereas the stochastic non logarithmic Laundau form of the equation is the following
[TABLE]
Convergence to equilibrium is usually determined through measures of the entropy production. We define the relative entropy for all positive functions as follows
[TABLE]
we have [54]
[TABLE]
where the dissipation functional is defined as
[TABLE]
In the classical setting and , where is the temperature, the steady state is given by the Maxwellian density (1.4) with and relation (2.7) coupled with the log–Sobolev inequality
[TABLE]
leads to the exponential decay of the relative entropy as proved in the following result [83].
Theorem 1**.**
Let be the solution to the initial value problem
[TABLE]
with the initial condition with finite entropy. Then converges for all to given by (1.4) and
[TABLE]
For more general diffusion functions analogous log–Sobolev inequality are not available. A strategy to study the convergence to equilibrium is to investigate the relation of relative entropy with the relative weighted Fisher information, see [26, 54, 69, 83] for more details.
2.2 Micro–Macro formulation
In this paragraph we describe the Micro–Macro approach to kinetic equations of the form (2.1). The approach is based on the classical Micro–Macro decomposition originally developed by Liu and Yu in [68] for the fluid limit of the Boltzmann equation. This method has been fruitfully employed for the development of numerical methods by several authors (see [13, 36, 46, 37, 66, 90] and the references therein). These techniques has been also recently developed in [42, 53, 72] to construct spectral methods for the collisional operator of the Boltzmann equation that preserves exactly the Maxwellian steady state of the system. Since under suitable regularity assumptions on the initial distribution the Fokker-Planck equation admits a unique steady state solution , the Micro–Macro formulation is obtained decomposing the solution of the differential problem into the equilibrium part and the non–equilibrium part as follows
[TABLE]
where is a distribution function such that
[TABLE]
for some moments . The above decomposition (2.10) applied to the Fokker-Planck problem (2.1)–(2.2) yields the following result.
Proposition 1**.**
If the nonlinear Fokker-Planck equations (2.1)-(2.2) admits the unique equilibrium state , the differential operator defined in (2.2) with given by (1.2) may be rewritten as
[TABLE]
where is a linear operator defined as
[TABLE]
The only admissible steady state solution of the problem
[TABLE]
is given by .
The proof is an immediate consequence of the fact that at the steady state we have . Note that the steady state solution of the reformulated problem (2.12) is therefore independent of the uncertainty.
Remark 1**.**
Under suitable assumptions, see Theorem 1, one can show that exponentially decays to the equilibrium solution. As a consequence, the non–equilibrium part of the Micro–Macro approximation exponentially decays to for all .
3 Collocation methods
One of the most popular computational approaches for UQ relies on the class of collocation methods [87, 88]. These methods are non intrusive and permit to couple existing solvers for the PDEs without random inputs with techniques for the quantification of the uncertainty. Moreover, since the structure of the solution remains unchanged, numerical analysis of collocation methods is a straightforward consequence of the results obtained for the underlying method used for solving the original equation.
In the following, since the linear transport part in (LABEL:eq:MF_general) can be discretized using standard approaches, see for instance [43], we concentrate on homogeneous Fokker-Planck problem of the form (2.1)-(2.2).
Collocation methods consist in solving the problem in a finite set of nodes of the random field. In this class of methods belongs the usual Monte Carlo sampling (MC) which will be treated in Section 4. If the distribution of the random input is known, an efficient way to treat the uncertainty is to select the nodes in the random space according to Gaussian quadrature rules related to such distribution. This is straightforward in the univariate case, whereas becomes more challenging in the multivariate case [49].
For each we obtain a totally deterministic and decoupled problem since the value of the random variable is fixed. Therefore, solving this system of equations poses no difficulty provided one has a well–established deterministic algorithm. The result is an ensemble of deterministic solutions which can be post–processed to recover the statistical values of interest. For example, in the univariate case if are the Gaussian weights on corresponding to we can use the approximations
[TABLE]
where and denote the mean and the variance respectively. In the following we concentrates on the construction of numerical schemes which preserve the structural properties of the solution, like non–negativity, entropy dissipation and accurate asymptotic behavior [74, 75]. These properties are essential for a correct description of the underlying physical problem.
3.1 Structure preserving methods
In the one-dimensional case for all the Fokker-Planck equation (2.1)-(2.2) may be written as
[TABLE]
where now
[TABLE]
using the compact notation . Typically, when is a finite size set the problem is complemented with no-flux boundary conditions at the extremal points. In the sequel we assume in the internal points of .
We introduce an uniform spatial grid , such that . We denote as usual and consider a conservative discretization of (3.3)
[TABLE]
where for each is the numerical flux function characterizing the discretization.
Let us set and adopt the notations , . We will consider a general flux function which is combination of the grid points and
[TABLE]
where
[TABLE]
For example, the standard approach based on central difference is obtained by considering for all the quantities
[TABLE]
It is well-known, however, that such a discretization method is subject to restrictive conditions over the mesh size in order to keep non negativity of the solution.
Here, we aim at deriving suitable expressions for the family of weight functions and for in such a way that the method yields nonnegative solutions, without restriction on , and preserves the steady state of the system with arbitrary order of accuracy.
First, observe that at the steady state the numerical flux should vanish. From (3.6) we get
[TABLE]
Similarly, if we consider the analytical flux imposing , we have
[TABLE]
which is in general not solvable, except in some special cases due to the nonlinearity on the right hand side. We may overcome this difficulty in the quasi steady-state approximation integrating equation (3.9) on the cell
[TABLE]
which gives
[TABLE]
Now, by equating the ratio and in (3.8)–(3.11) for the numerical and exact flux respectively, and setting
[TABLE]
we recover
[TABLE]
where
[TABLE]
We have the following result [75]
Proposition 2**.**
The numerical flux function (3.6)-(3.7) with and defined by (3.12) and (3.13)-(3.14) vanishes when the corresponding flux (3.4) is equal to zero over the cell . Moreover the nonlinear weight functions defined by (3.13)-(3.14) are such that .
By discretizing (3.14) through the midpoint rule
[TABLE]
we obtain the second order method defined by
[TABLE]
and
[TABLE]
Higher order accuracy of the steady state solution can be obtained using suitable higher order quadrature formulas for the integral (3.12). We will refer to this type of schemes as structure preserving Chang-Cooper (SP-CC) type schemes.
Some remarks are in order.
Remark 2**.**
- •
If we consider the limit case in (3.16)-(3.17) we obtain the weights
[TABLE]
and the scheme reduces to a first order upwind scheme for the corresponding aggregation equation.
- •
For linear problems of the form the exact stationary state can be directly computed from the solution of
[TABLE]
together with the boundary conditions. Explicit examples of stationary states will be reported in the last section. Using the knowledge of the stationary state we have
[TABLE]
therefore
[TABLE]
and
[TABLE]
In this case, the numerical scheme preserves the steady state exactly.
- •
The cases of higher dimension may be derived similarly using dimensional splitting (see **[74]** for details).
3.1.1 Main properties
In the following we recall some results on the preservation of the structural properties, like non negativity and entropy dissipation.
Non negativity
Concerning non negativity, first we report a result for an explicit time discretization scheme [74]. We introduce a time discretization with and and consider the simple forward Euler method
[TABLE]
for all .
Proposition 3**.**
Under the time step restriction
[TABLE]
the explicit scheme (3.22) with flux defined by (3.13)-(3.14) preserves nonnegativity for all , i.e if , , .
Higher order strong stability preserving (SSP) methods [57] are obtained by considering a convex combination of forward Euler methods. Therefore, the non negativity result can be extended to general SSP methods.
In practical applications, it is desirable to avoid the parabolic restriction of explicit schemes. Unfortunately, fully implicit methods originate a nonlinear system of equations due to the nonlinearity of and the dependence of the weights from the solution. However, we have the following nonnegativity result for the semi-implicit case
[TABLE]
where
[TABLE]
We have
Proposition 4**.**
Under the time step restriction
[TABLE]
the semi-implicit scheme (3.24) preserves nonnegativity, i.e if , for all .
We refer to [74] for a detailed proof. Higher order semi-implicit approximations can be constructed following [16].
Entropy property
In order to discuss the entropy property we consider the prototype equation for all
[TABLE]
with equipped with deterministic initial distribution , and boundary conditions
[TABLE]
It can be shown that the introduced structure preserving scheme dissipates the numerical entropy [74]
Theorem 2**.**
Let us consider as in equation (3.27). The numerical flux (3.6)-(3.7) with and given by (3.12)-(3.13) satisfies the discrete entropy dissipation for all
[TABLE]
where
[TABLE]
and is the positive discrete dissipation function
[TABLE]
For more general equations the above approach does not permit to prove the entropy dissipation, see [74]. In the following, we introduce a different class of structure preserving schemes that, in addition to preservation of the steady state of the problem, ensure the entropy dissipation.
3.2 Entropic average schemes
Let us consider the general class of nonlinear Fokker-Planck equation with gradient flow structure [10, 21, 25]
[TABLE]
with the collocation nodes of the random field, and no-flux boundary conditions, where
[TABLE]
with an uncertain interaction potential. A stochastic free energy functional is defined as follows
[TABLE]
which is dissipated along solutions as
[TABLE]
where is the entropy dissipation function. The corresponding discrete free energy is given by
[TABLE]
In this case it is not possible to show that the discrete entropy functional (LABEL:eq:discrete_entropy) is dissipated by the SP–CC type schemes developed in the previous sections, see [74]. For this reason we introduce the new entropic family of flux function
[TABLE]
for all . We will refer to the above approximation of the solution at the grid point as entropic average of the grid points and . In the general case of the flux function (3.4) with non constant diffusion the resulting numerical flux reads
[TABLE]
Concerning the stationary state, we obtain immediately by imposing the numerical flux equal to zero
[TABLE]
and therefore we get
[TABLE]
By equating the above ratio with the quasi-stationary approximation (3.11) we get the same expression for for all as in (3.12)
[TABLE]
A fundamental result concerning the entropic average (3.36) is the following
Lemma 1**.**
The entropy average defined in (3.36) may be written as a convex combination with nonlinear weights
[TABLE]
where
[TABLE]
Remark 3**.**
As a consequence the Chang-Cooper type average (3.7) and the entropic average (3.36) define the same quantity at the steady state when .
On the contrary to the Chang-Cooper average the restrictions for the non negativity property of the solution are stronger. Therefore, similar to central differences, we have a restriction on the mesh size which becomes prohibitive for small values of the diffusion function . It is possible to show that the same condition is necessary also for the non negativity of semi-implicit approximations.
Concerning the entropy dissipation we can summarize the main results in the following [74]
Theorem 3**.**
The numerical flux (LABEL:eq:F_E)-(3.36) for a constant diffusion satisfies the discrete entropy dissipation
[TABLE]
where is given by (LABEL:eq:discrete_entropy) and is the discrete entropy dissipation function
[TABLE]
with the discrete version of (3.33).
Further, we can state the following entropy dissipation results for problem (3.27) in the nonlogarithmic Landau form (2.5).
Theorem 4**.**
Let us consider as in equation (3.27). The numerical flux (LABEL:eq:F_E)-(3.36) with given by (3.12) satisfies the discrete entropy dissipation
[TABLE]
where is given by (3.30) and is the positive discrete dissipation function
[TABLE]
3.3 Numerical results
We report a numerical example obtained with the collocation approach in combination with the structure–preserving numerical methods. We consider a stochastic Fokker–Planck equation with uncertainty in the initial distribution, i.e.
[TABLE]
for all with
[TABLE]
and , , . In (3.46) the diffusion coefficient is the temperature
[TABLE]
It is well–known that the steady–state solution of this problem is the Maxwellian distribution (1.4).
In the previous paragraphs we showed how an essential aspect for the accurate description of the stochastic steady state relies in the approximation of the family of integrals , , see (3.13)–(3.14). In this case, however, since the steady state is known we can evaluate exactly these weight functions as in (3.20)-(3.21). We postpone to the last section of the present contribution the discussion on numerical results obtained with more general weight functions for which no exact formulation are given. In Figure 1 (right) we report the relative error for the expectation of the solution in time. As expected the schemes are capable to capture the stochastic steady state exactly. Next in Figure 2 the evolution of the expectation for the numerical entropy is given.
4 Variance reduction Monte Carlo methods
Among the different type of techniques used in UQ, certainly the Monte Carlo methods represent one of the most popular and important class [20, 44, 55, 73]. They show all their potential when the dimension of the uncertainty space becomes very large. In addition, Monte Carlo methods are effective when the probability distribution of the random inputs is not known analytically or lacks of regularity since other approaches based on orthogonal stochastic polynomials may be impossible to use of may produce poor results.
In this section, we first describe a standard Monte Carlo approach which deals with random initial data and then we describe a modification of this algorithm which permits to strongly decrease the computational costs and increase the accuracy close to the equilibrium steady state.
4.1 The standard Monte Carlo method
We describe the method when applied to the solution of a Vlasov-Fokker-Planck type equation of the type (LABEL:eq:MF_general) with deterministic parameters and and random initial data . First we assume that the kinetic equation has been discretized by a deterministic solver in the variables , and . In this setting, the simplest Monte Carlo (MC) method for UQ in kinetic equation is based on the following steps.
Algorithm 1** (Standard Monte Carlo (MC) method).**
Sampling*: Sample independent identically distributed (i.i.d.) initial data , from the random field and approximate these over the grid (for example by piece-wise constant cell averages).* 2. 2.
Solving*: For each realization the underlying kinetic equation (LABEL:eq:MF_general) is solved numerically by the deterministic solver. We denote the solutions at time by , , where and characterizes the discretizations in and .* 3. 3.
Estimating*: Estimate the expected value of the random solution field with the sample mean of the approximate solution*
[TABLE]
The above algorithm is straightforward to implement in any existing code for the Vlasov-Fokker-Planck equations. Furthermore, the only (data) interaction between different samples is in step , when ensemble averages are computed. Thus, the MC algorithms for UQ are non-intrusive and easily parallelizable as well.
The typical error estimate that one obtains using such an approach is of the type
[TABLE]
where is a suitable norm, , , and are positive constants depending only on the second moments of the initial data and the interaction term, and , and characterize the accuracy of the discretizations in the phase-space. Clearly, it is possible to equilibrate the discretization and the sampling errors in the a-priori estimate taking , and . This means that in order to have comparable errors the number of samples should be extremely large, especially when dealing with high order deterministic discretizations. This may make the collocation Monte Carlo approach very expensive in practical applications. In order to address the slow convergence of Monte Carlo methods, we discuss in the next paragrph the development of variance reduction Monte Carlo methods.
4.2 The Micro-Macro Monte Carlo method
In order to improve the performances of standard MC methods, we introduce a novel class of variance reduction Monte Carlo methods [44]. The key idea on which they rely is to take advantage of the knowledge of the steady state solution of the kinetic equation in order to reduce both the variance and the computational cost of the Monte Carlo estimate. The method is essentially a control variate strategy based on a suitable microscopic-macroscopic decomposition of the distribution function.
We describe the method in the space homogeneous case, an example of such technique in the non homogeneous case is reported in Section 6, while we refer to [44] for a detailed discussion and extensions of such method to more general kinetic equations. Following Section 2.2 we introduce the Micro–Macro decomposition
[TABLE]
where is the steady state solution of the problem considered. Then, the idea consists in using the Monte Carlo estimation procedure only to the non equilibrium part solution of (2.12).
The crucial aspect is that the equilibrium state is zero and therefore, independent from . More precisely, we can decompose the expected value of the distribution function in an equilibrium and non equilibrium part
[TABLE]
and then exploit the fact that is known to have an estimate of the error committed by the Monte Carlo integration of type
[TABLE]
instead of
[TABLE]
where and are the variances of respectively the perturbation and the distribution function and where we have supposed for simplicity that expected value of the equilibrium part is computed with a negligible error. Now, since it is known that the perturbation goes to zero in time exponentially fast, then also its variance goes to zero, which means that in the steady state solution limit the Monte Carlo integration becomes only dependent on the way in which the expected value of the equilibrium part is computed.
The simplest version of the algorithm consists of the following steps:
Algorithm 2** (Micro-Macro Monte Carlo (M3C) method).**
Small scale sampling*: Sample independent identically distributed (i.i.d.) initial data , from the random field . For each sample compute the corresponding equilibrium state from its moments evaluated through suitable quadrature rules in based on the discretization parameter .* 2. 2.
Large scale sampling*: Select samples , and compute and approximate these over the grid (for example by piece-wise constant cell averages).* 3. 3.
Solving*: For each realization the underlying kinetic equation (2.12) is solved numerically by the deterministic solver. We denote the solutions at time by , .* 4. 4.
Estimating*: We estimate the expected value of the random solution field*
[TABLE]
with the sample mean of the approximate solution
[TABLE]
Using such an approach one obtains an error estimate of the type
[TABLE]
where now the constant depends on time and on the second moment of the solution which will vanish for large times. In fact, independently of we have that as . Therefore, the method reduces the variance of the estimator in time and asymptotically, since as , depends only on the fine scale sampling which does not affect the overall computational cost.
The efficiency of the M3C can be further improved in the case of monotonic convergence to equilibrium of the distribution function by introducing a strategy of sampling reduction at each time step. The resulting algorithm is the following
Algorithm 3** (Fast Micro-Macro Monte Carlo (FM3C) method).**
Small scale sampling*: Sample independent identically distributed (i.i.d.) initial data , from the random field . For each sample compute the corresponding equilibrium state from its moments evaluated through suitable quadrature rules in based on the discretization parameter .* 2. 2.
Large scale sampling*: Select samples , and compute and approximate these over the grid (for example by piece-wise constant cell averages).* 3. 3.
Solving*: For each realization the underlying kinetic equation (2.12) is solved numerically by the deterministic solver. This is realized at each time step as follows.*
- (a)
Advance in time*: Starting from , compute the solution with one time step of the deterministic solver.* 2. (b)
Discard samples*: At each time step we compute the variance of as*
[TABLE]
Set where denotes the integer part and discard uniformly samples. 4. 4.
Estimating*: We estimate the expected value of the random solution field*
[TABLE]
with the sample mean of the approximate solution
[TABLE]
The algorithm preserves the advantages of the simple M3C method but with a greater computational efficiency since the number of samples, and therefore the number of deterministic equations that we have to solve, decreases in time and asymptotically vanishes.
Remark 4**.**
- •
In the case the underlying uncertainty probability density function is known, the M3C method can be applied without any small scale sampling since the estimate of the expected value reduces to
[TABLE]
In this case M3C methods achieves arbitrary accuracy for large times.
- •
In contrast with Multi Level Monte Carlo (MLMC) methods **[55]**, which can produce non monotone estimators, the estimators produced by the M3C method are monotonic, i.e. mean estimator of positive quantities (such as density) is also positive, the same holds true for the entropy property.
- •
The extension of the M3C method to the non homogeneous case is straightforward, whereas the advantages of the FM3C method depend on the type of problem considered. Applications are reported in Section 6. For more general cases we refer to **[44]** where e detailed discussion is done.
4.3 Numerical results
In this section we show some results concerning the Micro-Macro Monte Carlo methods by comparing them to the standard MC method for UQ. In particular, we study the behaviors of our approach in solving the stochastic Fokker–Planck equation with uncertainty in the initial distribution (3.46)-(3.47).
In Figure 3 we report the norm of the error for the expected solution with a standard MC method. Left image shows the error for an increasing number of samples for different times, while right image shows the trend of the error over time for a different number of random inputs. The final time is set to , the number of cells in velocity is , while the stability condition gives . The maximum number of samples which furnishes the set of initial conditions is , while the solution is averaged over 100 realizations. One can clearly see the slope for the error in the left picture.
In Figure 4, the norm of the error is reported in the same setting of the MC case for the M3C method both for the left and the right images. The same number of averages and stochastic initial condition is employed. We can see how the error decreases as a function of time in an exponential fashion at the contrary of the MC case for which the error is almost independent on time.
Finally, in Figure 5 we show the behavior of the fast M3C method. The number of samples for which the time evolution of the perturbation is considered is reported on the right and it diminishes with time. The corresponding norm of the error is shown on the left. For this case, we increased the initial number of random nodes to to highlight the behavior of the fast approach. The number of evolutions of the distribution function computed diminishes exponentially.
5 Stochastic Galerkin methods
Among the various methods for UQ in PDEs, stochastic Galerkin (SG) methods based on generalized polynomial chaos (gPC) expansions are very attractive thanks to the spectral convergence property with respect to the random input [32, 59, 63, 67, 77, 78, 87, 89]. On the other hand, their intrusive nature forces a complete reformulation of the problem and standard schemes for the corresponding deterministic problem cannot be used in a straightforward way.
In particular, it is well known that, this intrusive formulation may lead to the loss of important structural properties of the original problem, like hyperbolicity, positivity and preservation of large time behavior [30, 39, 76].
In this section we analyze gPC-SG methods for the numerical approximation of stochastic Vlasov-Fokker-Planck equations in the form (LABEL:eq:MF_general). In particular, using the Micro–Macro approach in the gPC-SG setting we show how it is possible to construct methods which preserve the asymptotic behavior of the solution [45]. We mention here related approaches for kinetic equations developed in [60, 91].
We recall first some basic notions on Galerkin approximation techniques for stochastic computations.
5.1 Preliminaries on gPC-SG techniques
Let us consider the function , in the random variable , solution of the differential problem
[TABLE]
with a given differential operator. The present setup of the problem may be naturally extended to a dimensional vector of random variables .
We consider the space of polynomials of degree up to , generated by a family of orthogonal polynomials with respect to the probability density function of the random variable , namely . They form an orthogonal basis of , i.e.
[TABLE]
where is the Kronecker delta function. Let us assume that has finite second order moment, we can represent the function through the complete polynomial chaos expansion as follows
[TABLE]
where is given by
[TABLE]
The generalized polynomial chaos expansion approximates the solution of (5.1) with its -th order truncation and considers the Galerkin projections of differential problem for each
[TABLE]
Thanks to the orthogonality of the polynomial basis of the space we obtain a coupled system of purely deterministic equations
[TABLE]
These subproblems must then be solved through suitable numerical techniques. The approximation of the statistical quantities of interest are defined in terms of the introduced projections. From (5.4) being we have
[TABLE]
and thanks to the orthogonality it is possible to show that
[TABLE]
5.1.1 gPC-SG methods for Vlasov–Fokker–Planck equations
Let us consider the stochastic Vlasov–Fokker–Planck equation (LABEL:eq:MF_general) with a nonlocal drift of the form (1.2).
The gPC-SG approximation is given by the following system of deterministic differential equations
[TABLE]
where
[TABLE]
Note that, due to the nonlinearity of Fokker–Planck problems, we obtain a coupled system of deterministic Vlasov–Fokker–Planck equations describing the evolution of each projection. In vector notations we have
[TABLE]
where and the component of the matrix are given by (5.10).
In a similar way, we can derive the gPC-SG formulation of stochastic Vlasov–Fokker–Planck equations with uncertain diffusion terms.
Remark 5**.**
In case the uncertainty is present only in the initial data, and therefore , the matrix is diagonal and we need to solve the decoupled system of Vlasov type equations
[TABLE]
. Hence, a structure preserving approach as in Section 3.1 may be introduced in order to preserve the large time behavior of the collision step of each projection by defining a family of weight functions
[TABLE]
In this setting the scheme capture with arbitrary accuracy the steady state and the expected value of the numerical solution is kept nonnegative. However, for more general nonlocal type operators this approach cannot be applied for the construction of a stochastic Galerkin expansion which preserves the steady state solution and nonnegativity of the mean.
5.2 A Micro–Macro gPC approach
We discussed in the previous section how the gPC-SG method for stochastic Fokker–Planck equations generates a coupled system of partial differential equations. Although gPC-SG guarantees spectral convergence on the random field under suitable regularity assumptions, its accuracy in describing the long–time solutions of the problems is limited and depends on the particular scheme for solving the coupled system.
Let us consider suitable regularity assumptions on the initial distribution such that the stochastic Fokker–Planck problem admits the unique steady state solution . With the aim of preserving the steady states of the problem in the Galerkin setting we introduce a Micro–Macro gPC-SG scheme. Thanks to the formalism introduced in Section 5.1 and by analogy with (2.10) the Micro–Macro gPC decomposition for all reads [45]
[TABLE]
where
[TABLE]
Being equation (5.14) equivalent to require , we can reformulate the original problem in terms of . Equation (5.9) may be reformulated for all in terms of the nonequilibrium part of the Micro–Macro gPC decomposition as follows
[TABLE]
where the operator is the Galerkin projection of the quadratic operators of the collisional type defined in (2.2) and is a linear operator defined as
[TABLE]
Now, the equilibrium state of each gPC projection is and any consistent schemes for the numerical approximations of the differential terms in (5.16) admits as equilibrium state for all . For example, we can use a standard central difference approximation scheme for the differential terms in (5.16) to achieve second order accuracy for transient times and exact preservation of the steady state asymptotically.
5.3 Numerical results
We consider the evolution of the Fokker–Planck equation (3.46) with the uncertain initial condition (3.47). Following the set–up introduced in the previous section, we obtain the SG system of equations
[TABLE]
with
[TABLE]
In order to build the Micro–Macro gPC decomposition of the SG system we take advantage of the analytical solution given by the Maxwellian distribution (1.4), which can be approximated by its –order truncation as in Section 5.2. Therefore we aim at solving the modified problem for all
[TABLE]
In all our numerical examples we use second order central difference approximations of the derivatives in . In Figure 6 we compare the numerical long time solution obtained through a standard SG system (5.17) and the Micro–Macro SG system (MM). We can observe how the Micro–Macro gPC-SG method gives an accurate description of the expected steady state of the problem, on the contrary the error of the standard gPC-SG method saturates at the accuracy obtained with the central differences.
6 Other applications
In this section we present several numerical examples of stochastic Fokker–Planck and Vlasov–Fokker-Planck equations solved with the schemes introduced in the previous sections. In particular we focus on some recent models in socio–economic and life sciences as discussed in the Introduction.
6.1 Example 1: Opinion model with uncertain interactions
Let us consider a distribution function describing the density of agents with opinion whose evolution is given in terms of a stochastic Fokker-Planck equation characterized by the nonlocal term (1.5) with uncertain compromise propensity function . In the following we will solve the problem both in the collocation and in the Galerkin setting.
We consider as deterministic initial distribution for all , with
[TABLE]
with a normalization constant and let the mean opinion. We choose a uniformly distributed random input and a random interaction function of the form .
We discretize the random variable by considering the first Gauss–Legendre collocation nodes. In Figure 7 we compute the relative error for mean and variance with respect to the exact steady state (1.6) using points for the scheme with various quadrature rules adopted for the evaluation of the weights function in (3.14). Singularities at the boundaries in the integration of (3.14) can be avoided using open Newton–Cotes methods. In the sequel, we will adopt the notation , , to denote the structure preserving schemes with Chang–Cooper flux when (3.14) is approximated with second, fourth, sixth order open Newton–Cotes or Gaussian quadrature, respectively.
In Figure 8 the time evolution of the expected solution and variance are given. We can observe from the estimation of the variance the regions of higher variability of the expected solution due to uncertain interactions. The evolutions of the statistical quantities have been computed through a collocation method with quadrature points for the evaluation of (3.14).
Finally, as in Section 5.2 we consider a Micro-Macro gPC Galerkin setting based on the knowledge of the stationary solution (1.6).
In Figure 9 we present the behavior numerical error for large time where the differential terms in are solved by central differences. We report also the large time behavior for the expected solution in both schemes, where it is possible to observe how the Micro–Macro gPC is able to capture with high accuracy the steady state of the problem.
6.2 Example 2: Wealth evolution with uncertain diffusion
We consider the Fokker-Planck equation defined by (1.7) where now with represent the wealth of the agents and the uncertainty acts on the diffusion parameter. We consider the deterministic initial distribution for all with
[TABLE]
where is a normalization constant. To deal with the truncation of the computational domain in the interval , following [74], after introducing grid points we consider the quasi stationary boundary condition in order to evaluate , i.e.
[TABLE]
for all . In Figure 10 we report in a semilog scale the relative error for mean and variance with respect to the exact steady state introduced in (1.8) of the semi–implicit SP–CC scheme for several integration methods with points over with , and an increasing number of collocation nodes . The time step is chosen in such a way that the CFL condition for the positivity of the semi–implicit scheme is satisfied, i.e. see Section 3.1.1. For the tests we considered , where . We can observe how the error decays exponentially for an increasing number of collocation nodes.
Next we consider the SG-gPC formulation of the equation for the wealth evolution. Since in this case the uncertainty enters in the definition of the diffusion variable taking the analytical steady state solution of the problem is given in (1.8) and we can consider the Micro–Macro gPC scheme as in Section 5.2.
In Figure 12 we compare the error for a standard gPC approximation and the Micro–Macro gPC. In both cases central differences have been used for the differential terms in . We can see how computed at time , close to the stationary solution, decreases in relation to the number of terms of the gPC approximation whereas the standard gPC show a limited accuracy given by the error in approximating the large time behavior of the problem.
6.3 Example 3: Swarming model with uncertainties
Finally, the last example is devoted to a Vlasov-Fokker-Planck equation describing the swarming behavior of large group of agents. It is worth to observe how for this problem one steady state solution is provided by the global Maxwellian, which is a locally stable pattern, see [24, 48]. We compare the numerical solution of the problem making use of MC and M3C scheme analyzed in Section 4.
We consider an uncertain self–propelled swarming model described by the Vlasov-Fokker-Planck equation (LABEL:eq:MF_general) characterized by (1.9). This describes the time evolution of a distribution function which represents the density of individuals in position having velocity at time . The initial data consists in a bivariate normal distribution of the form
[TABLE]
where
[TABLE]
and
[TABLE]
with , , , and is a normalization constant. The uncertainty is present in the diffusion coefficient, i.e. , and it is distributed accordingly to .
We compute the solution by using the structure preserving scheme discussed in Section 3.1 for solving the homogeneous Fokker-Planck equation and we combine this method with a WENO scheme for the linear transport part. A second order time splitting approach joins the two discretization in space and velocity space. More in details, we compare a Monte Carlo collocation with the Micro-Macro collocation discussed in Section 3. The number of cells in space is fixed to , in velocity space to while the number of random inputs is fixed to . The solution is averaged over different realization and the final time is fixed to . The size of the domain is with in space and in velocity space.
In Figure 13 the time evolution of the expected distribution with respect to the uncertain variable is reported for different times computed by the MC approach together with the time evolution of the expected perturbation from the steady state solution computed with the M3C method. In Figure 14, the variance of the distribution over time and the variance of the perturbation over time are reported, the firsts computed by the MC method, the seconds with the M3C one. Finally, in Figure 15, the norm of the error for the MC and the M3C methods are reported as a function of time for different number of random inputs. The gain in computational accuracy of the M3C method is clearly evident for large times. The reference solution has been computed by a collocation method which employs the Gauss nodes as quadrature nodes with random inputs.
Acknowledgements
The research that led to the present survey was partially supported by the research grant Numerical methods for uncertainty quantification in hyperbolic and kinetic equations of the group GNCS of INdAM. MZ acknowledges support from GNCS and ”Compagnia di San Paolo” (Torino, Italy).
The reference list from the paper itself. Each links out to its DOI / PubMed record.
- 1[1] S. M. Ahn, S. Y. Ha. Stochastic flocking dynamics of the Cucker–Smale model with multiplicative white noises. Journal of Mathematical Physics , 51(10): 103–301, 2010.
- 2[2] G. Ajmone Marsan, N. Bellomo, M. Egidi. Towards a mathematical theory of complex socio–economical systems by functional subsystems representation. Kinetic and Related Models , 1(2): 249–278, 2008.
- 3[3] G. Albi, M. Herty, L. Pareschi. Kinetic description of optimal control problems and applications to opinion consensus. Communications in Mathematical Sciences , 13(6): 1407–1429, 2015.
- 4[4] G. Albi, L. Pareschi. Binary interaction algorithms for the simulation of flocking and swarming dynamics. Multiscale Model. Simul. , 11(1):1–29, 2013.
- 5[5] G. Albi, L. Pareschi, G. Toscani, M. Zanella. Recent advances in opinion modeling: control and social influence. In Active Particles Vol.1: Theory, Methods, and Applications , N. Bellomo, P. Degond, and E. Tadmor Eds., Birkhäuser–Springer: 49–98, 2017
- 6[6] G. Albi, L. Pareschi, and M. Zanella, Uncertainty quantification in control problems for flocking models, Math. Probl. Eng. 2015 (2015), 1–14.
- 7[7] G. Albi, L. Pareschi, M. Zanella. Opinion dynamics over complex networks: Kinetic modelling and numerical methods. Kinetic and Related Models , 10(1): 1–32, 2017.
- 8[8] M. Ballerini, N. Cabibbo, R. Candelier, A. Cavagna, E. Cisbani, I. Giardina, A. Orlandi, G. Parisi, A. Procaccini, M. Viale, V. Zdravkovic. Empirical investigation of starling flocks: a benchmark study in collective animal behavior. Animal Behavior , 76(1): 201–215, 2008.