A Primal-Dual Weak Galerkin Finite Element Method for Fokker-Planck Type Equations

# A Primal-Dual Weak Galerkin Finite Element Method for Fokker-Planck Type Equations

Chunmei Wang Department of Mathematics, Texas State University, San Marcos, TX 78666, USA. The research of Chunmei Wang was partially supported by National Science Foundation Award DMS-1522586 and DMS-1648171.    Junping Wang Division of Mathematical Sciences, National Science Foundation, Arlington, VA 22230 (jwang@nsf.gov). The research of Junping Wang was supported by the NSF IR/D program, while working at National Science Foundation. However, any opinion, finding, and conclusions or recommendations expressed in this material are those of the author and do not necessarily reflect the views of the National Science Foundation.
###### Abstract

This paper presents a primal-dual weak Galerkin (PD-WG) finite element method for a class of second order elliptic equations of Fokker-Planck type. The method is based on a variational form where all the derivatives are applied to the test functions so that no regularity is necessary for the exact solution of the model equation. The numerical scheme is designed by using locally constructed weak second order partial derivatives and the weak gradient commonly used in the weak Galerkin context. Optimal order of convergence is derived for the resulting numerical solutions. Numerical results are reported to demonstrate the performance of the numerical scheme.

p

rimal-dual, weak Galerkin, finite element methods, Fokker-Planck equation, weak Hessian, weak gradient, polytopal partitions.

{AMS}

Primary, 65N30, 65N15, 65N12, 74N20; Secondary, 35B45, 35J50, 35J35

## 1 Introduction

The Fokker-Planck equation plays a critical role in statistical physics and in the study of fluctuations in physical and biological systems [11, 23, 24, 27, 12]. In statistical physics, it is a second order partial differential equation that describes the time evolution of the probability density function of the velocity of a particle under the influence of drag forces and random forces resulting from Gaussian white noise. The general setting of the Fokker-Planck equation is as follows. Given an open domain (the -dimensional Euclidean space) and a terminal time , we seek for a time-dependent density function satisfying

 (1.1) ∂tp+∇⋅(μp)−12d∑i,j=1∂2ij(aijp)= 0,t∈(0,T), x∈Ω,p(x,0)= p0(x),x∈Ω,

where is the second order partial derivative in the directions and , is the diffusion tensor, is the drift vector, and is the initial profile of the density function. Two common boundary conditions for (1.1) can be imposed: Dirichlet for the density function and Neumann condition for the flux. A homogeneous Dirichlet boundary data corresponds to the case where particles exit once they reach the boundary, and a prescribed flow or Neumann boundary condition represents a known current of particles crossing the boundary in the normal direction.

Numerical methods for the Fokker-Planck equation have several challenges involving various difficulties of different nature. Among them are the high dimensionality and lack of solution regularity for the probability density function. For example, in classical statistical mechanics, the Fokker-Planck equation characterizes a joint probability density function in many phase variables so that the dimension might be a big number to deal with. It can also be seen that for non-smooth diffusion tensor , the resulting probability density function exhibits a shock-like discontinuity that needs to be resolved numerically. In addition, some conservation properties such as mass conservation and solution non-negativity property must be retained by the numerical solutions.

Various finite element methods have been designed for the Fokker-Planck equation for a numerical computation of the probability density function, see [3, 17, 18, 2, 26, 21, 16] and the references cited therein. In [3, 17], the stationary Fokker-Planck equation was discretized by using Galerkin finite element methods based on a weak form obtained from the usual integration by parts. In [18], another Galerkin finite element method was used to solve the Fokker-Planck equation in combination with a generalized Lagrange multiplier method to handle the associated integral constraint. In [2], the authors applied the usual finite element method to the Fokker-Planck system subject to both additive and multiplicative white noise excitations. In [21], the authors developed a framework for multi-scale finite element methods for the solution of the multi-dimensional Fokker-Planck equation in stochastic structural dynamics. All the aforementioned finite element methods assumed smooth or constant diffusion tensor for the Fokker-Planck equation so that a regular weak form can be derived for the system. In [16], a finite element scheme was described and numerically tested for the transient Fokker-Planck equation without any smoothness assumption on ; but no theory of convergence was developed for the numerical method.

