Structure preserving stochastic Galerkin methods for Fokker-Planck equations with background interactions

TL;DR

This paper develops structure-preserving stochastic Galerkin methods for Fokker-Planck equations with uncertainties and background interactions, ensuring physical properties and high accuracy in transient and long-term regimes.

Contribution

It introduces novel stochastic Galerkin schemes that preserve key physical properties for Fokker-Planck equations with background interactions, applicable to collective behavior models.

Findings

Methods preserve nonnegativity, entropy dissipation, and asymptotic behavior.

Achieve second order accuracy in transient regimes and high order for large times.

Successfully applied to models with fixed and dynamic background distributions.

Abstract

This paper is devoted to the construction of structure preserving stochastic Galerkin schemes for Fokker-Planck type equations with uncertainties and interacting with an external distribution, that we refer to as a background distribution. The proposed methods are capable to preserve physical properties in the approximation of statistical moments of the problem like nonnegativity, entropy dissipation and asymptotic behaviour of the expected solution. The introduced methods are second order accurate in the transient regimes and high order for large times. We present applications of the developed schemes to the case of fixed and dynamic background distribution for models of collective behaviour.

Tables

Table 1. Table 1: Example 1. Estimation of the order of convergence toward the reference stationary state for SP–CC scheme with RK4 method. Rates have been computed using N = 21 , 41 , 81 𝑁 21 41 81 N=21,41,81 , σ 2 / 2 = 0.1 superscript 𝜎 2 2 0.1 \sigma^{2}/2=0.1 , Δ ​ t = Δ ​ w 2 / σ 2 Δ 𝑡 Δ superscript 𝑤 2 superscript 𝜎 2 \Delta t=\Delta w^{2}/\sigma^{2} .

$𝔼\left[f\right]$	$S{P}_k$
Time	2	4	6	G
1	2.0785	1.9989	2.0025	2.0026
5	1.9949	4.2572	2.2868	2.3361
10	1.9953	3.9141	6.4698	7.3367
$\text{Var}\left(f\right)$	$S{P}_k$
Time	2	4	6	G
1	2.0870	2.0001	2.0030	2.0031
5	1.9978	4.4192	2.2398	2.2789
10	1.9982	3.9309	6.6929	7.3405

Table 2. Table 2: Example 2. Estimation of the order of convergence of the scheme in the case of dynamic background g ​ ( v , t ) 𝑔 𝑣 𝑡 g(v,t) for second order semi-implicit method. The evolution of the background distribution follows an advection equation with α = 0.05 𝛼 0.05 \alpha=0.05 . The rates have been computed using N = 21 , 41 , 81 𝑁 21 41 81 N=21,41,81 , σ 2 / 2 = 0.1 superscript 𝜎 2 2 0.1 \sigma^{2}/2=0.1 , Δ ​ t = CFL ​ Δ ​ v Δ 𝑡 CFL Δ 𝑣 \Delta t=\textrm{CFL}\Delta v with CFL = 0.5 CFL 0.5 \textrm{CFL}=0.5 .

$𝔼\left[f\right]$	$S{P}_k$
Time	2	4	6	G
1	1.8976	2.0834	2.0972	2.0975
5	2.4162	4.7225	4.7940	4.7955
10	2.6446	4.5139	4.5082	4.5082
20	2.0685	2.5829	2.6271	2.6273
$\text{Var}\left(f\right)$	$S{P}_k$
Time	2	4	6	G
1	1.8834	2.1088	1.8565	1.8568
5	2.4162	4.3191	4.3937	4.3953
10	2.5793	4.5869	4.5809	4.5809
20	2.2616	2.4154	2.4678	2.4679

Equations

$${\partial }_{t}f(\theta ,v,t)={\nabla }_v\cdot \left[\mathcal{B}[g](v,t)f(\theta ,v,t)+{\nabla }_v(D(v)f(\theta ,v,t))\right],$$

$$\mathcal{B}[g](v,t)={\int }_{V}P(v,v_{\ast })(v-v_{\ast })g(v_{\ast },t)d{v}_{\ast },$$

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

$$\mathbb{E}[{\Phi }_h(\theta ){\Phi }_k(\theta )]={\int }_{I_{\Theta }}{\Phi }_h(\theta ){\Phi }_k(\theta )\Psi (\theta )d\theta =\Vert {\Phi }_h^2(\theta ){\Vert }_{L^2(\Omega )}{\delta }_{hk},$$

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

$${\hat{f}}_k(v,t)=\mathbb{E}[f(\theta ,v,t){\Phi }_k(\theta )],k=0,\dots ,M.$$

$${\partial }_t{f}^M(\theta ,v,t)={\partial }_v\left[\mathcal{B}[g](v,t)f^M(\theta ,v,t)+{\partial }_v(D(v)f^M(\theta ,v,t))\right].$$

$${\partial }_t{\hat{f}}_h(v,t)={\partial }_v\left[\mathcal{B}[g](v,t){\hat{f}}_h(v,t)+{\partial }_v(D(v){\hat{f}}_h(v,t))\right],$$

$${\hat{f}}_h(v,0)=\mathbb{E}[f(\theta ,v,0){\Phi }_h(\theta )].$$

$$\mathbb{E}[f(\theta ,v,t)]\approx {\hat{f}}_0(v,t),$$

$$\mathbb{E}[f(\theta ,v,t{)}^2]-\mathbb{E}[f(\theta ,v,t){]}^2\approx \mathbb{E}[(f^M(\theta ,v,t){)}^2]-\mathbb{E}[f^M(\theta ,v,t){]}^2,$$

$$\mathbb{E}\left[\sum _{k=0}^M{\hat{f}}_k^2(v,t){\Phi }_k^2(\theta )+2\sum _{k=0}^M\sum _{h=0}^{k-1}{\hat{f}}_k(v,t){\hat{f}}_h(v,t){\Phi }_k(\theta ){\Phi }_h(\theta )\right]-{\hat{f}}_0^2(v,t).$$

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

$$\mathcal{B}[g](v,t)+{\partial }_{v}D(v)=0,v\in \partial V.$$

