1 Introduction↩︎

Arising from the study of stochastic processes, the Fokker-Planck equation provides an indispensable mathematical apparatus for analyzing statistical fluctuations across physical systems, engineering applications, and biological complexes [1]–[5]. In the context of statistical mechanics, this second-order partial differential equation governs the temporal evolution of a probability density function for a particle moving under the coupled effects of deterministic drag forces and stochastic disturbances modeled as Gaussian white noise. Structurally, given an open spatial domain \(\Omega\subset{\mathbb{R}}^d\) (\(d\)-dimensional Euclidean space) and a terminal time horizon \(T\), the objective is to find a time-dependent distribution function \(p=p(x,t): \;\Omega\times [0,T]\to {\mathbb{R}}\) that satisfies \[\label{EQ:FPE} \begin{align} \partial_t p +\nabla\cdot ({\boldsymbol{\mu}}p)-\frac{1}{2}\sum_{i,j=1}^d \partial_{ij}^2(a_{ij}p) & = \;0,\qquad t\in (0, T), \;x\in \Omega,\\ p(x, 0) & = \;p_0(x),\qquad x\in \Omega, \end{align}\tag{1}\] where \(\partial_{ij}^2=\frac{\partial}{\partial x_j}\frac{\partial}{\partial x_i}\) signifies the second-order partial differentiation with respect to coordinates \(x_i\) and \(x_j\). Within this system, \(a(x)=\{a_{ij}(x)\}_{d\times d}\) denotes the diffusion tensor, \({\boldsymbol{\mu}}=(\mu_1,\cdots, \mu_d)\) captures the drift vector field, and \(p_0=p_0(x)\) indicates the initial density profile. To ensure well-posedness, 1 is typically closed using one of two standard boundary conditions: a Dirichlet constraint on the density or a Neumann condition on the total flux. A homogeneous Dirichlet boundary condition represents an absorbing barrier where particles immediately exit the system upon impact, whereas a prescribed Neumann flux dictates a known net current of particles moving across the boundary along the outward normal direction.

Designing accurate numerical schemes for the Fokker-Planck equation involves overcoming several fundamental mathematical and computational challenges. Chief among these issues are the curse of high dimensionality and the low regularity inherent to the target probability density function. For example, in multi-particle statistical mechanics, the Fokker-Planck framework models a joint probability density spanning a vast phase space, making the spatial dimension \(d\) prohibitively large. Concurrently, when the diffusion coefficients \(a(x)\) are non-smooth, the resulting density profile \(p(x,t)\) often develops steep, shock-like discontinuities that are notoriously difficult to resolve on discrete meshes. Beyond capturing these sharp interfaces, any viable numerical algorithm must also strictly respect physical invariants, including total mass conservation and the pointwise non-negativity of the computed density field.

To simulate these probability fields, a variety of finite element techniques have been proposed in the literature; see [6]–[12] and the references therein. For instance, in [6], [7], the stationary version of the equation was treated using standard Galerkin finite element methods based on a weak formulation derived via classical integration by parts. To enforce the global integral constraint associated with probability preservation, the authors in [8] augmented a Galerkin scheme with a generalized Lagrange multiplier approach. In [9], conventional continuous \(C^0\) elements were applied to Fokker-Planck systems driven by both additive and multiplicative white noise processes, while a multiscale finite element architecture was introduced in [11] to analyze stochastic structural dynamics in high dimensions. A major limitation of these existing methodologies is their strict reliance on a highly smooth or constant diffusion tensor \(a(x)\), which is necessary to ensure the validity of the underlying weak forms. While a \(C^0\) continuous finite element framework was implemented and tested for the transient equation under non-smooth diffusion coefficients in [12], their work was purely empirical and lacked a formal mathematical convergence analysis.

When the diffusion coefficients possess sufficient regularity, the second-order differential operator can be expanded using the product rule to yield the following equivalent identity: \[\label{EQ:Reformulation} \frac{1}{2} \partial_{ij}^2(a_{ij}p) = \frac{1}{2} \partial_j (a_{ij}(x) \partial_i p) + \frac{1}{2} \partial_j ( (\partial_i a_{ij}) p).\tag{2}\] This algebraic restructuring effectively transforms 1 into a standard time-dependent convection-diffusion equation in divergence form. Setting aside high-dimensional bottlenecks, this reformulated system represents a relatively straightforward problem for conventional numerical methods. However, when the diffusion tensor is non-smooth, the expansion 2 breaks down completely due to the non-differentiability of the coefficients. In this rough regime, the exact solution \(p(x,t)\) exhibits unmapped internal discontinuities, making classical finite element approximations highly unstable or inapplicable. The primary objective of this paper is to construct a novel finite element framework that successfully circumvents this structural limitation, providing a stable, accurate, and rigorous approach specifically designed to handle non-smooth diffusion tensors in the Fokker-Planck equation.