For a smooth diffusion tensor, the second order differential part of the Fokker-Plank equation can be reformulated as

 (1.2) 12∂2ij(aijp)=12∂j(aij(x)∂ip)+12∂j((∂iaij)p).

Therefore, the equation (1.1) can be viewed as a time-evolving convection-diffusion equation in divergence form. Aside from the high dimensionality issue, the corresponding Fokker-Planck equation is then considered as a relatively less challenging problem to solve numerically. But for non-smooth diffusion tensor, the formulation (1.2) no longer holds true, and the exact profile of the density function possesses discontinuities that are not known a priori so that the existing finite element methods have difficulty to apply. The goal of this paper is to develop a new finite element method that addresses the numerical challenges arising from the non-smoothness nature of the diffusion tensor in the Fokker-Planck equation.

For simplicity, we consider a Fokker-Planck type model equation with homogeneous Dirichlet boundary condition. The model problem seeks an unknown function satisfying

 (1.3) ∇⋅(μu)−12d∑i,j=1∂2ij(aiju)=f,in Ω,u=0,on ∂Ω,

where is an open bounded domain in with Lipschitz continuous boundary and is a given function. Two of our main motivations for this selection of the model problem are: (1) a complete understanding of the Fokker-Planck equation strongly depends on the numerical properties for (1.3) as it offers a projection operator that is extremely useful in the mathematical study of (1.1), and (2) the problem (1.3) itself is a poorly understood system from numerical aspects when the diffusion tensor is discontinuous. Therefore, the model problem (1.3) deserves a study as a research topic in numerical partial differential equations.

Throughout this paper, we assume that the diffusion tensor is symmetric, uniformly bounded and positive definite in , and that the drift vector . We will follow the usual notation for Sobolev spaces and norms [9, 13, 15, 14, 4]. For any open bounded domain with Lipschitz continuous boundary, we use and to denote the norm and seminorms in the Sobolev space for any , respectively. The inner product in is denoted by . The space coincides with , for which the norm and the inner product are denoted by and , respectively. When , we shall drop the subscript in the norm and inner product notation. For convenience, we use “ ” to denote “less than or equal to up to a general constant independent of the mesh size or functions appearing in the inequality”.

By a weak solution of (1.3) we mean a function satisfying

 (1.4) (u,Lv)=−(f,v),∀v∈H2(Ω)∩H10(Ω),

where is a differential operator given by . The differential operator is assumed to satisfy the -regularity property in the sense that for any given , there exists a unique strong solution satisfying

 (1.5) LΦ=χ,∥Φ∥2≲∥χ∥.

In [25], it was shown that the regularity assumption (1.5) holds true on bounded convex domain if and the diffusion tensor satisfies the following Cordès condition

 (1.6) ∑di,j=1a2ij(∑di=1aii)2≤1d−1+εin Ω,

for a parameter . The Cordès condition (1.6) is automatically satisfied in 2D for diffusion tensor that is bounded, symmetric, and uniformly positive definite in the domain, see [30] for a verification.

Our numerical scheme for the model problem (1.3) is based on the weak formulation (1.4) through a weak Galerkin approach that combines the primal variable with its dual. The dual problem for the primal equation (1.4) is given by

 (1.7) (w,Lρ)=0,∀w∈L2(Ω),

where is the dual variable. Under the -regularity assumption (1.5), the solution to the dual problem (1.7) is clearly trivial; i.e., . Note that the primal and the dual equations are formally uncorrelated to each other in the continuous case; but this changes significantly in the context of weak Galerkin finite element methods. In the weak Galerkin approach, the differential operator is discretized as

 Lw(v):=μ⋅∇wv+12d∑i,j=1aij∂2ji,wv,

where is a discrete weak gradient [31, 19, 33, 32] and is a discrete weak second order partial derivative [20, 28] (also see Section 2 for their definition). The corresponding primal and dual equation then become to be

and