$$\Vert \hat{\text{f}}(v,t){\Vert }_{L^2}={\left[{\int }_V\left(\sum _{h=0}^M{\hat{f}}_h^2(v,t)\right)dv\right]}^{1/2},$$

$$\Vert f^M(\theta ,v,t){\Vert }_{L^2(\Omega )}=\Vert \hat{\text{f}}(v,t){\Vert }_{L^2},$$

$$\Vert f^M(\theta ,v,t){\Vert }_{L^2(\Omega )}={\left({\int }_{I_{\Theta }}{\int }_V{\left(\sum _{h=0}^M{\hat{f}}_h(v,t){\Phi }_h(\theta )\right)}^{2}dv\Psi (\theta )d\theta \right)}^{1/2}.$$

$${\partial }_t\hat{\text{f}}={\partial }_v\left[\text{B}\hat{\text{f}}+\text{D}{\partial }_v\hat{\text{f}}\right],$$

$$\begin{array}{rl}{\text{B}}_{i,i}& =\mathcal{B}[g](v,t)+{\partial }_{v}D(v),{\text{B}}_{i,j}=0\\ {\text{D}}_{i,i}& =D(v),{\text{D}}_{i,j}=0.\end{array}$$

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

$$\hat{f}^{2}_{h}(v,t)\mathcal{B}[g](v,t)\Big{|}_{v\in\partial V}=0,\quad\hat{f}_{h}(v,t)\partial_{v}(D(v)\hat{f}_{h}(v,t))\Big{|}_{v\in\partial V}=0,\quad\hat{f}_{h}(v,t)D(v)\partial_{v}\hat{f}_{h}(v,t)\Big{|}_{v\in\partial V}=0$$

$$\begin{array}{rl}{\int }_V{\partial }_t\left(\frac{1}{2}{\hat{f}}_h^2(v,t)\right)dv& ={\int }_V{\hat{f}}_h(v,t){\partial }_v\left[\mathcal{B}[g](v,t){\hat{f}}_h(v,t)+{\partial }_v(D(v){\hat{f}}_h(v,t))\right]dv\\ & =\underset{A}{\underbrace{{\int }_V{\hat{f}}_h(v,t){\partial }_v\left(\mathcal{B}[g](v,t){\hat{f}}_h(v,t)\right)dv}}+\underset{B}{\underbrace{{\int }_V{\hat{f}}_h(v,t){\partial }_v^2\left(D(v){\hat{f}}_h(v,t)\right)dv}}\end{array}$$

$${\int }_V{\hat{f}}_h(v,t){\partial }_v\left(\mathcal{B}[g](v,t){\hat{f}}_h(v,t)\right)dv={\int }_V{\hat{f}}_h(v,t)\mathcal{B}[g](v,t){\partial }_v{\hat{f}}_h(v,t)dv+{\int }_V{\hat{f}}_h^2(v,t){\partial }_v\mathcal{B}[g](v,t)dv$$

$$\begin{array}{rl}& {\int }_V{\hat{f}}_h(v,t){\partial }_v\left(\mathcal{B}[g](v,t){\hat{f}}_h(v,t)\right)dv=\\ & -{\int }_V{\hat{f}}_h(v,t){\partial }_v\left(\mathcal{B}[g](v,t){\hat{f}}_h(v,t)\right)dv+{\int }_V{\hat{f}}_h^2(v,t){\partial }_v\mathcal{B}[g](v,t)dv\end{array}$$

$${\int }_V{\hat{f}}_h(v,t){\partial }_v\left(\mathcal{B}[g](v,t){\hat{f}}_h(v,t)\right)dv=\frac{1}{2}{\int }_V{\hat{f}}_h^2(v,t){\partial }_v\mathcal{B}[g](v,t)dv.$$

$$\sum _{h=0}^M{\int }_V{\hat{f}}_h(v,t){\partial }_v\left(\mathcal{B}[g](v,t){\hat{f}}_h(v,t)\right)dv=\frac{1}{2}\sum _{h=0}^M{\int }_V{\hat{f}}_h^2(v,t){\partial }_v\mathcal{B}[g](v,t)dv\le \frac{C_{\mathcal{B}}}{2}\Vert \hat{\text{f}}(v,t){\Vert }_{L^2}^2.$$

$$\begin{array}{rl}{\int }_V{\hat{f}}_h(v,t){\partial }_v^2\left(D(v){\hat{f}}_h(v,t)\right)dv& ={\int }_V\left({\partial }_v^2{\hat{f}}_h(v,t)\right)D(v){\hat{f}}_h(v,t)dv\\ & \le -C_D{\int }_V{\left({\partial }_v{\hat{f}}_h(v,t)\right)}^2,\end{array}$$

$$\begin{array}{rl}\frac{1}{2}\Vert \hat{\text{f}}(v,t){\Vert }_{L^2}^2& \le \frac{C_{\mathcal{B}}}{2}\Vert \hat{\text{f}}(v,t){\Vert }_{L^2}^2-\Vert {\partial }_v\hat{\text{f}}(v,t){\Vert }_{L^2}^2\\ & \le \left(\frac{C_{\mathcal{B}}}{2}+C_D\right)\Vert \hat{\text{f}}{\Vert }_{L^2},\end{array}$$

$$\mathcal{B}[g^{\infty }](v)f^{\infty }(\theta ,v)+{\partial }_v(D(v)f^{\infty }(\theta ,v))=0,$$

$$\frac{{\partial }_v{f}^{\infty }(\theta ,v)}{f^{\infty }(\theta ,v)}=-\frac{\mathcal{B}[g^{\infty }](v)+D^{\prime }(v)}{D(v)},$$

Full text

**Keywords: uncertainty quantification, stochastic Galerkin, Fokker-Planck equations, collective behaviour.

**

MSC: 35Q70,35Q83,65M70.

1 Introduction

Uncertainty quantification (UQ) for partial differential equations describing real world phenomena gained an increased interest in recent years [7, 11, 13, 16, 17, 18, 34, 35]. One of the main advantages of UQ methods relies in its capability to provide a sound mathematical framework to replicate realistic experiments. The introduction of stochastic parameters reflects our incomplete information on the initial configuration of a system, on its inner interactions forces and on the modelling parameters as well. Methods and ideas of UQ had a deep impact on the applied community in terms of effectivity for a robust conception, modelling and sensitivity analysis of the problem of interest.

