B. Leimkuhler,A. Sharma,M.V. Tretyakov
Abstract:We derive and analyze numerical methods for weak approximation of underdamped (kinetic) Langevin dynamics in bounded domains. First-order methods are based on an Euler-type scheme interlaced with collisions with the boundary. To achieve second order, composition schemes are derived based on decomposition of the generator into collisional drift, impulse, and stochastic momentum evolution. In a deterministic setting, this approach would typically lead to first-order approximation, even in symmetric compositions, but we find that the stochastic method can provide second-order weak approximation with a single gradient evaluation, both at finite times and in the ergodic limit. We provide theoretical and numerical justification for this observation using model problems and compare and contrast the numerical performance of different choices of the ordering of the terms in the splitting scheme.
What problem does this paper attempt to address?
This paper attempts to solve the problem of numerical simulation of restricted Langevin dynamics (i.e., stochastic differential equations with reflective boundaries) in a bounded domain. Specifically, the article mainly focuses on the following issues:
1. **Developing numerical methods**: In order to weakly approximate the weakly damped (or kinetic energy) Langevin dynamics in a bounded domain, the author proposes first - order and second - order numerical methods. These methods are especially suitable for situations where boundary collisions need to be considered.
2. **Handling boundary collisions**: For elastic collisions that occur when particles reach the boundary, the author designs a special numerical scheme to handle this phenomenon, ensuring that the orthogonal component of momentum is reversed while the tangential component remains unchanged.
3. **Improving computational accuracy**: By decomposing the generation operator into parts such as collision drift, impulse, and stochastic momentum evolution, the author proposes a combined scheme to achieve a second - order weak approximation. This method is effective in both finite - time and ergodic limits and requires only one gradient evaluation.
4. **Theoretical and numerical verification**: The author not only provides theoretical proofs but also conducts numerical verification through model problems, comparing the influence of different term - order selections in different splitting schemes on numerical performance.
### Mathematical formula summary
- **Constrained Langevin dynamics (CLD)**:
\[
Q(t)=Q(0)+\int_{0}^{t} P(s) \, ds,
\]
\[
P(t)=P(0)+\int_{0}^{t} b(P(s), Q(s)) \, ds+\sigma W(t)+R(t),
\]
where \(R(t)=-\sum_{0 < s \leq t} 2 (P(s^-) \cdot n_G(Q(s))) n_G(Q(s)) I_{\partial G}(Q(s)), \quad t \geq 0\), and \(W(t)\) is a \(d\)-dimensional standard Wiener process.
- **ECLD in specific cases**:
\[
Q(t)=Q(0)+\int_{0}^{t} P(s) \, ds,
\]
\[
P(t)=P(0)-\int_{0}^{t} \nabla U(Q(s)) \, ds-\gamma \int_{0}^{t} P(s) \, ds+\sqrt{\frac{2\gamma}{\beta}} W(t)+R(t),
\]
where \(R(t)=-\sum_{0 < s \leq t} 2 (P(s^-) \cdot n_G(Q(s))) n_G(Q(s)) I_{\partial G}(Q(s))\), \(\beta > 0\) represents the inverse temperature in molecular dynamics, and \(\gamma > 0\) is the friction parameter.
- **Ergodicity**:
The equation (ECLD) is ergodic and has an invariant density:
\[
\rho(q, p)=\frac{1}{Z} e^{-\beta(U(q)+|p|^2/2)}, \quad (q, p) \in \bar{G} \times \mathbb{R}^d.
\]
### Application background
The main motivation of this research comes from the requirements of multiple application fields, such as:
- **Sampling**: Sampling from a given distribution is an important task in computational statistics and has wide applications especially in machine learning.
- **Molecular dynamics**: Studying the molecular properties of substances at a constant temperature.
- **Computational fluid mechanics**: Simulating the behavior of discrete particles in turbulent flow in a confined domain.
- **Optimization**: Minimizing the loss function under certain constraint conditions.
- **Solving mirror boundary value problems**: Providing a stochastic representation of linear parabolic partial differential equations through the Feynman - Kac formula.
In summary, this paper aims to provide efficient and accurate numerical simulation methods for constrained Langevin dynamics, especially in cases involving boundary collisions.