where and are two weak finite element spaces used to approximate and respectively. While neither (1.8) nor (1.9) makes any computationally feasible schemes, their combination through the use of a suitably-defined stabilizer does provide numerical methods that are efficient, accurate, and stable for several model problems [7, 8, 6, 30]. A formal description of the scheme reads as follows: Find and such that

 (1.10) s(ρh,v)+(uh,Lwv)=−(f,v),∀v∈V0h,k,(w,Lwρh)=0,∀w∈Wh,s,

where is a bilinear form in known as stabilizer or smoother that enforces a certain weak continuity for the approximation . Numerical schemes in the form of (1.10) were named primal-dual weak Galerkin finite element methods in [30], but were broadly called stabilized finite element methods in [7, 8, 6].

In the rest of the paper, we will provide all the technical details for the numerical scheme (1.10), including the construction of the finite element spaces and , representation of the stabilizer or smoother , mathematical convergence for the corresponding numerical approximations, and some numerical results that demonstrate the performance of the method. One of the distinguished features of this approach lies on ultra weak regularity assumptions for the primal variable in the mathematical convergence theory. The method essentially assumes no regularity on the primal variable so that solutions with discontinuity can be well approximated by our primal-dual finite element method. This work is a non-trivial extension of [30] in both theory and algorithmic development.

The paper is organized as follows. In Section 2, we shall briefly discuss the computation of weak gradients and weak second order partial derivatives. In Section 3, we will present a detailed description of the primal-dual weak Galerkin finite element method for the Fokker-Planck type model problem (1.3) based on the weak formulation (1.4). In Section 4, we will study the solution properties for our numerical method. In particular, we shall derive an inf-sup condition, and then establish a result on the solution existence and uniqueness. In Section 5, we will derive an error equation for the numerical solutions. Then in Section 6, we will establish an error estimate for the primal variable that is of optimal order in . Section 7 is devoted to a presentation of error estimates for the usual projections. Finally in Section 8, we report some numerical results to demonstrate the performance of the primal-dual weak Galerkin finite element method.

## 2 Weak Partial Derivatives

The goal of this section is to brief the definition and computation of the discrete weak partial derivatives introduced in [28, 29, 32]. To this end, let be a polygonal or polyhedral region with boundary . By a weak function on we mean a triplet in which , and . The first and second components and are intended for the value of in the interior and on the boundary of , respectively. The third component is used to represent the gradient of on . In general, and are not required to be consistent with the trace of and on . Denote by the space of all weak functions on :

 (2.1) W(T)={v={v0,vb,vg}:v0∈L2(T),vb∈L2(∂T),vg∈[L2(∂T)]d}.

For any , the weak second order partial derivative , denoted as , is defined as a linear functional in the dual space of satisfying

 (2.2) ((∂2ij,wv,φ))T=(v0,∂2jiφ)T−⟨vbni,∂jφ⟩∂T+⟨vgi,φnj⟩∂T,

for all . Here, stands for the action of at , is the unit outward normal direction to , is the usual inner product in , and is the inner product in .

The weak gradient of , denoted by , is defined as a linear functional in the dual space of such that

 ((∇wv,ψ))T=−(v0,∇⋅ψ)T+⟨vb,ψ⋅n⟩∂T,

for all . Note that the weak gradient makes no use of the third component of the weak function .

Denote by the set of polynomials on with degree no more than . A discrete version of for , denoted by , is defined as the unique polynomial in satisfying

 (2.3) (∂2ij,w,r,Tv,φ)T=(v0,∂2jiφ)T−⟨vbni,∂jφ⟩∂T+⟨vgi,φnj⟩∂T,

for all , which, by using integration by parts, yields

 (2.4) (∂2ij,w,r,Tv,φ)T=(∂2ijv0,φ)T+⟨(v0−vb)ni,∂jφ⟩∂T−⟨∂iv0−vgi,φnj⟩∂T,

for all .

A discrete form of for , denoted by , is defined as the unique polynomial vector in satisfying

 (2.5) (∇w,r,Tv,ψ)T=−(v0,∇⋅ψ)T+⟨vb,ψ⋅n⟩∂T,∀ψ∈[Pr(T)]d,

which, from integration by parts, gives

 (2.6) (∇w,r,Tv,ψ)T=(∇v0,ψ)T−⟨v0−vb,ψ⋅n⟩∂T,∀ψ∈[Pr(T)]d,