In the context of kinetic equations, this issue can be translated on a general uncertainty affecting a distribution function of particles/agents, whose evolution is influenced by the presence of a random variable $\theta$, taking value in the set $I_{\Theta}\subseteq\mathbb{R}$, and with known probability distribution function $\Psi(\theta):I_{\Theta}\rightarrow\mathbb{R}^{+}$. In particular, in the present manuscript we are interested in Fokker-Planck type equations for the evolution of the distribution $f=f(\theta,v,t)$, $v\in V\subseteq\mathbb{R}^{d_{v}}$, $\theta\in I_{\Theta}$ and $t\geq 0$ is the time. The introduced distribution represents the proportion of particles/agents in $[v,v+dv]$ at time $t\geq 0$ and for given value of uncertainty $\theta\in I_{\Theta}$. In more details, we consider the partial differential equation

[TABLE]

where $v\in V\subseteq\mathbb{R}^{d_{v}}$ and $\mathcal{B}[\cdot]$ is the operator

[TABLE]

where $g=g(v,t)$ is a background distribution, whose dynamics do not incorporate the presence of the uncertain quantity $\theta\in I_{\Theta}$. In applications, problems with background interactions are very often considered to mimic the influence of environmental factors on the agents’ dynamics especially in socio-economic and life sciences. For example, the process of knowledge formation depends on social factors that determine the progress in competence acquisition of individuals, see [24, 25] and the references therein. Similarly, in soft-matter physics biological particles like cells undergo various heterogeneous stimuli forcing their observable motion [29]. Other examples have been studied in opinion dynamics, economic processes for the formation of wealth distributions, and urban growth theory, see [15] for a review.

We consider for (1) an uncertain initial distribution $f(\theta,v,0)$, no-flux boundary conditions are considered on the boundaries of the domain to enforce conservation of the total mass of the system. A clear understanding on the global behavior of the system governed by (1)-(2) is obtained in terms of expected statistical quantities whose accurate and physically admissible description is therefore of paramount importance.

Due to the increased dimensionality of the problem induced by the presence of uncertainties, the issue of developing fast converging numerical methods for the approximation of statistical quantities is of the highest importance. Among the most popular numerical methods for the UQ, stochastic Galerkin (SG) methods gained in recent years increasing interest since they provide spectral convergence in the random space under suitable regularity assumptions [1, 20, 35, 36, 37]. Despite their accuracy, direct application of SG techniques to high dimensional random vectors often suffers of the so-called curse of dimensionality and require ad-hoc solutions, see e.g. [31]. Similarly to classical spectral methods SG methods generally require a strong modification of the original problem and can lead to the loss of structural properties like positivity of the solution, entropy dissipation and hyperbolicity, when applied to hyperbolic and kinetic equations, see [11, 19]. The loss of structural properties of the solution induces an evident gap in its true physical meaning. To overcome this problem, recently has been proposed a novel methods that combines both Monte Carlo and SG generalized polynomial chaos methods (MCgPC) and which preserves spectral accuracy in the random space. In particular, MCgPC methods can mitigate the curse of dimensionality induced by SG-gPC. We refer to [6, 7] for a detailed discussion on these methods.

In the present manuscript we construct structure preserving methods for the SG formulation of the problem in the case of background interactions. In order to do that we will take advantages of structure preserving (SP) methods [26, 27], that have been designed to preserve the mentioned structural properties of the solution of nonlinear Fokker-Planck equations without restriction on the mesh size. We consider applications of the developed schemes both in case of fixed and dynamic background.

The rest of the paper is organized as follows. In Section 2 we briefly introduce stochastic Galerkin methods for the problem of interest where the interactions take place with respect to a deterministic background, stability results are proved and discussed together with the analysis of trends to asymptotic states. In Section 3 we derive structure preserving methods in the Galerkin setting, positivity conditions for explicit and semi-implicit schemes are discussed and we prove entropy inequality for a class of one dimensional Fokker-Planck models. Several applications of the schemes are finally considered in Section 4 for several problems arising in the description of collective phenomena in socio-economic and life-sciences. Some conclusions are reported at the end of the manuscript.

2 Stochastic Galerkin methods for kinetic equations

For simplicity of presentation we consider the case $d_{v}=1$. We focus on real-valued distributions depending on a one dimensional random input. Let $(\Omega,F,P)$ be a probability space where as usual $\Omega$ is the sample space, $F$ is a $\sigma-$algebra and $P$ a probability measure, and let us defined a random variable

[TABLE]

with $I_{\Theta}\subset\mathbb{R}$ and $\mathcal{B}_{\mathbb{R}}$ is the Borel set. We focus on real-valued distributions of the form $f(\theta,v,t):\Omega\times V\times[0,T]\rightarrow\mathbb{R}^{d}$. In the present section we derive a stochastic Galerkin approximation for Fokker-Planck equation with uncertain initial distribution and background interactions (1).

Let us consider the linear space $\mathbb{P}_{M}$ of polynomials of degree up to $M$ generated by a family of orthogonal polynomials $\{\Phi_{h}(\theta)\}_{h=0}^{M}$ such that

[TABLE]

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

[TABLE]

being $\hat{f}_{k}(v,t)$ the projection of $f$ into the polynomial space of degree $k$, i.e.

[TABLE]

Plugging $f^{M}$ into (1) we obtain

[TABLE]

Hence, by multiplying (4) by $\Phi_{h}(\theta)$ for all $h=0,\dots,M$ and after projection in each polynomial space we obtain the following system of $M+1$ deterministic kinetic-type PDEs

[TABLE]

with the initial conditions

[TABLE]

The related deterministic subproblems can be tackled through suitable numerical methods and the approximation of statistical quantities of interest are defined in terms of the projections. In particular we have

[TABLE]

whose evolution is given by (5) in the case $h=0$. Thanks to the orthogonality in $L^{2}(\Omega)$ of the polynomials $\{\Phi_{h}\}_{h=0}^{M}$ we have

[TABLE]

and from (3) we have

[TABLE]

Therefore the variance of the solution is approximated in terms of the projections as follows

[TABLE]

