Structure preserving schemes for nonlinear Fokker-Planck equations and applications

TL;DR

This paper develops second-order accurate numerical schemes for nonlinear Fokker-Planck equations that preserve key structural properties, enabling accurate long-term behavior simulation and applications in socio-economic and life sciences.

Contribution

It introduces structure-preserving, mesh-unrestricted numerical schemes for nonlinear Fokker-Planck equations with nonlocal terms, improving accuracy and physical fidelity.

Findings

Schemes accurately capture asymptotic steady states.

Methods preserve non-negativity and entropy dissipation.

Applications demonstrate relevance to socio-economic models.

Abstract

In this paper we focus on the construction of numerical schemes for nonlinear Fokker-Planck equations that preserve the structural properties, like non negativity of the solution, entropy dissipation and large time behavior. The methods here developed are second order accurate, they do not require any restriction on the mesh size and are capable to capture the asymptotic steady states with arbitrary accuracy. These properties are essential for a correct description of the underlying physical problem. Applications of the schemes to several nonlinear Fokker-Planck equations with nonlocal terms describing emerging collective behavior in socio-economic and life sciences are presented.

Tables

Table 1. Table 1: Different expressions of the weights in ( 9 )

Scheme	${\delta }_{i+1/2}$	${\lambda }_{i+1/2}$
SP-CC	$\frac{1}{{\lambda }_{i+1/2}}+\frac{1}{1-\exp \left({\lambda }_{i+1/2}\right)}$	${\int }_{w_i}^{w_{i+1}}\frac{ℬ\left[f\right](w,t)+D^{\prime }\left(w\right)}{D\left(w\right)}𝑑w$
SP-CC2 (midpoint)	$\frac{1}{{\lambda }_{i+1/2}}+\frac{1}{1-\exp \left({\lambda }_{i+1/2}\right)}$	$\frac{\Delta w\left({ℬ}_{i+1/2}+D_{i+1/2}^{\prime }\right)}{D_{i+1/2}}$
SP-CCE (exact)	$\frac{1}{\log \left(f_i^{\infty }\right)-\log \left(f_{i+1}^{\infty }\right)}+\frac{f_{i+1}^{\infty }}{f_{i+1}^{\infty }-f_i^{\infty }}$	$\log \left(\frac{f_i^{\infty }}{f_{i+1}^{\infty }}\right)$

Table 2. Table 2: Example 1. Estimation of the order of convergence toward the reference stationary state for SP–CC scheme with explicit Euler (left) and RK4 (right) methods. Rates have been computed using N = 40 , 80 𝑁 40 80 N=40,80 and a reference solution with N = 640 𝑁 640 N=640 , σ 2 / 2 = 0.1 superscript 𝜎 2 2 0.1 \sigma^{2}/2=0.1 , Δ ​ t = Δ ​ w 2 / 2 ​ σ 2 Δ 𝑡 Δ superscript 𝑤 2 2 superscript 𝜎 2 \Delta t=\Delta w^{2}/2\sigma^{2} .

	$SP-C{C}_k$				$SP-C{C}_k$
Time	2	4	6	G	2	4	6	G
1	1.9456	1.9751	1.9740	1.9740	1.9470	1.9773	1.9762	1.9762
5	1.9700	3.2328	2.3690	2.3487	1.9700	3.2323	2.3724	2.3522
10	1.9695	3.9156	6.8498	7.3299	1.9695	3.9156	6.8517	7.3252
15	1.9695	3.9156	6.8715	7.3304	1.9695	3.9156	6.8761	7.3223

Table 3. Table 3: Example 1. Estimation of the order of convergence toward the reference stationary state for SP–CC scheme with first (left) and second order (right) semi–implicit methods. Rates have been computed using N = 40 , 80 𝑁 40 80 N=40,80 and a reference solution with N = 640 𝑁 640 N=640 , σ 2 / 2 = 0.1 superscript 𝜎 2 2 0.1 \sigma^{2}/2=0.1 , Δ ​ t = Δ ​ w / 2 ​ σ 2 Δ 𝑡 Δ 𝑤 2 superscript 𝜎 2 \Delta t=\Delta w/2\sigma^{2} .

	$SP-C{C}_k$				$SP-C{C}_k$
Time	2	4	6	G	2	4	6	G
1	1.0681	1.0648	1.0648	1.0648	2.0991	2.0931	2.0932	2.0932
5	2.0093	1.9760	1.9648	1.9648	2.0700	2.8755	2.6082	2.6092
10	2.0155	3.9714	3.4966	3.2983	2.0700	4.0006	6.0768	7.4096
15	2.0155	3.9776	5.3045	7.3283	2.0700	3.9982	5.8780	9.0173

Table 4. Table 4: Example 2. Estimation of the order of convergence toward the reference stationary state for the semi–implicit SP-CC scheme, N = 51 , 101 𝑁 51 101 N=51,101 , reference solution computed with N = 1601 𝑁 1601 N=1601 , σ 2 / 2 = 0.1 superscript 𝜎 2 2 0.1 \sigma^{2}/2=0.1 .

	$SP-C{C}_k$
Time	2	4	6	G
1	1.3047	1.5010	1.5021	1.5021
10	1.9893	4.0634	2.8122	2.8682
20	1.9894	3.9842	6.0784	10.0422

Table 5. Table 5: Example 3. Estimation of the order of convergence for the one-dimensional swarming model for the explicit SP–CC and SP–EA over the domain [ − L , L ] 𝐿 𝐿 [-L,L] with L = 5 𝐿 5 L=5 , N = 21 , 41 , 81 𝑁 21 41 81 N=21,41,81 , D = 0.4 𝐷 0.4 D=0.4 , Δ ​ t = Δ ​ w 2 / L 2 Δ 𝑡 Δ superscript 𝑤 2 superscript 𝐿 2 \Delta t=\Delta w^{2}/L^{2} .