To streamline the presentation, we consider a Fokker-Planck type model equation subject to a homogeneous Dirichlet boundary condition. The model problem seeks an unknown function \(u=u(x)\) satisfying \[\label{model} \begin{align} \nabla \cdot ({\boldsymbol{\mu}}u)-\frac{1}{2}\sum_{i,j=1}^d \partial^2_{ij}(a_{ij}u)=&f,\quad \text{in}\;\Omega,\\ u =& 0,\quad \text{on}\;\partial\Omega, \end{align}\tag{3}\] where \(\Omega\) is an open, bounded domain in \(\mathbb{R}^d (d=2, 3)\) with a Lipschitz continuous boundary \(\partial\Omega\), and \(f\in L^2(\Omega)\) is a prescribed source term. We assume 1) the diffusion tensor \(a(x)=\{a_{ij}(x)\}_{d\times d}\in (W_2^\infty (\Omega))^{d\times d}\) is symmetric, positive definite, and uniformly bounded across the domain \(\Omega\); 2) the drift vector satisfies the regularity condition \({\boldsymbol{\mu}}\in (W_1^\infty (\Omega))^d\).

The Weak Galerkin (WG) finite element method [13]–[38] provides a highly adaptable alternative to conventional continuous finite element frameworks. The core mechanism of the WG approach relies on the deployment of discrete weak derivatives, coupled with targeted stabilization forms that enforce weak continuity along element interfaces. This structural flexibility allows the WG methodology to seamlessly operate on arbitrary polygonal and polyhedral mesh configurations.

Building upon these foundational principles, the Primal-Dual Weak Galerkin (PDWG) method [39]–[53] was designed to cast numerical approximations as constrained optimization problems. The PDWG architecture delivers enhanced numerical stability, particularly when applied to non-self-adjoint systems such as linear transport models [54]. However, this theoretical advantage comes with a computational trade-off: the necessary inclusion of dual variables inherently expands the total number of global unknowns within the discrete system.

To address this limitation, the proposed Least-Squares Weak Galerkin (LS-WG) method offers several transformative advantages for the Fokker-Planck model problem 3 . The least-squares formulation naturally yields a symmetric positive-definite (SPD) discrete system. This property is a significant computational asset, permitting the direct application of high-performance iterative solvers such as the Conjugate Gradient (CG) method. By marrying the inherent stability of least-squares minimization with the geometric versatility of weak derivatives, the LS-WG approach provides a robust and scalable numerical platform.

In this work, we provide a rigorous theoretical foundation for the LS-WG scheme, establish the well-posedness and uniqueness of the discrete solution, and derive optimal-order error estimates in a discrete energy norm. Extensive numerical experiments are subsequently presented to validate these theoretical findings and demonstrate the overall robustness of the method.

We will follow the usual notation for Sobolev spaces and norms. For any open bounded domain \(D\subset \mathbb{R}^d\) with Lipschitz continuous boundary, we use \(\|\cdot\|_{s,D}\) and \(|\cdot|_{s,D}\) to denote the norm and seminorms in the Sobolev space \(H^s(D)\) for any \(s\ge 0\), respectively. The inner product in \(H^s(D)\) is denoted by \((\cdot,\cdot)_{s,D}\). The space \(H^0(D)\) coincides with \(L^2(D)\), for which the norm and the inner product are denoted by \(\|\cdot \|_{D}\) and \((\cdot,\cdot)_{D}\), respectively. When \(D=\Omega\), we shall drop the subscript \(D\) in the norm and inner product notation.

The paper is organized as follows. In Section 2, we shall discuss the computation of weak divergence and weak second order partial derivatives. In Section 3, we will present a detailed description of the least squares weak Galerkin finite element method for the Fokker-Planck type model problem and study the solution existence and uniqueness for our numerical method. In Section 4, we will establish an error estimate in a discrete norm. Finally, in Section 5, we report some numerical results to demonstrate the performance of the LS-WG finite element method.

2 Weak Differential Operators↩︎