We observe that the initial mass defined by $\int_{V}\hat{f}_{h}(v,0)dv$ is conserved in time assuming no-flux boundary conditions, i.e.

[TABLE]

Let us introduce the vector $\hat{\textbf{f}}(v,t)=\left(\hat{f}_{0}(v,t),\dots,\hat{f}_{M}(v,t)\right)$. If we define as $\|\hat{\textbf{f}}(v,t)\|_{L^{2}}$ the standard $L^{2}$ norm of the vector $\hat{\textbf{f}}(v,t)$

[TABLE]

then from the orthonormality of the introduced basis $\{\Phi_{h}\}_{h=0}^{M}$ in $L^{2}(\Omega)$ we have that

[TABLE]

where

[TABLE]

We can reformulate the problem (5) in a more compact form as follows

[TABLE]

where $\textbf{B}=\{B_{ij}\}_{i,j=1}^{M+1}$ and $\textbf{D}=\{D_{i,j}\}_{i,j=1}^{N}$ are diagonal matrices with components

[TABLE]

The following stability result can be established

Theorem 1**.**

If $\|\partial_{v}\mathcal{B}[g](v,t)\|_{L^{\infty}}\leq C_{\mathcal{B}}$, with $C_{\mathcal{B}}>0$, and if $D\leq C_{D}$ we have for all times $t>0$

[TABLE]

provided that the boundary terms vanish for all $h=0,\dots,M$

[TABLE]

Proof.

We multiply (5) by $\hat{f}_{h}(v,t)$ and integrate over $V\subseteq\mathbb{R}$

[TABLE]

The integral $A$ may be rewritten as follows by direct computation

[TABLE]

and, integrating by parts, it is equivalent to

[TABLE]

being \hat{f}^{2}_{h}(v,t)\mathcal{B}[g](v,t)\Big{|}_{v\in\partial V}=0. Therefore, the integral $A$ may be written as

[TABLE]

Hence, the following estimate holds

[TABLE]

Furthermore, for the integral $B$ in (8) we have

[TABLE]

provided \hat{f}_{h}(v,t)\partial_{v}(D(v)\hat{f}_{h}(v,t))\Big{|}_{v\in\partial V}=0, and \hat{f}_{h}(v,t)D(v)\partial_{v}\hat{f}_{h}(v,t)\Big{|}_{v\in\partial V}=0.

Finally, after summation on $h=0,\dots,M$ of the obtained bounds, we get

[TABLE]

and thanks to the Gronwall’s theorem we can conclude. ∎

Remark 1**.**

The background distribution $g(v,t)$ is in general ruled by an additional PDE that does not depend on the stochastic density function $f(\theta,v,t)$ and does not incorporate additional uncertainties. In the case of evolving background we need to couple to (1) its dynamics.

2.1 Asymptotic behaviour

Under suitable smoothness assumptions the introduced Fokker-Planck equation has a unique smooth solution, see e.g. [28, 30]. In the present section we concentrate on the large time solution of the introduced class of problems, known in the literature as equilibrium solutions or steady states. Assuming that the dynamics of the background $g(v,t)$ admit a unique stationary state the asymptotic distribution of (1) is solution of the differential equation

[TABLE]

which gives

[TABLE]

and therefore the analytical stationary distribution of the original problem reads

[TABLE]

being $C(\theta)>0$ a normalization constant depending only on the initial uncertainties of the problem.

On the other hand, the asymptotic solutions $f^{\infty}_{h}(v)$ of (5) in each polynomial space of degree $h=0,\dots,M$ are defined by solving the following set of differential equations

[TABLE]

whose stationary states are

[TABLE]

being $C_{h}$ such that

[TABLE]

We can observe how, if the initial state has deterministic mass $\int_{V}f(\theta,v,0)dv=\bar{\rho}>0$, the asymptotic state of the problem given by (10) does not incorporate any uncertainty. In particular, this is induced by the fact that the normalization constant does not depend anymore on the uncertainty of the problem, i.e. $C(\theta)=\bar{C}$ for all $\theta\in I_{\Theta}$. This fact reflects on the asymptotic state of each projection $\hat{f}_{h}^{\infty}(v)$, $h=0,\dots,M$, since $\mathbb{E}[\bar{C}\Phi_{h}(\theta)]=0$ for $h>0$. Therefore, in the case of deterministic initial mass we obtain

[TABLE]

and the variance of $f(\theta,v,t)$ vanishes asymptotically. In the general case of uncertain initial mass the asymptotic state still depends on $\theta\in I_{\Theta}$.

In the following we explicit the trend to equilibrium defined by stochastic background interaction models following the ideas in [14].

2.1.1 Constant background

Let us assume that the background is fixed, i.e. $\mathcal{B}[g](v,t)=\mathcal{B}[g](v)$. In particular, the stationary distribution of (1) is solution of the following differential equation

[TABLE]

Taking advantage of the condition (13) we can observe that (1) admits several equivalent formulations. Indeed we have for all $t\geq 0$

[TABLE]

In particular from (13) it follows that the Fokker-Planck equation (1) with constant background can be rewritten in Landau form

[TABLE]

or equivalently in the non-logarithmic Landau form

[TABLE]

From these reformulations we can obtain the evolution for $F(\theta,v,t)=\frac{f(\theta,v,t)}{f^{\infty}(\theta,v)}$ since

[TABLE]

Therefore, from (13) we get

[TABLE]

with no-flux boundary conditions

[TABLE]

In the following we summarise several findings on the entropy production of Fokker-Planck problems. Their importance to detected the correct trends to equilibrium is essential in kinetic theory, we point the interested reader to [32] for more details in the classic Fokker-Planck setting and to [22] for results in the nonconstant diffusions setting.

In the case of interaction with a constant background distribution the following result holds

Theorem 2**.**

Let the smooth function $\Phi(x)$, $x\in\mathbb{R}_{+}$ be convex. Then, if $F(\theta,v,t)$ is the solution to (14) in $V\subseteq\mathbb{R}$ and $F(\theta,v,t)$ is bounded for all $\theta\in I_{\Theta}$ the functional

[TABLE]

is monotonically decreasing in time and its evolution is given by

[TABLE]