$\alpha =0$	$SP-C{C}_k$				$SP-E{A}_k$
Time	2	4	6	G	2	4	6	G
1	2.1387	2.1387	2.1387	2.1387	2.4142	2.4142	2.4142	2.4142
5	6.9430	6.9430	6.9430	6.9430	10.0712	10.0712	10.0712	10.0712
10	20.0127	20.0127	20.0127	20.0127	23.9838	23.9838	23.9838	23.9838
$\alpha =1$	$SP-C{C}_k$				$SP-E{A}_k$
Time	2	4	6	G	2	4	6	G
1	2.5310	2.5310	2.5310	2.5310	2.2614	2.2892	2.2892	2.2892
5	2.0498	7.6659	7.6659	7.6659	2.0635	10.9818	10.9818	10.9818
10	2.0503	18.7697	18.7697	18.7697	2.0613	14.8321	14.8321	14.8321

Equations

$$\begin{cases}\partial_{t}f(w,t)=\nabla_{w}\cdot\Big{[}\mathcal{B}[f](w,t)f(w,t)+\nabla_{w}(D(w)f(w,t))\Big{]},\\ f(w,0)=f_{0}(w),\end{cases}$$

$$\mathcal{B}[f](w,t)=S(w)+{\int }_{{\mathbb{R}}^d}P(w,w_{\ast })(w-w_{\ast })f(w_{\ast },t)d{w}_{\ast },$$

$${\partial }_{t}f(w,t)={\nabla }_w\cdot [(\mathcal{B}[f](w,t)+{\nabla }_{w}D(w))f(w,t)+D(w){\nabla }_{w}f(w,t)].$$

$$\mathcal{F}[f](w,t)=(\mathcal{B}[f](w,t)+{\nabla }_{w}D(w))f(w,t)+D(w){\nabla }_{w}f(w,t),$$

$${\partial }_{t}f(w,t)={\nabla }_w\cdot \mathcal{F}(w,t).$$

$${\partial }_{t}f(w,t)={\partial }_w\mathcal{F}[f](w,t),w\in \Omega \subseteq \mathbb{R},$$

$$\mathcal{F}[f](w,t)=(\mathcal{B}[f](w,t)+D^{\prime }(w))f(w,t)+D(w){\partial }_{w}f(w,t)$$

$$\frac{d}{dt}{f}_i(t)=\frac{{\mathcal{F}}_{i+1/2}(t)-{\mathcal{F}}_{i-1/2}(t)}{\Delta w},$$

$$\begin{array}{r}{\mathcal{F}}_{i+1/2}={\tilde{\mathcal{C}}}_{i+1/2}{\tilde{f}}_{i+1/2}+D_{i+1/2}\frac{f_{i+1}-f_i}{\Delta w},\end{array}$$

$${\tilde{f}}_{i+1/2}=(1-{\delta }_{i+1/2})f_{i+1}+{\delta }_{i+1/2}{f}_i.$$

$${\delta }_{i+1/2}=1/2,{\tilde{\mathcal{C}}}_{i+1/2}=\mathcal{C}[f](w_{i+1/2},t).$$

$$\frac{f_{i+1}}{f_i}=\frac{-{\delta }_{i+1/2}{\tilde{\mathcal{C}}}_{i+1/2}+\frac{D_{i+1/2}}{\Delta w}}{(1-{\delta }_{i+1/2}){\tilde{\mathcal{C}}}_{i+1/2}+\frac{D_{i+1/2}}{\Delta w}}.$$

$$D(w){\partial }_{w}f(w,t)=-(\mathcal{B}[f](w,t)+D^{\prime }(w))f(w,t),$$

$${\int }_{w_i}^{w_{i+1}}\frac{1}{f(w,t)}{\partial }_{w}f(w,t)dw=-{\int }_{w_i}^{w_{i+1}}\frac{1}{D(w)}(\mathcal{B}[f](w,t)+D^{\prime }(w))dw,$$

$$\frac{f(w_{i+1},t)}{f(w_i,t)}=\exp \left\{-{\int }_{w_i}^{w_{i+1}}\frac{1}{D(w)}(\mathcal{B}[f](w,t)+D^{\prime }(w))dw\right\}.$$

$${\tilde{\mathcal{C}}}_{i+1/2}=\frac{D_{i+1/2}}{\Delta w}{\int }_{w_i}^{w_{i+1}}\frac{\mathcal{B}[f](w,t)+D^{\prime }(w)}{D(w)}dw$$

$${\delta }_{i+1/2}=\frac{1}{{\lambda }_{i+1/2}}+\frac{1}{1-\exp ({\lambda }_{i+1/2})},$$

$${\lambda }_{i+1/2}={\int }_{w_i}^{w_{i+1}}\frac{\mathcal{B}[f](w,t)+D^{\prime }(w)}{D(w)}dw=\frac{\Delta w{\tilde{\mathcal{C}}}_{i+1/2}}{D_{i+1/2}}.$$

$${\int }_{w_i}^{w_{i+1}}\frac{\mathcal{B}[f](w,t)+D^{\prime }(w)}{D(w)}dw=\frac{\Delta w({\mathcal{B}}_{i+1/2}+D_{i+1/2}^{\prime })}{D_{i+1/2}}+O(\Delta w^3),$$

$${\lambda }_{i+1/2}^{\text{mid}}=\frac{\Delta w({\mathcal{B}}_{i+1/2}+D_{i+1/2}^{\prime })}{D_{i+1/2}}$$

$${\delta }_{i+1/2}^{\text{mid}}=\frac{D_{i+1/2}}{\Delta w({\mathcal{B}}_{i+1/2}+D_{i+1/2}^{\prime })}+\frac{1}{1-\exp ({\lambda }_{i+1/2}^{\text{mid}})}.$$