The goal of this section is to formally introduce the continuous and discrete weak differential operators. To this end, let \(T\) be a polygonal or polyhedral element with boundary \(\partial T\). By a weak function on \(T\), we mean a triplet \(v=\{v_0,v_b,\mathbf{v}_g\}\) in which \(v_0\in L^2(T)\) represents the value in the interior of \(T\), \(v_b\in L^{2}(\partial T)\) represents the value on the boundary of \(T\), and \(\mathbf{v}_g=(v_{g1},\cdots,v_{gd})\in [L^{2}(\partial T)]^d\) represents the gradient of \(v\) on \(\partial T\). In general, the boundary components \(v_b\) and \(\mathbf{v}_g\) are not required to be equal to the traditional traces of \(v_0\) and \(\nabla v_0\) on \(\partial T\). We denote the space of all well-defined weak functions on \(T\) by \({\mathcal{W}}(T)\): \[\label{2461} {\mathcal{W}}(T)=\{v=\{v_0,v_b,\mathbf{v}_g\}: v_0\in L^2(T), v_b\in L^{2}(\partial T), \mathbf{v}_g\in [L^{2}(\partial T)]^d\}.\tag{4}\]

For any \(v\in {\mathcal{W}}(T)\), the weak second-order partial derivative of \(a_{ij}v\), denoted by \(\partial^2_{ij,w} (a_{ij}v)\), is defined as a linear functional in the dual space of \(H^2(T)\) satisfying \[\label{2463} (\partial^2_{ij,w}(a_{ij}v),\varphi)_T=(a_{ij}v_0,\partial^2_{ji}\varphi)_T- \langle a_{ij}v_b n_i,\partial_j\varphi\rangle_{\partial T}+ \langle a_{ij} v_{gi},\varphi n_j\rangle_{\partial T},\tag{5}\] for all \(\varphi\in H^2(T)\). Here, \(\mathbf{n}=(n_1,\cdots,n_d)\) denotes the unit outward normal vector to \(\partial T\), \((\cdot,\cdot)_T\) represents the standard \(L^2\) inner product in \(L^2(T)\), and \(\langle \cdot, \cdot\rangle_{\partial T}\) is the \(L^2\) inner product on the boundary \(\partial T\).

Analogously, for any \(v\in {\mathcal{W}}(T)\), the weak divergence of \({\boldsymbol{\mu}}v\), denoted by \(\nabla_w \cdot({\boldsymbol{\mu}}v)\), is defined as a linear functional in the dual space of \(H^1(T)\) satisfying \[\label{246495continuous} (\nabla_w \cdot({\boldsymbol{\mu}}v), \psi)_T = -({\boldsymbol{\mu}}v_0, \nabla \psi)_T + \langle ({\boldsymbol{\mu}}\cdot \mathbf{n}) v_b, \psi \rangle_{\partial T},\tag{6}\] for all \(\psi\in H^1(T)\).

Let \(P_r(T)\) denote the space of polynomials on \(T\) of degree no greater than \(r\). A discrete version of the weak second-order partial derivative \(\partial^2_{ij,w} (a_{ij}v)\) for \(v\in {\mathcal{W}}(T)\), denoted by \(\partial^2_{ij,w,r,T} (a_{ij}v)\), is defined as the unique polynomial in \(P_r(T)\) satisfying \[\label{2464} (\partial^2_{ij,w,r,T}(a_{ij}v),\varphi)_T=(a_{ij}v_0,\partial^2_{ji}\varphi)_T-\langle a_{ij}v_b n_i,\partial_j\varphi\rangle_{\partial T} +\langle a_{ij}v_{gi},\varphi n_j\rangle_{\partial T},\tag{7}\] for all \(\varphi \in P_r(T)\). Employing standard integration by parts, this definition yields the equivalent representation: \[\label{dissec} \begin{align} (\partial^2_{ij,w,r,T}(a_{ij}v),\varphi)_T&=(\partial^2 _{ ij}(a_{ij}v_0),\varphi)_T+\langle a_{ij} (v_0-v_b) n_i,\partial_j\varphi\rangle_{\partial T} \\&-\langle a_{ij}(\partial_i v_0-v_{gi}),\varphi n_j\rangle_{\partial T}, \end{align}\tag{8}\] for all \(\varphi \in P_r(T)\).