provided that is sufficiently regular.

## 3 Numerical Schemes

The goal of this section is to present a finite element method for the variational problem (1.4). To this end, let be a finite element partition of the domain into polygons in 2D or polyhedra in 3D which is shape regular in the sense as described in [32]. For three dimensional domains, all the polyhedral elements are assumed to have flat faces. Denote by the set of all edges or faces in and the set of all interior edges or faces. Denote by the meshsize of and the meshsize of the partition .

Let be a given integer. Denote by the discrete local weak function space on given by

 Vk(T)={{v0,vb,vg}:  v0∈Pk(T), vb∈Pk(e), vg∈[Pk−1(e)]d, e⊂∂T}.

By patching over all the elements through a common value and on the interior edges/faces, we obtain a global weak finite element space :

 Vh,k={{v0,vb,vg}:  {v0,vb,vg}|T∈Vk(T), T∈Th}.

Denote by the subspace of with vanishing boundary value for on , i.e.,

 V0h,k={{v0,vb,vg}∈Vh,k:  vb|e=0, e⊂∂Ω}.

For any given integer , denote by the usual finite element space consisting of piecewise polynomials of degree ; i.e.,

 Wh,s={w:  w|T∈Ps(T), T∈Th}.

For application in the approximation of (1.4), the integer will be chosen as either or . In the case of , the only viable option for this integer would be .

For simplicity of notation, denote by the discrete weak gradient computed by using (2.5) on each element with :

 (∇wσ)|T=∇w,k−1,T(σ|T),σ∈Vh,k.

Analogously, we use to denote the discrete weak second order partial derivative computed by using (2.3) on each element with :

 (∂2ij,wσ)|T=∂2ij,w,s,T(σ|T),   σ∈Vh,k.

The corresponding weak differential operator is defined by using weak partial derivatives as follows

 Lw(σ)=μ⋅∇wσ+12d∑i,j=1aij∂2ji,wσ,

for any .

Let us introduce the following bilinear forms

 s(ρ,σ)= ∑T∈ThsT(ρ,σ),ρ,σ∈Vh,k, b(σ,v)= ∑T∈ThbT(σ,v),σ∈Vh,k, v∈Wh,s,

where

 (3.1) sT(ρ,σ)=h−3T∫∂T(ρ0−ρb)(σ0−σb)ds+h−1T∫∂T(∇ρ0−ρg)(∇σ0−σg)ds+δ∫TL(ρ0)L(σ0)dT

and

 bT(σ,v)=(v,Lw(σ))T.

Here, is a parameter independent of the meshsize and the functions involved.

We are now in a position to state our primal-dual weak Galerkin finite element scheme for the model variational problem (1.4).

###### Primal-Dual Weak Galerkin Algorithm 3.1

Let be a given integer and be another integer. A numerical approximation for the solution of (1.4) is the component in satisfying

 (3.2) s(ρh,σ)+b(σ,uh) = −(f,σ0),∀σ∈V0h,k, (3.3) b(ρh,v) = 0,∀v∈Wh,s.

## 4 Solution Existence, Uniqueness, and Stability

The goal of this section is to study the solution for the numerical scheme (3.2)-(3.3). In particular, we shall prove the existence and uniqueness for the numerical solution under certain assumptions on the finite element partition and the differential operator .

On each element , denote by the projection onto , . Similarly, on each edge or face , denote by and the projections onto and , respectively. For any , define the projection so that on each element one has

 (4.1) Qhw={Q0w,Qbw,Qg(∇w)}.

Denote by the projection onto - the space of piecewise polynomials of degree . In the rest of this paper, the integer will be taken as either or .

{lemma}

[28, 29, 32] The aforementioned projection operators satisfy the following commutative properties: For any , one has

 (4.2) ∂2ij,w(Qhw)= Q(s)h(∂2ijw),i,j=1,…,d, (4.3) ∇w(Qhw)= Q(k−1)h(∇w).
{proof}

