Numerical solution of the Fokker–Planck equation using physics-based mixture models

https://doi.org/10.1016/j.cma.2022.115424Get rights and content

Abstract

The Fokker–Planck equation governs the uncertainty propagation of dynamical systems driven by stochastic processes. The solution of the Fokker–Planck equation is a time-varying Probability Density Function (PDF) that is usually of high dimension with unbounded support. There are also properties that must be conserved in time for the joint and any marginal solution PDF (i.e., a time-varying PDF is a non-negative function that must integrate to unity at any given time.) Satisfying these properties poses a significant challenge to the numerical solution of the Fokker–Planck equation. Another challenge is to capture the tail behavior of the solution PDF with unbounded support, required to predict the probability of rare events like a system failure. Satisfying the conditions of a proper PDF with unbounded support limits the application of traditional grid-based numerical methods, like finite difference and finite element methods. This paper develops a novel numerical method based on physics-based mixture models for the transient and steady-state solutions of the Fokker–Planck equation. In the proposed numerical method, the trial function space consists of the convex combinations of some parametric PDFs (i.e., mixture models) that trivially satisfy the necessary conditions of the solution. Estimating the unknown parameters of the mixture model is via Bayesian inference while considering the constraints on model parameters. Bayesian inference facilitates integrating data on the responses of dynamical systems (e.g., from simulations) with the governing Fokker–Planck equation to estimate the unknown parameters. Since the solution PDF is not observable, combining it with observable data on system responses is far from trivial based on current numerical methods. The formulation of Bayesian inference also introduces a weighting function to reduce the estimation error in specific regions of interest like the tail of the solution PDF. To reduce the computational demand, the paper develops an importance sampling algorithm that generates a small set of collocation points at which the residual of the Fokker–Planck equation is evaluated. The performance of the proposed numerical method is demonstrated with several benchmark problems. The obtained results are compared with analytical solutions when available and otherwise with simulations.

Introduction

The Fokker–Planck equation governs the time-varying response probability density function (PDF) of dynamical systems driven by stochastic processes [1]. It arises in many areas of science and engineering, ranging from statistical physics [2], [3], [4] and chemistry [5], [6] to system biology [7], [8], mathematical finance [9], and structural dynamics [10], [11]. The Fokker–Planck equation of a dynamical system with n state variables is a parabolic partial differential equation (PDE) of dimension n+1 (i.e., n spatial dimensions plus time). Finding the transient and steady-state solutions of the Fokker–Planck equation is a daunting task. In particular, the Fokker–Planck equation of realistic dynamical systems usually involves many state variables leading to a high-dimensional PDE. There are also properties of the solution that must be conserved in time. The solution of the Fokker–Planck equation is a time-varying PDF, a non-negative function that must integrate to unity at any given time. These properties apply to the joint and any marginal PDF of the system responses. Another significant challenge is to capture the tail behavior of the response PDF with unbounded support, required to predict the probability of rare events like a system failure. The analytical solution of the Fokker–Planck equation is available only for a few idealized dynamical systems under restrictive conditions [1], [12].

Since there is no general analytical solution to the Fokker–Planck equation, various numerical methods and simulation techniques have been developed. Notable examples of numerical methods include the variants of finite element [13], [14], [15], finite difference [16], [17], [18], and path integration [19], [20], [21] methods. The general idea of these numerical methods is to discretize the computational domain into lattices or grids and develop an interpolation function that approximates the solution at grid points. However, the computational demand of these numerical methods makes them impractical for dynamical systems with dimensions larger than three [16] and capturing the tail of the response PDF with unbounded support. This is because their computational demand is directly related to the grid size, which rapidly grows with the extent and dimension of the computational domain. Recent accounts of different numerical methods to solve kinetic equations, including the Fokker–Planck equation, can be found in [22], [23]. More recently, physics-informed neural networks [24], [25], [26] have been used to solve various PDEs, including the Fokker–Planck equation. The general idea is to approximate the solution of the PDE by a neural network and find its parameters via optimization. The objective function enforces the neural network to match observable data and satisfy the governing equation and its initial and boundary conditions at a set of prescribed collocation points in the computational domain. The numerical solution based on neural networks can address the scaling issues of traditional PDE solvers, like finite element and finite difference. It has been used to solve specific parabolic PDEs on domains of dimension O102 [24], [27]. However, satisfying the necessary conditions of a proper PDF, discussed earlier, and capturing its tail behavior in the Fokker–Planck equation remains a significant challenge [28], [29], [30]. Alternatively, the variants of Monte Carlo simulation technique generate sample paths of the stochastic response process from which various response statistics can be estimated [31], [32], [33]. Though straightforward to implement, simulation techniques are greatly limited by their computational complexity (see, for example, [34] for more discussion on computational complexity.)

In this paper, we develop a novel numerical method that approximates the transient and steady-state solutions of the Fokker–Planck equation by a mixture model via Bayesian inference. Since the mixture model is the convex combination of some parametric PDFs, it trivially satisfies the necessary conditions of the solution. Applying the time integration scheme to the strong form of the Fokker–Planck equation, we formulate a weighted regression problem to estimate the unknown parameters of the mixture model. The graph-based automatic differentiation technique [35], [36] greatly facilitates evaluating the partial derivatives in the strong form of the Fokker–Planck equation. Using Bayesian inference, we then estimate the unknown parameters of the mixture model while enforcing required constraints on their values. Examples of such constraints include the convexity of the mixture model and the symmetry and positive-definiteness of covariance matrices in the Gaussian mixture model. Bayesian inference also allows us to incorporate any other information like simulated responses of the dynamical system in estimating the unknown model parameters. We introduce an importance sampling algorithm to generate the collocation points at which we evaluate the residual of the Fokker–Planck equation. The proposed adaptive approach for generating the collocation points at each time step can significantly reduce the computational demand of Bayesian inference in high-dimensional problems. The specific design of the importance sampling density also allows us to improve estimating the tail of the response PDF [32], [37], [38]. As a proof of concept, we illustrate the capabilities of the proposed physics-based mixture model by solving the Fokker–Planck equations of several dynamical systems.