$$\delta_{i+1/2}=\left\{\begin{array}[]{cc}0,&\mathcal{B}_{i+1/2}>0,\\ 1,&\mathcal{B}_{i+1/2}<0\\ \end{array}\right.$$

$$D(w){\partial }_w{f}^{\infty }(w)=-(\mathcal{B}(w)+D^{\prime }(w))f^{\infty }(w),$$

$$\frac{f_{i+1}^{\infty }}{f_i^{\infty }}=\exp \left\{-{\int }_{w_i}^{w_{i+1}}\frac{1}{D(w)}(\mathcal{B}(w)+D^{\prime }(w))dw\right\}=\exp \left(-{\lambda }_{i+1/2}^{\infty }\right),$$

$${\lambda }_{i+1/2}^{\infty }=\log \left(\frac{f_i^{\infty }}{f_{i+1}^{\infty }}\right)$$

$${\delta }_{i+1/2}^{\infty }=\frac{1}{\log (f_i^{\infty })-\log (f_{i+1}^{\infty })}+\frac{f_{i+1}^{\infty }}{f_{i+1}^{\infty }-f_i^{\infty }}.$$

$$\sum _{i=0}^N\frac{d}{dt}{f}_i(t)=0,\forall t>0.$$

$$f_i^{n+1}=f_i^n+\Delta t\frac{{\mathcal{F}}_{i+1/2}^n-{\mathcal{F}}_{i-1/2}^n}{\Delta w}.$$

$$\Delta t\le \frac{\Delta w^2}{2(M\Delta w+D)},M=\underset{i}{\max }|{\tilde{\mathcal{C}}}_{i+1/2}^n|,D=\underset{i}{\max }{D}_{i+1/2},$$

$$\begin{split}&f_{i}^{n+1}=f_{i}^{n}+\dfrac{\Delta t}{\Delta w}\Bigg{[}\left((1-\delta_{i+1/2}^{n})\tilde{\mathcal{C}}_{i+1/2}^{n}+\dfrac{D_{i+1/2}}{\Delta w}\right)f_{i+1}^{n}\\ &+\left(\tilde{\mathcal{C}}_{i+1/2}^{n}\delta_{i+1/2}^{n}-\tilde{\mathcal{C}}_{i-1/2}^{n}(1-\delta_{i-1/2}^{n})-\dfrac{1}{\Delta w}(D_{i+1/2}+D_{i-1/2})\right)f_{i}^{n}\\ &-\left(\tilde{\mathcal{C}}^{n}_{i-1/2}\delta_{i-1/2}^{n}-\dfrac{D_{i-1/2}}{\Delta w}\right)f_{i-1}^{n}\Bigg{]}.\end{split}$$

Full text

**Keywords: structure preserving methods, finite difference schemes, Fokker-Planck equations, emerging collective behavior.

**

1 Introduction

In this paper we construct and discuss a steady-state preserving method for a wide class of nonlinear Fokker-Planck equations of the form

[TABLE]

where $t\geq 0$, $w\in\Omega\subseteq\mathbb{R}^{d}$, $d\geq 1$, $f(w,t)\geq 0$ is the unknown distribution function, $\mathcal{B}[\cdot]$ is a bounded operator which describes aggregation dynamics and $D(\cdot)\geq 0$ is a diffusion function.

A typical example is given by mean-field models of collective behavior where the nonlocal operator $\mathcal{B}[\cdot]$ has the form

[TABLE]

with $P:\mathbb{R}^{d\times d}\rightarrow\mathbb{R}^{+}$ and $S:\mathbb{R}^{d}\rightarrow\mathbb{R}^{d}$. With the choice (2) equation (1) describes typical features of the collective behavior in multiagent systems with nonlocal type interactions. These models of collective behavior has been extensively discussed in the last decades at the particle, kinetic and hydrodynamic level [2, 3, 4, 5, 11, 12, 13, 17, 19, 30]. In particular, many heterogeneous phenomena like swarming behaviors, human crowds motion and formation of wealth distributions are described by these type of PDEs under special assumptions. We refer to [26, 27], and the references therein, for a recent overview of such models.

In the following, we focus on the construction of numerical methods for such problems which are able to preserve the structural properties of the PDE, like non negativity of the solution, entropy dissipation and large time behavior. The methods here developed are second order accurate, they do not require any restriction on the mesh size and are capable to capture the asymptotic steady states with arbitrary accuracy. These properties are essential for a correct description of the underlying physical problem.

The derivation of the schemes follows the main lines of the seminal work of Chang–Cooper for the linear Fokker-Planck equation [9, 16, 25, 29]. However, in the nonlinear case, the exact stationary solution is unknown and a more advanced treatment is needed in order to find a good approximation to the problem. Similar approaches for nonlinear Fokker-Planck equations were previously derived in [8, 24]. Related methods for the case of nonlinear degenerate diffusions equations were proposed in [6, 15] and with nonlocal terms in [10, 11]. We refer also to [1] for the development of methods based on stochastic approximations and to [21] for a recent survey on schemes which preserve steady states of balance laws and related problems.

Although we derive the schemes in the case of Fokker-Planck equations, the methods can be easily applied to more general problems where the solution depends on additional parameters and the PDE is of Vlasov-Fokker-Planck type. In this case, preservation of the steady states is of fundamental importance in order to develop asymptotic-preserving methods [18].

The rest of the paper is organized as follows. In the next Section we first derive the Chang-Cooper type schemes in one-dimension with a particular attention to the steady state preserving properties. We then prove non negativity of solutions for explicit and semi-implicit schemes and entropy inequality for a class of one dimensional Fokker-Planck models. In Section 3 we introduce a modification of the schemes based on a more general entropy dissipation principle. We show that these entropic schemes preserve stationary solutions and derive sufficient conditions for non negativity of explicit and semi-implicit schemes. Several applications of the schemes are finally presented in Section 4 for various nonlinear Fokker-Planck problems describing collective behaviors in socio-economic and life sciences. Some conclusions are reported at the end of the manuscript. Extension to high order semi-implicit methods and the multi-dimensional case are considered in separate appendices.

2 Chang-Cooper type schemes

In the following we focus on the design of numerical schemes for (1) which we rewrite in divergence form as

[TABLE]

We can define the $d$-dimensional flux

[TABLE]

therefore (3) reads

[TABLE]

2.1 Derivation of the schemes

In the one-dimensional case $d=1$ equation (4) becomes

[TABLE]

where now

[TABLE]

using the compact notation $D^{\prime}(w)=\partial_{w}D(w)$. Typically, when ${\Omega}$ is a finite size set the problem is complemented with no-flux boundary conditions at the extremal points. In the sequel we assume $D(w)\geq\alpha>0$ in the internal points of $\Omega$.

We introduce an uniform spatial grid $w_{i}\in\Omega$, such that $w_{i+1}-w_{i}=\Delta w$. We denote as usual $w_{i\pm 1/2}=w_{i}\pm\Delta w/2$ and consider the conservative discretization

[TABLE]

where for each $t\geq 0$, $f_{i}(t)$ is an approximation of $f(w_{i},t)$ and $\mathcal{F}_{i\pm 1/2}[f](t)$ is the numerical flux function characterizing the discretization.

Let us set ${\mathcal{C}}[f](w,t)=\mathcal{B}[f](w,t)+D^{\prime}(w)$ and adopt the notations $D_{i+1/2}=D(w_{i+1/2})$, $D^{\prime}_{i+1/2}=D^{\prime}(w_{i+1/2})$. We will consider a general flux function which is combination of the grid points $i+1$ and $i$

[TABLE]

where

[TABLE]

Here, we aim at deriving suitable expressions for $\delta_{i+1/2}$ and $\tilde{\mathcal{C}}_{i+1/2}$ in such a way that the method yields nonnegative solutions, without restriction on $\Delta w$, and preserves the steady state of the system with arbitrary order of accuracy.

For example, the standard approach based on central difference is obtained by considering for all $i$ the quantities

[TABLE]

It is well-known, however, that such a discretization method is subject to restrictive conditions over the mesh size $\Delta w$ in order to keep non negativity of the solution.

First, observe that when the numerical flux vanishes from (8) we get

[TABLE]

Similarly, if we consider the analytical flux imposing ${\cal F}[f](w,t)\equiv 0$, 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 (10) on the cell $[w_{i},w_{i+1}]$

[TABLE]

which gives

[TABLE]

The terminology quasi steady-state was introduced originally by Chang and Cooper in [16], the ”quasi” coming from the fact that, with general time-dependent coefficients, no time-stabilization can be expected.

Now, by equating the ratio $f_{i+1}/f_{i}$ and $f(w_{i+1},t)/f(w_{i},t)$ of the numerical and exact flux, and setting

[TABLE]

we recover

[TABLE]

where

[TABLE]

We can state the following

Proposition 1**.**

The numerical flux function (8)-(9) with $\tilde{\mathcal{C}}_{i+1/2}$ and $\delta_{i+1/2}$ defined by (12) and (13)-(14) vanishes when the corresponding flux (6) is equal to zero over the cell $[w_{i},w_{i+1}]$. Moreover the nonlinear weight functions $\delta_{i+1/2}$ defined by (13)-(14) are such that $\delta_{i+1/2}\in(0,1)$.

The latter result follows from the simple inequality $\exp(x)\geq 1+x$. We refer to this type of schemes as structure preserving Chang-Cooper (SP-CC) type schemes.

By discretizing (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 (12). We refer to Section 4 for examples and more details. For linear problems of the form $\mathcal{B}[f](w,t)=\mathcal{B}(w)$ with constant diffusion $D^{\prime}=0$, the above scheme (15)-(16) is usually referred to as the Chang-Cooper method [16, 25]. In particular, if $\mathcal{B}(w)$ is a first order polynomial in $w$ as in [9] the midpoint rule is equivalent to the exact evaluation of the integral (12).

Remark 1*.*

• •

If we consider the limit case $D_{i+1/2}\to 0$, $D^{\prime}_{i+1/2}\to 0$ in (15)-(16) we obtain the weights

[TABLE]

and the scheme locally reduces to a first order upwind scheme for the corresponding continuity equation.

• •

For linear problems of the form $\mathcal{B}[f](w,t)=\mathcal{B}(w)$ the exact stationary state $f^{\infty}(w)$ can be directly computed from the solution of

[TABLE]

together with the boundary conditions. Explicit examples of stationary states will be reported in Section 4.

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. Finally, in Table 1 we summarize the different expressions of the weight functions.

2.2 Main properties

In order to study the structural properties of the numerical schemes, like conservations, non negativity and entropy property, we restrict to the one-dimensional case. To start with we consider the following simple result.

Lemma 1**.**

Let us consider the scheme (7)-(8) for $i=0,\ldots,N$ with no flux boundary conditions $\mathcal{F}_{N+1/2}=\mathcal{F}_{-1/2}=0$. We have

[TABLE]

The proof is a simple consequence of the telescopic summation property and the no flux boundary conditions.

2.2.1 Positivity preservation

Concerning non negativity, first we prove a result for the explicit scheme. We introduce a time discretization $t^{n}=n\Delta t$ with $\Delta t>0$ and $n=0,\dots,T$ and consider the simple forward Euler method

[TABLE]

Proposition 2**.**

Under the time step restriction

[TABLE]

the explicit scheme (21) with flux defined by (13)-(14) preserves nonnegativity, i.e $f^{n+1}_{i}\geq 0$ if $f^{n}_{i}\geq 0$, $i=0,\dots,N$.

Proof.

The scheme reads

[TABLE]

From (LABEL:eq:prop1_f) we have a convex combination if the coefficients of $f_{i+1}^{n}$ and $f_{i-1}^{n}$ satisfy

[TABLE]

or equivalently

[TABLE]

which holds true thanks to the properties of the exponential function. In order to ensure the non negativity of the scheme the time step should satisfy the restriction $\Delta t\leq{\Delta w}/{\nu}$, with

[TABLE]

Being $M$ and $D$ defined in (22), and $0\leq\delta_{i\pm 1/2}\leq 1$, we obtain the prescribed bound. ∎

Remark 2*.*

Higher order SSP methods [22] 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 $\Delta t=O(\Delta w^{2})$ of explicit schemes. Unfortunately, fully implicit methods originate a nonlinear system of equations due to the nonlinearity of $\mathcal{B}[f]$ and the dependence of the weights $\delta_{i\pm 1/2}$ from the solution. However, we can prove that nonnegativity of the solution holds true also for the semi-implicit case

[TABLE]

where

[TABLE]

We have

Proposition 3**.**

Under the time step restriction

[TABLE]

the semi-implicit scheme (24) preserves nonnegativity, i.e $f^{n+1}_{i}\geq 0$ if $f^{n}_{i}\geq 0$, $i=0,\dots,N$.

Proof.

Equation (24) corresponds to

[TABLE]

thanks to the definition of the flux function introduced in (8)-(9). Using the indentity $\lambda_{i+1/2}^{n}=\Delta w{\tilde{\mathcal{C}}_{i+1/2}^{n}}/{D_{i+1/2}}$ we obtain

[TABLE]

Let us denote $\alpha_{i+1/2}^{n}=\dfrac{\lambda_{i+1/2}^{n}}{\exp(\lambda_{i+1/2}^{n})-1}\geq 0$ and

[TABLE]

we can write

[TABLE]

If we introduce the matrix

[TABLE]

with $R_{i}^{n}>0$, $Q_{i}^{n}\geq 0$, $P_{i}^{n}\geq 0$ defined in (LABEL:eq:RQP) the semi-implicit scheme may be expressed in matrix form as follows

[TABLE]

with $\textbf{f}^{n}=\left(f_{0}^{n},\dots,f_{N}^{n}\right)$. Since $\textbf{f}^{n}\geq 0$, in order to prove that $\textbf{f}^{n+1}\geq 0$ it is sufficient to show $\mathcal{A}[f^{n}]^{-1}\geq 0$. Now, $\mathcal{A}[\cdot]$ is a tridiagonal matrix with positive diagonal elements and if $\mathcal{A}$ is strictly diagonally dominant we can conclude that $\mathcal{A}^{-1}\geq 0$.

The matrix $\mathcal{A}$ is strictly diagonally dominant if and only if

[TABLE]

condition which holds true if

[TABLE]

∎

Remark 3*.*

• •

Higher order semi-implicit approximations can be constructed following [7]. An example of second order semi-implicit IMEX scheme is given in Appendix A. We mention also [25] where a second order semi-implicit BDF method has been considered. Note, however, that the determination of nonnegative semi-implicit schemes with optimal stability regions is an open problem which goes beyond the purpose of the present manuscript.

• •

Although fully implicit schemes originate a nonlinear system of equations, we remark that the same argument used in Proposition 3 permits to prove nonnegativity of the scheme even with the fully implicit fluxes

[TABLE]

with

[TABLE]

In fact, we obtain the nonlinear system

[TABLE]

where the matrix $\mathcal{A}[{\bf f}^{n+1}]$ has the same structure (29) with the entries evaluated at time $n+1$. The above system can be solved iteratively at each time step

[TABLE]

where now each iteration step is explicit and can be made non negative under a stability restriction analogous to (26). Therefore, if ${\bf{f}}_{k}^{n+1}\rightarrow{\bf{f}}^{n+1}$ as ${k\rightarrow+\infty}$ we can infer the nonnegativity of the scheme under the condition (32), being $\mathcal{A}[{\bf{f}}^{n+1}]\geq 0$ strictly diagonally dominant and then $\mathcal{A}[{\bf{f}}^{n+1}]^{-1}\geq 0$. In general, the convergence properties of the above iterative method depend on the nonlinear flux function $\mathcal{B}[f]$ which defines the coefficients in $\mathcal{A}[{\bf{f}}]$. For example, since by nonnegativity and mass conservation $\|f_{k+1}^{n+1}\|_{1}=\|f_{k}^{n+1}\|_{1}$, $k\geq 0$, where $\|\cdot\|_{1}$ denotes the discrete $L_{1}$ norm, convergence is guaranteed for any choice of the initial value $f^{n+1}_{0}$ if $\mathcal{A}^{-1}[{\bf{\cdot}}]{\bf{f}}^{n}$ is a contraction [23].

2.2.2 Relative entropy dissipation

In order to discuss the entropy property we consider the prototype equation

[TABLE]

with $-1<u<1$ a given constant and boundary conditions

[TABLE]

If the stationary state $f^{\infty}$ exists equation (35) may be written in the Landau form as

[TABLE]

or in the non logarithmic Landau form as

[TABLE]

We define the relative entropy for all positive functions $f(w,t),g(w,t)$ as follows

[TABLE]

we have [20]

[TABLE]

where the dissipation functional $\mathcal{I}_{D}(\cdot,\cdot)$ is defined as

[TABLE]

Of course we might consider other entropies like the $L^{2}$ entropy which is defined as

[TABLE]

with

[TABLE]

see [20] for further examples.

Note that, since the definition of relative entropy implies the existence of a steady state $f^{\infty}$ then the considerations of the second part of Remark 1 apply. Therefore, the weights in the Chang-Cooper type method can be evaluated exactly and are given by (19)-(20). This property is used to prove the following results.

Lemma 2**.**

In the case $\mathcal{B}[f](w,t)=\mathcal{B}(w)$ the numerical flux function (8)-(9) with $\tilde{\mathcal{C}}_{i+1/2}$ and $\delta_{i+1/2}$ given by (12)-(13) can be written in the form (38) and reads

[TABLE]

with

[TABLE]

Proof.

In the hypothesis $\mathcal{B}[f](w,t)=\mathcal{B}(w)$ the definition of $\lambda_{i+1/2}$ does not depends on time, i.e. $\lambda_{i+1/2}=\lambda_{i+1/2}^{\infty}$ and if a steady state exists we may write

[TABLE]

Furthermore, the flux function $\mathcal{F}^{n}_{i+1/2}$ assumes the following form

[TABLE]

where

[TABLE]

Hence we have

[TABLE]

which gives (44). ∎

Theorem 1**.**

Let us consider $\mathcal{B}[f](w,t)=w-u$ as in equation (35). The numerical flux (8)-(9) with $\tilde{\mathcal{C}}_{i+1/2}$ and $\delta_{i+1/2}$ given by (12)-(13) satisfies the discrete entropy dissipation

[TABLE]

where

[TABLE]

and $I_{\Delta}$ is the positive discrete dissipation function

[TABLE]

Proof.

From the definition of relative entropy we have

[TABLE]

and after summation by parts we get

[TABLE]

Thanks to the identity of Lemma 2 we may conclude since the function $(x-y)\log(x/y)$ is non-negative for all $x,y\geq 0$. ∎

3 Entropic schemes for gradient flow problems

In this section we introduce a second class of structure preserving numerical scheme based on the entropy dissipation principle. To this aim, let us consider the general class of nonlinear Fokker-Planck equation with gradient flow structure [5, 11, 14]

[TABLE]

and no-flux boundary conditions. In the case of equation (1) with constant diffusion $D$ we have

[TABLE]

We focus on the following prototype of function $\xi(w,t)$, $w\in\mathbb{R}^{d}$

[TABLE]

which in our case corresponds to

[TABLE]

with $U=U(|w|)$ an interaction potential.

The corresponding free energy is given by

[TABLE]

Using the fact that $U$ is an even function we can write

[TABLE]

Hence, upon integration by parts we obtain the dissipation of the free energy $\mathcal{E}(t)$ along solutions

[TABLE]

where $\mathcal{I}(\cdot)$ is the entropy dissipation function.

3.1 Derivation of the schemes

In the one-dimensional case $d=1$ equation (52) reads

[TABLE]

where

[TABLE]

and we assume ${\mathcal{B}}[f](w,t)=\partial_{w}(U\ast f)(w,t)$.

In order to derive schemes which satisfy the entropy dissipation property (57) we consider the semi-discrete version of the entropy (55)

[TABLE]

Therefore, we have

[TABLE]

Now using the general discrete conservative formulation (7) and the fact that $\xi_{i}=U*f_{i}+D\log f_{i}$ we get

[TABLE]

Furthermore, after summation by parts we can write the last term as follows

[TABLE]

Now, integrating (53) over the cell $[w_{i},w_{i+1}]$ we obtain

[TABLE]

where

[TABLE]

Let us now consider a general scheme in the Chang-Cooper form (8), which in our case can be rewritten as

[TABLE]

with

[TABLE]

Therefore, we have

[TABLE]

Thus we cannot prove that the discrete entropy functional (60) is dissipated by the Chang-Cooper type scheme developed in the previous sections, unless $K_{j+1/2}\equiv 1$. This latter requirement is satisfied if we consider the new entropic flux function

[TABLE]

We will refer to the above approximation of the solution at the grid point $i+1/2$ as entropic average of the grid points $i$ and $i+1$. In the general case of the flux function (6) with non constant diffusion the resulting numerical flux reads

[TABLE]

Finally, concerning the stationary state, we obtain immediately imposing the numerical flux equal to zero

[TABLE]

and therefore we get

[TABLE]

By equating the above ratio with the quasi-stationary approximation (11) we get the same expression for $\tilde{\mathcal{C}}_{i+1/2}$ as in (12)

[TABLE]

We can state the following

Proposition 4**.**

The numerical flux function (67) with $\tilde{\mathcal{C}}_{i+1/2}$ defined by (68) vanishes when the corresponding flux (59) is equal to zero over the cell $[w_{i},w_{i+1}]$.

3.2 Main properties

A fundamental result concerning the entropic average (66) is the following Lemma.

Lemma 3**.**

The entropy average defined in (66) may be written as a convex combination with nonlinear weights

[TABLE]

where

[TABLE]

Proof.

From (70) we have

[TABLE]

that is (67). It is a easy computation to verify that $\delta^{E}_{i+1/2}$ lies in the interval $(0,1)$. ∎

Remark 4*.*

As a consequence the Chang-Cooper type average (9) and the entropic average (66) define the same quantity at the steady state when $f_{i}=f_{i}^{\infty}$. In fact, in this case (70) are the same as (20).

We can summarize our findings of Section 3.1 as follows.

Theorem 2**.**

The numerical flux (67)-(66) for a constant diffusion $D$ satisfies the discrete entropy dissipation

[TABLE]

where $\mathcal{E}_{\Delta}$ is given by (60) and $I_{\Delta}$ is the discrete entropy dissipation function

[TABLE]

with $\xi_{j+1}-\xi_{j}$ defined as in (62).

Remark 5*.*

On the contrary to the Chang-Cooper average the restrictions for the non negativity property of the solution are stronger. In fact, by the same arguments we used in the previous section, non negativity of the explicit scheme requires

[TABLE]

However, the weights do not possess any special structure that permits to avoid a constraint of the mesh size $\Delta w$ which now must satisfy

[TABLE]

Therefore, similar to central differences, we have a restriction on the mesh size which becomes prohibitive for small values of the diffusion function $D(w)$. It is easy to verify that the same condition is necessary also for the non negativity of semi-implicit approximations.

3.2.1 The case ${\mathcal{B}}[f](w,t)=B(w)$

Next we consider the case of linear flux ${\mathcal{B}}[f](w,t)=B(w)$. The following Lemma holds true.

Lemma 4**.**

In the case $\mathcal{B}[f](w,t)=B(w)$ the numerical flux (67)-(66) corresponds to the form (37) and reads

[TABLE]

Proof.

If a stationary $f^{\infty}(w)$ state exists it nullify the flux and we have

[TABLE]

From the definition of the entropic flux (67) we obtain

[TABLE]

from which we conclude. ∎

We can now state the following entropy dissipation results for problem (35) in the nonlogarithmic Landau form (38).

Theorem 3**.**

Let us consider $\mathcal{B}[f](w,t)=w-u$ as in equation (35). The numerical flux (67)-(66) with $\tilde{\mathcal{C}}_{i+1/2}$ given by (12) satisfies the discrete entropy dissipation

[TABLE]

where $\mathcal{H}_{\Delta w}(f,f^{\infty})$ is given by (49) and $I^{E}_{\Delta}$ is the positive discrete dissipation function

[TABLE]

Proof.

From

[TABLE]

and being

[TABLE]

we have

[TABLE]

∎

4 Applications

In this section we present several numerical examples of Fokker-Planck equations solved with the structure-preserving schemes here introduced. An essential aspect for the accurate description of the steady state is the approximation of the integral defining the quasi-stationary solution

[TABLE]

Except for simple linear cases, a suitable quadrature formula is required. In the following numerical examples we consider open Newton-Cotes formulas up to order 6 and Gauss-Legendre quadrature with $6$ points. In the sequel, we will adopt the notation SP–CCk and SP–EAk, $k=2,4,6,G$ to denote the structure preserving schemes with Chang-Cooper (CC) and entropic average (EA) flux and when (78) is approximated with second, fourth, sixth order Newton–Cotes quadrature or Gaussian quadrature, respectively. Singularities at the boundaries in the integration of (78) can be avoided using open Newton–Cotes rules.

4.1 Example 1: Opinion dynamics in bounded domains

Let us consider the evolution of a distribution function described by (1), with $w\in I$, where $I=[-1,1]$, and

[TABLE]

The model describes the evolution of the distribution functions of agents having opinion $w$ at time $t$ (see [27, 30] for more details).

In the simplified case $P\equiv 1$ the corresponding stationary distribution reads

[TABLE]

with $\sigma\in\mathbb{R}$ a given parameter, $C>0$ is a normalization constant and $u=\int_{I}wf(w,t)dw,$.

We consider as initial distribution

[TABLE]

with $\beta>0$ a normalization constant. Since diffusion vanishes at the boundaries we present results for the Chang-Cooper type numerical schemes SP–CC only.

In Figure 1 we compute the relative $L^{1}$ error of the numerical solution with respect to the exact (80) stationary state using $N=80$ points for the SP–CC scheme with various quadrature rules. It is possible to observe how the different integration methods capture the steady state with different accuracy. In particular low order quadrature rules achieve the numerical steady state faster than high order quadratures, with the Gaussian quadrature that essentially reach machine precision. In the same figure we illustrate how SP–CC scheme dissipates the relative entropy (49) in the case of two coarse grids with $N=10$ and $N=20$ points.

In Table 2 we estimate the overall order of convergence of explicit SP–CC schemes for several integration methods. Here we used $N=40,80$, with reference solutions computed with $N=640$ points. For comparison, the time integration has been performed with the explicit Euler and RK4 methods and the time step chosen in such a way that the CFL condition for the positivity of the scheme is satisfied, i.e. $\Delta t=O(\Delta w^{2})$. As expected the two methods are essentially equivalent and both schemes are second order accurate in the transient regimes and assume the order of the quadrature method close to the steady state.

In Table 3 we estimate the order of convergence of SP–CC schemes with first and second order semi–implicit methods (see Appendix A for a detailed description of the methods). In this case the CFL condition is $\Delta t=O(\Delta w)$ and a second order method is necessary to achieve second order accuracy in the transient regime, whereas both schemes achieve higher order accordingly to the quadrature used for large times.

In the general case $P(w,w_{*})\neq 1$ and it is not possible to give an analytical formulation of the steady state solution $f^{\infty}(w,t)$. In Figure 2 we represent a typical evolution of an aggregation model in the bounded confidence case [27]

[TABLE]

where $\chi(\cdot)$ is the indicator function, for $\Delta=0.4$, $\Delta=0.8$. Here, the evolution has been computed through a SP–CC with Gauss quadrature, the integral $\mathcal{B}[f](w,t)$ has been evaluated through a trapezoidal method.

4.2 Example 2: Wealth evolution in unbounded domains

Let us consider equation (1) with $w\in\mathbb{R}^{+}$ and

[TABLE]

With the above choice, the Fokker-Planck equation describes the evolution of the wealth distribution $w$ at time $t$ in a large set of interacting economic agents (see [26, 27] for details).

In the case of constant interaction $a(w,w_{*})\equiv 1$ the steady state of the equation is analytically computable

[TABLE]

where $\mu=1+2/\sigma^{2}$ is the so-called Pareto exponent. In the numerical test we consider the initial distribution

[TABLE]

with $\beta>0$ a normalization constant.

Again, due to degeneracy of the diffusion on the left boundary we report results only for SP–CC schemes. In Figure 3 we present the solution with $u=1$ in the domain $[0,L]$, $L=10$. In both figures $a(\cdot,\cdot)=1$ whereas the diffusion constant assumes different values. We report the evolution of the solution and the relative $L^{1}$ error with respect to the stationary state using $N=201$ points for the semi–implicit SP–CC scheme (SISP–CC). We observe how the introduced methods describe the stationary state with different levels of accuracy. Note that, at the right boundary we must introduce an artificial boundary condition in order to truncate the computational domain. In our numerical results we impose the quasi stationary condition (11) in order to evaluate $f_{N+1}(t)$, that is

[TABLE]

In Table 4 we estimate the overall order of convergence of the semi–implicit SP–CC scheme for several integration methods with $N=51,101$ for the domain $[0,L]$, $L=10$, with reference solutions computed with $N=1601$ gridpoints. The time step is chosen in such a way that the CFL condition for the positivity of the scheme is satisfied, i.e. $\Delta t=O(\Delta w)$. We can observe that for short times the order of accuracy is limited by the semi–implicit method, which is first order accurate, whereas as we approach to the stationary solution the order depends on the quadrature formula used.

4.3 Example 3: 2D model of swarming

Let us consider a self-propelled swarming model of Cucker-Smale type [5] with diffusion. In this model the evolving distribution $f(x,w,t)$ represents the density of individuals (birds, fishes, $\ldots$) in position $x\in\mathbb{R}^{d}$ having velocity $w\in\mathbb{R}^{d}$ at time $t>0$. We have the following dynamics

[TABLE]

with

[TABLE]

and $K(w,y)>0$ a localization kernel, $\alpha>0$ a self-propulsion term and $D>0$ a constant noise intensity.

The space homogeneous version of the model (86) may be formulated in terms of the nonlinear Fokker–Planck equation (1) with

[TABLE]

with $\alpha$ a positive constant and $P(w,w_{*})\equiv 1$. The above equation can be written as a gradient flow. In fact, if we define

[TABLE]

with $U(w)$ a Coloumb potential and $\Phi(w)$ a confining potential given by

[TABLE]

the equation reads

[TABLE]

A free energy functional which dissipates along solutions is defined by

[TABLE]

with

[TABLE]

Stationary solutions should satisfy the identity $\nabla_{w}\xi=0$ and have the form

[TABLE]

with $C>0$ a normalization constant. It is possible to prove the following result (see [5] for more details).

Theorem 4**.**

Let us consider equation (86) in the space-homogeneous case, i.e. (1) with $\mathcal{B}[f](w,t)$ and diffusion as in (88), exhibits a phase transition in the following sense

• i)

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

• ii)

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