For any and , from (2.3), the usual property of projections, and the integration by parts, we have

 (∂2ij,w(Qhw),φ)T=(Q0w,∂2jiφ)T−⟨Qbw,∂jφ⋅ni⟩∂T+⟨Qgi(∂iw)⋅nj,φ⟩∂T=(w,∂2jiφ)T−⟨w,∂jφ⋅ni⟩∂T+⟨∂iw⋅nj,φ⟩∂T=(∂2ijw,φ)T=(Q(s)h∂2ijw,φ)T,

which completes the proof of (4.2). The other identity (4.3) can be derived in a similar fashion, and details can be found in [28, 29, 32].

The stabilizer defined through (3.1) naturally induces a semi-norm in the weak finite element space as follows:

 (4.4) |||ρ|||w=s(ρ,ρ)12,ρ∈Vh,k.
{lemma}

(inf-sup condition) Assume that the drift term and the coefficient tensor is uniformly piecewise continuous with respect to the finite element partition . Then, there exists a constant such that for any , there exists a weak function satisfying

 (4.5) b(v,σ) ≥ 12∥v∥2, (4.6) ≤ β∥v∥,

provided that meshsize satisfies for a small, but fixed parameter value .

{proof}

Let be the solution of the following auxiliary problem:

 (4.7) LΦ = v,in Ω, (4.8) Φ = 0,on ∂Ω.

From the assumption (1.5), the problem (4.7)-(4.8) has the following -regularity estimate

 (4.9) ∥Φ∥2≤C∥v∥.

With , we have from Lemma 4 that

 (4.10) b(v,σ)=∑T∈Th(v,Lw(QhΦ))T=∑T∈Th(μv,∇wQhΦ)T+12d∑i,j=1(aijv,∂2ji,wQhΦ)T=∑T∈Th(μv,Q(k−1)h(∇Φ))T+12d∑i,j=1(aijv,Q(s)h∂2jiΦ)T=∑T∈Th(μv,∇Φ)T+12d∑i,j=1(aijv,∂2jiΦ)T+(μv,(Q(k−1)h−I)∇Φ)T+12d∑i,j=1(aijv,(Q(s)h−I)∂2jiΦ)T=∑T∈Th(v,v)T+∑T∈Th(μv,(Q(k−1)h−I)∇Φ)T+12d∑i,j=1((aij−¯aij)v,(Q(s)h−I)∂2ijΦ)T,

where stands for the average of on each element . As the coefficient tensor is uniformly piecewise continuous in , there exists a small parameter depending on the meshsize and the continuity of on each element such that

 ∣∣ ∣∣∑T∈Th(μv,(Q(k−1)h−I)∇Φ)T∣∣ ∣∣≤Ch∥Φ∥2∥v∥∣∣ ∣∣d∑i,j=1((aij−¯aij)v,(Q(s)h−I)∂2ijΦ)T∣∣ ∣∣≤Cε(h)∥Φ∥2∥v∥.

Substituting the above estimates into (4.10) yields

 (4.11) b(v,σ)≥∥v∥2−C(h+ε(h))∥Φ∥2∥v∥≥(1−Ch−Cε(h))∥v∥2,

where we have used the regularity estimate (4.9). In particular, for , there exists a parameter value such that when . It follows that holds true for , and hence

 (4.12) b(v,σ)≥12∥v∥2,

which verifies the inequality (4.5).

It remains to establish the estimate (4.6) for . To this end, from the usual trace inequality we have

 (4.13) ∑T∈Thh−3T∫∂T|σ0−σb|2ds=∑T∈Thh−3T∫∂T|Q0Φ−QbΦ|2ds≤∑T∈Thh−3T∫∂T|Q0Φ−Φ|2ds≲∑T∈Thh−4T∫T|Q0Φ−Φ|2dT+h−2T∫T|∇Q0Φ−∇Φ|2dT≲∥Φ∥22≲∥v∥2.

A similar analysis can be applied to yield the following estimate:

 (4.14) ∑T∈Thh−1T∫∂T|∇σ0−σg|2ds≲∥v∥2.

Furthermore, we have

 (4.15) ∑T∈Thδ∫T|L(σ0)|2dT≲∑T∈Th∥σ0∥22,T≲∑T∈Th∥Q0Φ∥22,T≲∥Φ∥22≲∥v∥2.