The rest of the paper consists of four sections. Section 2 provides a formal definition of the problem. Section 3 discusses the proposed physics-based mixture model for the numerical solution of the Fokker–Planck equation. Section 4 explains the proposed formulation through numerical examples. Finally, the last section summarizes the paper and draws some conclusions.

Section snippets

The Fokker–Planck equation and its solution

Consider a general (nonlinear) dynamical system defined by the stochastic differential equation dyt=ayt,tdt+Byt,tdwt,with initial condition y0YRn, where yt is the state vector; a, is the drift vector; B, is an Rn×m-valued matrix that yields the diffusion matrix D=BB; and wt is the Rm-valued standard Wiener process. The transition probability density function (PDF) fy,t of the Markov process yt satisfies the Fokker–Planck equation fy,tt=ay,tfy,t+12Dy,tfy,t=LFPfy,t,with boundary

Proposed physics-based mixture model to find the solution of the Fokker–Planck equation

Physics-based models are a class of statistical models that directly satisfy the governing physics in addition to learning model parameters from data. This section presents the proposed physics-based mixture models that satisfy the Fokker–Planck equation of a given dynamical system and can learn model parameters from any data on system responses. For illustration purposes, we first explain the general idea by considering the numerical solution of the Fokker–Planck equation for an equivalent

Numerical examples

This section illustrates the proposed physics-based mixture model to solve the Fokker–Planck equations of three dynamical systems. The selected examples consist of two nonlinear systems (one elastic and the other hysteretic) driven by Gaussian white-noise processes and a nonlinear stochastic climate model. These examples illustrate the feasibility of the proposed approach to solve the Fokker–Planck equation of dynamical systems with diverse nonlinear behavior. We use Gaussian mixture models in

Conclusions

This paper developed a novel numerical method for the solution of the Fokker–Planck equation. The main idea is to design the trial function space that 1) facilitates satisfying the necessary conditions of the solution function and 2) can be expanded to ensure convergence to the actual solution. Since the solution of the Fokker–Planck equation is a time-varying probability density function (PDF), the trial functions must be non-negative and integrate to unity at any given time. To meet the

Declaration of Competing Interest

The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.

References (59)

  • PeskovN.

    Finite element solution of the Fokker-Planck equation for single domain particles

    Physica B

    (2020)
  • QianY. et al.

    A conservative, free energy dissipating, and positivity preserving finite difference scheme for multi-dimensional nonlocal Fokker-Planck equation

    J. Comput. Phys.

    (2019)
  • UreñaF. et al.

    Non-linear Fokker-Planck equation solved with generalized finite differences in 2D and 3D

    Appl. Math. Comput.

    (2020)
  • YuJ. et al.

    A new path integration procedure based on Gauss-Legendre scheme

    Int. J. Non-Linear Mech.

    (1997)
  • NaessA. et al.

    Efficient path integration methods for nonlinear dynamic systems

    Probab. Eng. Mech.

    (2000)
  • ChoH. et al.

    Numerical methods for high-dimensional probability density function equations

    J. Comput. Phys.

    (2016)
  • SirignanoJ. et al.

    DGM: A deep learning algorithm for solving partial differential equations

    J. Comput. Phys.

    (2018)
  • RaissiM. et al.

    Physics-informed neural networks: A deep learning framework for solving forward and inverse problems involving nonlinear partial differential equations

    J. Comput. Phys.

    (2019)
  • MengX. et al.

    PPINN: Parareal physics-informed neural network for time-dependent PDEs

    Comput. Methods Appl. Mech. Engrg.

    (2020)
  • MackeM. et al.

    Importance sampling for randomly excited dynamical systems

    J. Sound Vib.

    (2003)
  • NatarajanH. et al.

    A high-order semi-Lagrangian method for the consistent Monte-Carlo solution of stochastic Lagrangian drift–diffusion models coupled with Eulerian discontinuous spectral element method

    Comput. Methods Appl. Mech. Engrg.

    (2021)
  • PeherstorferB. et al.

    Multifidelity importance sampling

    Comput. Methods Appl. Mech. Engrg.

    (2016)
  • TabandehA. et al.

    A review and assessment of importance sampling methods for reliability analysis

    Struct. Saf.

    (2022)
  • CaugheyT. et al.

    The exact steady-state solution of a class of non-linear stochastic systems

    Int. J. Non-Linear Mech.

    (1982)
  • HornikK. et al.

    Multilayer feedforward networks are universal approximators

    Neural Netw.

    (1989)
  • Di PaolaM. et al.

    Differential moment equations of FE modelled structures with geometrical non-linearities

    Int. J. Non-Linear Mech.

    (1990)
  • TabandehA. et al.

    Nonlinear random vibration analysis: A Bayesian nonparametric approach

    Probab. Eng. Mech.

    (2021)
  • ChenN. et al.

    Efficient statistically accurate algorithms for the Fokker-Planck equation in large dimensions

    J. Comput. Phys.

    (2018)
  • RiskenH. et al.

    The Fokker-Planck Equation: Methods of Solution and Applications

    (1996)
  • Cited by (8)

    • SPARSE–R: A point-cloud tracer with random forcing

      2024, International Journal of Multiphase Flow
    View all citing articles on Scopus
    View full text