Since diffusion is constant, we compute the solution both using SP–CC type schemes and the entropic average schemes SP–EA. In Table 5 we estimate the order of convergence of the SP–CC and SP–EA schemes in the 1D case for several integration methods. We can observe how each method reach spectral accuracy in the case $\alpha=0$, in this case, in fact, all quadrature methods become exact being the quantity $(\mathcal{B}[f]+D^{\prime})/D$ a first order polynomial in $w$.

In Figure 4 we show that, as expected, on a coarse grid the SP–EA method becomes unstable for vanishing diffusions, whereas the SP–CC scheme remains stable and reduces to first order upwinding. In this case the solution becomes close to a Dirac delta in the velocity space. Finally, in Figure 5 we present the resulting 2D nonlinear Fokker–Planck equation for swarming with $\mathcal{B}[f](w,t)$ and $D(w)$ in (88), for several values of the diffusion coefficient $D=0.1,0.3,0.5$, and of the self-propulsion $\alpha=0,2,4$. The initial distribution is a bivariate normal distribution of the form

[TABLE]

with $\mu_{w}=\mu_{v}=2$ and $\sigma_{w}^{2}=\sigma_{v}^{2}=0.5$. The generalization of the schemes to the multidimensional case is done dimension by dimension and is summarized in Appendix B. The semi-implicit numerical scheme has been used, with a $6$th order open Newton–Cotes quadrature method. It is possible to observe the threshold phenomenon occurring for an increasing diffusion prescribed by Theorem 4. The results obtained with the two different schemes are essentially equivalent in this case.

