A Rosenau-type approach to the approximation of the linear Fokker--Planck equation

{The numerical approximation of the solution of the Fokker--Planck equation is a challenging problem that has been extensively investigated starting from the pioneering paper of Chang and Cooper in 1970. We revisit this problem at the light of the approximation of the solution to the heat equation proposed by Rosenau in 1992. Further, by means of the same idea, we address the problem of a consistent approximation to higher-order linear diffusion equations.


1.
Introduction. The Fokker-Planck equation is a partial differential equation describing the time evolution of a density function f (v, t), where v ∈ IR d , d ≥ 1 and t ≥ 0, departing from a nonnegative initial density ϕ(v). The standard assumptions on ϕ(v) is that it possesses finite mass ρ, mean velocity u and temperature θ, where for any given density g(v) is the mass density, is the mean velocity, and θ is the temperature defined by The Fokker-Planck equation is a fundamental model in kinetic theories and statistical mechanics. Its general form reads (1.4) The one-particle friction constant γ is usually assumed to be a function of ρ, u, θ. Equation (1.4) has a stationary solution of given mass ρ, mean velocity u and temperature θ given by the Maxwellian density function which is such that J F P (M ρ,u,θ ) = 0. The Fokker-Planck equation appears in many different contexts. It was originally derived for the distribution function of a Brownian particle in a fluid [7], and is applicable in a more general form to a plasma [9]. A detailed investigation of this model has been performed by Frisch, Helfand, and Lebowitz [13] in connection with the kinetic theory of liquids. As shown more recently [27] (cf. also [6]), it provides also a good description of the grazing collisions in a one-dimensional gas. The Fokker-Planck operator J F P has the usual conservation properties of mass, mean velocity, and temperature, and logf J F P (f ) dv < 0, which guarantees the increasing in time of the (Shannon) entropy It is interesting to remark that, if the friction γ is taken to be proportional to the pressure p = ρθ , then J F P (f ) has the same kind of nonlinearity (quadratic) as the true Boltzmann equation.
For the purpose of accurate numerical simulations, a discretized Fokker-Planck equation must guarantee most of the conservation laws of the original equation, starting from mass conservation. Furthermore, since the solution of the Fokker-Planck equation represents a density function, any numerical scheme that approximates equation (1.4) is required to guarantee the positivity of the solution. In addition, it would be desirable that an approximation scheme must be accurate and stable.
The seminal paper for the approximation to equation (1.4) is due to Chang and Cooper [8]. Other classical references are the paper by Larsen, Levermore Pomraning and Sanderson [16], and the well-known book by Risken [24]. Various aspects of the numerical approximation of Fokker-Planck equation were subsequently dealt with by a number of authors [2,3,10,11,12,20]. Also in recent times, this problem has attracted the interest of research [19,22].
The aim of this paper is to present a discretized version of equation (1.4) which maintains most of the physical properties of the original equation. These properties include conservation of mass and positivity of the discrete solution, same evolution for the mean velocity and temperature, monotonicity in time of the discrete Shannon entropy, and the existence of an explicit discrete equilibrium density. In addition, the problem of the large-time behavior of the approximation and the convergence to the corresponding equilibrium density has been dealt with in the one-dimensional situation.
This discrete version is largely inspired by a recent paper [23], in which the kinetic meaning of the approximation to the heat equation proposed by Rosenau in [25] has been deeply investigated.
The last part of the paper is devoted to show how this idea could be fruitfully applied to construct a numerical approximation to one-dimensional linear diffusion equations of higher order. In particular, it is shown that starting from this approximation one can easily obtain an explicitly computable formula for the central difference approximation of a derivative of any even order.
2. Main properties of the Fokker-Planck equation. Given a nonnegative initial value ϕ(v) with finite mean velocity u and temperature θ, easy computations show that the mass, mean velocity, and temperature of the solution to the Fokker-Planck equation (1.4) do not change with time. It is convenient to normalize f to be a probability density instead of a mass density, and change equation (1.4) to a dimensionless form. To do this, one introduces the dimensionless variablesv,t, and the dimensionless functionsφ,f defined by the formulas (2.1) Substituting (2.1) into (1.4), carrying out elementary calculations, and then omitting the bars, we obtain that the function f (v, t) will now satisfy the equation with the initial condition ϕ(v) and consequently f (v, t) satisfying the following simple normalization conditions The normalization (2.3) corresponds to the equilibrium Maxwellian density Let ϕ be any probability density on IR d with finite second moment. Let X be any random variable with this density, and let W be any independent Gaussian random variable with density M given by (2.4). For every t > 0 define Then, the random variable Z t has a density f (v, t) at each t ≥ 0, and it is well-known that f (t) is evolved from ϕ under the action of the adjoint Ornstein-Uhlenbeck semigroup. Therefore f (v, t) satisfies equation (2.2), which can of course be checked directly from the definition. Mean velocity and temperature of the solution at any time t ≥ 0 can be obtained directly from expression (2.5). We obtain and A direct computation then shows that the following laws of evolution hold Of course, in the case in which ϕ satisfies (2.3), (2.8) imply conservation of both mean velocity and temperature. It is well-known that the solution to the Fokker-Planck equation (2.2) converges exponentially in time to zero in relative entropy [26,28], which implies exponential convergence to equilibrium in L 1 -norm. A slightly less known result is that exponential convergence to equilibrium for non regular initial data can be directly shown to hold also in weaker norms by resorting to the Fourier transform. Given a probability density f (v), v ∈ IR d , we define its Fourier transform f (ξ), ξ ∈ IR d by Let us consider a family of metrics that has been introduced in the paper [15] to study the trend to equilibrium of solutions to the space homogeneous Boltzmann equation for Maxwell molecules, and subsequently applied to a variety of problems related to kinetic models of Maxwell type. For a more detailed description, we address the interested reader to the lecture notes [5].
Given s > 0 and two random variables X, Y with probability distributions f (respectively g), their Fourier based distance d s (X, Y ) is given by the quantity The distance is finite, provided that X and Y have the same moments up to order [s], where, if s / ∈ N, [s] denotes the entire part of s, or up to order s − 1 if s ∈ N. Moreover d s is an ideal metric. Its main properties are the following 1. Let X 1 , X 2 , X 3 , with X 3 independent of the pair X 1 , X 2 be random variables with probability distributions f 1 , f 2 , f 3 . Then where the symbol * denotes convolution; 2. Define, for a given nonnegative constant a, the dilation of a function Then, given two random variables X, Y with probability distributions f and g, for any nonnegative constant a d s (aX, aY ) = d s (f a , g a ) ≤ a s d s (f, g) = a s d s (X, Y ).
Consider that, in view of the representation formula (2.5), given the initial values ϕ andφ, the two solutions f (t) andf (t) are the probability distributions of the random variables Z t andZ t expressed by (2.5). Hence if for some s > 0 the distance (2.9) The first inequality follows from property 1. Then the dilation property 2 is applied to conclude. The same result [4] can be easily obtained also resorting to the Fourier transform version of the Fokker-Planck equation (2.10), that reads By choosing φ(v) = e −iξv we obtain the (one-dimensional) Fourier transform version (2.10) of the Fokker-Planck equation (2.2). The advantage of working with a weak version of the equation, is that we can allow the initial value ϕ(v) to be a measure on IR.
3.1. A Rosenau-type approximation. Rosenau [25] proposed a regularized version of the Chapman-Enskog expansion of hydrodynamics. This regularized expansion resembles the usual Navier-Stokes viscosity terms at law wave-numbers, but unlike the latter, it has the advantage of being a bounded macroscopic approximation to the linearized collision operator. The model originally considered by Rosenau is given by the scalar equation The operator on the right side looks like the usual viscosity term εf vv at low wave-numbers ξ, while for higher wave numbers it is intended to model a bounded approximation of a linearized collision operator, thereby avoiding the artificial instabilities that occur when the Chapman-Enskog expansion for such an operator is truncated after a finite number of terms.
The right side of (3.12) can be written as where * denotes convolution and is a non-negative symmetric function satisfying M γ L 1 = 1. In (3.14) with the linear kinetic equation A detailed study of the properties of the kinetic equation (3.16), as well as its connections with the diffusion equation (3.15) was recently given in [23] (cf. also [18,21]). Also, a similar approximation was used in connection with the one-dimensional fractional diffusion equation [14]. These studies showed that the approximation maintains most of the properties of the original diffusion equations. Moreover, in the case of the heat equation, the Rosenau-type approximation gives a new physical inside into the numerical approximation of (3.15). In particular the moments at the first two order of the approximate solutions follow the same evolution of the original diffusion equation. However, it is important to notice that (3.13) is only a physically relevant way to write the right-hand side of (3.12). Indeed, by making use of standard properties of the convolution operation one has This shows that Rosenau approximation is obtained by smoothing out the righthand side of the linear diffusion equation (3.15) by means of its convolution with M ε . Hence, the linear kinetic equation (3.16) can be alternatively written as It is tempting to use the same approximation for the one-dimensional Fokker-Planck equation (2.2). In this case one considers the equation In Fourier transform, equation (3.19) reads It is clear that the diffusion term is given by a linear kinetic equation of type (3.13). Then, the Fourier transform of the drift term can be written as ∂ξ .
Hence we obtain the identity Note that, while M γ (v) is a symmetric probability density,M γ is antisymmetric, and it is obtained from M γ (v) by changing the sign on the domain v ≥ 0. Finally, the approximated Fokker-Planck equation (3.19) can be equivalently written as

The discrete Fokker-Planck equation.
Having in mind the previous discussion, given a small positive parameter ε 1 we consider the approximation to the operator J F P given by acting on probability densities f (v) satisfying the normalization conditions (2.3). In (3.24) where δ(v) denotes as usual the Dirac delta function concentrating on v = 0. Note that, in analogy with (3.23), P ε (v) is a symmetric probability measure, andP ε is antisymmetric, and obtained from P ε by changing its sign when v > 0.
Using the symmetry properties of P ε andP ε , one shows that the weak form of the evolution equation Alternatively, by choosing Φ(v) = e −iξv one obtains the evolution equation for the The Fourier description clearly shows that as ε → 0, the right-hand side of equation (3.28) converges pointwise to the right-hand side of equation (2.10). Let us examine in details the properties of the solution to equation (3.27). In view of definition (3.25), it follows that for 0 ≤ n ∈ N A further fundamental property of the solution to equation (3.27) follows by considering as initial datum the law of a random variable X ε that takes values only on a discrete number of points. To be more precise, for a given positive number N ∈ N, N 1, let us set ε = 1/N . Then we define where the nonnegative constants ϕ j satisfy Let D be the space of functions of type (3.30), subject to condition (3.31). Let us consider a nonnegative measure g(v) in D. Owing to definition (3.25), it is immediate to verify that T (g) = ε 2 J ε F P (g)/2 + g belongs to D. Indeed (3.32) Therefore T (g)(v) is a linear combination of the values of g in the points g(v + ε) and g(v − ε), and whenever v belongs to the interval (−2/ε, 2/ε) the coefficients of the linear combination are nonnegative. Moreover, a direct inspection shows that Last, condition (3.31) remains verified. This property can be verified by owing to the mass conservation property of J ε F P , or by direct inspection. Indeed, using (3.33) we have Consequently, T (g) is a linear mapping of D into D. Now, given the nonnegative measure ϕ ε ∈ D, consider that the initial value problem for the discrete Fokker-Planck equation (3.26) can be rewritten as The (unique) solution to (3.34) can be explicitly expressed in the form of a Wild sum [21,29] f ε (v, t) = e −2t/ε 2 is the initial value, and the nonnegative coefficients f . Since at any time t ≥ 0 the solution (3.35) is a convex combination of the time-independent nonnegative coefficients f (i) ε , the solution to the initial value problem (3.34) with ϕ ε ∈ D belongs to D at any subsequent time t ≥ 0. In addition, the solution is nonnegative for all t ≥ 0.
Therefore, integrating both sides from 0 to ξ, and assuming that the solution has mass equal to one, so that f ∞ (0) = 1 one obtains Hence, the function (3.42) is the characteristic function of the random variable X/N . Consequently, if ε = 1/N , the function (3.41) is the characteristic function of the random variable where the X j are independent and identically distributed copies of X, defined as in (3.43). By construction, the random variable S N takes values only on the discrete set of points εj, where |j| ≤ 2N 2 . Consequently, S N has a probability distribution that belongs to D. Considering now that X has zero mean, and variance 1/2, it follows easily by the central limit theorem that the law of S N is an approximation of the Gaussian distribution, that is an approximation of the stationary solution of the original Fokker-Planck equation. In addition, the law of S N satisfies the normalization conditions (2.3).

3.4.
Large-time behavior. The results of Section 3.2 showed that the solution to the discrete Fokker-Planck equation (3.26) in D maintains most of the properties of the original Fokker-Planck equation, like preservation of positivity and mass, same evolution of moments up to order two, and entropy monotonicity. A further property of the original Fokker-Planck equation is the exponential convergence of its solution towards the Maxwellian equilibrium. As discussed in Section 2, exponential convergence to equilibrium can be shown also in Fourier metric. In what follows, we will investigate about the large-time behavior of the solution to (3.26) and its (eventual) convergence to equilibrium in terms of the metric d s . Let εξ ∈ [mπ, (m + 1)π), where |m| ∈ N. Then, since the function sin εξ does not change sign in this interval, the differential equation can be solved uniquely as soon as εξ(t) ∈ [mπ, (m + 1)π) to give the relationship tan εξ(t) 2 = e t tan εξ 2 . (3.46) Then, since both εξ(t)/2 and εξ/2 belong to the interval [m π 2 , (m + 1) π 2 ) identity (3.46) can be used to relate in a unique way ξ to ξ(t) (or vice-versa)   We remark that, if the initial value ϕ ε (v) ∈ D, its Fourier transform is given by where the nonnegative constants ϕ j satisfy condition (3.31). Then, for any valuē ξ = 2mπ/ε, with m ∈ N, ϕ ε (ξ) = 1. By letting t → ∞ in (3.50), for any ξ ∈ IR such that εξ ∈ [mπ, (m + 1)π) Since m is arbitrary, pointwise convergence to the stationary solution (3.41) follows for all ξ ∈ IR.
A stronger result about convergence to the steady state follow by restricting the allowed set of values of ξ.
Indeed, since the function (3.41) is a stationary solution to the Fokker-Planck equation, for every t ≥ 0 it satisfies the identity (3.53) Therefore, considering that for t ≥ 0 1 + e −2t (tan(εξ/2)) 2 1 + (tan(εξ/2)) if the initial value ϕ ε (v) has zero mean, we have the inequality arctan (e −t tan(εξ/2)) |ξ| 2 , and this inequality, on the set |εξ| ≤ π/2, clearly implies arctan (e −t tan(εξ/2)) arctan (e −t tan(εξ/2)) This implies exponential convergence of the solution to the stationary state (on the set |εξ| ≤ π/2) at the same rate of the Fokker-Planck equation (2.2). On the other hand, on the set |εξ| > π/2, since both | f |(ξ, t) ≤ 1 and | f ∞ |(ξ) ≤ 1, we have the bound (3.55) 3.5. Stability of the approximation. Let us consider a Taylor expansion of the trigonometric functions on the right-hand side of equation (3.28) up to the second order. We obtain where the remainder term has the expression andξ belongs to the interval (0, ξ). Let the initial value ϕ of the Fokker-Planck equation (2.2) possess finite moments up to the order three, and let us consider an approximation ϕ ε ∈ D with the same moments of ϕ up to the second order. Then, since the mean velocity and the temperature of the Fokker-Planck equation and of its approximation follow the laws (2.8), it is immediate to conclude that 58) where u 0 is the initial mean velocity defined by (1.2). In addition, if the Fourier metric d 3 (f, f ε )(t) is initially bounded, it remains bounded at any subsequent time.
Let h(ξ, t) be defined as Since for ξ = 0 we have the identity by considering the difference between the Fourier transforms of one-dimensional Fokker-Planck equation (2.10) and (3.56), we conclude that h satisfies Integrating along characteristics, (3.59) is equivalent to Hence, by using (3.58) we obtain In conclusion (3.60) Hence, by choosing initial data for the discrete Fokker-Planck equation such that d 3 (ϕ, ϕ ε ) ≤ Cε, we obtain that in the Fourier distance d 3 the uniform in time estimate 4. The general case. The discretization of the one-dimensional Fokker-Planck equation introduced in Section 3 can be easily extended to any dimension d > 1.
To this aim, it is enough to outline that the Fourier transformed equation (2.10) can be rewritten as The Rosenau approximation in this case is given by that in the physical space reads As in Section 3, one can easily verify that the mean velocity and temperature of the solution to (3.23) follow the same laws of evolution of the original Fokker-Planck equation (2.2). Also, using the functions (3.25) one obtains the evolution equation for the Fourier transform f ε (ξ, t) In the limit ε → 0, the right-hand side of equation (3.28) converges pointwise to the right-hand side of equation (4.1). As in the one-dimensional case, let us consider as initial datum the law ϕ ε (v) = ϕ(v 1 , v 2 , . . . , v d ) of a random vector (X (1) ε , X where the nonnegative constants ϕ j1,...,j d satisfy Let D d be the space of functions of type (3.30), subject to condition (3.31). Proceeding as in Section 3, we can prove that, starting from a nonnegative initial value in D d , the solution to equation (4.3) belongs to D d for each time t ≥ 0. Moreover, the solution is nonnegative. Also, mass is preserved, and the laws of evolution of the mean velocity and temperature follow the same laws of evolution (2.8) of the continuous equation. Last, Shannon entropy is monotonically increasing. Concerning the equilibrium distribution, we assume that (4.7) Then we have is the Fourier transform of the joint distribution function of the random vector X = (X 1 , . . . , X d ) where the random variables X k , for k = 1, . . . , d are independent each other and distributed according to (3.43). Consequently, if ε = 1/N , the function (4.7) is the characteristic function of the random vector where the X j are independent and identically distributed copies of the random vector X Considering now that X has zero mean, and variance d/2, it follows by the central limit theorem that the law of S N is an approximation of the d-dimensional Gaussian distribution, that is an approximation of the stationary solution of the original Fokker-Planck equation in IR d . Unlikely, the analysis of Section 3.4 is no more valid in dimension d > 1, and the study of the large-time behavior of the solution to the discretized Fokker-Planck equation requires further efforts. We leave this problem open for future research.
Following the idea of Rosenau [25], for a given ε 1 we consider the approximation given by This approximation is consistent with the analogous one introduced for the Fokker-Planck equation. Note that, since with M ε defined as in (3.14), the approximation (5.3) corresponds to modify equation (5.1) by taking the convolution of the right-hand side with M ε * M ε . Hence the approximation in the physical space reads Since

equation (5.4) can be rewritten as a linear kinetic equation of the form
where the function has integral equal to one, but, at difference with the case of the heat equation, is no more a probability density function. Indeed, it becomes negative on part of the real line. This fact is in agreement with the theory of higher-order diffusions, and it is connected with the absence of a maximum principle for the solution to these equations. It is remarkable that, in view of the expression (5.4), the moments up to order four of the solution to the approximated equation follow the same evolution of the original equation.
The same property is maintained by considering in the definition of G ε a probability density different from M ε . In particular, we can resort to P ε , as defined in (3.25). Consequently, we consider the approximation where the constant 4 in front of the interaction operator is chosen to preserve the evolution of the fourth-order moment. By resorting to the definition of P ε (v), it is immediate to conclude that Hence, for a given (smooth) function h = h(v) (5.8) Expression (5.8) coincides with one of the central differences approximation of the fourth order derivative of a function (cf. [17] and the references therein). The previous reasoning allows to conclude that, at difference with other approximations, (5.8) is well-adapted, in view of its properties about preservation of moments evolution, to approximate the diffusion equation (5.1).
By construction, the solution to the kinetic equation (5.13) is such that, for all k ≤ 2n − 1 Likewise, given a (smooth) function h(v), the expression where as usual P ε (v) is given by (3.25), gives an explicitly computable central difference approximation of order 2n of h with a number of good properties with respect to the higher-order diffusion equation (5.10). Indeed, the solution to the approximated diffusion equation ∂g(v, t) ∂t = (2n)! 1 ε 2n [G n (P ε ) * g(v, t) − g(v, t)] is such that, in agreement with the solution to the original diffusion equation, all moments up to the order 2n − 1 remain constant in time, while the moment of order 2n is linearly increasing with the same rate. 6. Conclusions. We introduced and discussed a discretized version of the Fokker-Planck equation which maintains most of the properties of the continuous version. The basic idea was to use a suitable modification of the Fourier transform of the equation, similar to the one considered by Rosenau [25] for the linear heat equation. This approach has been subsequently applied to higher-order linear diffusion operators, to obtain an easy-to-handle way to recover explicitly a central difference approximation to derivatives of any even order. A main problem, however, remains open. It is not clear wether the solution to the discrete Fokker-Planck equation converges towards the stationary discrete solution or not. A partial result in this direction has been derived in Section 3.4. Also, it would be interesting to know if, in the set D, and for random variables with a law that satisfies the normalization condition (2.3), the Shannon entropy (3.39) attains the maximum value in correspondence to the law of the stationary distribution S N defined in (3.44).