where with $\mathcal{I}$ we denote the nonnegative quantity

[TABLE]

Proof.

The proof of this result follows the strategy adopted in [14] for all $\theta\in I_{\Theta}$. ∎

Now in the case $\Phi(x)=x\log(x)$ we obtain the relative Shannon entropy $\mathcal{H}(f,f^{\infty})(\theta,t)$ which is a functional depending on the uncertainties of the model. From the above result it follows that this quantity is dissipated with the rate given for all $\theta\in I_{\Theta}$ by

[TABLE]

and we have

[TABLE]

In the stochastic Galerkin approximation the relative Shannon entropy for $f^{M}(\theta,v,t,)$ in (3) reads

[TABLE]

from which approximated statistical moments can be obtained by projection in the space defined by the polynomial basis

[TABLE]

being

[TABLE]

We observe that, due to the nonlinearities in the definition of the convex functional $\mathcal{H}(f,f^{\infty})$, a coupled system of differential equations must be solved to estimate the expected trends to equilibrium provided by the relative entropy functional. Nevertheless, at the Galerkin level we have no guarantee that the weighted Fisher information defines a positive quantity for the obtained truncated distribution and, hence, that the entropy monotonically decreases.

On the other hand, following the same computations presented in the beginning of this section, the system of $M+1$ projections defined in (5) can be rewritten for all $h=0,\dots,M$ in the case of fixed background as follows

[TABLE]

Therefore, by introducing the ratio $F_{h}=\frac{\hat{f}_{h}(v,t)}{\hat{f}^{\infty}_{h}(v)}>0$ we have

[TABLE]

complemented with no-flux boundary conditions. Hence, in analogy with what we discussed above, the following result holds.

Theorem 3**.**

Let the smooth function $\Phi(x)$, $x\in\mathbb{R}_{+}$ be convex. Then, if $F_{h}(v,t)$ is the solution to (16) in $V\subseteq\mathbb{R}$ and $F_{h}(v,t)$ is bounded the functional

[TABLE]

is monotonically decreasing in time and its evolution is given by

[TABLE]

where with ${\mathcal{I}}$ we denote the nonnegative quantity

[TABLE]

Now, in the case of relative Shannon entropy $\Phi(x)=x\log x$ we obtain in each polynomial space

[TABLE]

Therefore, each projection of $f(\theta,v,t)$ in the linear space of arbitrary degree $h=0,\dots,M$ converges monotonically in time to its equilibrium $\hat{f}^{\infty}_{h}(v)$. In particular this is true for the expected quantities of the problem.

3 Structure preserving methods

In this section we introduce the class of so-called structure preserving (SP) numerical methods for the solution of Fokker-Planck equations with nonlocal terms. These methods preserve the fundamental structural properties of the problem like nonnegativity of the solution, entropy dissipation and capture the steady state of each problem with arbitrarily accuracy, see [8, 13, 26, 27]. The applications of the SP methods is here particularly appropriate since, thanks to background interactions, the system of $M+1$ equations (5) is decoupled.

In the following we summarise the construction ideas at the basis of SP methods in dimension $d=1$, extension to general dimension can be found in [26]. The case of nonlocal Fokker-Planck equations with anisotropic diffusion has been studied in [21].

3.1 Derivation of the SP method

For all $h=0,\dots,M$ we may rewrite (5) in flux form as follows

[TABLE]

where

[TABLE]

and $\mathcal{C}[g](v,t)=\mathcal{B}[g](v,t)+\partial_{v}D(v)$. Let us introduce a uniform grid $v_{i}\in V$, such that $v_{i+1}-v_{i}=\Delta v>0$ and let $v_{i\pm/2}=v_{i}\pm\Delta v/2$. We consider the conservative discretization

[TABLE]

where $\hat{f}_{h,i}(t)$ is an approximation of $\hat{f}_{h}(v_{i},t)$ and where $\mathcal{F}_{h,i\pm 1/2}$ a numerical flux having the form

[TABLE]

where

[TABLE]

Hence, we aim at finding the weight functions $\delta_{i+1/2}$ and $\tilde{\mathcal{C}}[g]_{i+1/2}$ such that the scheme produces nonnegative solutions without restrictions on the mesh size $\Delta v$, and is able to capture with arbitrary accuracy the steady state of the (5) for all $h=0\dots,M$.

We observe that for a vanishing numerical flux from (18) we obtain

[TABLE]

At the analytical level we obtained from (11) in Section 2.1 that

[TABLE]

which admits the quasi-steady state approximation in the cell $[v_{i},v_{i+1}]$ for all $h=0,\dots,M$

[TABLE]

that is

[TABLE]

Now equating $\hat{f}_{h}(v_{i+1},t)/\hat{f}_{h}(v_{i},t)$ in (20) and $\hat{f}_{h,i+1}(t)/\hat{f}_{h,i}(t)$ in (19) and setting

[TABLE]

we can determine weight functions

[TABLE]

where

[TABLE]

It is worth pointing out that, by construction, the numerical flux of the SP scheme vanishes when the analytical flux is equal to zero. Furthermore, the long time behavior of (11) is described with the accuracy with which we evaluate the weights (21). At the numerical level this fact reflects on the considered quadrature methods as we will observe in Section 4. In the following, we will show that suitable restrictions on the time discretization can be defined to guarantee positivity preservation of the SP scheme, without any restriction on the mesh for $v\in V$. Moreover, we will show that the scheme dissipates the numerical entropy with a rate which is coherent with what we observed in Section 2.1.

Remark 2**.**

The obtained weights do not depend on the degree of the linear space since they are equal for all $h=0,\dots,M$. Furthermore, in the case of interaction with a constant background, i.e. $\mathcal{B}[g](v,t)=\mathcal{B}[g](v)$, a stationary state $\hat{f}^{\infty}_{h}(v)$ is defined for all $h=0\dots,M$, see equation (11) together with boundary conditions. Hence, thanks to the knowledge of the stationary state in each polynomial space we have

[TABLE]

Which leads to

[TABLE]

and

[TABLE]