Conclusion

The construction of structure–preserving schemes for nonlinear Fokker–Planck equations has been studied. Two different types of schemes have been constructed. The first type represents a natural extension of the so–called Chang–Cooper scheme to the nonlinear case. The second type of schemes represents a modification with better entropy dissipation properties. Both methods are second order accurate and capable to preserve the stationary state with arbitrary accuracy. However, non negativity restrictions are more severe for the second type of schemes. Even if the analysis is performed in the one-dimensional case, extensions to multidimensional situations are also considered. Several applications to linear and nonlinear Fokker-Planck equations arising in socio-economic sciences are presented and show the generality of the present approach. Extensions of the schemes to include nonlinear diffusion terms and higher order schemes in the limiting of vanishing diffusion are actually under study and will be presented elsewhere.

Acknowledgements

The research that led to the present paper 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 “Compagnia di San Paolo”.

Appendix A High order semi-implicit methods

Here we follow the approach in [7]. We write the semi-discrete scheme (7) in the equivalent form

[TABLE]

where ${\bf f}=(f_{0},\ldots,f_{N})$, ${\bf g}=(g_{0},\ldots,g_{N})$,

[TABLE]

and

[TABLE]

with initial conditions $f_{i}(0)=f_{0}(w_{i})$, $g_{i}(0)=f_{0}(w_{i})$. In the above equations we used the notation $[\cdot]$ to denote the functional dependence.