Correspondingly, the discrete weak divergence operator \(\nabla_{w,r,T} \cdot ({\boldsymbol{\mu}}v)\) for any \(v\in {\mathcal{W}}(T)\) is defined as the unique polynomial in \(P_r(T)\) that satisfies \[\label{disgradient} (\nabla_{w,r,T} \cdot ({\boldsymbol{\mu}}v), \psi)_T = -({\boldsymbol{\mu}}v_0, \nabla \psi)_T + \langle ({\boldsymbol{\mu}}\cdot \mathbf{n}) v_b, \psi \rangle_{\partial T}, \quad \forall \psi \in P_r(T).\tag{9}\] Applying integration by parts to the volume term in 9 provides the alternative formulation: \[\label{disdiv} (\nabla_{w,r,T} \cdot ({\boldsymbol{\mu}}v), \psi)_T = (\nabla \cdot ({\boldsymbol{\mu}}v_0), \psi)_T + \langle ({\boldsymbol{\mu}}\cdot \mathbf{n})(v_b - v_0), \psi \rangle_{\partial T}, \quad \forall \psi \in P_r(T).\tag{10}\]

3 Least Squares Weak Galerkin Scheme↩︎

The primary objective of this section is to formulate a least-squares weak Galerkin finite element method for the model problem 3 . To this end, let \(\mathcal{T}_h\) be a finite element partition of the domain \(\Omega\) consisting of polygons in \(\mathbb{R}^2\) or polyhedra in \(\mathbb{R}^3\), which is shape-regular in the sense described in [36]. We denote the set of all edges (in 2D) or faces (in 3D) in \(\mathcal{T}_h\) by \(\mathcal{E}_h\), and define the set of all interior edges or faces as \(\mathcal{E}_h^0=\mathcal{E}_h\setminus\partial\Omega\). For each element \(T\in\mathcal{T}_h\), \(h_T\) denotes its diameter, and the global mesh size of the partition is given by \(h=\max_{T\in\mathcal{T}_h}h_T\).

Let \(k\ge 1\) be a specified integer. On each element \(T\in\mathcal{T}_h\), we define the local discrete weak function space \(V(T)\) as: \[\begin{align} \label{V-T} V(T)= \{\{v_0, v_b, \mathbf{v}_g\}:\;v_0\in P_k(T),\;v_b\in P_k(e),\;\mathbf{v}_g\in [P_{k-1}(e)]^d,\;e\subset \partial T \}. \end{align}\tag{11}\] By patching \(V(T)\) over all elements \(T\in\mathcal{T}_h\) via common value of \(v_b\) on the interior edges or faces, we construct the global weak finite element space \(V_h\): \[V_h= \{\{v_0, v_b, \mathbf{v}_g\}:\;\{v_0, v_b, \mathbf{v}_g\}|_T\in V(T),\;\forall T\in\mathcal{T}_h \}.\] To incorporate the homogeneous Dirichlet boundary condition, we define \(V_h^0\) as the subspace of \(V_h\) with vanishing boundary value for \(v_b\) on \(\partial\Omega\): \[V_h^0= \{\{v_0, v_b, \mathbf{v}_g\}\in V_h:\;v_b|_e=0,\;\forall e\subset \partial\Omega \}.\]

For simplicity of notation, we denote the discrete weak divergence \(\nabla_{w,k-1,T}\cdot({\boldsymbol{\mu}}\sigma)\) computed using 9 on each element \(T\) with \(r=k-1\) simply by \(\nabla_w\cdot({\boldsymbol{\mu}}\sigma)\): \[(\nabla_w\cdot({\boldsymbol{\mu}}\sigma))|_T=\nabla_{w,k-1,T}\cdot({\boldsymbol{\mu}}\sigma|_T),\qquad \forall \sigma\in V_h.\] Analogously, \(\partial^2_{ij, w}(a_{ij}\sigma)\) denotes the discrete weak second-order partial derivative \(\partial^2_{ij,w,k-1,T}(a_{ij}\sigma)\) computed via 7 on each element \(T\) with \(r=s\) where \(s=k-1\) or \(s=k-2\) for \(k\geq 2\) and \(s=k-1\) for \(k=1\): \[(\partial^2_{ij, w}(a_{ij}\sigma))|_T=\partial^2_{ij,w,s,T}(a_{ij}\sigma|_T),\qquad \forall \sigma\in V_h.\]