We highlight how in this case the dependence on $h=0,\dots,M$ is only apparent since for each times $t\geq 0$ the ratio $\hat{f}_{h,i+1}/\hat{f}_{h,i}$ in (20) does not depend on the specific projection thanks to background interactions. In addition, note that if the steady state is analytically known, for some special form of the operator $\mathcal{B}[\cdot]$ and of the diffusion function, the SP scheme does not introduce any additional source of errors on the steady state distribution, see [8, 26].

3.2 Positivity of statistical moments

In general positivity of the solution, or of its statistical moments, is not achievable once we apply stochastic Galerkin methods and the solution of the system $f^{M}(\theta,v,t)$ looses a genuine physical meaning. In this section we provide explicit conditions to preserve nonnegativity of projections $\hat{f}_{h}(v,t)$ and, therefore, of the statistical moments of $f^{M}(\theta,v,t)$, that have been obtained in Section 2 from direct inspection of the Galerkin projections (6)-(7). In particular, we will show how in the background interactions case we are able to provide reliable conditions, without restriction on $\Delta v$, for positivity preservation.

In recent works [6, 7] a particle scheme has been proposed to enforce positivity of statistical quantities for uncertainty quantification of kinetic models. The core idea of the approach presented in the cited works is to approximate the expected solution of a mean-field type model by its Monte Carlo (MC) formulation in the phase space, which is then expanded through a SG generalized polynomial chaos (SG-gPC) method. The expected solution is then reconstructed from expected positions and velocities of the microscopic system, which is considered in the gPC setting. The authors of these works refer to this method as MCgPC. The solution of the MCgPC approach is still spectrally accurate in the random space whereas in the phase space it assumes to accuracy of the Monte Carlo method. The approach presented in the present manuscript for the linear case provide instead global high accuracy.

Let us introduce the time discretization $t^{n}=n\Delta t$, $\Delta t>0$ and $n=0,\dots,T$ and consider the following forward Euler method for all $h=0,\dots,M$

[TABLE]

where now $\hat{f}_{h,i}^{n}$ is an approximation of $\hat{f}_{h}(v_{i},t^{n})$ and the flux has the form introduced in (18). We can prove the following result

Theorem 4**.**

Under the time step restriction

[TABLE]

the explicit scheme (22) preserves nonnegativity, i.e. $\hat{f}_{h,i}^{n+1}\geq 0$ provided $\hat{f}_{h,i}^{n}\geq 0$.

Proof.

The proof of this result is analogous for all $h=0,\dots,M$ to the result for explicit scheme obtained in [26]. ∎

We observe that no explicit dependence on the expansion degree $h=0,\dots,M$ appears in the derived restriction thanks to the background-type interactions. Furthermore, the restriction on $\Delta t$ in Theorem 4 ensures nonnegativity without additional bounds on the spatial grid as for example happen for central type schemes. The derived condition automatically holds for higher order strong stability preserving (SSP) methods like Runge-Kutta and multistep methods since these are convex combinations of the forward Euler integration. The proved nonnegativity of the scheme is extended straightforwardly to each SSP-type time integration.

We highlight how the derived parabolic restriction to enforce nonnegativity of explicit schemes can be quite heavy for practical applications. A convenient strategy to lighten this burden resorts to the technology of semi-implicit methods, see [3] for an introduction. Indeed, we can prove nonnegativity of the numerical approximations of solutions $\{\hat{f}_{h}\}_{h=0}^{M}$ by considering the set of modified fluxes

[TABLE]

The scheme is semi-implicit since we compute the background dependent $\tilde{\mathcal{C}}_{i+1/2}(t)$ and weight functions $\delta_{i+1/2}(t)$ at time $t^{n}$. As a consequence, it is easily seen how in the case of a fixed background the scheme is coherent with a fully implicit method.

The following result holds

Proposition 1**.**

Let us consider a semi-implicit method for all $h=0,\dots,M$

[TABLE]

with fluxes defined in (23). Under the time step restriction

[TABLE]

the semi-implicit scheme preserves nonneagivity, i.e. $\hat{f}_{h,i}^{n+1}\geq 0$ if $\hat{f}_{h,i}^{n}\geq 0$ for all $i=1,\dots,N$ and $h=0,\dots,M$.

Proof.

The proof of this result is analogous for all $h=0,\dots,M$ to the result for semi-implicit scheme obtained in [26]. ∎

Extensions to higher order semi-implicit schemes have been obtained in [3].

3.3 Entropy dissipation

In Section 2.1 we have observe how an entropy can be defined for the introduced class of problems. In details, a Shannon entropy operator has been defined and we proved that it is dissipated in time, leading to interpret the steady state as the zero entropy state of the system. At the numerical level, a remarkable property of a scheme relies on the correct description of the trends toward equilibrium. In the following we will show how the introduced methods are entropic.

We concentrate on the case of fixed background. In Section 2.1 we have seen how the Fokker-Planck problems of interest can be rewritten in Landau form (15). In particular, it can be proven how the numerical flux for this reformulation is given by the following equivalent form

[TABLE]

with

[TABLE]

since for all $h=0,\dots,M$ we are looking at the constant background case we have $\lambda_{i+1/2}^{\infty}=\log\hat{f}_{h,i}^{\infty}-\log f_{h,i+1}^{\infty}$ and the weight functions are rewritten as

[TABLE]

see Remark 2. We can prove the following result

Theorem 5**.**

Let us consider the conservative discretization (17) for all $t\geq 0$ and $h=0\dots,M$. The numerical flux (18) satisfies the discrete entropy dissipation

[TABLE]

where

[TABLE]

and

[TABLE]

Proof.

From the definition of relative entropy for all $h=0,\dots,M$ we have

[TABLE]

After summation by parts we have

[TABLE]

and from the reformulation of the flux in (24) we may conclude since $(x-y)\log(x/y)\geq 0$ for all $x,y\geq 0$. ∎

4 Numerical tests

In the present section, we present several tests for Fokker-Planck equations with background interactions and uncertain initial distribution. In particular, we numerically show the performance of the discussedstructure preserving stochastic Galerkin (SP-SG) method. Both the scenarios with fixed and evolving backgrounds will be considered. As discussed in Section 3 the essential aspect for the accurate computation of the large time distribution of the problem (1) lies in the numerical approximation of the integral

[TABLE]