System (93) is then solved by an implicit-explicit (IMEX) Runge-Kutta method [28] where the variables ${f}_{i}$ are treated implicitly and the variables ${g}_{i}$ are treated explicitly. More precisely, using standard notations we can write an implicit-explicit Runge-Kutta scheme for (93) as follows. First we set $f_{i}^{n}=g_{i}^{n}$ and compute for $h=1,\ldots,s$

[TABLE]

where ${\bf F}^{k}=(F_{0}^{k},\ldots,F_{N}^{k})$, ${\bf G}^{k}=(G_{0}^{k},\ldots,G_{N}^{k})$ and next we update the numerical solution

[TABLE]

In particular, the IMEX scheme is chosen such that $b_{k}={\tilde{b}}_{k}$, $k=1,\ldots,s$ so that $f^{n+1}=g^{n+1}$ and therefore the duplication of the system is only apparent since there is only one set of numerical solutions. In our numerical tests we coupled the structure preserving discretizations with the second order semi-implicit scheme obtained as a combination of Heun method (explicit) and Crank-Nicolson (implicit) characterized by $s=2$ and

[TABLE]

As observed in [7] this represents a natural choice when dealing with convection-diffusion type equations, since the Heun method is an SSP explicit RK method, and Crank-Nicolson is an A-stable method, widely used for diffusion problems. Higher order methods can be found in [7, 28].