Next, we introduce the following global bilinear forms on \(V_h\times V_h\): \[\begin{align} a(u, v) &= \sum_{T\in\mathcal{T}_h} a_T(u, v), \\ s(u, v) &= \sum_{T\in\mathcal{T}_h} s_T(u, v), \end{align}\] where the local bilinear forms are defined on each element \(T\) as \[a_T(u, v)= (\nabla_w\cdot({\boldsymbol{\mu}}u)-\frac{1}{2}\sum_{i,j=1}^d \partial^2_{ij, w}(a_{ij}u),\;\nabla_w\cdot({\boldsymbol{\mu}}v)-\frac{1}{2}\sum_{i,j=1}^d \partial^2_{ij, w}(a_{ij}v) )_T,\] and \[s_T(u, v)=h_T^{-3}\langle u_0-u_b, v_0-v_b\rangle_{\partial T} + h_T^{-1}\langle\nabla u_0-\mathbf{u}_g, \nabla v_0-\mathbf{v}_g\rangle_{\partial T}.\]

We are now in a position to state the least-squares weak Galerkin finite element scheme for the model problem 3 .

The well-posedness of this numerical scheme is established by the following lemma.

Finally, we introduce a semi-norm on the space \(V_h\) defined by: \[{|||} v{|||}^2=a(v, v) + s(v, v).\] Following the exact same algebraic arguments utilized in the proof of Lemma 1, one can easily verify that \({|||}\cdot{|||}\) defines a rigorous norm on the constrained subspace \(V_h^0\).

4 Error Estimates↩︎

In this section, we establish optimal-order error estimates for the proposed LS-WG finite element method. Throughout the analysis, the mesh partition \(\mathcal{T}_h\) is assumed to be shape-regular. For any element \(T \in \mathcal{T}_h\) and any function \(\phi \in H^1(T)\), we rely on the following standard continuous trace inequality [36]: \[\label{tracein} \|\phi\|^2_{\partial T} \leq C ( h_T^{-1}\|\phi\|_T^2 + h_T \|\nabla \phi\|_T^2 ).\tag{15}\]

On each element \(T \in \mathcal{T}_h\), let \(Q_0\) denote the standard \(L^2\) projection onto the polynomial space \(P_k(T)\) for \(k \geq 1\). Similarly, on each edge or face \(e \subset \partial T\), let \(Q_b\) and \(\mathbf{Q}_g := (Q_{g1}, \cdots, Q_{gd})\) denote the \(L^2\) projections onto \(P_k(e)\) and \([P_{k-1}(e)]^d\), respectively. For any sufficiently smooth function \(w \in H^2(\Omega)\), we define the global projection operator \(Q_h w \in V_h\) such that its restriction to each element \(T\) is given by: \[Q_h w = \{Q_0 w, Q_b w, \mathbf{Q}_g(\nabla w)\}.\]

For simplicity, we assume that the diffusion tensor \(a=\{a_{ij}\}_{d\times d}\) and the convection coefficient \({\boldsymbol{\mu}}\) are piecewise constant. The analysis readily extends to the case where \(a\) and \({\boldsymbol{\mu}}\) are non-smooth coefficients.

5 Numerical experiment↩︎

In the numerical test, we solve the model Fokker-Planck problem 3 , where \(\Omega=(0,1)\times(0,1)\), and \[\begin{align} \label{sol1} \begin{aligned} \;{\boldsymbol{\mu}}&=\begin{pmatrix}0\\0\end{pmatrix}, \;\;(a_{ij})=\begin{pmatrix}4&0 \\ 0& 4\end{pmatrix}, \;\;u = 2^8(x-x^2)^2(y-y^2)^2, \end{aligned} \end{align}\tag{20}\] or \[\begin{align} \label{sol2} \begin{aligned} \;{\boldsymbol{\mu}}&=\begin{pmatrix}1\\1\end{pmatrix}, \;\;(a_{ij})=\begin{pmatrix}4&2 \\ 2& 4\end{pmatrix}, \;\;u = 2^8(x-x^2)^2(y-y^2)^2. \end{aligned} \end{align}\tag{21}\]

The computations are done on four types of grids, ordered by their level of difficulty: the square grids shown in Figure 2, the triangular grids shown in Figure 3, the quadrilateral grids shown in Figure 4, and the non-convex pentagonal grids shown in Figure 5. We apply the least squares weak Galerkin method [ls] using \(P_{k}\) polynomial spaces (\(k=2, 3, 4\)) in 11 to solve the two problems. The results are listed in Tables 1-8, which show that the optimal orders of convergence are achieved. It appears that the solutions for 20 are consistently more accurate than those for 21 .

References↩︎