which defines the quasi-stationary states of each projection. In general a high order quadrature method is needed. In the following numerical examples we will consider open Newton-Cotes quadrature methods up to the $6th$ order and the Gauss-Legendre quadrature, see e.g. [10]. Through the text we will refer to these methods as $SP_{k}$, $k=2,4,6,G$, where the index $k$ indicates the order of the adopted quadrature method with $G$ referring to the Gauss-Legendre case. To highlight the advantages of this approach, a nonconstant diffusion function is considered for bounded domains. In all the tests we considered suitable restrictions on the time discretization to guarantee positivity of the expected solution of the problems both in the explicit and semi-implicit integration. Extension to the multidimensional case is considered at the end of this section.

4.1 Test 1: Stationary background distribution

Let us consider the evolution of a distribution function $f(\theta,v,t)$ in the presence of uncertainty that follows (1), with $v\in[-1,1]$, and interacting with a given background distribution $g(v,t)=g(v)$ for all $t\geq 0$ of the form

[TABLE]

with $\beta>0$ a constant such that $\int_{-1}^{1}g(v)dv=1$. We consider in this test a nonconstant diffusion $D(v)=\frac{\sigma^{2}}{2}(1-v^{2})^{2}$ with given $\sigma^{2}$ that will be specified later on. Furthermore, the nonlocal operator in (2) is defined in terms of the interaction function

[TABLE]

where $\Delta>0$ is a constant measuring the maximal distance under which interactions may occur. The introduced function $P(\cdot,\cdot)$ is usually defined as bounded confidence function. This model has been proposed in the literature to describe the evolution of the distribution of agents having opinion $v$ at time $t\geq 0$, see [23, 33]. In particular, the presence of background interactions is generally considered to take into account the influence of external actors in opinion dynamics like the case of media [4] or the action of possible control strategies [2]. Extensions to the case of uncertain interactions have been proposed in [34].

In this first test we consider as initial distribution

[TABLE]

with $C(\theta)$ such that $\int_{V}f(\theta,v,0)=\rho(\theta)>0$ for all $\theta\in I_{\Theta}$ and $u_{i}(\theta)$, $i=1,2$ given by

[TABLE]

being $\theta\sim\mathcal{U}([-1,1])$. In the case $\Delta=2$ it follows that $P\equiv 1$ and we can compute the explicit stationary distribution

[TABLE]

The stochastic Galerkin decomposition of the resulting problem can be performed by considering a Legendre polynomial basis $\{\Phi_{h}\}_{h=0}^{M}$ being $\Psi(\theta)=\frac{1}{2}\chi(\theta\in[-1,1])$. The resulting system of equations have the form (5) whose asymptotic solution for all $h=0,\dots,M$ reads

[TABLE]

being $C_{h}=\frac{1}{2}\int_{I_{\Theta}}C(\theta)\Phi_{h}(\theta)d\theta$. In Figure 1 we present the evolution of the $L^{1}$ relative error computed with respect to the exact stationary state for the $SP_{k}$, $k=2,4,6,G$, schemes for various quadrature methods. To exemplify the advantages, we consider two projections $h=0$ (left) and $h=1$ (right). For each $SP_{k}$ we considered $N=41$ gridpoints for the discretization of the state variable.

We can observe how we achieve different accuracy in terms of the steady states of the problem in relation to the considered quadrature rules for both $h=0,1$. Further, with low order quadrature we approach to the numerical steady state of the method faster than with high order rules. We observe that with a Gauss-Legendre method we essentially reach machine precision in finite time for each projection. In the same figure, we show the dissipation of the relative entropy functional $\mathcal{H}(\hat{f}_{h},\hat{f}_{h}^{\infty})$ discussed in Section 2.1 with $h=0,1$ obtained with the structure preserving method. We present the case of two coarse grids obtained with $N=11$, $N=21$ gridpoints compare with the exact dissipation of the relative entropy.

The high accuracy of the scheme in the description of the large time behavior, in each polynomial space, reflects in a high accuracy in the approximation of statistical moments of the solution of the problem, see Figure 2. Here, we considered the schemes $SP_{2}$ and $SP_{G}$, that is the structure preserving schemes with approximation of $\lambda_{i+1/2}$ with a $2$nd order and Gauss-Legendre method respectively. We highlight how the approximated expectation $\hat{f}_{0}$ is positive thanks to the properties of the scheme.

In Table 1 we estimate the order of convergence of the SP-SG method in terms of accuracy of the expected quantities $\mathbb{E}[f]$, and $\textrm{Var}(f)$ in their stochastic Galerkin approximation (6)-(7). It is easily observed how for the approximation of the variance it is required the solutions of the whole set of projections $h=0,\dots,M$. In the present test we considered $M=10$. Furthermore, we used $N=21,41,81$ and the order of convergence of the explicit structure preserving schemes is measured as $\log_{2}\frac{e_{1}(t)}{e_{2}(t)}$, where $e_{1}(t)$ is the relative error at time $t\geq 0$ of the expected solution and its variance computed with $N=21$ gridpoints with respect to that computed with $N=41$ gridpoints and, likewise, $e_{2}(t)$ is the relative error at time $t\geq 0$ computed with $N=41$ with respect to that computed with $N=81$ gridpoints. The time integration has been performed with RK4 at each time step chosen in such a way that the restriction for positivity of the scheme in Theorem 4 is satisfied, i.e. $\Delta t=O(\Delta v^{2})$. We can observe how the $SP_{k}$ schemes are second order accurate in the transient regimes and assume the order of the quadrature method near the expected steady state and its related variance.

More generally, $\Delta<2$ and therefore $P(v,v_{*})\neq 1$, so that we have no analytical insight on the large time solution $\hat{f}^{\infty}_{h}(v)$ in each polynomial space. In Figure 3 we consider the case of bounded confidence type interactions (26) with $\Delta=1.0$ and a fixed background distribution $g(v)$ of the form

[TABLE]

with $u_{g}=\frac{1}{2}$ and $\sigma_{g}^{2}=10^{-2}$. We considered the uncertain initial density (27) with deterministic initial mass $\rho(\theta)=1$ and uncertainty in $u_{1}(\theta)$, $u_{2}(\theta)$ so that