Appendix B The multi-dimensional case

In this section we report for the sake of completeness the details of the numerical schemes for multi-dimensional situations. We consider the case of Chang-Cooper type fluxes, and to keep notations simple we restrict to two dimensional problems $d=2$. We introduce a uniform mesh $(w_{i},v_{j})\in\Omega\subseteq\mathbb{R}^{2}$, with $\Delta w=w_{i+1}-w_{i}$ and $\Delta v=v_{j+1}-v_{j}$. We denote by $w_{i+1/2}=w_{i}+\Delta w/2$ and $v_{j+1/2}=v_{j}+\Delta v/2$. Let $f_{i,j}(t)$ be an approximation of the solution $f(w_{i},v_{j},t)$ and consider the following discretization of the nonlinear Fokker-Planck equation (4)

[TABLE]

being $\mathcal{F}_{i\pm 1/2,j}[f]$, $\mathcal{F}_{i,j\pm 1/2}[f]$ flux functions characterizing the numerical discretization. The quasi-stationary approximations over the cell $[w_{i},w_{i+1}]\times[v_{i},v_{i+1}]$ of the two dimensional problem now read

[TABLE]

Therefore, setting

[TABLE]

and by considering the natural generalization of the one-dimensional numerical flux

[TABLE]

we define $\delta_{i+1/2,j}$ and $\delta_{i,j+1/2}$ in such a way that we preserve the steady state solution for each dimension. The Chang-Cooper type structure preserving methods are then given by

[TABLE]

and

[TABLE]

The cases of higher dimension $d\geq 3$ and entropic average fluxes may be derived in a similar way.