Finally, by combining the estimates (4.13)-(4.15) and the definition of we obtain

 |||σ|||2w≤β2∥v∥2

for some constant . This completes the proof of the lemma.

We are now in a position to state the main result on solution existence and uniqueness.

{theorem}

Assume that the drift vector and the coefficient tensor are uniformly piecewise continuous in with respect to the finite element partition . Then, there exists a fixed such that the primal-dual weak Galerkin finite element algorithm (3.2)-(3.3) has one and only one solution if the meshsize satisfies .

{proof}

It suffices to show that zero is the only solution to the problem (3.2)-(3.3) with homogeneous data and . To this end, assume and in (3.2)-(3.3). By choosing and , the difference of (3.3) and (3.2) gives , which implies and on each , and hence . Moreover, we have

 Lρ0=0.

Thus, from the solution existence and uniqueness for the differential operator with homogeneous Dirichlet boundary condition we obtain , and hence .

With , the equation (3.2) now becomes

 (4.16) b(uh,τ)=0,∀τ∈V0h,k.

From Lemma 4, there exists a such that . It then follows from (4.16) that . This completes the proof of the theorem.

## 5 Error Equations

In this section we shall derive an error equation for the numerical solution arising from the primal-dual weak Galerkin finite element algorithm (3.2)-(3.3). To this end, let and be the solution of (1.4) and (3.2)-(3.3), respectively. Note that is supposed to approximate the trivial function .

{lemma}

For any and , the following identity holds true:

 (5.1) (Lwσ,v)T=(Lσ0,v)T+RT(σ,v),

where

 (5.2) RT(σ,v)=12d∑i,j=1⟨σ0−σb,nj∂i(Q(s)h(aijv))⟩∂T−12d∑i,j=1⟨∂jσ0−σgj,niQ(s)h(aijv)⟩∂T−⟨σ0−σb,Q(k−1)h(μv)⋅n⟩∂T.
{proof}

From the formula (2.4) and (2.6) for the weak derivatives, we have

 (5.3) (Lw(σ),v)T=(μ⋅∇wσ,v)T+12d∑i,j=1(aij∂2ji,wσ,v)T=(∇wσ,Q(k−1)h(μv))T+12d∑i,j=1(∂2ji,wσ,Q(s)h(aijv))T=(∇σ0,μv)T+12d∑i,j=1(∂2jiσ0,aijv)T+RT(σ,v)=(Lσ0,v)T+RT(σ,v),

where is given by (5.2).

By error functions we mean the difference between the numerical solution and the interpolation of the exact solution; namely,

 (5.4) eh =uh−Q(s)hu, (5.5) ϵh =ρh−Qhρ=ρh.
{lemma}

Let and be the solutions arising from (1.3) and (3.2)-(3.3), respectively. Then, the error functions and satisfy the following equations

 (5.6) s(ϵh,σ)+b(σ,eh) = ℓu(σ),∀σ∈V0h,k, (5.7) b(ϵh,v) = 0,∀v∈Wh,s,

where is given by

 (5.8) ℓu(σ)=∑T∈Th⟨σb−σ0,(μu−Q(k−1)h(μQ(s)hu))⋅n⟩∂T+12∑T∈Thd∑i,j=1⟨σ0−σb,nj∂i(aiju−Q(s)h(aijQ(s)hu))⟩∂T−12∑T∈Thd∑i,j=1⟨(∂jσ0−σgj)ni,aiju−Q(s)h(aijQ(s)hu)⟩∂T+∑T∈Th(L(σ0),Q(s)hu−u)T.
{proof}

First of all, from (5.5) and (3.3) we have

 b(ϵh,v)=b(ρh,v)=0,∀v∈Wh,s,

which gives (5.7).

Next, notice that . Thus, by using (3.2) we arrive at

 (5.9)

The rest of the proof shall deal with the term . To this end, we use Lemma 5 to obtain

 (5.10) b(σ,Q(s)hu)=∑T∈Th(Lwσ,Q(s)hu)T=∑T∈Th(Lσ0,Q(s)hu)T+RT(σ,Q(s)hu)=∑T∈Th(Lσ0,u)T+(Lσ0,Q(s)hu−u)T+RT(σ,Q