[TABLE]

with $\theta\sim\mathcal{U}([-1,1])$. The integral $\mathcal{B}[g](v)$ has been evaluated through a trapezoidal rule. As observed in Section 2.1 the large time solution for all $h=0,\dots,M$ does not depend on the uncertainties of the initial distribution. Indeed, the variance annihilates as we can observe in 3(d) and the asymptotic state coincides with $\mathbb{E}[f]$.

4.2 Example 2: Evolving background distribution

In this section we test the performance of the introduced structure preserving stochastic Galerkin scheme in the case of an evolving background distribution. To exemplify a dynamic background distribution we consider the solution of a linear advection equation

[TABLE]

which is coupled to the original stochastic Fokker-Planck equation in (1) through the operator $\mathcal{B}[g](v,t)$. The initial background is considered of the form (25), with $u_{g}=-\frac{1}{2}$, we consider periodic boundary conditions for (29) and $\alpha=0.05$. The advection equation is solved numerically with a Lax-Wendroff scheme for each time $t\geq 0$. In the following we consider as uncertain initial distribution (27) with $\bar{u}=0.5$, $\kappa=0.25$, and the mass is $\rho(\theta)=1+\frac{1}{2}\theta$ with uniform perturbation $\theta\sim\mathcal{U}([-1,1])$.

In Table 2 we estimate the order of convergence of the stochastic Galerkin structure preserving scheme in terms of the two statistical quantities $\mathbb{E}(f)$, $\textrm{Var}(f)$ in the case of a background distribution evolving as (29). The numerical integration is a second-order semi-implicit method, see [3, 26]. The approximation of the variance we considered $M=5$ projections. We can observe that the dynamic background distribution prevents the formation of steady state solution in the original Fokker-Planck problem. Indeed, for each $SP_{k}$, $k=2,4,6,G$ the scheme initially increases its order according to the quadrature method and for large times it is reduced to the initial second-order.

In Figure 4 we can observe the behavior of (1) in the bounded domain $v\in[-1,1]$ and interacting through a bounded confidence type $P(v,v_{*})$ with $\Delta=1$, where the evolving background follows the advection (29).

4.3 Example 3: 2D model of swarming

We consider a kinetic swarming model with self-propulsion and diffusion with uncertain initial distribution. In a deterministic framework we refer to [5] where a similar model has been studied in the nonlocal setting. Here the authors proved a sharp phase transition between self-propulsion forces and diffusion, that discriminate the minimal amount of noise needed to obtain symmetric distribution with zero mean. The study of possible uncertain quantities in the dynamics is here of paramount importance, since coefficients like the noise intensity and the self-propulsion strength, are commonly based on field observations and empirical evidence. We refer to [6] for a more detailed analysis of the influence of uncertain quantities in problems with phase transition. In the following, we concentrate on the case of uncertain initial distribution.

We consider a model for the evolution of the density of individuals $f=f(\theta,v,t)$ having velocity $v\in\mathbb{R}^{2}$ at time $t\geq 0$ and uncertain initial condition $f(\theta,v,0)$ having mass $\rho(\theta)$. In details the model reads

[TABLE]

being $\alpha>0$ the self-propulsion strength and $D>0$ a constant noise intensity. At a difference with the original nonlinear case here the agents interact with a background distribution $g(v)$ through its mean velocity $u_{g}=\int_{V}vg(v)dv$. It may be shown that a free energy functional is defined which dissipates along solutions. Further, stationary solutions have the form

[TABLE]

with $C(\theta)>0$ a normalization constant. In particular, we focus on the 2D case and we consider the fixed background distribution

[TABLE]

The extension of the presented structure preserving methods to the multidimensional case has been established in [26, 27]. The idea is to apply a structure preserving scheme to each dimension of the stochastic Galerkin projections

[TABLE]

with initial distribution $\hat{f}_{h}(v,0)=\mathbb{E}[f(\theta,v,0)\Phi_{h}(\theta)]$. In Figure 5 we present the large time distributions for the choices of diffusion $D=0.2$, $D=0.8$ and three configuration of the background distribution. We applied the $SP_{G}$ scheme to (32) for $h=0,\dots,5$ to obtain the approximation of expected distribution and of the variance. The initial distribution is here such that $\int_{V}f(\theta,v,0)dv=1+\frac{1}{2}\theta$, $\theta\sim\mathcal{U}([-1,1])$.

A uniform grid over $[-4,4]\times[-4,4]$ with $N=51$ gridpoints in each direction has been considered. The integration over the time interval $[0,100]$ has been performed taking advantage of the second order semi-implicit method with $\Delta t=O(\Delta v)$. The surfaces represent the expected solution whereas the red grids represent the upper confidence band that may be computed as $\mathbb{E}[f]+\sqrt{\textrm{Var}(f)}$. The provided confidence bands gives an indication on the regions where $f(\theta,v,t)$ is more affected by the action of the uncertainty $\theta\in I_{\Theta}$.

We can clearly observe the influence of the background in shaping the large time distribution of the problem, which is steered towards the background mean. The computed confidence bands, furthermore, make clear how the behavior is stable under the action of initial uncertainties.

Conclusion

We studied the application of structure preserving type schemes to the stochastic Galerkin approximation of Fokker-Planck equations with uncertain initial distribution and background interactions. The developed methods are capable to preserve the stationary state of the problem with arbitrary accuracy and define nonnegative expected solutions under suitable time step restrictions. Both explicit and semi-implicit type time integrations have been taken into account. Furthermore, we have proven discrete relative entropy dissipation property for the derived scheme for each projection of the original model. Several applications to prototype problems in socio-economic and life sciences have been proposed both in case of fixed and evolving background distribution together with the extension of the method to the multidimensional case. Extensions of the scheme to the uncertain background case and to the case of vanishing diffusion are under study both in the deterministic and uncertain setting.

Acknowledgments

This work has been written within the activities of GNFM (Gruppo Nazionale per la Fisica Matematica) of INdAM (Istituto Nazionale di Alta Matematica), Italy. The author acknowledges partial support from the Excellence Project CUP: E11G18000350001 of the Department of Mathematical Sciences “G. L. Lagrange”, Politecnico di Torino, Italy.