# Numerical integrators for confined Langevin dynamics

B. Leimkuhler School of Mathematics and the Maxwell Institute for Mathematical
Sciences, University of Edinburgh, UK<EMAIL_ADDRESS>A. Sharma
Department of Mathematical Sciences, Chalmers University of Technology and the
University of Gothenburg, Sweden<EMAIL_ADDRESS>M.V. Tretyakov School of
Mathematical Sciences, University of Nottingham, UK;
<EMAIL_ADDRESS>
###### 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.

Keywords:  stochastic differential equations with reflection, weak
approximation, sampling with constraints, computing ergodic limits, splitting
schemes.

AMS Classification:  65C30, 60H35, 60H10, 37H10.

##  1 Introduction

There is great current demand for reliable numerical methods for solving
stochastic differential equations (SDEs). In part this can be attributed to
their role in large scale statistical computation, where long stochastic paths
are often used to calculate ergodic averages. Many standard algorithms for
SDEs are formulated on unbounded domains, e.g. Euclidean space, but in
practice are subject to inequality constraints (e.g. positivity relations) or
are employed with more complicated domain restrictions, often implemented in
an ad hoc manner. In this article, we explore the rational design of numerical
algorithms for accurate long-term simulation of Langevin SDEs with domain
restriction, including sampling from stationary measures with compact support.

Let $G\subset\mathbb{R}^{d}$ be a bounded domain with sufficiently smooth
boundary $\partial G$, $n_{G}(q)$ be the outward normal unit vector at
$q\in\partial G$, and $I_{\partial G}(q)$ be the indicator function of
$q\in\partial G$. Let $\sigma$ be a positive constant,
$b:\bar{G}\times\mathbb{R}^{d}\rightarrow\mathbb{R}^{d}$ satisfying certain
conditions, and ($\cdot$) denote scalar product. Consider the confined
Langevin dynamics (CLD) governed by the stochastic differential equations
(SDEs):

$\displaystyle Q(t)=Q(0)+$
$\displaystyle\int_{0}^{t}P(s)ds,\;\;\;\;P(t)=P(0)+\int_{0}^{t}b(P(s),Q(s))ds+\sigma
W(t)+\mathcal{R}(t),$ (CLD) $\displaystyle\mathcal{R}(t)$
$\displaystyle=-\sum\limits_{0<s\leq t}2(P({s^{-}})\cdot
n_{G}(Q(s)))n_{G}(Q(s))I_{\partial G}(Q(s)),\;\;\;t\geq 0,$

where $W(t)$ is a $d$-dimensional standard Wiener process.

The process $Q(t)$ does not leave the set $\bar{G}$ and it elastically
reflects from the boundary $\partial G$ [44, 5, 6] which is achieved due to
the cádlág process $\mathcal{R}(t)$. At the collision of the position $Q(t)$
with $\partial G$, $\mathcal{R}(t)$ acts on the momentum $P(t)$ to reverse the
direction of its component orthogonal to $\partial G$ keeping the tangential
component unchanged.

If $b(q,p):=-\nabla U(q)-\gamma p$ and $\sigma=\sqrt{\frac{2\gamma}{\beta}}$,
i.e., consider the particular case of confined Langevin dynamics (CLD) defined
in the set $\bar{G}\times\mathbb{R}^{d}$:

$\displaystyle Q(t)$ $\displaystyle=Q(0)+\int_{0}^{t}P(s)ds,$ (ECLD)
$\displaystyle P(t)$ $\displaystyle=P(0)-\int_{0}^{t}\nabla
U(Q(s))ds-\gamma\int_{0}^{t}P(s)ds+\sqrt{\frac{2\gamma}{\beta}}W(t)+\mathcal{R}(t),$
$\displaystyle\mathcal{R}(t)$ $\displaystyle=-\sum\limits_{0<s\leq
t}2\big{(}P(s^{-})\cdot n_{G}(Q(s))\big{)}n_{G}(Q(s))I_{\partial G}(Q(s)),$

where $\beta>0$ is interpreted as the inverse temperature in molecular
dynamics and $\gamma>0$ is a friction parameter. Eq. (ECLD) is ergodic with
the invariant density $\rho(q,p)=\frac{1}{Z}e^{-\beta(U(q)+|p|^{2}/2)}$,
$(q,p)\in\bar{G}\times\mathbb{R}^{d}$. The central theme of this paper is the
efficient weak-sense numerical approximation of (CLD) including the
computation of ergodic limits associated with (ECLD) which is motivated by
applications such as those discussed below.

Sampling. Sampling from a desired distribution is of prime importance for
computational statistics and it has found a large number of machine learning
applications. Many scenarios involve sampling in constrained spaces, for
example, truncated data problems arising in survival time studies [23],
ordinal data models [21], constrained Lasso and Ridge regression [11], Latent
Dirichlet Allocation [2], non-negative matrix factorization [42], and
stochastic optimal control [19]. In [27] (see also references therein) the
task of sampling from a given distribution with compact support is addressed
using reflected Brownian dynamics (in other words, overdamped Langevin
dynamics with reflection), where the dynamics are confined within $\bar{G}$
thanks to making use of the local time. In the case of the whole space,
advantages of underdamped Langevin equations vs overdamped Langevin dynamics
are well studied (see e.g. [26, 14]) which are also expected for the
considered confined-space case.

Molecular dynamics. To study molecular properties of a given substance under
constant temperature $1/\beta$, Langevin dynamics are often employed in which
$N$ particles (atoms or molecules) with positions $q$ move under the action of
a potential energy function $U(q)$. These properties can be expressed as
ergodic limits associated with Langevin dynamics (see [3, 29, 25] for more
details). To perform such simulations in the case of the dynamics being
constrained to a confined domain, one typically uses a steep (often singular)
barrier in the potential $U(q)$ preventing particles from leaving the confined
space. The barrier leads to the drift $-\nabla U(q)$ in the Langevin equations
increasing rapidly near the boundary of the confined space which can
negatively affect stability of numerical methods, or require the use of a
variable stepsize [30]. Confined Langevin dynamics (ECLD) serves the purpose
of simulating molecular configurations in a restricted region, with atoms
undergoing perfectly elastic collisions at the walls, without causing
stability problems or requiring stepsize restriction, in contrast to the case
of using a steep prohibiting barrier. Furthermore, confined Langevin dynamics
offers a natural physical behavior of particles (elastic collisions) at the
boundary. In the kinetic theory of gases, for example, specular reflection
models are used to describe reflection of gas particles at a totally elastic
solid surface [12].

Computational fluid dynamics. The SDEs (CLD) find applications in modelling
discrete particles immersed in turbulent flows inside a confined domain, where
colloidal particles are carried by the fluid flow, experiencing a drag force
and collisions with the domain boundary. Numerical experiments using models of
(CLD) type allow to study, e.g., convection, polydispersity etc. related to
mesoscopic particles suspended in turbulent flow (see e.g., [15, 32, 10] for
more details).

Optimization. Minimizing a loss/error function is at the core of machine
learning. In many cases, there are certain constraints which are to be
satisfied as a requirement of the model, e.g. positivity constraints. Also,
constraining the parameter set can result in improvement of the performance of
the model due to regularization effects of constraints, e.g. [43, 28]. One can
exploit (ECLD) to minimize a function $U(q)$, where $q$ satisfies constraints
represented by domain $G$. Moreover, one can also consider replacing $\nabla
U$ with a noisy gradient $\tilde{\nabla}U$ to introduce batching of data for
machine learning specific tasks.

Solving specular boundary value problem. Via the Feynman-Kac formula [6, 7],
the SDEs (CLD) provide a probabilistic representation of the specular boundary
value problem for linear parabolic partial differential equations (PDEs) (see
(2.2)-(2.4)). Solving such PDEs in high dimensions using deterministic methods
is computationally infeasible. We can utilise the proposed here numerical
methods for (CLD) together with the Monte-Carlo technique to solve such high-
dimensional problems.

The above applications give strong motivation to construct and study numerical
methods for (CLD). Despite this, very limited attention has been paid so far
to numerical integration of confined Langevin dynamics. In this respect we
note that in [7] the authors investigated a numerical scheme when
$G=\mathbb{R}^{d-1}\times(0,\infty)$ and proved its first order weak
convergence.

In this paper (Section 3), we derive first-order weak methods to numerically
approximate (CLD) as well as second-order schemes in the case of (ECLD). To
construct the splitting schemes, we decompose the generator of the Markov
process $(Q(t),P(t))$ into components representing collisional drift, impulse,
and stochastic momentum evolution. In Section 3 we also state theoretical
results on error bounds for the methods, both when the SDEs are integrated
over a finite time interval and when associated ergodic limits are
approximated. Proofs of the error bounds are given in Section 4 with necessary
background material for (CLD) and (ECLD) provided in Section 2. The main
difficulty in proving convergence arises due to the need to quantify the
number of steps with collision. We note the following remarkable fact about
second-order symmetric integrators for (ECLD). In a deterministic setting, the
discussed splitting approach would typically lead to first-order
approximation, even in symmetric compositions [18, 25]. Then stochastic
numerics experience [39] suggests that finite-time convergence of a
generalization of such deterministic methods to (ECLD) should be of weak order
one as well. However, it turns out that the stochastic method with a single
gradient evaluation has second-order convergence. Simply put, this
counterintuitive and remarkable result is due to the average collision time
(within a step of the numerical method in which a collision is encountered)
being at the midpoint of the step up to terms of higher order. In Section 5 we
present results of numerical tests which confirm our theoretical predictions
and compare the methods from Section 3. In Section 6 we summarize this work
and suggest topics of interest for further research related to confined
Langevin dynamics.

##  2 Preliminaries

Notation.  For the boundary sets we write: $\Sigma^{+}=\\{(q,p)\in\partial
G\times\mathbb{R}^{d};(p\cdot n_{G}(q))>0\\}$ and
$\Sigma_{t}^{+}=(0,t)\times\Sigma^{+}$. The space
$C^{l_{1},l_{2},l_{3}}([0,T]\times A\times\mathbb{R}^{d})$ denotes the space
of functions whose continuous derivatives exist with respect to $t\in[0,T]$,
$q\in A$ and $p\in\mathbb{R}^{d}$ up to order $l_{1}$, $l_{2}$ and $l_{3}$,
respectively; $C^{l_{1},l_{2}}(A\times\mathbb{R}^{d})$ has the same meaning
but with respect to $q$ and $p$; $C_{b}^{l_{1},l_{2}}(A\times\mathbb{R}^{d})$
is a space of functions in $C^{l_{1},l_{2}}(A\times\mathbb{R}^{d})$ with
bounded derivatives up to order $l_{1}$ in $q$ and $l_{2}$ in $p$; and
$C_{c}^{l_{1},l_{2}}(A\times\mathbb{R}^{d})$ are functions in
$C^{l_{1},l_{2}}(A\times\mathbb{R}^{d})$ with compact support. As usual,
$|\cdot|$ denotes the $L^{2}$-norm of a vector in $\mathbb{R}^{d}$. For a
differentiable vector-valued function $F(x)$ from
$\bar{G}\subset\mathbb{R}^{d}$ or $\mathbb{R}^{d}$ to $\mathbb{R}^{d}$,
$\textbf{J}_{F}(x)$ denotes the Jacobian matrix of $F$ evaluated at $x$.

To prove the first-order weak-sense convergence of the schemes proposed in
Section 3.1, we need assumptions on the boundary $\partial G$ of $G$ and the
drift coefficient $b(q,p)$.

######  Assumption 2.1.

We assume the boundary $\partial G$ of set $G$ belongs to $C^{4}$.

######  Assumption 2.2.

The coefficient $b(q,p)$ is a $C^{2,2}(\bar{G}\times\mathbb{R}^{d})$ function
which grows at most linearly in $p$ at infinity, i.e. there exists a constant
$K>0$ such that,

$\displaystyle|b(q,p)|\leq K(1+|p|),$ (2.1)

for all $p\in\mathbb{R}^{d}$ and $q\in\bar{G}$.

In [44], existence of weak solutions and uniqueness in law of (CLD) in
$G=(0,\infty)\times\mathbb{R}^{d-1}$ with multiplicative noise and
$b(q,p)=b(q)-\gamma p$, $\gamma>0$, were established under boundedness
assumptions on $b(q)$ and the diffusion coefficient $\sigma(q)$. In [5],
existence of weak solution of (CLD) in $G=(0,\infty)\times\mathbb{R}^{d-1}$
with bounded drift coefficient $b$ was shown. This well-posedness result was
extended in [6] to the case of $\partial G\in C^{3}$ being a compact
submanifold of $\mathbb{R}^{d}$, by establishing the existence of a weak
solution and its pathwise uniqueness for $b(q,p)\equiv 0$, and then, by using
Girsanov’s theorem, for any bounded $b(q,p)$. We note that Girsanov’s theorem
can be applied if $b(q,p)$ has linear growth in $p$ (see [22, Proposition
5.3.6]) as done in [44]. Pathwise uniqueness of the solution of (CLD) with
$b(q,p)\equiv 0$ is proved in [47] when $G$ is a convex polytope.

###  2.1 Backward Kolmogorov equation for (CLD)

Consider the backward Kolmogorov equation

$\displaystyle\frac{\partial u}{\partial
t}+\frac{\sigma^{2}}{2}\Delta_{p}u+(p\cdot\nabla_{q}u)+(b(q,p)\cdot\nabla_{p}u)=0,\;\;\;\;(t,q,p)\in(0,T)\times
G\times\mathbb{R}^{d},$ (2.2)

with terminal condition

$u(T,q,p)=\varphi(q,p),\;\;\;\;(q,p)\in G\times\mathbb{R}^{d},$ (2.3)

and specular boundary condition

$u(t,q,p)=u(t,q,p-2(p\cdot n_{G})n_{G}),\;\;\;\;(t,q,p)\in\Sigma_{T}^{+},$
(2.4)

where $\Delta_{p}$ is the Laplacian with respect to $p$, $\nabla_{q}$ and
$\nabla_{p}$ are gradients with respect to $q$ and $p$, respectively, and
$n_{G}:=n_{G}(q)$.

The solution of (2.2)-(2.4) has the probabilistic representation [6, 7]:

$u(t_{0},q_{0},p_{0})=\mathbb{E}\varphi(Q_{t_{0},q_{0},p_{0}}(T),P_{t_{0},q_{0},p_{0}}(T)),$
(2.5)

where $(Q_{t_{0},q_{0},p_{0}}(t),P_{t_{0},q_{0},p_{0}}(t))$, $t\geq t_{0}$,
solves (CLD) with the starting point $(t_{0},q_{0},p_{0})$.

Below we give three examples of the specular boundary value problem
(2.2)-(2.4), for which we can obtain solutions in an explicit form.

######  Example 2.1.

Take $G=(0,\infty)$ and $\partial G=\\{0\\}$,
$b(q,p)=\Big{(}\frac{1}{q+c}+p\Big{)}$, $\sigma=1$,
$\varphi(q,p)=(q+c)^{2}e^{-p^{2}}$, and boundary condition
$u(t,q,p)=u(t,q,-p)$ at $\Sigma_{T}^{+}$, where $c\geq 0$. The solution is
given by $u(t,q,p)=(q+c)^{2}e^{-p^{2}}e^{-(T-t)}$.

######  Example 2.2.

Let $G\subset\mathbb{R}^{d}$, $U(q)\in C^{1}(\bar{G})$ and $\alpha,\;\beta>0$.
Take $b(q,p)=-\nabla_{q}U(q)+\alpha p$, $\sigma=\sqrt{2\alpha/\beta}$,
$\varphi=e^{-\beta\big{(}|p|^{2}/2+U(q)\big{)}}$ and boundary condition is
$u(t,q,p)=u(t,q,p-2(p\cdot n_{G})n_{G}),(t,q,p)\in\Sigma_{T}^{+}$. The
solution is $u(t,q,p)=\exp{\Big{(}-\beta\big{(}|p|^{2}/2+U(q)\big{)}-\alpha
d(T-t)\Big{)}}$.

######  Example 2.3.

Let $G\subset\mathbb{R}^{d}$,
$b(q,p)=-\frac{q\cos{|q|^{2}}+2(d+2)p}{2|p|^{2}+1}$, $\sigma=\sqrt{2}$, and
$\varphi(q,p)=|p|^{4}+|p|^{2}+\sin{|q|^{2}}$. The boundary condition is
$u(t,q,p)=u(t,q,p-2(p\cdot n_{G})n_{G}),\;(t,q,p)\in\Sigma_{T}^{+}$. Then, the
solution is $u(t,q,p)=|p|^{4}+|p|^{2}+\sin{|q|^{2}}+2d(T-t)$.

We make the following assumptions on the terminal function $\varphi$ and the
solution of (2.2)-(2.4).

######  Assumption 2.3.

$\varphi(q,p)\in C^{4,4}(\bar{G}\times\mathbb{R}^{d})$.

######  Assumption 2.4.

(i) The solution of (2.2)-(2.4) exists and belongs to $C^{2,2,4}((0,T)\times
G\times\mathbb{R}^{d})$ and
$C^{0,1,0}([0,T]\times\bar{G}\times\mathbb{R}^{d})$, and Hessian derivative of
$u(t,q,p)$ with respect to $q$ exists on
$[0,T]\times\bar{G}\times\mathbb{R}^{d}$. (ii) The solution $u(t,q,p)$ and all
its existing derivatives have growth at most polynomially in $p$.

We note that there is a lack of regularity results for solutions of the
specular boundary value problem (2.2)-(2.4) available in the literature. The
above three examples demonstrate that there is a sufficiently large class of
specular boundary value problems whose solutions satisfy Assumption 2.4. Also,
it was shown in [7] that if $\varphi\in C_{c}^{1,1}(G\times\mathbb{R}^{d})$
and $b\in C_{b}^{1,1}(\mathbb{R}^{d},\mathbb{R}^{d})$, then $u(t,q,p)\in
C([0,T]\times\bar{G}\times\mathbb{R}^{d})$ and $\nabla_{q}u$ and $\nabla_{p}u$
exist and belong to $C((0,T)\times\bar{G}\times\mathbb{R}^{d})$.

###  2.2 Ergodicity

Consider (ECLD) which is a particular case of (CLD).

######  Assumption 2.5.

$U(q)\in C^{\infty}(\bar{G})$ and $U(q)>1$ for all $q\in\bar{G}$.

Exploiting the approach used in proving Theorem 3.2 in [34, Section 3] (see
also [45]) one can establish under Assumptions 2.1 and 2.5 (note that $\nabla
U(q)$ is bounded thanks to $G$ being bounded) that (ECLD) is exponentially
ergodic, i.e. there exists a unique invariant measure $\mu(dq,dp)$ defined on
$G\times\mathbb{R}^{d}$ and the following inequality holds for any
$\varphi(q,p)$, $(q,p)\in\bar{G}\times\mathbb{R}^{d}$, bounded in $q$ and with
not faster than polynomial growth in $p$ and for any initial point
$(q_{0},p_{0})$ and all $T>0$:

$|\mathbb{E}(\varphi(Q(T),P(T)))-\bar{\varphi}|\leq Ke^{-\lambda T},$ (2.6)

where

$\bar{\varphi}=\int_{\mathbb{R}^{d}}\int_{G}\varphi(q,p)\mu(dq,dp).$

Moreover, the invariant measure $\mu$ has a probability density
$\rho_{\infty}$ with respect to Lebesgue measure. In [41, Section 9.1.5],
geometric ergodicity is proved in the case of the semi-space
$G=(0,\infty)\times\mathbb{R}^{d-1}$.

The stationary Fokker-Plank equation for the density $\rho_{\infty}$ has the
form

$\frac{\gamma}{\beta}\Delta_{p}\rho_{\infty}-(p\cdot\nabla_{q}\rho_{\infty})+(\nabla_{q}U(q)\cdot\nabla_{p}\rho_{\infty})+\gamma(\nabla_{p}\cdot(p\rho_{\infty}))=0,\;\;\;\;(q,p)\in
G\times\mathbb{R}^{d},$ (2.7)

with boundary condition

$\rho_{\infty}(q,p)=\rho_{\infty}(q,p-2(n_{G}\cdot
p)n_{G}),\;\;\;\;(q,p)\in\Sigma^{+}.$ (2.8)

As one may check, the Gibbs density

$\rho_{\infty}(q,p)=\frac{1}{\mathbb{Z}}\exp{\Big{(}-\beta\big{(}|p|^{2}/2+U(q)\big{)}\Big{)}}$
(2.9)

indeed satisfies the stationary Fokker-Plank equation (2.7) and the boundary
condition (2.8). Here $\mathbb{Z}$ is the normalization constant so that
$\int_{\mathbb{R}^{d}}\int_{\bar{G}}\rho_{\infty}(q,p)dqdp=1$.

In [41, Section 9.1.5], the stationary Fokker-Plank equation (2.7)-(2.8) in
the case of $G=(0,\infty)\times\mathbb{R}^{d-1}$ was considered. In [6],
existence of the solution (in the weak sense) of a McKean-Vlasov Fokker-Plank
equation with specular boundary condition was shown.

The above discussion of geometric ergodicity together with numerical methods
for (ECLD) introduced in Section 3 equip us with the tools needed to
approximately sample from a targeted density $\rho(q)$, $q\in\bar{G}$, by
choosing $U(q)=-\frac{1}{\beta}\ln(\rho(q))$, $\beta>0$.

To prove error estimates for computing ergodic limits, we need the following
assumption on solutions of the backward Kolmogorov equation with specular
boundary condition (2.2)-(2.4) related to (ECLD), i.e., of (2.2)-(2.4) with
$b(q,p)=-\nabla U(q)-\gamma p$.

######  Assumption 2.6.

The solution of (2.2)-(2.4) with $b(q,p)=-\nabla U(q)-\gamma p$ satisfies
Assumption 2.4 and the condition that there are positive constants $C$,
$\lambda$ and $m$ independent of $T$ such that the following bound holds for
some $m\geq 1$:

$\displaystyle\sum\limits_{l=1}^{2}\sum\limits_{i+|j|=l}$
$\displaystyle|D_{t}^{i}D_{q}^{j}u(t,q,p)|+\sum\limits_{l=1}^{4}\sum\limits_{|j|=l}|D_{p}^{j}u(t,q,p)|\leq
C(1+|p|^{2m})e^{-\lambda(T-t)},$ (2.10)

where $D_{x}^{j}u=\frac{\partial^{|j|}u}{\partial x^{j_{1}},\dots,\partial
x^{j_{d}}}$, $x=q\in\bar{G}\;\text{or}\;p\in\mathbb{R}^{d}$ and $j$ is a
multi-index, and $C>0$ is independent of $T$ and $p$.

We note that an estimate of the form (2.10) in the case of Langevin equations
in the whole space is proved in [45].

Let us discuss the relationship between the confined Langevin dynamics (ECLD)
and the reflected gradient SDE (Brownian dynamics with reflection or, in other
words, the overdamped reflected Langevin dynamics) [31, 20, 16, 27]:

$\displaystyle dX(t)$ $\displaystyle=-\nabla
U(X(t))dt+\sqrt{\frac{2}{\beta}}dW(t)+\nu(X(t))I_{\partial G}(X(t))dL(t),$
(2.11)

where $L(t)$ is the local time of the process $X(t)$ on the boundary $\partial
G$. The solution $X(t)$ of (2.11) belongs to $\bar{G}$ and it is
instantaneously reflected along the internal normal $\nu$ to $\partial G$ when
it hits the boundary $\partial G$. Under the same conditions as we have
imposed on the potential $U$ and the boundary $\partial G$ above, the SDE
(2.11) is ergodic with the invariant measure $\rho_{{\rm
RSDE}}(x)\propto\exp{(-\beta U(x)})$ [1], which coincides with
$\rho_{\infty}(q)=\int_{\mathbb{R}^{d}}\rho_{\infty}(q,p)dp$, i.e., both
(ECLD) and (2.11) can be used for sampling from the Gibbs distribution in the
position space. Moreover, analogously to the relationship between Langevin
equation and their overdamped limit in the whole $\mathbb{R}^{d}$ space [40]
(see also [38, 25]), the solution $Q(t)$ of (ECLD) tends in distribution to
the solution $X(t)$ of (2.11) when $\gamma\rightarrow\infty$ [13, 44, 32]. In
[27] we proposed and studied numerical algorithms for sampling from
$\rho_{{\rm RSDE}}(x)$, $x\in\bar{G}$, based on (2.11). Benefits of using
Langevin equations vs overdamped Langevin dynamics in the case of the whole
space are well established (see e.g. [25, 14]). Such a comparison in the case
of confined spaces will be investigated in a future work.

##  3 Numerical methods

In Section 3.1, we present first-order weak methods for the SDEs (CLD). In
Section 3.2, we propose splitting schemes for the SDEs (ECLD) which
demonstrate (in numerical studies) second-order weak convergence. In Section
3.3, we compare known numerical methods for (CLD) with the schemes of Sections
3.1 and 3.2, and we also discuss multiple collisions. We start this section by
describing the ingredients needed for approximating collision dynamics.

Consider a uniform discretization of the time interval $[0,T]$ with time step
$h=T/N$, i.e. $t_{n+1}-t_{n}=h$, $n=0,\dots,N-1$, and $t_{0}=0$. We take
$Q_{0}\in G$ and $P_{0}\in\mathbb{R}^{d}\backslash\\{0_{d}\\}$.

Approximating $P(t)$, $Q(t)$ when the position $Q$ is inside the domain can be
done by any suitable standard method of weak approximation of SDEs in the
whole space (see e.g. [39, 25]). However, to construct numerical methods for
(CLD) and (ECLD), we also need to approximate elastic reflection. Hence the
critical part of proposing such methods is to consider steps of a numerical
scheme when collisions with $\partial G$ happen.

One step with possible collision. Let the current state of a Markov chain
approximating the solution of (CLD) be $(q,p)\in\bar{G}\times\mathbb{R}^{d}$.
We construct the next state $P_{1}$, $Q_{1}$ of the chain as follows.

If $q+ph\in G$,

$A_{c}(q,p;h)=(q+ph,p),$

else (i.e., if $q+ph\in G^{c}$) $A_{c}(q,p;h)=(Q_{1},P_{1}),$ obtained from

$\displaystyle Q_{c}=q+p\tau\in\partial G,$ $\displaystyle
P_{1}=p-(n(Q_{c})\cdot p)n(Q_{c}),\,\,Q_{1}=Q_{c}+(h-\tau)P_{1}.$

The random collision time $\tau$ is found so that if $q+ph\notin G$ then
$\tau$ ensures that the intermediate position $Q_{c}\in\partial G$ (note that
$\big{(}p\cdot n_{G}(Q_{c})\big{)}>0$). Thereafter, the numerical method
follows the step $P_{1}=p-(n(Q_{c})\cdot p)n(Q_{c})$ which results in specular
reflection (also known as elastic reflection) of the direction of momentum,
and then completes the remaining step in position using the reversed momentum.

If there is more than one solution of $q+p\tau\in\partial G$, we take the
smallest positive value of $\tau$. Such a choice is possible due to Assumption
2.1 (see [4]). For illustration, consider an annular domain in
$\mathbb{R}^{2}$:
$G=\\{(q_{1},q_{2})\;|\;r_{1}^{2}<q_{1}^{2}+q_{2}^{2}<r_{2}^{2}\\}$ and
$\partial
G=\\{(q_{1},q_{2})\;|\;q_{1}^{2}+q_{2}^{2}=r_{1}^{2}\\}\;\cup\;\\{(q_{1},q_{2})\;|\;q_{1}^{2}+q_{2}^{2}=r_{2}^{2}\\}$,
where $r_{1}$ and $r_{2}$ are positive constants such that $r_{1}<r_{2}$. Let
$(q_{1},q_{2})=((r_{1}+r_{2})/2,0)$ and $p=(-1,0)$, then four possible
solutions of $q+\tau p\in\partial G$ are $\tau_{1}=(r_{2}-r_{1})/2$,
$\tau_{2}=(3r_{1}+r_{2})/2$, $\tau_{3}=(r_{1}+3r_{2})/2$ and
$\tau_{4}=(r_{1}-r_{2})/2$. Here, $\tau_{4}$ is negative, so we take
$\tau=\min\\{\tau_{1},\tau_{2},\tau_{3}\\}=\tau_{1}=(r_{2}-r_{1})/2$.

###  3.1 First-order schemes

In this section we propose two methods, [PAc] and [AcP]. In the [PAc] scheme
the Markov chain first takes an auxiliary step in momentum $P$ and then
updates position $Q$ via an Ac step; the order of updates is reversed in
[AcP]. We introduce the following notation:

$\delta(q,p\;;h,\xi)=hb(q,p)+h^{1/2}\sigma\xi,$ (3.1)

where $\xi=(\xi^{1},\dots,\xi^{d})^{\intercal}$ is an $\mathbb{R}^{d}$-valued
random variable such that $\xi^{i}$, $i=1,\dots,d$, are i.i.d. real-valued
random variables with bounded all moments and satisfying

$\mathbb{E}\xi^{i}=\mathbb{E}(\xi^{i})^{3}=0,\;\mathbb{E}(\xi^{i})^{2}=1.$
(3.2)

The conditions (3.2) are, for instance, satisfied if the $\\{\xi^{i}\\}$ are
drawn from the standard normal distribution or as i.i.d. discrete random
variables with the distribution:

$P(\xi^{i}=\pm 1)=\frac{1}{2}.$ (3.3)

In Fig. 3.1 we present the proposed schemes. For the sake of clarity, we
describe [PAc] scheme in detail. As already mentioned, the aim is to construct
a Markov chain $(t_{k},Q_{k},P_{k})_{(k\geq 0)}$ to approximate (CLD) in the
weak sense. At every iteration moving from $(t_{k},Q_{k},P_{k})$ to
$(t_{k+1},Q_{k+1},P_{k+1})$, the Markov chain first takes auxiliary steps in
velocity and in position, denoted by $\hat{P}_{k+1}$ and $\hat{Q}_{k+1}$
respectively:

$\displaystyle\hat{P}_{k+1}=P_{k}+\delta(Q_{k},P_{k}\;;h,\xi_{k+1}),\;\;\;\;\hat{Q}_{k+1}=Q_{k}+h\hat{P}_{k+1},$
(3.4)

where $\xi_{k+1}\sim\xi$ and they are mutually independent. The next state of
the Markov chain, $Q_{k+1}$ and $P_{k+1}$, is evaluated according to
Ac($Q_{k}$, $\hat{P}_{k+1};h)$, which means that if $\hat{Q}_{k+1}\in G$, then
$P_{k+1}=\hat{P}_{k+1}$ and $Q_{k+1}=\hat{Q}_{k+1}$, else
$P_{k+1}=\hat{P}_{k+1}-(n_{G}(Q_{c,k+1})\cdot\hat{P}_{k+1})n_{G}(Q_{c,k+1})$
and $Q_{k+1}=Q_{c,k+1}+(h-\tau)P_{k+1}$, where
$Q_{c,k+1}=Q_{k}+\hat{P}_{k+1}\tau\in\partial G$.

[PAc] $\displaystyle\hat{P}_{k+1}=P_{k}+\delta(Q_{k},P_{k};h,\xi_{k+1})$
$\displaystyle(Q_{k+1},P_{k+1})=A_{c}(Q_{k},\hat{P}_{k+1};h)$

[AcP] $\displaystyle(Q_{k+1},\hat{P}_{k+1})=A_{c}(Q_{k},P_{k};h)$
$\displaystyle
P_{k+1}=\hat{P}_{k+1}+\delta(Q_{k+1},\hat{P}_{k+1};h,\xi_{k+1})$

Figure 3.1: Description of [PAc] and [AcP] schemes.

Schemes [PAc] and [AcP] have finite-time convergence with first weak order
when they are applied to (CLD). By way of illustration, we state below Theorem
3.1 for scheme [PAc] which is proved in Section 4.1. When these schemes are
applied to the ergodic SDE (ECLD), they approximate the ergodic limit with
error $\mathcal{O}(h+\exp(-\lambda T))$. We state below the corresponding
Theorem 3.2 for scheme [PAc] which is proved in Section 4.2.

######  Theorem 3.1.

Under Assumptions 2.1-2.4, the order of weak convergence of scheme [PAc] for
(CLD) is $\mathcal{O}(h)$, i.e. the following inequality holds, for
sufficiently small $h>0$:

$|\mathbb{E}\varphi(Q_{N},P_{N})-\mathbb{E}\varphi\big{(}Q(T),P(T)\big{)}|\leq
Ch,$ (3.5)

where $(Q(t),P(t))$ is from (CLD) and $C>0$ is independent of $h$.

######  Theorem 3.2.

Let Assumptions 2.1, 2.3, 2.5-2.6 hold. For sufficiently small $h>0$ the
following error estimate holds for scheme [PAc] applied to (ECLD):

$|\mathbb{E}\varphi(Q_{N},P_{N})-\bar{\varphi}|\leq C(h+\exp(-\lambda T)),$
(3.6)

where $C>0$ and $\lambda>0$ are constants independent of $T$ and $h$.

###  3.2 [Ac, B, O] splitting schemes

In the whole space (i.e., without a reflective boundary) splitting schemes of
second weak order for the Langevin equations, which rely on a concatenation of
the exact integration for the drift and impulse and of the Ornstein-Uhlenbeck
process at each time step, are known for their excellent performance in
molecular dynamics and statistics applications (see, e.g. [35, 9, 24, 25, 39]
and references therein). Hence, it is natural to extend those methods to the
confined Langevin dynamics (ECLD). The [Ac, B, O] splitting schemes proposed
in this subsection degenerate to the standard deterministic Störmer-Verlet
schemes with collision step Ac when (ECLD) degenerates to the (deterministic)
Hamiltonian system with elastic collisions. Such deterministic schemes have
been shown to be first order accurate due to the collision step Ac having
local error of order ${\cal O}(h)$ [18, 25]. Because of this, stochastic
numerics experience [39] suggests that finite-time weak convergence of the
[Ac, B, O] schemes for (ECLD) should be first order, too. However, the [Ac, B,
O] schemes for (ECLD) unexpectedly turned out to be of weak order $2$, which
we demonstrate in several carefully constructed and performed numerical
experiments (Section 5) and theoretically justify when $G$ is a half-space,
i.e. $G=(0,\infty)\times\mathbb{R}^{d-1}$, in Section 4.3.

Let $\xi=(\xi^{1},\dots,\xi^{d})^{\intercal}$ be an $\mathbb{R}^{d}$ valued
random variable such that $\xi^{i}$, $i=1,\dots,d$, are i.i.d. real-valued
random variables with bounded all moments and satisfying

$\mathbb{E}\xi^{i}=\mathbb{E}(\xi^{i})^{3}=\mathbb{E}(\xi^{i})^{5}=0,\;\;\mathbb{E}(\xi^{i})^{2}=1,\;\;\mathbb{E}(\xi^{i})^{4}=3.$
(3.7)

The conditions (3.2) are, for instance, satisfied for $\xi^{i}$ having the
standard normal distribution or being i.i.d. discrete random variables with
the distribution:

$P(\xi^{i}=\pm\sqrt{3})=\frac{1}{6},\,\,P(\xi^{i}=0)=\frac{2}{3}.$ (3.8)

In Fig. 3.2, $\xi_{k}\sim\xi$ and $\zeta_{k}\sim\xi$ and they all are mutually
independent. Introduce the notation

$O(p;h,\xi,\gamma):=pe^{-\gamma
h}+\sqrt{\frac{\sigma^{2}}{2\gamma}(1-\mathrm{e}^{-2\gamma h})}\xi.$

The three proposed methods are presented in Fig. 3.2. It is straightforward to
write down the other three possible symmetric concatenations, namely [BOAcOB],
[AcBOBAc], [AcOBOAc], which we test in Experiment 5.2.

[OBAcBO] $\displaystyle P_{k+1/2}=O(P_{k};h/2,\xi_{k+1},\gamma)$
$\displaystyle\hat{P}_{k+1/2}=P_{k+1/2}-\frac{h}{2}\nabla U(Q_{k})$
$\displaystyle(Q_{k+1},P_{k+1/2}^{r})=A_{c}(Q_{k},\hat{P}_{k+1/2};h)\;\;\;$
$\displaystyle\hat{P}_{k+1}=\hat{P}_{k+1/2}^{r}-\frac{h}{2}\nabla U(Q_{k+1})$
$\displaystyle P_{k+1}=O(\hat{P}_{k+1};h/2,\zeta_{k+1},\gamma)$

[BAcOAcB] $\displaystyle P_{k+1/2}=P_{k}-\frac{h}{2}\nabla U(Q_{k})$
$\displaystyle(Q_{k+1/2},P_{k+1/2}^{r})=A_{c}(Q_{k},P_{k+1/2};h/2)\;\;\;$
$\displaystyle\hat{P}_{k+1/2}=O(P_{k+1/2}^{r};h,\xi_{k+1},\gamma)$
$\displaystyle(Q_{k+1},P_{k+1}^{r})=A_{c}(Q_{k+1/2},\hat{P}_{k+1/2};h/2)$
$\displaystyle P_{k+1}=P_{k+1}^{r}-\frac{h}{2}\nabla U(Q_{k+1})$

[OAcBAcO] $\displaystyle P_{k+1/2}=O(P_{k};h/2,\xi_{k+1},\gamma)$
$\displaystyle(Q_{k+1/2},P_{k+1/2}^{r})=A_{c}(Q_{k},P_{k+1/2};h/2)$
$\displaystyle\hat{P}_{k+1/2}=P_{k+1/2}^{r}-h\nabla U(Q_{k+1/2})$
$\displaystyle(Q_{k+1},P_{k+1}^{r})=A_{c}(Q_{k+1/2},\hat{P}_{k+1/2};h/2)$
$\displaystyle P_{k+1}=O(P_{k+1}^{r};h/2,\zeta_{k+1},\gamma)$

Figure 3.2: Description of [OBAcBO], [BAcOAcB], [OAcBAcO] schemes for (ECLD).

###  3.3 Comparison with other numerical integrators for (CLD)

Here we present a comparison of previously known schemes for (CLD) with the
numerical integrators proposed in the previous two subsections.

Let $G=(0,\infty)$. In [7] (see also [8, 32]), the authors studied a
symmetrized scheme in which, given $(q,p)$ at time $t_{0}$, the approximate
position is always updated as

$\displaystyle\bar{Q}(t_{0}+h)=|q+ph|.$ (3.9)

At every step, $\theta=t_{0}-\frac{q}{p}$ is calculated. A collision happens
if $\theta\in(t_{0},t_{0}+h)$. In that case

$\displaystyle\bar{Q}(\theta)=0,\,\,\,\bar{P}(\theta-)=p+b(q,p)(\theta-
t_{0})+\sigma(W(\theta)-W(t_{0})),\,\,\,\bar{P}(\theta)=-\bar{P}(\theta-),$
$\displaystyle\bar{P}(t+h)=p+b(\bar{Q}(\theta),\bar{P}(\theta))(t_{0}+h-\theta)+\sigma(W(t_{0}+h)-W(\theta)).$

If no collision occurs then
$\bar{P}(t_{0}+h)=p+b(q,p)h+\sigma(W(t_{0}+h)-W(t_{0}))$ (i.e., the
straightforward Euler update). This scheme was generalized to half-space in
higher dimensions, i.e. $G=(0,\infty)\times\mathbb{R}^{d-1}$, in [8, 7, 32].

The main difference between the above scheme and our schemes lies in the
update of momentum. In the scheme above, it is possible for the momentum to
get reversed already in the interval $(t_{0},\theta)$ due to a realized large
Wiener increment. In such a case the update
$\bar{P}(\theta)=-\bar{P}(\theta-)$ results in a direction towards the
negative axis, which leads to multiple collisions during the time interval
$(t_{0},t_{0}+h)$. However, this pathological situation does not arise, by
construction, in the schemes [AcP], [PAc] and the splitting schemes of Section
3.2 if implemented for $G=(0,\infty)\times\mathbb{R}^{d-1}$. We also observe
that if one extends the symmetrized reflection of (3.9) around the boundary in
curved domains, then the symmetrized reflection around the boundary along the
normal does not imitate specular reflection of (CLD), unlike the schemes of
Sections 3.1 and 3.2.

Further, some penalized schemes for (CLD) are also proposed in [32] without
analysis. It is known (see [4, 39, 27] and references therein) that
symmetrized reflection and penalized schemes are good candidates to
approximate the reflected SDEs (2.11), in which the local time $L(t)$ is
responsible for keeping the trajectories $X(t)$ in $\bar{G}$. In contrast, the
confinement of position $Q(t)$ in $\bar{G}$ is achieved in (CLD) by the cádlág
process $\mathcal{R}(t)$ whose explicit expression is known via values of
$(P(t),Q(t))$. In comparison with the existing schemes, all our schemes are
constructed utilizing the expression for $\mathcal{R}(t)$, which, in
particular, ensures conservation of (conjugate) momenta in the event of
collision, i.e.,

$|P_{k}|^{2}=|\hat{P}_{k}|^{2}I_{\bar{G}}({\hat{Q}}_{k})+|\hat{P}_{k}-2(n^{\pi}\cdot\hat{P}_{k})n^{\pi}|^{2}I_{\bar{G}^{c}}({\hat{Q}}_{k})=|\hat{P}_{k}|^{2},$
(3.10)

where $n^{\pi}=n_{G}({Q}_{c,k})$. Since the collision is interlaced in the
one-step procedure of the discrete dynamics, estimates related to the average
number of collisions become crucial for the global error analysis (see Lemmas
4.3 and 4.6).

Although we present our results in a bounded domain $G$ with boundary
$\partial G$, the results can be extended to the half-space
$G=(0,\infty)\times\mathbb{R}^{d-1}$ or any other unbounded domain with
boundary satisfying Assumption 2.1. The only additional requirement is moment
bounds for $Q_{k}$ (moments of $Q(t)$ are typically bounded for most forces
$b$ and potentials $U$ of interest), which can be obtained under a global
Lipschitz assumption on $\nabla U$ or replaced by using the concept of
rejecting trajectories [36, 39] in the non-globally Lipschitz case.
Furthermore, the proposed schemes can be justified in the case of convex
domains with piecewise $C^{2}$ boundary, too.

In our schemes, multi-collisions in a single time-step can occur for a curved
bounded domain $G$ when $|p|$ is large and/or $(n(Q_{c})\cdot p)$ (see Ac
step) is small, i.e. when with sufficiently large speed, the particle collides
with the curved wall in the grazing manner. In this case $Q_{1}$ calculated in
Ac will not lie in $G$. We note that when $\partial G$ is a convex polytope,
we can always choose $h>0$ sufficiently small such that multi-collisions occur
with exponentially low probability as in this case multi-collisions can only
happen when $|p|$ is very large. It is also clear that there is no possibility
in our schemes for multi-collisions in half-space. But when $\partial G$ is
curved, geometry comes into the picture and it may potentially result in
multi-collisions in the discrete dynamics. Obviously, our schemes can be
adapted to this situation by applying the Ac step recursively. The question
that remains to be answered is how many times we need to utilize this
recursion. Let us again consider the one-step procedure of [PAc] scheme:
depending on the update of the momentum P, Ac is deterministic which means the
dynamics becomes billiard-like in high dimensions, i.e., a straight line
motion between two collisions.

For the sake of the discussion here, let us denote by $(\tilde{q},\tilde{p})$
the position and new velocity, respectively, after the P step in [PAc]. Also,
denote the continuous-time dynamics evolving according to the Ac step by
$(\tilde{Q}(t),\tilde{P}(t))$. As already highlighted, the dynamics of
$\tilde{Q}(t)$ follows straight paths, and indeed is the exact solution of the
following ODE system:

$\displaystyle d\tilde{Q}(t)$
$\displaystyle=\tilde{P}(t)dt,\;\;\;\;\tilde{Q}(t_{0})=\tilde{q},$
$\displaystyle d\tilde{P}(t)$
$\displaystyle=d\mathcal{R}(t),\;\;\;\;\tilde{P}(t_{0})=\tilde{p},$

in the interval $(t_{0},t_{0}+h)$. It is clear that the dynamics in the
duration of $h$ may have multiple collisions. It is proved in [33, p. 5], that
the Lebesgue measure of a set of configurations $(\tilde{q},\tilde{p})$ for
which one encounters infinite collisions in finite time (in the above case in
time from $t_{0}$ to $t_{0}+h$) and for which dynamics hit the boundary with
zero normal component of velocity is zero. Hence, we expect that the number of
trajectories corresponding to our schemes which experience multiple collisions
in a single step are negligible for sufficiently small $h>0$ and, indeed, this
was confirmed in our experiments (see Section 5).

We observe that all the schemes of Sections 3.1 and 3.2 can be implemented to
account for more than one collision in a step. In practice, for simplicity, we
suggest to reject trajectories for which more than one collision occurs in a
step of the Markov chain. In experiments for the ergodic case (i.e., during
long time simulations) on a circular domain of radius two in two dimensions,
we encountered multiple collisions in less than $10$ out of $10^{6}$ simulated
trajectories for $0.2<h<0.5$, and we observed no steps with multiple
collisions for $h<0.2$, i.e. the effect of multiple collisions in our schemes
is negligible, as expected. Similar results were seen in the case where we
chose $\nabla U$ to deliberately push the dynamics towards the boundary to
enhance the number of collisions (see Experiment 5.4).

##  4 Proofs

In this section, for the first-order scheme, we prove the finite-time
convergence Theorem 3.1 (Section 4.1) and Theorem 3.2 (Section 4.2) on the
error estimate for computing ergodic limits. We also provide a justification
(Section 4.3) for second-order weak convergence of the splitting schemes from
Section 3.2.

###  4.1 Proof of Theorem 3.1

The proof is based on two main ingredients: a one-step error estimate (Lemma
4.2) and an estimate related to the average number of collisions (Lemma 4.3).
According to the lemmas, inside the domain, where we use the usual Euler-type
walker, the one-step error is $\mathcal{O}(h^{2})$ and the walker makes
$\mathcal{O}(1/h)$ steps; the error in the collision steps is $\mathcal{O}(h)$
but the average number of collision steps is bounded and independent of the
time step $h$. The finite-time error is naturally $\mathcal{O}(h)$. We note
that the position $Q(t)$ in the continuous dynamics (CLD) and its discrete
approximation $Q_{k}$ is bounded as we consider (CLD) in a bounded domain
$\bar{G}$. At the same time, momentum $P(t)$ lives in $\mathbb{R}^{d}$ and
hence, to prove convergence, we need that moments of its discrete
approximation $P_{k}$ are bounded (see Lemma 4.1). The proof relies on the
fact that the drift in (CLD) is growing at most linearly in $P(t)$ and the
noise is additive.

######  Lemma 4.1.

Let Assumptions 2.1-2.2 hold, then the moments of the momentum $P_{k}$
constructed according to the scheme [PAc] applied to (CLD) are bounded, i.e.
for any $m\geq 1$:

$\mathbb{E}(|P_{k}|^{2m})\leq C(1+|P_{0}|^{2m}),$ (4.1)

where $C>0$ is a constant independent of $h$.

###### Proof.

Using (3.10), we get

$\displaystyle\mathbb{E}|P_{k+1}|^{2m}=\mathbb{E}|\hat{P}_{k+1}|^{2m}=\mathbb{E}|(\hat{P}_{k+1}-P_{k})+P_{k}|^{2m}$
$\displaystyle\leq\mathbb{E}|P_{k}|^{2m}+m\mathbb{E}|P_{k}|^{2m-2}[2(P_{k}\cdot(\hat{P}_{k+1}-P_{k}))+(2m-1)|\hat{P}_{k+1}-P_{k}|^{2}]$
$\displaystyle+K\sum_{l=3}^{2m}\mathbb{E}|P_{k}|^{2m-l}|\hat{P}_{k+1}-P_{k}|^{l},$

where $K>0$ depends on $m$ only. Thanks to (2.1), it is not difficult to
obtain

$|\mathbb{E}[(\hat{P}_{k+1}-P_{k})|P_{k}]|\leq
Kh|P_{k}|,\,\,\,\mathbb{E}[|\hat{P}_{k+1}-P_{k}|^{l}|P_{k}]\leq
Kh^{l/2}(1+|P_{k}|^{l}),\;l\geq 2,$

with $K>0$ dependent on $m$, $b$ and the moments of $\xi_{k}$. We observe that

$\displaystyle\mathbb{E}|P_{k+1}|^{2m}\leq\mathbb{E}|P_{k}|^{2m}+Kh\mathbb{E}|P_{k}|^{2m}+Kh$
(4.2)

and, using the Gronwall lemma, we get (4.1). ∎

Next we estimate the one-step weak error.

######  Lemma 4.2.

Under Assumptions 2.1-2.4, the one-step error of the scheme [PAc] is estimated
as

$\displaystyle\big{|}\mathbb{E}[u(t_{k+1},Q_{k+1},P_{k+1})-u(t_{k},Q_{k},P_{k})\mid
Q_{k},P_{k}]\big{|}$ $\displaystyle\leq
C\big{(}h^{2}(1+|P_{k}|^{2m})+h\mathbb{E}\big{[}(n_{G}(Q_{c,k+1})\cdot\hat{P}_{k+1})(1+|P_{k+1}|^{2m})I_{G^{c}}(\hat{Q}_{k+1})\big{|}Q_{k},P_{k}\big{]}\big{)},$

where $u(t,q,p)$ is the solution of (2.2)-(2.4) and $C>0$ is a constant
independent of $h$ and $P_{k}$.

###### Proof.

In the proof the letter $C$ is used for positive constants and $m$ for
positive integers which can change from line to line. We note that it follows
from (4.2) that

$\displaystyle\mathbb{E}[|P_{k+1}|^{2m}\mid Q_{k},P_{k}]\leq
C(h+|P_{k}|^{2m}),$

where $C>0$ is independent of $h$, $P_{k}$, $Q_{k}$.

We have

$\displaystyle u($ $\displaystyle
t_{k+1},Q_{k+1},P_{k+1})-u(t_{k},Q_{k},P_{k})=(u(t_{k+1},Q_{k+1},P_{k+1})-u(t_{k+1},Q_{k},\hat{P}_{k+1}))I_{G}(\hat{Q}_{k+1})$
(4.3)
$\displaystyle+(u(t_{k+1},Q_{k+1},P_{k+1})-u(t_{k+1},Q_{k},\hat{P}_{k+1}))I_{G^{c}}(\hat{Q}_{k+1})+(u(t_{k+1},Q_{k},\hat{P}_{k+1})-u(t_{k},Q_{k},P_{k}))$
$\displaystyle=:A_{1}+A_{2}+A_{3}.$

Consider the first term. Using Taylor’s expansion and Assumption 2.4, we
obtain

$\displaystyle A_{1}$
$\displaystyle:=(u(t_{k+1},Q_{k+1},P_{k+1})-u(t_{k+1},Q_{k},\hat{P}_{k+1}))I_{G}(\hat{Q}_{k+1})$
$\displaystyle=(u(t_{k+1},Q_{k+1},\hat{P}_{k+1})-u(t_{k+1},Q_{k},\hat{P}_{k+1}))I_{G}(\hat{Q}_{k+1})$
$\displaystyle=\Big{(}h\big{(}\hat{P}_{k+1}\cdot\nabla_{q}u(t_{k+1},Q_{k},\hat{P}_{k+1})\big{)}+r_{1}\Big{)}I_{G}(\hat{Q}_{k+1}),$

where the remainder $r_{1}$ satisfies

$|\mathbb{E}[r_{1}I_{G}(\hat{Q}_{k+1})\mid Q_{k},P_{k}]|\leq
Ch^{2}(1+|P_{k}|^{2m})$

with $C>0$ being a constant independent of $h$, $P_{k}$, $Q_{k}$.

Using Taylor’s expansion again for $u(t_{k+1},Q_{k+1},P_{k+1})$ around
$Q_{c,k+1}$ and $\nabla_{q}u(t_{k+1},Q_{c,k+1},P_{k+1})$ around $Q_{k}$, we
get for the second term in (4.3):

$\displaystyle
A_{2}:=(u(t_{k+1},Q_{k+1},P_{k+1})-u(t_{k+1},Q_{k},\hat{P}_{k+1}))I_{G^{c}}(\hat{Q}_{k+1})$
$\displaystyle=\Big{(}u(t_{k+1},Q_{c,k+1},P_{k+1})+(h-\tau)\big{(}P_{k+1}\cdot\nabla_{q}u(t_{k+1},Q_{c,k+1},P_{k+1})\big{)}$
$\displaystyle-u(t_{k+1},Q_{k},\hat{P}_{k+1})+r_{2}\Big{)}I_{G^{c}}(\hat{Q}_{k+1})$
$\displaystyle=\Big{(}u(t_{k+1},Q_{c,k+1},P_{k+1})+(h-\tau)\big{(}P_{k+1}\cdot\nabla_{q}u(t_{k+1},Q_{k},P_{k+1})\big{)}$
$\displaystyle\;\;\;-u(t_{k+1},Q_{k},\hat{P}_{k+1})+r_{3}\Big{)}I_{G^{c}}(\hat{Q}_{k+1})$
$\displaystyle=\Big{(}u(t_{k+1},Q_{c,k+1},\hat{P}_{k+1})+(h-\tau)\big{(}\hat{P}_{k+1}\cdot\nabla_{q}u(t_{k+1},Q_{k},P_{k+1})\big{)}$
$\displaystyle\;\;\;-2(h-\tau)(n_{G}(Q_{c,k+1})\cdot\hat{P}_{k+1})\big{(}n_{G}(Q_{c,k+1})\cdot\nabla_{q}u(t_{k+1},Q_{k},P_{k+1})\big{)}$
$\displaystyle\;\;\;-u(t_{k+1},Q_{k},\hat{P}_{k+1})+r_{3}\Big{)}I_{G^{c}}(\hat{Q}_{k+1}),$

where we have used
$P_{k+1}=\hat{P}_{k+1}-2(n_{G}(Q_{c,k+1})\cdot\hat{P}_{k+1})n_{G}(Q_{c,k+1})I_{G^{c}}(\hat{Q}_{k+1})$
and $u(t_{k+1},Q_{c,k+1},P_{k+1})=u(t_{k+1},Q_{c,k+1},\hat{P}_{k+1})$ from
(2.4). Here $r_{i}$, $i=2,3$, satisfy

$|\mathbb{E}[r_{i}I_{G^{c}}(\hat{Q}_{k+1})\mid Q_{k},P_{k}]|\leq
Ch^{2}(1+|P_{k}|^{2m})$ (4.4)

with $C>0$ being a constant independent of $h$, $P_{k}$, $Q_{k}$. Now,
applying Taylor’s formula to $u(t_{k+1},Q_{c,k+1},\hat{P}_{k+1})$ around
$Q_{k}$ and $\nabla_{q}u(t_{k+1},Q_{k},P_{k+1})$ around $\hat{P}_{k+1}$, we
obtain

$\displaystyle A_{2}$
$\displaystyle=\Big{(}\tau(\hat{P}_{k+1}\cdot\nabla_{q}u(t_{k+1},Q_{k},\hat{P}_{k+1}))+(h-\tau)\big{(}\hat{P}_{k+1}\cdot\nabla_{q}u(t_{k+1},Q_{k},P_{k+1})\big{)}$
$\displaystyle\;\;\;-2(h-\tau)(n_{G}(Q_{c,k+1})\cdot\hat{P}_{k+1})\big{(}n_{G}(Q_{c,k+1})\cdot\nabla_{q}u(t_{k+1},Q_{k},P_{k+1})\big{)}+r_{4}\Big{)}I_{G^{c}}(\hat{Q}_{k+1})$
$\displaystyle=\Big{(}\tau(\hat{P}_{k+1}\cdot\nabla_{q}u(t_{k+1},Q_{k},\hat{P}_{k+1}))+(h-\tau)\big{(}\hat{P}_{k+1}\cdot\nabla_{q}u(t_{k+1},Q_{k},\hat{P}_{k+1})\big{)}$
$\displaystyle\;\;\;-2(h-\tau)(\hat{P}_{k+1}\cdot\textbf{J}_{\nabla_{q}{u}}(\hat{P}_{k,\theta})n_{G}(Q_{c,k+1}))(\hat{P}_{k+1}\cdot
n_{G}(Q_{c,k+1}))$
$\displaystyle\;\;\;-2(h-\tau)(n_{G}(Q_{c,k+1})\cdot\hat{P}_{k+1})\big{(}n_{G}(Q_{c,k+1})\cdot\nabla_{q}u(t_{k+1},Q_{k},P_{k+1})\big{)}+r_{4}\Big{)}I_{G^{c}}(\hat{Q}_{k+1})$
$\displaystyle=\Big{(}h(\hat{P}_{k+1}\cdot\nabla_{q}u(t_{k+1},Q_{k},\hat{P}_{k+1}))-2(h-\tau)(\hat{P}_{k+1}\cdot\textbf{J}_{\nabla_{q}{u}}(\hat{P}_{k,\theta})n_{G}(Q_{c,k+1}))(\hat{P}_{k+1}\cdot
n_{G}(Q_{c,k+1}))$
$\displaystyle\;\;\;-2(h-\tau)(n_{G}(Q_{c,k+1})\cdot\hat{P}_{k+1})\big{(}n_{G}(Q_{c,k+1})\cdot\nabla_{q}u(t_{k+1},Q_{k},P_{k+1})\big{)}+r_{4}\Big{)}I_{G^{c}}(\hat{Q}_{k+1}),$

where $\textbf{J}_{\nabla_{q}{u}}(\hat{P}_{k,\theta})$ is the Jacobian of
$\nabla_{q}{u}$ with respect to $p$ evaluated at
$(t_{k+1},Q_{k},\hat{P}_{k,\theta})$ with

$\displaystyle\hat{P}_{k,\theta}=\hat{P}_{k+1}-2\theta(n_{G}(Q_{c,k+1})\cdot\hat{P}_{k+1})n_{G}(Q_{c,k+1}),\,\,\theta\in(0,1),$
(4.5)

and $r_{4}$ satisfies (4.4).

Combining the obtained expressions for $A_{1}$ and $A_{2}$, we ascertain that

$\displaystyle A_{1}+A_{2}$
$\displaystyle=h(\hat{P}_{k+1}\cdot\nabla_{q}u(t_{k+1},Q_{k},\hat{P}_{k+1}))$
$\displaystyle\;\;\;-2\Big{(}(h-\tau)(\hat{P}_{k+1}\cdot\textbf{J}_{\nabla_{q}{u}}(\hat{P}_{k,\theta})n_{G}(Q_{c,k+1}))(\hat{P}_{k+1}\cdot
n_{G}(Q_{c,k+1}))$
$\displaystyle\;\;\;+(h-\tau)(n_{G}(Q_{c,k+1})\cdot\hat{P}_{k+1})\big{(}n_{G}(Q_{c,k+1})\cdot\nabla_{q}u(t_{k+1},Q_{k},P_{k+1})\big{)}\Big{)}I_{G^{c}}(\hat{Q}_{k+1})+r_{5}$

with $r_{5}$ satisfying

$|\mathbb{E}[r_{5}\mid Q_{k},P_{k}]|\leq Ch^{2}(1+|P_{k}|^{2m})$ (4.6)

with $C>0$ being a constant independent of $h$, $P_{k}$, $Q_{k}$. Taking the
expectation conditioned on $(Q_{k},P_{k})$, we get

$\displaystyle\mathbb{E}(A_{1}+A_{2}|\;Q_{k},P_{k})=h\mathbb{E}\big{[}(\hat{P}_{k+1}\cdot\nabla_{q}u(t_{k+1},Q_{k},\hat{P}_{k+1}))|Q_{k},P_{k}\big{]}$
$\displaystyle\;\;\;-2\mathbb{E}\Big{[}\Big{(}(h-\tau)(\hat{P}_{k+1}\cdot\textbf{J}_{\nabla_{q}{u}}(\hat{P}_{k,\theta})n_{G}(Q_{c,k+1}))(\hat{P}_{k+1}\cdot
n_{G}(Q_{c,k+1}))$
$\displaystyle\;\;\;+(h-\tau)(n_{G}(Q_{c,k+1})\cdot\hat{P}_{k+1})\big{(}n_{G}(Q_{c,k+1})\cdot\nabla_{q}u(t_{k+1},Q_{k},P_{k+1})\big{)}\Big{)}I_{G^{c}}(\hat{Q}_{k+1})\big{|}Q_{k},P_{k}\Big{]}$
$\displaystyle\;\;\;+\mathbb{E}[r_{5}\mid Q_{k},P_{k}]$
$\displaystyle=h\mathbb{E}\big{[}(P_{k}\cdot\nabla_{q}u(t_{k},Q_{k},P_{k}))|Q_{k},P_{k}\big{]}+r_{6}$
$\displaystyle-2\mathbb{E}\Big{[}\Big{(}(h-\tau)(\hat{P}_{k+1}\cdot\textbf{J}_{\nabla_{q}{u}}(\hat{P}_{k,\theta})n_{G}(Q_{c,k+1}))(\hat{P}_{k+1}\cdot
n_{G}(Q_{c,k+1}))$
$\displaystyle+(h-\tau)(n_{G}(Q_{c,k+1})\cdot\hat{P}_{k+1})\big{(}n_{G}(Q_{c,k+1})\cdot\nabla_{q}u(t_{k+1},Q_{k},P_{k+1})\big{)}\Big{)}I_{G^{c}}(\hat{Q}_{k+1})\big{|}Q_{k},P_{k}\Big{]},$
(4.7)

where we have expanded $\nabla_{q}u(t_{k+1},Q_{k},\hat{P}_{k+1})$ around
$P_{k}$ and then $\nabla_{q}u(t_{k+1},Q_{k},P_{k})$ around $t_{k}$, and
$r_{6}$ satisfies

$|r_{6}|\leq Ch^{2}(1+|P_{k}|^{2m})$ (4.8)

with $C>0$ being a constant independent of $h$, $P_{k}$, $Q_{k}$.

With the usual Taylor expansion calculation for the Euler scheme for SDEs (see
e.g. [39, Chap. 7]), we can write

$\displaystyle\mathbb{E}(A_{3}|\;Q_{k},P_{k})$
$\displaystyle=h\frac{\partial}{\partial
t}u(t_{k},Q_{k},P_{k})+h(b_{k}\cdot\nabla_{p}u(t_{k},Q_{k},P_{k}))+h\frac{\sigma^{2}}{2}\Delta_{p}u(t_{k},Q_{k},P_{k})+r_{7},$
(4.9)

where $r_{7}$ satisfies (4.8).

Note that (see Assumption 2.4, (4.5) and (3.10)):

$\displaystyle|n_{G}(Q_{c,k+1})\cdot\nabla_{q}u(t_{k+1},Q_{k},P_{k+1})|$
$\displaystyle\leq C(1+|P_{k+1}|^{2m}),$
$\displaystyle|\hat{P}_{k,\theta}|\leq C|\hat{P}_{k+1}|=C|P_{k+1}|,$
$\displaystyle|(\hat{P}_{k+1}\cdot\textbf{J}_{\nabla_{q}u}(\hat{P}_{k,\theta})n_{G}(Q_{c,k+1}))|\leq
C(1+|P_{k+1}|^{2m}).$

Thus, using (4.7), (4.9) and (2.2), we infer

$\displaystyle|\mathbb{E}(u(t_{k+1},Q_{k+1},P_{k+1})-u(t_{k},Q_{k},P_{k})|\;Q_{k},P_{k})|$
$\displaystyle\leq
C\big{(}h^{2}(1+|P_{k}|^{2m})+h\mathbb{E}\big{[}(n_{G}({Q}_{c,k+1})\cdot\hat{P}_{k+1})(1+|P_{k+1}|^{2m})I_{G^{c}}(\hat{Q}_{k+1})\big{|}Q_{k},P_{k}\big{]}\big{)},$

where $C>0$ is independent of $h$ and $P_{k}$.

∎

The estimate in the above local error lemma contains the term of order
$\mathcal{O}(h)$ corresponding to the error of approximating collisions.
Consequently, the next lemma, related to the average number of collisions when
the Markov chain $(Q_{k},P_{k})$ encounters the boundary $\Sigma_{+}$, plays a
major role in global convergence analysis of the considered first-order
methods. Its proof rests on the boundary value problem associated with the
Markov chain $(t_{k},Q_{k},P_{k})$:

$\displaystyle\mathbb{T}V(t,q,p)-V(t,q,p)$ $\displaystyle=-g(t,q,p),$
$\displaystyle(t,q,p)\in[T_{0},T-h]\times\bar{G}\times\mathbb{R}^{d},$ (4.10)
$\displaystyle V(t,q,p)$ $\displaystyle=0,$
$\displaystyle(t,q,p)\in\\{T\\}\times\bar{G}\times\mathbb{R}^{d},$ (4.11)

where $\mathbb{T}$ is the one-step transition operator:

$\mathbb{T}V(t,q,p)=\mathbb{E}V(t_{1},Q_{1},P_{1}),\;\;(t_{0},Q_{0},P_{0})=(t,q,p),$

and $g(t,q,p)\geq 0$. The solution to this problem is given by [46, 37, 39]:

$V(t_{0},q,p)=\mathbb{E}\bigg{(}\sum\limits_{k=0}^{N-1}g(t_{k},Q_{k},P_{k})\Big{|}Q_{0}=q,P_{0}=p\bigg{)}.$
(4.12)

######  Lemma 4.3.

Let Assumptions 2.1-2.2 hold, then for any integer $m\geq 1$:

$\mathbb{E}\bigg{(}\sum\limits_{k=0}^{N-1}(n_{G}(Q_{c,k+1})\cdot\hat{P}_{k+1})(1+|P_{k+1}|^{2m})I_{G^{c}}\big{(}\hat{Q}_{k+1}\big{)}\bigg{)}\leq
C,$ (4.13)

where $\hat{Q}_{k+1}=Q_{k}+\hat{P}_{k+1}h$, $k=0,\dots,N-1$, are from [PAc]
scheme and $C>0$ is a constant independent of $h$.

###### Proof.

Let us take $g(t_{k},q,p)$ in (4.10) as

$g(t_{k},q,p)=\mathbb{E}\big{[}(n_{G}(Q_{c,k+1})\cdot\hat{P}_{k+1})(1+|P_{k+1}|^{2m})I_{G^{c}}(\hat{Q}_{k+1})|\;Q_{k}=q,P_{k}=p\big{]}$

then the solution of (4.10)-(4.11) is given by

$\displaystyle v$
$\displaystyle(t_{0},q,p)=\mathbb{E}\bigg{(}\sum_{k=0}^{N-1}g(t_{k},Q_{k},P_{k})\big{|}Q_{0}=q,P_{0}=p\bigg{)}$
$\displaystyle=\mathbb{E}\bigg{(}\sum\limits_{k=0}^{N-1}\mathbb{E}\Big{(}(n_{G}(Q_{c,k+1})\cdot\hat{P}_{k+1})(1+|P_{k+1}|^{2m})I_{G^{c}}\big{(}\hat{Q}_{k+1}\big{)}\big{|}Q_{k},P_{k}\Big{)}\Big{|}Q_{0}=q,P_{0}=p\bigg{)}$
$\displaystyle=\mathbb{E}\bigg{(}\sum\limits_{k=0}^{N-1}\big{(}n_{G}(Q_{c,k+1})\cdot\hat{P}_{k+1}\big{)}\big{(}1+|P_{k+1}|^{2m}\big{)}I_{G^{c}}\big{(}\hat{Q}_{k+1}\big{)}\Big{|}Q_{0}=q,P_{0}=p\bigg{)}.$

If we can find a function $V(t,q,p)$ with $g(t,q,p)$, satisfying
(4.10)-(4.11), such that

$g(t,q,p)\geq\mathbb{E}\big{[}(n_{G}(Q_{c})\cdot\hat{P}_{1})(1+|P_{1}|^{2m})I_{G^{c}}(\hat{\mathcal{Q}}_{1})\big{]},$
(4.14)

then the inequality, $v(t_{0},Q_{0},P_{0})\leq V(t_{0},Q_{0},P_{0})$ will
hold. Therefore, the aim is to construct a function $g=-(\mathbb{T}V-V)$
satisfying (4.14).

Consider a boundary zone $G_{R_{0}}\subset G$ defined as
$G_{R_{0}}:=\\{q\in\bar{G}\;;\;\operatorname{dist}(q,\partial G)<R_{0}\\}$.
Due to [4, Proposition 1], there is an $R_{0}$ such that the projection of $q$
on $\partial G$ is unique whenever $q\in G_{R_{0}}$. Introduce the function

$V(t,q,p)=\begin{cases}0,&(t,q,p)\in\\{T\\}\times\bar{G}\times\mathbb{R}^{d},\\\
(A+e^{L(T-t)})(1+|p|^{2m+2})+(B(q)\cdot
p)(1+|p|^{2m}),&(t,q,p)\in[0,T-h]\times\bar{G}\times\mathbb{R}^{d},\end{cases}$

where $B(q)\in C^{3}(\bar{G})$ is an $\mathbb{R}^{d}$ valued function such
that $B(q):=n_{G}(q^{\pi})$ when $q\in\bar{G}_{R_{0}}$; here $q^{\pi}$ is the
projection of $q$ on $\partial G$. Note that $n_{G}(q)\in C^{4}(\bar{G})$ due
to Assumption 2.1 (see [4, Proposition 1]). We extend the function $B(q)$ from
$G_{R_{0}}$ to $\bar{G}$ so that $B(q)\in C^{3}(\bar{G})$ and
$|B(q)|_{G}^{2,\alpha}\leq C|n(q^{\pi})|_{G_{R_{0}}}^{2,\alpha}$, where
$|\cdot|^{2,\alpha}$ is the Hölder norm and $C>0$ depends on $G$. Such an
extension is possible due to [17, Lemma 6.37]. We choose the constant $A>0$ so
that $V(t,q,p)$ is always positive. Note that $B(q)$ is a bounded function
(with bounded derivatives up to third order), hence it is possible to make
such a choice of $A$. The choice of the constant $L>0$ will be discussed later
in the proof.

The proof of Lemma 4.1 implies for any $\ell\geq 1$

$\displaystyle|P_{1}|^{2\ell}=|\hat{P}_{1}|^{2\ell}=|p|^{2\ell}+2\ell\sigma
h^{1/2}(\xi\cdot p)|p|^{2\ell-2}+r_{1},$ (4.15)

and

$\displaystyle\mathbb{E}[|P_{1}|^{2\ell}\;|q,p]=|p|^{2\ell}+r_{2},$ (4.16)

where $r_{1}$ and $r_{2}$ satisfy for any $\ell\geq 1$:

$(\mathbb{E}(|r_{1}|^{2\ell}))^{1/(2\ell)}\leq Ch(1+|p|^{2\ell})$

and

$|r_{2}|\leq Ch(1+|p|^{2\ell})$

with $C>0$ being independent of $h$ and $p$.

Given $(t_{k},Q_{k},P_{k})=(t,q,p)$, we now aim at calculating
$-(\mathbb{T}V(t,q,p)-V(t,q,p))$. Using Taylor’s theorem, we obtain

$\displaystyle B(Q_{1})$
$\displaystyle=\big{(}B(q)+h\textbf{J}_{B}(q+\epsilon_{1}h\hat{P}_{1})\hat{P}_{1}\big{)}I_{G}(\hat{Q}_{1})$
$\displaystyle\;\;\;+\big{(}B(Q_{c})+(h-\tau)\textbf{J}_{B}(Q_{c}+\epsilon_{2}(h-\tau)P_{1})P_{1}\big{)}I_{G^{c}}(\hat{Q}_{1})$
(4.17)

and

$\displaystyle B(Q_{c})$
$\displaystyle=B(q)+\tau\textbf{J}_{B}(q+\epsilon_{3}\tau\hat{P}_{1})\hat{P}_{1},$
(4.18)

where $\epsilon_{i}\in(0,1)$, $i=1,2,3$, and $\textbf{J}_{B}(x)$ denotes the
Jacobian matrix of $B$ evaluated at $x$. Therefore, we get

$\displaystyle(B(Q_{1})$ $\displaystyle\cdot
P_{1})=(B(Q_{1})\cdot\hat{P}_{1})I_{G}(\hat{Q}_{1})+(B(Q_{1})\cdot
P_{1})I_{G^{c}}(\hat{Q}_{1})$
$\displaystyle=\big{(}(B(q)\cdot\hat{P}_{1})+h(\textbf{J}_{B}(q+\epsilon_{1}h\hat{P}_{1})\hat{P}_{1}\cdot\hat{P}_{1})\big{)}I_{G}(\hat{Q}_{1})$
$\displaystyle\;\;\;+\big{(}(B(Q_{c})\cdot
P_{1})+(h-\tau)(\textbf{J}_{B}(Q_{c}+\epsilon_{2}(h-\tau)P_{1})P_{1}\cdot
P_{1})\big{)}I_{G^{c}}(\hat{Q}_{1}))$
$\displaystyle=\big{(}(B(q)\cdot\hat{P}_{1})+h(\textbf{J}_{B}(q+\epsilon_{1}h\hat{P}_{1})\hat{P}_{1}\cdot\hat{P}_{1})\big{)}I_{G}(\hat{Q}_{1})$
$\displaystyle\;\;\;+\big{(}(B(Q_{c})\cdot\hat{P}_{1})+(h-\tau)(\textbf{J}_{B}(Q_{c}+\epsilon_{2}(h-\tau)P_{1})P_{1}\cdot
P_{1})\big{)}I_{G^{c}}(\hat{Q}_{1}))$ $\displaystyle\;\;\;-2(\hat{P}_{1}\cdot
n_{G}(Q_{c}))(B(Q_{c})\cdot n_{G}(Q_{c}))I_{G^{c}}(\hat{Q}_{1}).$ (4.19)

Since $(B(Q_{c})\cdot n_{G}(Q_{c}))=1$ and using (4.18), we obtain

$\displaystyle(B(Q_{1})$ $\displaystyle\cdot
P_{1})=\big{(}(B(q)\cdot\hat{P}_{1})+h(\textbf{J}_{B}(q+\epsilon_{1}h\hat{P}_{1})\hat{P}_{1}\cdot\hat{P}_{1})\big{)}I_{G}(\hat{Q}_{1})$
$\displaystyle\;\;\;+\big{(}(B(Q_{c})\cdot\hat{P}_{1})+(h-\tau)(\textbf{J}_{B}(\mathcal{Q}_{c}+\epsilon_{2}(h-\tau)P_{1})P_{1}\cdot
P_{1})\big{)}I_{G^{c}}(\hat{Q}_{1}))$ $\displaystyle\;\;\;-2(\hat{P}_{1}\cdot
n_{G}(Q_{c}))I_{G^{c}}(\hat{Q}_{1})$
$\displaystyle=(B(q)\cdot\hat{P}_{1})+h(\textbf{J}_{B}(q+\epsilon_{1}h\hat{P}_{1})\hat{P}_{1}\cdot\hat{P}_{1})\big{)}I_{G}(\hat{Q}_{1})$
$\displaystyle\;\;\;+\big{(}\tau(\textbf{J}_{B}(q+\epsilon_{3}\tau\hat{P}_{1})\hat{P}_{1}\cdot\hat{P}_{1})+(h-\tau)(\textbf{J}_{B}(\mathcal{Q}_{c}+\epsilon_{2}(h-\tau)P_{1})P_{1}\cdot
P_{1})\big{)}I_{G^{c}}(\hat{Q}_{1})$ $\displaystyle\;\;\;-2(\hat{P}_{1}\cdot
n_{G}(Q_{c}))I_{G^{c}}(\hat{Q}_{1}).$ (4.20)

Due to the boundedness of first derivatives of $B$, we ascertain

$\displaystyle(\textbf{J}_{B}(q+\epsilon_{1}h\hat{P}_{1})\hat{P}_{1}\cdot\hat{P}_{1})$
$\displaystyle\leq C|\hat{P}_{1}|^{2},$
$\displaystyle(\textbf{J}_{B}(\mathcal{Q}_{c}+\epsilon_{2}(h-\tau)P_{1})P_{1}\cdot
P_{1})$ $\displaystyle\leq C|P_{1}|^{2}=C|\hat{P}_{1}|^{2},$
$\displaystyle(\textbf{J}_{B}(q+\epsilon_{3}\tau\hat{P}_{1})\hat{P}_{1}\cdot\hat{P}_{1})$
$\displaystyle\leq C|\hat{P}_{1}|^{2},$

where $C>0$ is a generic constant independent of $\hat{P}_{1}$. Using these
estimates, (4.20) becomes

$\displaystyle(B(Q_{1})$ $\displaystyle\cdot
P_{1})=(B(q)\cdot\hat{P}_{1})-2(\hat{P}_{1}\cdot
n_{G}(Q_{c}))I_{G^{c}}(\hat{Q}_{1})+r_{3},$ (4.21)

where $|r_{3}|\leq Ch|\hat{P}_{1}|^{2}$ with $C>0$ being independent of $h$
and $\hat{P}_{1}$. Using (4.15) and linear growth of $b$ in $p$, we obtain

$\displaystyle\mathbb{E}[(B(Q_{1})\cdot
P_{1})(1+|P_{1}|^{2m})]=\mathbb{E}\bigg{[}\Big{(}(B(q)\cdot\hat{P}_{1})-2(\hat{P}_{1}\cdot
n_{G}(Q_{c}))I_{G^{c}}(\hat{Q}_{1})+r_{3}\Big{)}\big{(}1+|P_{1}|^{2m}\big{)}\bigg{]}$
$\displaystyle=\mathbb{E}\big{[}(B(q)\cdot\hat{P}_{1})(1+|P_{1}|^{2m})\big{]}-2\mathbb{E}\big{[}(\hat{P}_{1}\cdot
n_{G}(Q_{c}))(1+|P_{1}|^{2m})I_{G^{c}}(\hat{Q}_{1})\big{]}+r_{4},$ (4.22)

where $|r_{4}|\leq Ch[1+|p|^{2m+2}]$ with $C>0$ being independent of $h$ and
$p$. Using (4.15) and (4.16), we get

$\displaystyle\mathbb{E}\big{[}(B(q)$
$\displaystyle\cdot\hat{P}_{1})(1+|P_{1}|^{2m})\big{]}=(B(q)\cdot
p)\mathbb{E}[(1+|P_{1}|^{2m})]+h(B(q)\cdot b)\mathbb{E}[(1+|P_{1}|^{2m})]$
$\displaystyle\;\;\;+h^{1/2}\mathbb{E}[(B(q)\cdot\xi)(1+|P_{1}|^{2m})]=(B(q)\cdot
p)(1+|p|^{2m})+r_{5},$ (4.23)

where $|r_{5}|\leq Ch(1+|p|^{2m+1})$ with $C>0$ being independent of $h$ and
$p$. From (4.22) and (4.23), we have

$\displaystyle\mathbb{E}[(B(Q_{1})\cdot P_{1})(1$
$\displaystyle+|P_{1}|^{2m})]=(B(q)\cdot
p)(1+|p|^{2m})-2\mathbb{E}\big{[}(\hat{P}_{1}\cdot
n_{G}(Q_{c}))(1+|P_{1}|^{2m})I_{G^{c}}(\hat{Q}_{1})\big{]}+r_{6},$ (4.24)

where $|r_{6}|\leq Ch(1+|p|^{2m+2})$ with $C>0$ being independent of $h$ and
$p$.

Using (4.16) with $\ell=m+1$, we have

$\displaystyle(A+e^{L(T-t-h)})(1+\mathbb{E}[|P_{1}|^{2m+2}])=\big{(}A+e^{L(T-t)}(1-Lh+L^{2}r_{7}\big{)}\big{(}1+|p|^{2m+2}+r_{2}\big{)}$
$\displaystyle=(A+e^{L(T-t)})(1+|p|^{2m+2}+r_{2})-Lhe^{L(T-t)}(1+|p|^{2m+2})+r_{8}e^{L(T-t)},$
(4.25)

where $r_{7}\leq Ch^{2}$ and $r_{8}\leq Ch^{2}(1+L^{2})(1+|p|^{2m+2})$ with
$C>0$ independent of $h$, $L$ and $p$. Recall that $r_{2}\leq
Ch(1+|p|^{2m+2})$, where $C>0$ is independent of $h$, $L$ and $p$. Therefore,
using (4.24) and (4.25), we get

$\displaystyle\mathbb{E}[$ $\displaystyle
V(t_{1},Q_{1},P_{1})]-V(t,q,p)=(A+e^{L(T-t-h)})(1+\mathbb{E}[|P_{1}|^{2m+2}])$
$\displaystyle\;\;\;\;+\mathbb{E}[(B(Q_{1})\cdot
P_{1})(1+|P_{1}|^{2m})]-(A+e^{L(T-t)})(1+|p|^{2m+2})-(B(q)\cdot
p)(1+|p|^{2m})$
$\displaystyle=(A+e^{L(T-t)})r_{2}-Lhe^{L(T-t)}(1+|p|^{2m+2})+r_{8}e^{L(T-t)}+r_{6}$
$\displaystyle\;\;\;-2\mathbb{E}\big{[}(\hat{P}_{1}\cdot
n_{G}(Q_{c}))(1+|P_{1}|^{2m})I_{G^{c}}(\hat{Q}_{1})\big{|}q,p\big{]}$ (4.26)

Hence, there are sufficiently large $L>0$ and sufficiently small $h>0$ so that
$g(t,q,p):=-(\mathbb{T}V-V)\geq\mathbb{E}\big{[}(\hat{P}_{1}\cdot
n_{G}(Q_{c}))(1+|P_{1}|^{2m})I_{G^{c}}(\hat{Q}_{1})|q,p\big{]}$. This proves
that (4.13) holds with a constant $C>0$ independent of $h$ (it depends on
$P_{0}$ and $T$).

∎

###### Proof of Theorem 3.1.

Using Lemma 4.2, we have

$\displaystyle|\mathbb{E}(\varphi(Q_{N},P_{N}))-\mathbb{E}(\varphi(Q(T),P(T)))|=\Big{|}\mathbb{E}\Big{(}\sum\limits_{k=0}^{N-1}\mathbb{E}\big{(}u(t_{k+1},Q_{k+1},P_{k+1})-u(t_{k},Q_{k},P_{k})|Q_{k},P_{k}\big{)}\Big{)}\Big{|}$
$\displaystyle\leq
Ch^{2}\mathbb{E}\Big{(}\sum\limits_{k=0}^{N-1}|P_{k}|^{2m}\Big{)}+Ch\mathbb{E}\Big{(}\sum\limits_{k=0}^{N-1}(n_{G}({Q}_{c,k+1})\cdot\hat{P}_{k+1})I_{G^{c}}(\hat{Q}_{k+1})(1+|P_{k+1}|^{2m})\Big{)},$

which on applying Lemma 4.1 and Lemma 4.3 gives the required result. ∎

###  4.2 Proof of Theorem 3.2

The proof of Theorem 3.2 has the same ingredients as the proof of Theorem 3.1
in the previous section with the major difference that the global error bound
in Theorem 3.2 should not depend on the integration time $T$ since here we are
interested in ergodic limits associated with (ECLD), while the error bound in
Theorem 3.1 for the approximation of (CLD) on the finite time interval $[0,T]$
can depend on $T$. We start (Lemma 4.4) by estimating moments of the momentum
$P(t)$ uniformly in $T$. Next, in Lemma 4.5, we consider the one-step (local)
error. An estimate related to the average number of collisions of the Markov
chain with the boundary is obtained in Lemma 4.6.

######  Lemma 4.4.

Let Assumptions 2.1 and 2.5 hold and $\gamma h<1$, then the moments of the
momentum $P_{k}$ constructed according to the scheme [PAc] applied to (ECLD)
are bounded, i.e. for any $m\geq 1$:

$\mathbb{E}(|P_{k}|^{2m})\leq C,$ (4.27)

where $C>0$ is independent of $h$ and $T$.

###### Proof.

The constant $K>0$ is changing from line to line. Using (3.10), we get (see
the proof of Lemma 4.1):

$\displaystyle\mathbb{E}|P_{k+1}|^{2m}=\mathbb{E}|\hat{P}_{k+1}|^{2m}\leq\mathbb{E}|P_{k}|^{2m}+m\mathbb{E}|P_{k}|^{2m-2}[2(P_{k}\cdot(\hat{P}_{k+1}-P_{k}))+(2m-1)|\hat{P}_{k+1}-P_{k}|^{2}]$
$\displaystyle+K\sum_{l=3}^{2m}\mathbb{E}|P_{k}|^{2m-l}|\hat{P}_{k+1}-P_{k}|^{l},$

where $K>0$ depends on $m$ only. It is not difficult to obtain

$\displaystyle\mathbb{E}[P_{k}\cdot(\hat{P}_{k+1}-P_{k})|P_{k}]$
$\displaystyle\leq$ $\displaystyle-\gamma
h|P_{k}|^{2}+Kh|P_{k}|,\,\,\,\mathbb{E}[|\hat{P}_{k+1}-P_{k}|^{2}|P_{k}]\leq
Kh^{2}|P_{k}|^{2}+\frac{2\gamma}{\beta}h,$
$\displaystyle\mathbb{E}[|\hat{P}_{k+1}-P_{k}|^{l}|P_{k}]$ $\displaystyle\leq$
$\displaystyle Kh^{l}|P_{k}|^{l}+Kh^{l/2},\;l\geq 3,$

with $K>0$ independent of $h$ (it depends on $\max_{q\in\bar{G}}|\nabla
U(q)|,\ \gamma,$ $\beta,$ $l$ and moments of $\xi)$. Using the generalized
Young’s inequality, i.e.
$ab\leq\epsilon^{i}\frac{a^{i}}{i}+\epsilon^{-j}\frac{b^{j}}{j}$; $a,b\geq 0$,
$\epsilon>0$; $i,j>1$; $\frac{1}{i}+\frac{1}{j}=1$, with $a=|P_{k}|^{2m-1}$,
$b=K$, $i=\frac{2m}{(2m-1)}$, $j=2m$ and $\epsilon=\gamma^{\frac{2m-1}{2m}}$,
we get

$\displaystyle
2m\mathbb{E}\left[|P_{k}|^{2m-2}(P_{k}\cdot(\hat{P}_{k+1}-P_{k}))\right]\leq-2m\gamma
h\mathbb{E}|P_{k}|^{2m}+2mhK\mathbb{E}(|P_{k}|^{2m-1})$
$\displaystyle\leq-2m\gamma h\mathbb{E}|P_{k}|^{2m}+(2m-1)\gamma
h\mathbb{E}|P_{k}|^{2m}+hK^{2m}\left(\frac{1}{\gamma}\right)^{2m-1}$
$\displaystyle\leq-\gamma
h\mathbb{E}|P_{k}|^{2m}+hK^{2m}\left(\frac{1}{\gamma}\right)^{2m-1}.$

In the similar manner, we have

$\displaystyle(2m-1)\mathbb{E}|P_{k}|^{2m-2}[|\hat{P}_{k+1}-P_{k}|^{2}]\leq
Kh^{2}\mathbb{E}|P_{k}|^{2m}+\frac{\gamma h}{2}\mathbb{E}|P_{k}|^{2m}+Kh,$

and similarly with some new $K>0$

$K\mathbb{E}|P_{k}|^{2m-l}|\hat{P}_{k+1}-P_{k}|^{l}\leq\frac{\gamma
h^{l/2}}{2m}\mathbb{E}|P_{k}|^{2m}+Kh^{l/2},\;l\geq 3.$

Next note that

$\mathbb{E}|P_{k+1}|^{2m}\leq\mathbb{E}|P_{k}|^{2m}-\frac{\gamma
h}{2}(1-\sqrt{h})\mathbb{E}|P_{k}|^{2m}+Kh$

and, using the Gronwall lemma, we get

$\mathbb{E}|P_{k}|^{2m}\leq\left(1-\frac{\gamma
h}{2}(1-\sqrt{h})\right)^{k}|P_{0}|^{2m}+\frac{2K}{\gamma(1-\sqrt{h})}\left[1-\left(1-\frac{\gamma
h}{2}(1-\sqrt{h})\right)^{k}\right].$

∎

######  Lemma 4.5.

Under Assumptions 2.1, 2.3, 2.5-2.6, the one-step error of numerical scheme
[PAc] applied to (ECLD) is estimated as

$\displaystyle\big{|}\mathbb{E}[u(t_{k+1},Q_{k+1},P_{k+1})-u(t_{k},Q_{k},P_{k})\mid
Q_{k},P_{k}]\big{|}\leq C\big{(}h^{2}(1+|P_{k}|^{2m})e^{-\lambda(T-t_{k})}$
$\displaystyle\;\;\;\;+he^{-\lambda(T-t_{k})}\mathbb{E}\big{[}(n_{G}(Q_{c,k+1})\cdot\hat{P}_{k+1})(1+|P_{k+1}|^{2m})I_{G^{c}}(\hat{Q}_{k+1})\big{|}Q_{k},P_{k}\big{]}\big{)},$

where $u(t,q,p)$ is the solution of (2.2)-(2.4) with $b(q,p)=-\nabla
U(q)-\gamma p$ and $C,\lambda>0$ are constants independent of $h$, $T$ and
$P_{k}$.

The proof follows exactly the same steps as in the proof of Lemma 4.2 with the
modification that Assumption 2.6 is used for derivatives’ decay to bound the
terms in the Taylor expansion.

######  Lemma 4.6.

Let Assumptions 2.1 and 2.5 hold. For sufficiently small $h>0$ the following
inequality holds for any integer $m\geq 1$:

$\mathbb{E}\bigg{(}\sum\limits_{k=0}^{N-1}e^{\lambda
t_{k}}(n_{G}({Q}_{c,k+1})\cdot\hat{P}_{k+1})(1+|P_{k+1}|^{2m})I_{G^{c}}\big{(}\hat{Q}_{k+1}\big{)}\bigg{)}\leq
Ce^{\lambda T},$ (4.28)

where $\hat{Q}_{k+1}=Q_{k}+\hat{P}_{k+1}h$ and ${Q}_{c,k+1}$, ${Q}_{k+1}$,
$\hat{P}_{k+1}$, $k=0,\dots,N-1$, are according to the Markov chain
constructed via scheme [PAc] and $C>0$ is a constant independent of $T$ and
$h$.

###### Proof.

If in (4.10) $g(t,q,p)=\mathbb{E}\big{[}e^{\lambda
t_{k}}(n_{G}(Q_{c})\cdot\hat{P}_{1})(1+|P_{1}|^{2m})I_{G^{c}}(\hat{Q}_{1})|\;q,p\big{]}$
then the solution of (4.10)-(4.11) is

$\displaystyle v(t_{0},q,p)$
$\displaystyle=\mathbb{E}\bigg{(}\sum_{k=0}^{N-1}g(t_{k},Q_{k},P_{k})\big{|}Q_{0}=q,P_{0}=p\bigg{)}$
$\displaystyle=\mathbb{E}\bigg{(}\sum\limits_{k=0}^{N-1}\mathbb{E}\Big{(}e^{\lambda
t_{k}}(n_{G}(Q_{c,k+1})\cdot\hat{P}_{k+1})(1+|P_{k+1}|^{2m})I_{G^{c}}\big{(}\hat{Q}_{k+1}\big{)}\big{|}Q_{k},P_{k}\Big{)}\Big{|}Q_{0}=q,P_{0}=p\bigg{)}$
$\displaystyle=\mathbb{E}\bigg{(}\sum\limits_{k=0}^{N-1}e^{\lambda
t_{k}}\Big{(}n_{G}(Q_{c,k+1})\cdot\hat{P}_{k+1}\Big{)}\big{(}1+|P_{k+1}|^{2m}\big{)}I_{G^{c}}\big{(}\hat{Q}_{k+1}\big{)}\Big{|}Q_{0}=q,P_{0}=p\bigg{)}.$

If we can find a function $V(t,q,p)$ with $g(t,q,p)$, satisfying
(4.10)-(4.11), such that

$g(t,q,p)\geq\mathbb{E}\big{[}e^{\lambda
t}(n_{G}(Q_{c})\cdot\hat{P}_{1})(1+|P_{1}|^{2m})I_{G^{c}}(\hat{Q}_{1})|q,p\big{]},$

then the inequality $v(t_{0},Q_{0},P_{0})\leq V(t_{0},Q_{0},P_{0})$ holds.

As in the proof of Lemma 4.3, we consider the boundary zone $G_{R_{0}}$.
Introduce the function

$V(t,q,p)=\begin{cases}0,&(t,q,p)\in\\{T\\}\times\bar{G}\times\mathbb{R}^{d},\\\
Ae^{\lambda t}|p|^{2\ell}+L(e^{\lambda T}-e^{\lambda t})&\\\ +(B(q)\cdot
p)(1+|p|^{2m})e^{\lambda
t},&(t,q,p)\in[0,T-h]\times\bar{G}\times\mathbb{R}^{d},\end{cases}$

where $A=\max_{q\in\bar{G}}|B(q)|$; integer $\ell\geq m+1$ will be chosen
later in the proof; $L>0$ is a constant, which value will also be chosen later
in the proof; and $B(q)$ is the same as in the proof of Lemma 4.3 (i.e.,
$B(q)=n(q^{\pi})$ for $q\in G_{R_{0}}$ and $q^{\pi}$ is projection of $q$ on
$\partial G$). Note that $V(t,q,p)\geq 0$.

Using Taylor’s expansion, we get for any $l\geq 1$

$\displaystyle\mathbb{E}|P_{1}|^{2l}$
$\displaystyle=\mathbb{E}|\hat{P}_{1}|^{2l}=|p|^{2l}-2lh\gamma(p\cdot
p)|p|^{2l-2}-2lh(\nabla U(q)\cdot p)|p|^{2l-2}$
$\displaystyle+\frac{2\gamma}{\beta}lh(2(l-1)+d)|p|^{2l-2}+r_{1},$

where $r_{1}$ satisfies

$|r_{1}|\leq Cl^{4}h^{2}(1+|p|^{2l})$ (4.29)

with $C>0$ being independent of $h,$ $p$ and $l$. Expanding
$e^{\lambda(t+h)}$, using the generalised Young inequality twice and then
simplifying the expression, we obtain

$\displaystyle\mathbb{E}e^{\lambda(t+h)}|P_{1}|^{2\ell}-e^{\lambda
t}|p|^{2\ell}$ (4.30) $\displaystyle\leq e^{\lambda t}\left[-\ell
h\gamma|p|^{2\ell}+\ell
h\gamma\left[\frac{2\ell}{\beta(\ell-1)}(2(\ell-1)+d)\right]^{\ell}+(\lambda+1)h|p|^{2\ell}+h(2\ell|\nabla
U(q)|)^{2l}\right]+r_{2},$

where $r_{2}$ satisfies (4.29) with $l=\ell$.

Applying similar arguments used to obtain (4.24), we get

$\mathbb{E}[(B(Q_{1})\cdot P_{1})(1+|P_{1}|^{2m})]=(B(q)\cdot
p)(1+|p|^{2m})-2\mathbb{E}\big{[}(\hat{P}_{1}\cdot
n_{G}(Q_{c}))(1+|P_{1}|^{2m})I_{G^{c}}(\hat{Q}_{1})\big{]}+r_{3},$ (4.31)

where $r_{3}$ satisfies

$|r_{3}|\leq Ch(1+|p|^{2m})$ (4.32)

with $C>0$ being independent of $h$ and $p$ (it depends on $m)$.

From (4.30) and (4.31) and using $e^{\lambda(t+h)}\geq e^{\lambda t}+\lambda
he^{\lambda t},$ it follows that

$\displaystyle\mathbb{T}V-V=\mathbb{E}Ae^{\lambda(t+h)}|P_{1}|^{2\ell}+L(e^{\lambda
T}-e^{\lambda(t+h)})+\mathbb{E}[(B(Q_{1})\cdot
P_{1})(1+|P_{1}|^{2m})]e^{\lambda(t+h)}$ $\displaystyle-Ae^{\lambda
t}|p|^{2\ell}-L(e^{\lambda T}-e^{\lambda t})-(B(q)\cdot p)(1+|p|^{2m})$
$\displaystyle\leq$ $\displaystyle Ae^{\lambda t}\left[-\ell
h\gamma|p|^{2\ell}+\ell
h\gamma\left[\frac{2\ell}{\beta(\ell-1)}(2(\ell-1)+d)\right]^{\ell}+(\lambda+1)h|p|^{2\ell}+h(2\ell|\nabla
U(q)|)^{2l}\right]$ $\displaystyle-Le^{\lambda t}\lambda
h-2\mathbb{E}\big{[}(\hat{P}_{1}\cdot
n_{G}(Q_{c}))(1+|P_{1}|^{2m})I_{G^{c}}(\hat{Q}_{1})\big{]}e^{\lambda
t}+r_{2}+r_{4},$

where $r_{4}$ satisfies (4.32). Choosing a sufficiently large $\ell$ and,
depending on the choice of $\ell,$ choosing a sufficiently small $h>0$ and
large $L$, we get $g(t,q,p)=-(\mathbb{T}V-V)\geq\mathbb{E}\big{[}e^{\lambda
t}(n_{G}(Q_{c})\cdot\hat{P}_{1})(1+|P_{1}|^{2m})I_{G^{c}}(\hat{Q}_{1})\big{]}$.
Hence, the lemma is proved.

∎

###### Proof of Theorem 3.2.

Using Assumption 2.6 and Lemma 4.5, we have

$\displaystyle|\mathbb{E}\varphi(Q_{N},P_{N})-\bar{\varphi}|\leq|\mathbb{E}(\varphi(Q(T),P(T)))-\bar{\varphi}|+|\mathbb{E}(\varphi(Q_{N},P_{N}))-\mathbb{E}(\varphi(Q(T),P(T)))|$
$\displaystyle\leq C\exp(-\lambda
T)+\Big{|}\mathbb{E}\Big{(}\sum\limits_{k=0}^{N-1}\mathbb{E}\big{(}u(t_{k+1},Q_{k+1},P_{k+1})-u(t_{k},Q_{k},P_{k})|Q_{k},P_{k}\big{)}\Big{)}\Big{|}$
$\displaystyle\leq C\exp(-\lambda
T)+Ch^{2}\mathbb{E}(|P_{k}|^{2m})\sum\limits_{k=0}^{N-1}e^{-\lambda(T-t_{k})}$
$\displaystyle+Che^{-\lambda
T}\mathbb{E}\Big{(}\sum\limits_{k=0}^{N-1}e^{\lambda
t_{k}}(n_{G}({Q}_{c,k+1})\cdot\hat{P}_{k+1})(1+|P_{k+1}|^{2m})I_{G^{c}}(\hat{Q}_{k+1})\Big{)},$

which together with Lemmas 4.4 and 4.6 gives the desired result. ∎

###  4.3 One-step lemma for scheme [OBAcBO]

Consider $G=(0,\infty)\times R^{d-1}.$ In this case the one step of scheme
[OBAcBO] from Section 3.2 takes the form:

$\displaystyle\mathcal{P}_{1}=P_{0}e^{-\gamma
h/2}+\sqrt{\frac{\sigma^{2}}{2\gamma}(1-\mathrm{e}^{-\gamma
h})}\xi,\,\,\mathcal{P}_{2}=\mathcal{P}_{1}+\frac{h}{2}f(Q_{0}),\,\,\mathcal{Q}=Q_{0}+h\mathcal{P}_{2};$
$\displaystyle\textbf{If }\mathcal{Q}^{1}<0,\ \
\tau=-\frac{Q_{0}^{1}}{\mathcal{P}_{2}^{1}},\
\mathcal{P}_{3}^{1}=-\mathcal{P}_{2}^{1},\
\mathcal{Q}_{c}^{1}=0=Q_{0}^{1}-\tau\mathcal{P}_{3}^{1}=Q_{0}^{1}+\tau\mathcal{P}_{2}^{1},\
$
$\displaystyle\mathcal{Q}_{c}^{i}=Q_{0}^{i}+\tau\mathcal{P}_{3}^{i}=Q_{0}^{i}+\tau\mathcal{P}_{2}^{i}\
,\ i\neq
1;\,\,Q^{1}=(h-\tau)\mathcal{P}_{3}^{1}=\left(h-\frac{Q_{0}^{1}}{\mathcal{P}_{3}^{1}}\right)\mathcal{P}_{3}^{1},\
$ $\displaystyle
Q^{i}=\mathcal{Q}_{c}^{i}\mathcal{+}(h-\tau)\mathcal{P}_{3}^{i}=Q_{0}^{i}+h\mathcal{P}_{2}^{i},\
i\neq 1;$ $\displaystyle\textbf{else }Q=\mathcal{Q},\
\mathcal{P}_{3}=\mathcal{P}_{2};$
$\displaystyle\mathcal{P}_{4}=\mathcal{P}_{3}+\frac{h}{2}f(Q),\;\;P=\mathcal{P}_{4}e^{-\gamma
h/2}+\sqrt{\frac{\sigma^{2}}{2\gamma}(1-\mathrm{e}^{-\gamma h})}\zeta.$

Here the components $\xi^{i}$ of $\xi$ and $\zeta^{i}$ of $\zeta$ are i.i.d.
with the standard normal distribution and $\mathcal{Q}^{i}$, $i=1,\dots,d$,
denotes $i-$th component of $\mathcal{Q}$ and likewise for $q$ and
$\mathcal{P}_{j}$, $j=2,3$.

In this section we need the following assumptions. The potential $U(q)$
satisfies Assumption 2.5, and $\nabla U(q)$ is globally Lipschitz and
derivatives of $\nabla U(q)$ are bounded. To ensure ergodicity (see [34] and
[41][Section 9.1.5]), we need dissipativity of the extended generator of the
Markov process $(Q(t),P(t))$ which can be guaranteed by assuming that there
exists an $\alpha_{1}>0$ and $0<\alpha_{2}<1$ such that

$\frac{1}{2}(\nabla
U(q),q)\geq\alpha_{2}U(q)+\gamma^{2}\frac{\alpha_{2}(2-\alpha_{2})}{8(1-\alpha_{2})}|q|^{2}+\alpha_{1}.$
(4.33)

Further, we assume that the function $\varphi(q,p)\in
C^{6,6}(\bar{G}\times\mathbb{R}^{d})$ and all its existing derivatives have
growth at most polynomially in $q$ and $p$. As usual for proofs of higher
order numerical methods, we also require that the solution of the backward
Kolmogorov equation (2.2)-(2.4) with $b(q,p)=-\nabla U(q)-\gamma p$ is
sufficiently smooth and, in particular, satisfies the condition that there are
positive constants $C$ and $\lambda$ independent of $T$ such that the
following bound holds for some $m_{1},m_{2}\geq 1$ (cf. (2.10) in Assumption
2.6):

$\sum\limits_{l=1}^{3}\sum\limits_{i+|j|=l}|D_{t}^{i}D_{q}^{j}u(t,q,p)|+\sum\limits_{l=1}^{6}\sum\limits_{|j|=l}|D_{p}^{j}u(t,q,p)|\leq
C(1+|q|^{2m_{1}}+|p|^{2m_{2}})e^{-\lambda(T-t)},$

where $C>0$ is independent of $T,$ $q$ and $p$.

By standard means it is not difficult to establish that under the stated
conditions the moments of $Q_{k}$ and $P_{k}$ from scheme [OBAcBO] are bounded
uniformly in time (cf. Lemma 4.4).

######  Lemma 4.7.

Assume that the above conditions hold and that $(Q_{0},P_{0})$ is random with
a density so that all the required moments with respect to this density exist,
then the one-step error of scheme [OBAcBO] has the following estimate

$|\mathbb{E}\left[u(t+h,Q,P)-u(t,Q_{0},P_{0})\right]|\leq
Ce^{-\lambda(T-t)}\left(h^{3}+h^{2}\mathbb{E}\left[I_{G^{c}}(\mathcal{Q})\right]\right),$
(4.34)

where $C>0$ is independent of $T,$ $t,$ and $h$.

Discussion The number of steps is $T/h$ while we expect (cf. Lemmas 4.3 and
4.6) that the average number of steps with collisions is independent of $h$.
This together with (4.34) leads to the claim that [OBAcBO] (and analogously
other schemes of Section 3.2) has second weak convergence order in $h$, both
at finite time and for the ergodic limit. We recall again that the
deterministic counterpart of [OBAcBO] has first order convergence and second
order is achieved here thanks to an average (see details in the proof of this
lemma) of the collision time $\tau$ being equal to $h/2$ up to terms of order
$\mathcal{O}(h^{2})$.

###### Proof.

We write the one step error as

$\displaystyle R:=\mathbb{E}\left[u(t+h,Q,P)-u(t,Q_{0},P_{0})\right]$
$\displaystyle=\mathbb{E}\left[u(t+h,Q,P)-u(t+h/2,Q,\mathcal{P}_{3})\right]+\mathbb{E}\left[u(t+h/2,Q,\mathcal{P}_{3})-u(t+h/2,\mathcal{Q}_{c},\mathcal{P}_{3})\right]I_{G^{c}}(\mathcal{Q})$
$\displaystyle+\mathbb{E}\left[u(t+h/2,\mathcal{Q}_{c},\mathcal{P}_{2})-u(t+h/2,Q_{0},\mathcal{P}_{2})\right]I_{G^{c}}(\mathcal{Q})$
$\displaystyle+\mathbb{E}\left[u(t+h/2,Q,\mathcal{P}_{2})-u(t+h/2,Q_{0},\mathcal{P}_{2})\right]I_{G}(\mathcal{Q})+\mathbb{E}\left[u(t+h/2,Q_{0},\mathcal{P}_{2})-u(t+h/2,Q_{0},\mathcal{P}_{1})\right]$
$\displaystyle+\mathbb{E}\left[u(t+h/2,Q_{0},\mathcal{P}_{1})-u(t+h/2,Q_{0},P_{0})\right]+\left[u(t+h/2,Q_{0},P_{0})-u(t,Q_{0},P_{0})\right].$

Appropriately applying Taylor’s expansion to each of the brackets in the above
re-writing of $R$, we arrive at

$\displaystyle R$ $\displaystyle=$
$\displaystyle\mathbb{E}\left[(\frac{h}{2}-\tau)\mathcal{P}_{3}\cdot\nabla_{q}u(t+h/2,\mathcal{Q}_{c},\mathcal{P}_{3})I_{G^{c}}(\mathcal{Q})\right]+\mathbb{E}\left[(\tau-\frac{h}{2})\mathcal{P}_{2}\cdot\nabla_{q}u(t+h/2,Q_{0},\mathcal{P}_{2})I_{G^{c}}(\mathcal{Q})\right]$
$\displaystyle+r_{1}+r_{2},$

where $|r_{1}|\leq Ce^{-\lambda(T-t)}h^{3}$ and $|r_{2}|\leq
Ce^{-\lambda(T-t)}h^{2}\mathbb{E}\left[I_{G^{c}}(\mathcal{Q})\right]$.
Introduce $\pi$ so that $\pi\mathcal{P}_{3}=\mathcal{P}_{2},$ i.e.,
$\pi\mathcal{P}_{3}=(-\mathcal{P}_{3}^{1},\mathcal{P}_{3}^{2},\ldots,\mathcal{P}_{3}^{d}).$
We have

$\displaystyle R$ $\displaystyle=$
$\displaystyle\mathbb{E}\left[(\frac{h}{2}-\tau)\left(\mathcal{P}_{3}\cdot\nabla_{q}u(t+h/2,\mathcal{Q}_{c},\mathcal{P}_{3})-\pi\mathcal{P}_{3}\cdot\nabla_{q}u(t+h/2,\mathcal{Q}_{c},\pi\mathcal{P}_{3})\right)I_{G^{c}}(\mathcal{Q})\right]$
$\displaystyle+r_{1}+r_{3},$

where $r_{3}$ satisfies the same estimate as $r_{2}$ above.

Let us first consider the following component of (4.3). Using the expression
for $\tau,$ we get

$\displaystyle R_{1}$ $\displaystyle=$
$\displaystyle\mathbb{E}\tau\mathcal{P}_{3}^{1}\left[\frac{\partial}{\partial
q^{1}}u(t+h/2,0,Q_{0}^{2},\ldots,Q_{0}^{d},\mathcal{P}_{3})+\frac{\partial}{\partial
q^{1}}u(t+h/2,0,Q_{0}^{2},\ldots,Q_{0}^{d},\pi\mathcal{P}_{3})\right]I_{G^{c}}(\mathcal{Q})$
$\displaystyle=\mathbb{E}\left[Q_{0}^{1}v_{1}(Q_{0}^{2},\ldots,Q_{0}^{d},\mathcal{P}_{3})I_{G^{c}}(\mathcal{Q})\right],$

where

$v_{1}(q^{2},\ldots,q^{d},p):=\frac{\partial}{\partial
q^{1}}u(t+h/2,0,q^{2},\ldots,q^{d},p)+\frac{\partial}{\partial
q^{1}}u(t+h/2,0,q^{2},\ldots,q^{d},\pi p).$

Let $\rho(q,p)$ be the joint density for $Q_{0},$ $\mathcal{P}_{3}.$ Assume
that all the required moments with respect to this density exist. We have

$R_{1}=\frac{\int_{\mathbb{R}^{2d-2}}\int_{0}^{\infty}\int_{0}^{p^{1}h}q^{1}v_{1}(q^{2},\ldots,q^{d},p)\rho(q,p)dq^{1}dp^{1}dq^{2}\cdots
dq^{d}dp^{2}\cdots
dp^{d}}{\int_{R^{2d-2}}\int_{0}^{\infty}\int_{0}^{p^{1}h}\rho(q,p)dq^{1}dp^{1}dq^{2}\cdots
dq^{d}dp^{2}\cdots dp^{d}}.$

Note that

$\int_{0}^{p^{1}h}q^{1}\rho(q,p)dq^{1}=\frac{1}{2}\left(p^{1}\right)^{2}h^{2}\rho(0,q^{2},\ldots,q^{d},p)+{\cal
O}(h^{3}),$
$\int_{0}^{p^{1}h}\rho(q,p)dq^{1}=p^{1}h\rho(0,q^{2},\ldots,q^{d},p)+{\cal
O}(h^{2}).$

Hence

$R_{1}=\frac{h}{2}\frac{\int_{\mathbb{R}^{2d-2}}\int_{0}^{\infty}\left(p^{1}\right)^{2}\rho(p,0,q^{2},\ldots,q^{d})U_{1}(p,q^{2},\ldots,q^{d})dp^{1}dq^{2}\cdots
dq^{d}dp^{2}\cdots dp^{d}+{\cal
O}(h)}{\int_{R^{2d-2}}\int_{0}^{\infty}p^{1}\rho(p,0,q^{2},\ldots,q^{d})dp^{1}dq^{2}\cdots
dq^{d}dp^{2}\cdots dp^{d}+{\cal O}(h)}.$

Next consider the following component of (4.3):

$\displaystyle R_{2}$ $\displaystyle=$
$\displaystyle\frac{h}{2}\mathbb{E}\bigg{(}\mathcal{P}_{3}^{1}\left[\frac{\partial}{\partial
q^{1}}u(t+h/2,0,Q_{0}^{2},\ldots,Q_{0}^{d},\mathcal{P}_{3})+\frac{\partial}{\partial
q^{1}}u(t+h/2,0,Q_{0}^{2},\ldots,Q_{0}^{d},\pi\mathcal{P}_{3})\right]I_{G^{c}}(\mathcal{Q})\bigg{)}$
$\displaystyle=$
$\displaystyle\frac{h}{2}\frac{\int_{\mathbb{R}^{2d-2}}\int_{0}^{\infty}\int_{0}^{p^{1}h}p^{1}v_{1}(q^{2},\ldots,q^{d},p)\rho(q,p)dq^{1}dp^{1}dq^{2}\cdots
dq^{d}dp^{2}\cdots
dp^{d}}{\int_{R^{2d-2}}\int_{0}^{\infty}\int_{0}^{p^{1}h}\rho(q,p)dq^{1}dp^{1}dq^{2}\cdots
dq^{d}dp^{2}\cdots dp^{d}}$ $\displaystyle=$
$\displaystyle\frac{h}{2}\frac{\int_{\mathbb{R}^{2d-2}}\int_{0}^{\infty}\left(p^{1}\right)^{2}\rho(0,q^{2},\ldots,q^{d},p)v_{1}(q^{2},\ldots,q^{d},p)dp^{1}dq^{2}\cdots
dq^{d}dp^{2}\cdots dp^{d}+{\cal
O}(h)}{\int_{R^{2d-2}}\int_{0}^{\infty}p^{1}\rho(0,q^{2},\ldots,q^{d},p)dp^{1}dq^{2}\cdots
dq^{d}dp^{2}\cdots dp^{d}+{\cal O}(h)}.$

Thus $|R_{1}-R_{2}|\leq
Ch^{2}e^{-\lambda(T-t)}\mathbb{E}\left[I_{G^{c}}(\mathcal{Q})\right]$. The
terms for $i\neq 1$

$\displaystyle A_{i}$
$\displaystyle=\mathbb{E}\left(\tau-h/2\right)\mathcal{P}_{3}^{i}$
$\displaystyle\times\left[\frac{\partial}{\partial
q^{i}}u(t+h/2,0,Q_{0}^{2},\ldots,Q_{0}^{d},\mathcal{P}_{3})-\frac{\partial}{\partial
q^{i}}u(t+h/2,0,Q_{0}^{2},\ldots,Q_{0}^{d},\pi\mathcal{P}_{3})\right]I_{G^{c}}(\mathcal{Q})$

are analyzed analogously. ∎

##  5 Numerical tests

In this section we confirm theoretical results for the numerical integrators
from Section 3 in both the finite-time and ergodic case.

The following notation is used in this section. The Monte Carlo estimator for
$\mathbb{E}(\varphi(Q_{N},P_{N}))$ is
$\check{\varphi}_{M}=\frac{1}{M}\sum_{k=1}^{M}\varphi(Q_{N}^{(k)},P_{N}^{(k)})$,
where $(Q_{N}^{(k)},P_{N}^{(k)})$ are independent realizations of
$(Q_{N},P_{N})$; $\chi$ denotes the number of collisions along a trajectory
and $\check{\chi}_{M}$ is the Monte-Carlo estimator for the average number of
collisions; and ${D}_{M}$ denotes the usual sample variance of Monte Carlo
estimators.

######  Experiment 5.1.

In this experiment we consider finite-time weak convergence. Consider Example
2.2 with $G:=\\{q^{2}_{1}+q^{2}_{2}<4\\}]$ and
$U(q_{1},q_{2})=-q_{1}^{2}-q_{2}^{2}$, $\beta=1$, $\alpha=0.25$ and the final
time $T=4$. Hence, the solution is $u(0,q,p)=\exp(-|p|^{2}/2+|q|^{2}-2)$. The
point at which $u(0,q,p)$ is evaluated in the experiments is
$(q_{1},q_{2})=(1,1)$ and $(p_{1},p_{2})=(-0.1,-0.1)$. To compute the global
error $e=\check{\varphi}_{M}-u(0,(1,1),(-0.1,-0.1))$, we calculate the exact
solution $u(0,(1,1),(-0.1,-0.1))=0.99005$ (5 d.p.). We test the first-order
scheme [PAc] and two second-order schemes [OBAcBO] and [BAcOAcB]. Tables 5.2
and 5.3 and Fig. 5.1 confirm second-order accuracy of [OBAcBO] and [BAcOAcB].
The first-order method [PAc] had slightly higher accuracy than first order but
clearly lower than second order. In this experiment $O$ step in [OBAcBO] and
[BAcOAcB] schemes is

$O(p;h,\xi,\alpha):=pe^{\alpha
h}+\sqrt{\frac{\sigma^{2}}{2\alpha}(\mathrm{e}^{2\alpha h}-1)}\xi.$

Table 5.1: Experiment 5.1. The global error $e$ and estimated average number
of collisions $\check{\chi}_{M}$ in interval $[0,4]$ for [PAc] scheme.

$h$ | $M$ | $\check{\varphi}_{M}\pm 2\sqrt{D_{M}}$ | $e$ | $\check{\chi}_{M}\pm 2\sqrt{D_{M}}$ |
---|---|---|---|---|---
0.08 | $4\times 10^{8}$ | $1.13004\pm 3.5\times 10^{-4}$ | 0.139994 | $3.739\pm 1.3\times 10^{-4}$ |
0.04 | $4\times 10^{8}$ | $1.03974\pm 3.2\times 10^{-4}$ | $0.049692$ | $3.679\pm 1.2\times 10^{-4}$ |
0.02 | $4\times 10^{8}$ | $1.00826\pm 3.1\times 10^{-4}$ | 0.018207 | $3.660\pm 1.2\times 10^{-4}$ |
0.01 | $4\times 10^{8}$ | $0.99762\pm 3.1\times 10^{-4}$ | 0.007568 | $3.654\pm 1.2\times 10^{-4}$ |

Table 5.2: Experiment 5.1. The global error $e$ and the estimated average
number of collisions $\check{\chi}_{M}$ in interval $[0,4]$ for [OBAcBO]
scheme.

$h$ | $M$ | $\check{\varphi}_{M}\pm 2\sqrt{D_{M}}$ | $e$ | $\check{\chi}_{M}\pm 2\sqrt{D_{M}}$ |
---|---|---|---|---|---
0.08 | $4\times 10^{8}$ | $1.08582\pm 3.4\times 10^{-4}$ | 0.095775 | $3.714\pm 1.2\times 10^{-4}$ |
0.04 | $4\times 10^{8}$ | $1.01649\pm 3.2\times 10^{-4}$ | $0.026441$ | $3.665\pm 1.2\times 10^{-4}$ |
0.02 | $4\times 10^{8}$ | $0.99683\pm 3.1\times 10^{-4}$ | 0.006776 | $3.653\pm 1.2\times 10^{-4}$ |
0.01 | $4\times 10^{8}$ | $0.99164\pm 3.1\times 10^{-4}$ | 0.001591 | $3.650\pm 1.2\times 10^{-4}$ |

Table 5.3: Experiment 5.1. The global error $e$ and the estimated average
number of collisions $\check{\chi}_{M}$ in interval $[0,4]$ for [BAcOAcB]
scheme.

$h$ | $M$ | $\check{\varphi}_{M}\pm 2\sqrt{D_{M}}$ | $e$ | $\check{\chi}_{M}\pm 2\sqrt{D_{M}}$ |
---|---|---|---|---|---
0.08 | $4\times 10^{8}$ | $1.07393\pm 3.4\times 10^{-4}$ | 0.083880 | $3.710\pm 1.2\times 10^{-4}$ |
0.04 | $4\times 10^{8}$ | $1.01308\pm 3.2\times 10^{-4}$ | $0.023029$ | $3.664\pm 1.2\times 10^{-4}$ |
0.02 | $4\times 10^{8}$ | $0.99586\pm 3.1\times 10^{-4}$ | 0.005807 | $3.653\pm 1.2\times 10^{-4}$ |
0.01 | $4\times 10^{8}$ | $0.99148\pm 3.1\times 10^{-4}$ | 0.001431 | $3.650\pm 1.2\times 10^{-4}$ |

Figure 5.1: Experiment 5.1. The global error $e$ for the first-order method
[PAc] and for [OBAcBO] and [BAcOAcB]. Error bars indicate the Monte Carlo
error.

Table 5.4: Experiment 5.2. Dependence of $\check{\varphi}_{M}$ on $h$ for the
[AcBOBAc] and [AcOBOAc] schemes.

|  | [AcBOBAc] | [AcOBOAc]
---|---|---|---
$h$ | $M$ | $\check{\varphi}_{M}\pm 2\sqrt{D_{M}}$ | $\check{\varphi}_{M}\pm 2\sqrt{D_{M}}$
0.5 | $5\times 10^{6}$ | $1.32450\pm 7.8852\times 10^{-4}$ | $1.34389\pm 8.0906\times 10^{-4}$
0.4 | $5\times 10^{6}$ | $1.30362\pm 7.7085\times 10^{-4}$ | $1.31691\pm 7.8433\times 10^{-4}$
0.2 | $5\times 10^{6}$ | $1.27357\pm 7.4647\times 10^{-4}$ | $1.27655\pm 7.4981\times 10^{-4}$
0.1 | $5\times 10^{7}$ | $1.26542\pm 2.3408\times 10^{-4}$ | $1.26605\pm 2.3426\times 10^{-4}$
0.08 | $5\times 10^{7}$ | $1.26437\pm 2.3383\times 10^{-4}$ | $1.26497\pm 2.3400\times 10^{-4}$

Table 5.5: Experiment 5.2. Dependence of $\check{\varphi}_{M}$ on $h$ for
[BAcOAcB] and [BOAcOB] schemes.

| [BAcOAcB] |  | [BOAcOB] |
---|---|---|---|---
$h$ | $M$ | $\check{\varphi}_{M}\pm 2\sqrt{D_{M}}$ | $M$ | $\check{\varphi}_{M}\pm 2\sqrt{D_{M}}$
0.5 | $5\times 10^{6}$ | $1.30502\pm 7.5645\times 10^{-4}$ | $5\times 10^{6}$ | $1.42944\pm 8.6166\times 10^{-4}$
0.4 | $5\times 10^{6}$ | $1.28807\pm 7.4837\times 10^{-4}$ | $5\times 10^{6}$ | $1.36950\pm 8.1676\times 10^{-4}$
0.2 | $5\times 10^{6}$ | $1.26871\pm 7.4072\times 10^{-4}$ | $5\times 10^{6}$ | $1.28988\pm 7.5830\times 10^{-4}$
0.1 | $5\times 10^{8}$ | $1.26403\pm 7.3866\times 10^{-5}$ | $5\times 10^{7}$ | $1.26922\pm 2.3492\times 10^{-4}$
0.08 | $5\times 10^{8}$ | $1.26352\pm 7.3854\times 10^{-5}$ | $5\times 10^{7}$ | $1.26685\pm 2.3440\times 10^{-4}$

Table 5.6: Experiment 5.2. Dependence of $\check{\varphi}_{M}$ on $h$ for
[OAcBAcO] and [OBAcBO] schemes.

| [OAcBAcO] |  | [OBAcBO] |
---|---|---|---|---
$h$ | $M$ | $\check{\varphi}_{M}\pm 2\sqrt{D_{M}}$ | $M$ | $\check{\varphi}_{M}\pm 2\sqrt{D_{M}}$
0.5 | $5\times 10^{6}$ | $1.22339\pm 7.0554\times 10^{-4}$ | $5\times 10^{6}$ | $1.42944\pm 8.4456\times 10^{-4}$
0.4 | $5\times 10^{6}$ | $1.23709\pm 7.1800\times 10^{-4}$ | $5\times 10^{6}$ | $1.3695\pm 8.4456\times 10^{-4}$
0.2 | $5\times 10^{6}$ | $1.25579\pm 7.3281\times 10^{-4}$ | $5\times 10^{6}$ | $1.28988\pm 7.5295\times 10^{-4}$
0.1 | $5\times 10^{8}$ | $1.26094\pm 7.3689\times 10^{-5}$ | $5\times 10^{7}$ | $1.26922\pm 2.3457\times 10^{-4}$
0.08 | $5\times 10^{8}$ | $1.26151\pm 7.3736\times 10^{-5}$ | $5\times 10^{7}$ | $1.26685\pm 2.3417\times 10^{-4}$

Figure 5.2: Experiment 5.2. Dependence of the global error $e$ on $h$ for
[AcBOBAc], [AcOBOAc], [BAcOAcB], [BOAcOB], [OAcBAcO] and [OBAcBO]. Error bars
indicate the Monte Carlo error.

######  Experiment 5.2.

In this experiment we test convergence of second-order integators [AcBOBAc],
[AcOBOAc], [BAcOAcB], [BOAcOB], [OAcBAcO] and [OBAcBO] to ergodic limits.

Take $G:=(1,\infty)$ and $\partial G:=\\{1\\}$. Consider (ECLD) with the
potential

$\displaystyle U(q)=q^{2}/2$

and $\beta=1$, $\gamma=1$. For numerical simulation, we choose $T=20$ as a
sufficiently large time to achieve accurate values of the ergodic limit of
$\varphi(q)=q^{2}/2$ with respect to this choice of (ECLD). The exact value of
$\bar{\varphi}$ is $1.2626$ (4 d.p.). The starting point of Markov chain is
$Q_{0}=2$ and $P_{0}=-0.1$. We denote
$e:=|\check{\varphi}_{M}-\bar{\varphi}|$. The results are presented in Tables
5.4-5.6 and Fig. 5.2, which clearly show second order convergence for all the
integrators. As expected [24, 25], [BAcOAcB] and [OAcBAcO] are the most
accurate schemes.

Figure 5.3: Experiment 5.3. Dependence of the global error $e$ on $h$ for
[OBAcBO] and [BAcOAcB]. Error bars indicate the Monte Carlo error.

Table 5.7: Experiment 5.3. The estimate $\check{\varphi}_{M}$ for [OBAcBO]
scheme.

$h$ | $M$ | $\check{\varphi}_{M}\pm 2\sqrt{D_{M}}$
---|---|---
0.1 | $5\times 10^{6}$ | $-4.09714\pm 1.2002\times 10^{-3}$
0.075 | $5\times 10^{7}$ | $-4.13173\pm 3.7841\times 10^{-4}$
0.05 | $5\times 10^{7}$ | $-4.15843\pm 3.7753\times 10^{-4}$
0.04 | $5\times 10^{7}$ | $-4.16630\pm 3.7727\times 10^{-4}$

Table 5.8: Experiment 5.3. The estimate $\check{\varphi}_{M}$ for [BAcOAcB]
scheme.

$h$ | $M$ | $\check{\varphi}_{M}\pm 2\sqrt{D_{M}}$
---|---|---
0.1 | $5\times 10^{6}$ | $-4.14494\pm 1.1933\times 10^{-3}$
0.075 | $5\times 10^{7}$ | $-4.16074\pm 3.7709\times 10^{-4}$
0.05 | $5\times 10^{7}$ | $-4.17151\pm 3.7709\times 10^{-4}$
0.04 | $5\times 10^{8}$ | $-4.17511\pm 1.1920\times 10^{-4}$

######  Experiment 5.3.

In this experiment we test convergence of second-order integators [BAcOAcB]
and [OBAcBO] to ergodic limits. Take G:= $\\{q_{1}^{2}+q_{2}^{2}<4\\}$.
Consider (ECLD) with the potential

$U(q_{1},q_{2})=\frac{1}{2}(q_{1}-q_{2})^{2}+\frac{q_{1}^{2}(q_{1}^{2}-12)}{12}+\frac{q_{2}^{2}(q_{2}^{2}-24)}{12}$

and $\beta=1$, $\gamma=4$. We choose

$\varphi(q_{1},q_{2})=U(q_{1},q_{2}).$

We note that the potential $U$ is not symmetric around zero and that its local
minima lie outside $G$. We are interested in approximating $\bar{\varphi}$. We
take $T=12$. We calculate the integral $\bar{\varphi}$ using integral2
function in Matlab with tolerance $10^{-10}$ which gives the reference value
$-4.18006$ (5 d.p.). The starting point of the Markov chain is
$(q_{1},q_{2})=(1,1)$ with $(p_{1},p_{2})=(-0.1,-0.1)$. We denote
$e:=|\check{\varphi}_{M}-\bar{\varphi}|$. The results are presented in Tables
5.7-5.8 and Fig. 5.3, which clearly show second-order convergence for both
integrators with [BAcOAcB] being more accurate.

Figure 5.4: Experiment 5.4. Dependence of the global error $e$ on $h$ for
[OBAcBO] and [BAcOAcB]. Error bars indicate the Monte Carlo error.

Table 5.9: Experiment 5.4. Dependence of global error $e$ on $h$ for [OBAcBO]
scheme with $T=20$.

$h$ | $M$ | $\check{\varphi}_{M}\pm 2\sqrt{D_{M}}$ | $e$ | $\check{\chi}_{M}\pm 2\sqrt{D_{M}}$
---|---|---|---|---
0.25 | $10^{4}$ | $1.697\pm 0.024$ | 2.102 | $32.52\pm 0.08$
0.2 | $10^{4}$ | $1.912\pm 0.023$ | 1.888 | $34.71\pm 0.10$
0.1 | $10^{4}$ | $2.860\pm 0.017$ | 0.940 | $53.01\pm 0.19$
0.08 | $10^{4}$ | $3.103\pm 0.014$ | $0.697$ | $61.16\pm 0.21$
0.05 | $10^{4}$ | $3.4514\pm 0.0092$ | 0.3486 | $78.58\pm 0.26$
0.04 | $10^{4}$ | $3.5569\pm 0.0077$ | 0.2431 | $86.92\pm 0.30$
0.025 | $10^{4}$ | $3.6859\pm 0.0058$ | 0.1141 | $102.86\pm 0.37$
0.02 | $10^{4}$ | $3.7238\pm 0.0053$ | 0.0762 | $108.89\pm 0.40$
0.01 | $10^{5}$ | $3.7785\pm 0.0014$ | 0.0215 | $120.97\pm 0.15$
0.008 | $10^{7}$ | $3.78550\pm 0.00013$ | 0.01450 | $123.12\pm 0.016$
0.004 | $10^{7}$ | $3.79630\pm 0.00013$ | 0.003670 | $126.33\pm 0.017$
0.002 | $10^{7}$ | $3.79907\pm 0.00013$ | 0.000928 | $127.24\pm 0.017$

Table 5.10: Experiment 5.4. Dependence of global error $e$ on $h$ for
[BAcOAcB] scheme with $T=20$.

$h$ | $M$ | $\check{\varphi}_{M}\pm 2\sqrt{D_{M}}$ | $e$ | $\check{\chi}_{M}\pm 2\sqrt{D_{M}}$
---|---|---|---|---
0.25 | $10^{4}$ | $1.712\pm 0.024$ | 2.088 | $33.52\pm 0.08$
0.2 | $10^{4}$ | $1.945\pm 0.023$ | 1.855 | $34.97\pm 0.10$
0.1 | $10^{4}$ | $2.868\pm 0.017$ | 0.932 | $53.68\pm 0.19$
0.08 | $10^{4}$ | $3.124\pm 0.014$ | 0.676 | $61.86\pm 0.21$
0.05 | $10^{4}$ | $3.4689\pm 0.0090$ | 0.331 | $79.92\pm 0.27$
0.04 | $10^{4}$ | $3.5690\pm 0.0076$ | 0.2310 | $88.07\pm 0.31$
0.025 | $10^{4}$ | $3.6859\pm 0.0058$ | 0.1141 | $103.70\pm 0.38$
0.02 | $10^{4}$ | $3.7292\pm 0.0051$ | 0.0708 | $109.94\pm 0.41$
0.01 | $10^{5}$ | $3.7798\pm 0.0014$ | 0.0202 | $121.56\pm 0.16$
0.008 | $10^{7}$ | $3.78663\pm 0.00013$ | 0.01337 | $123.538\pm 0.016$
0.004 | $10^{7}$ | $3.79644\pm 0.00013$ | 0.00356 | $129.494\pm 0.017$
0.002 | $10^{7}$ | $3.79910\pm 0.00013$ | 0.00090 | $127.279\pm 0.017$

Table 5.11: Experiment 5.4. Stability of [OBAcBO] and [BAcOAcB] schemes with
$T=20$ for large values of $h$.

| [BAcOAcB] | [OBAcBO] |
---|---|---|---
$h$ | $M$ | $\check{\varphi}_{M}\pm 2\sqrt{D_{M}}$ | $\check{\varphi}_{M}\pm 2\sqrt{D_{M}}$
0.25 | $10^{4}$ | $1.712\pm 0.024$ | $1.697\pm 0.024$
0.3 | $10^{4}$ | $1.588\pm 0.024$ | $--$
0.35 | $10^{4}$ | $1.496\pm 0.024$ | $--$
0.42 | $10^{4}$ | $1.469\pm 0.024$ | $--$
0.5 | $10^{4}$ | $--$ | $--$

######  Experiment 5.4.

We take $G=\\{q_{1}^{2}+q_{2}^{2}<4\\}$ and $\partial
G=\\{q_{1}^{2}+q_{2}^{2}=4\\}$. We sample from distribution given by

$\displaystyle\rho(q,p)=\frac{1}{\mathbb{Z}}e^{K|q|^{2}-|p|^{2}},\;\;\;\;q=(q_{1},q_{2})\in\bar{G},\;p=(p_{1},p_{2})\in\mathbb{R}^{d},$
(5.1)

and compute

$\displaystyle\bar{\varphi}=\int_{G}(q_{1}^{2}+q_{2}^{2})e^{K(q_{1}^{2}+q_{2}^{2})}dq_{1}dq_{2}$
(5.2)

via its Monte Carlo estimator based on two schemes [OBAcBO] and [BAcOAcB]
which are used to discretize (ECLD). We take $K=10$ in order to increase the
number of collisions with the boundary and the aim is to check if in this case
the two schemes [OBAcBO] and [BAcOAcB] maintain second order of convergence
asymptotically as the discretization step $h$ decreases. We take the initial
conditions $(q_{1},q_{2})=(0,0)$ and $(p_{1},p_{2})=(-0.1,-0.1)$.

We divide our analysis in two parts, in the first we are interested in this
experiment to investigate the asymptotic convergence of [OBAcBO] and [BAcOAcB]
as $h$ becomes smaller. Here it is crucial as our choice of $U(q)$ results in
large number of collisions in time period from $0$ to $T=20$. It is clear from
Tables 5.9 \- 5.10 and Fig. 5.4 that both numerical schemes [OBAcBO] and
[BAcOAcB] have $\mathcal{O}(h)$ convergence then $\mathcal{O}(h^{3/2})$
convergence and, as $h$ decreases, we clearly observe $\mathcal{O}(h^{2})$
convergence. In the second part, we study the large step size regime and
investigate which of the methods is more stable. It is clear from Table 5.11
that [BAcOAcB] is more stable than [OBAcBO]; it allows the use of step sizes
of around 0.42, whereas [OBAcBO] becomes unstable for a step size greater than
about 0.30.

##  6 Conclusions

We constructed and analysed a number of numerical methods for simulating
confined Langevin dynamics. The considered SDE model can be used in a variety
of applications from molecular dynamics to sampling from distributions with
compact support and optimization with constraints. We proved finite-time
convergence and convergence to ergodic limits for our first-order schemes. We
observed and justified (in the case of half-space) a surprising result that
the proposed [Ac, B, O] Störmer-Verlet-type splitting schemes are of second
weak order despite their deterministic counterparties are only of first order
due to the collision step Ac having local error of order ${\cal O}(h)$.

A number of extensions for this work are of interest. Those include further
testing of the proposed geometric integrators, e.g. applying them to train
neural networks with imposed constraints on weights. It is possible to extend
the developed approach to the case when particles are not only colliding with
the boundary of a domain but also with each other. Constructing efficient
numerical methods for Langevin equations with non-elastic collisions on
boundary is a possible topic for future work. Further theoretical analysis of
the constructed second-order integrators deserves a separate attention too.

## Acknowledgments

AS and MVT were supported by the Engineering and Physical Sciences Research
Council [grant number EP/X022617/1]. The authors express their gratitude to
the International Centre for Mathematical Sciences (ICMS) for its support via
the research-in-groups scheme. For the purpose of open access, the authors
have applied a CC BY public copyright licence to any Author Accepted
Manuscript version arising.

## References

  * [1] A. Benchérif-Madani and É. Pardoux.  A probabilistic formula for a Poisson equation with Neumann boundary condition.  Stoch. Anal. Appl., 27(4):739–746, 2009.
  * [2] D. M. Blei, A. Y. Ng, and M. I. Jordan.  Latent Dirichlet allocation.  J. Mach. Learn. Res., 3:993–1022, 2003.
  * [3] S. D. Bond and B. Leimkuhler.  Molecular dynamics and the accuracy of numerically computed averages.  Acta Numerica, 16:1, 2007.
  * [4] M. Bossy, E. Gobet, and D. Talay.  A symmetrized Euler scheme for an efficient approximation of reflected diffusions.  J. Appl. Prob., 41(3):877–889, 2004.
  * [5] M. Bossy and J.-F. Jabir.  On confined McKean Langevin processes satisfying the mean no-permeability boundary condition.  Stoch. Proc. Applic., 121(12):2751 – 2775, 2011.
  * [6] M. Bossy and J.-F. Jabir.  Lagrangian stochastic models with specular boundary condition.  J. Funct. Anal., 268(6):1309 – 1381, 2015.
  * [7] M. Bossy, J.-F. Jabir, and R. Maftei.  Convergence rate analysis of time discretization scheme for confined Lagrangian processes.  ⟨hal-01656716⟩, 2017.
  * [8] M. Bossy, R. Maftei, J.-P. Minier, and C. Profeta.  A stochastic approach to colloidal particle collision/agglomeration.  arXiv:1611.06248, 2016.
  * [9] N. Bou-Rabee and H. Owhadi.  Long-run accuracy of variational integrators in the stochastic context.  SIAM J. Numer. Anal., 48(1):278–297, 2010.
  * [10] L. Campana, M. Bossy, and C. Henry.  Lagrangian stochastic model for the orientation of inertialess spheroidal particles in turbulent flows: An efficient numerical method for CFD approach.  Computers and Fluids, 257:105870, 2023.
  * [11] G. Celeux, M. El Anbari, J.-M. Marin, and C. P. Robert.  Regularization in regression: Comparing Bayesian and frequentist methods in a poorly informative situation.  Bayesian Anal., 7(2):477–502, 06 2012.
  * [12] C. Cercignani.  The Boltzmann equation and its applications.  Springer, 1988.
  * [13] C. Costantini.  Diffusion approximation for a class of transport processes with physical reflection boundary conditions.  Ann. Prob., 19(3):1071–1101, 1991.
  * [14] R. L. Davidchack, T. E. Ouldridge, and M. V. Tretyakov.  New Langevin and gradient thermostats for rigid body dynamics.  J. Chem. Phys., 142(14):144114, 2015.
  * [15] M. Ebrahimian, R. S. Sanders, and S. Ghaemi.  Dynamics and wall collision of inertial particles in a solid–liquid turbulent channel flow.  J. Fluid Mech., 881:872–905, 2019.
  * [16] M. Freidlin.  Functional integration and partial differential equations.  Princeton Univ. Press, 1985.
  * [17] D. Gilbarg and N.S. Trudinger.  Elliptic partial differential equations of second order.  Springer, 1983.
  * [18] Y. A. Houndonougbo, B. B. Laird, and B. Leimkuhler.  A molecular dynamics algorithm for mixed hard-core/continuous potentials.  Molecular Physics, 98(5):309–316, 2000.
  * [19] V. A. Huynh, S. Karaman, and E. Frazzoli.  An incremental sampling-based algorithm for stochastic optimal control.  Intern. J. Robotics Res., 35(4):305–333, 2016.
  * [20] N. Ikeda and S. Watanabe.  Stochastic differential equations and diffusion processes.  North-Holland, 1989.
  * [21] V.E. Johnson and J.H. Albert.  Ordinal data modeling.  Springer, 2006.
  * [22] I. Karatzas and S.E. Shreve.  Brownian motion and stochastic calculus.  Springer, 1991.
  * [23] J.P. Klein and M.L. Moeschberger.  Survival analysis: techniques for censored and truncated data.  Springer, 2005.
  * [24] B. Leimkuhler and C. Matthews.  Rational construction of stochastic numerical methods for molecular sampling.  Appl. Math. Res. Express, 2013(1):34–56, 2013.
  * [25] B. Leimkuhler and C. Matthews.  Molecular dynamics with deterministic and stochastic numerical methods.  Springer, 2015.
  * [26] B. Leimkuhler, C. Matthews, and G. Stoltz.  The computation of averages from equilibrium and nonequilibrium Langevin molecular dynamics.  IMA J. Num. Anal., 36:13–79, 2016.
  * [27] B. Leimkuhler, A. Sharma, and M. V. Tretyakov.  Simplest random walk for approximating Robin boundary value problems and ergodic limits of reflected diffusions.  Ann. Appl. Prob., 33(3):1904–1960, 2023.
  * [28] B. Leimkuhler, T. J. Vlaar, T. Pouchon, and A. Storkey.  Better training using weight-constrained stochastic dynamics.  In Marina Meila and Tong Zhang, editors, Proceedings of the 38th International Conference on Machine Learning, volume 139, pages 6200–6211. PMLR, 2021.
  * [29] T. Lelièvre, M. Rousset, and G. Stoltz.  Free energy computations: A mathematical perspective.  Imperial College Press, 2010.
  * [30] A. Leroy, B. Leimkuhler, J. Latz, and D. Higham.  Adaptive stepsize algorithms for Langevin dynamics.  arXiv:2403.11993, 2024.
  * [31] P.L. Lions and A.S. Sznitman.  Stochastic differential equations with reflecting boundary conditions.  Comm. Pure Appl. Math., 37(4):511–537, 1984.
  * [32] R. Maftei.  Stochastic analysis for Lagrangian particles simulation : application to colloidal particle collision.  Theses, Université Côte d’Azur, December 2017.
  * [33] C. Marchioro, A. Pellegrinotti, E. Presutti, and M. Pulvirenti.  On the dynamics of particles in a bounded region: A measure theoretical approach.  J. Math. Phys., 17(5):647–652, 08 1976.
  * [34] J.C. Mattingly, A.M. Stuart, and D.J. Higham.  Ergodicity for SDEs and approximations: locally Lipschitz vector fields and degenerate noise.  Stoch. Proc. Applic., 101(2):185 – 232, 2002.
  * [35] G. N. Milstein and M. V. Tretyakov.  Numerical methods for Langevin type equations based on symplectic integrators.  IMA J. Numer. Anal., 23:593–626, 2003.
  * [36] G. N. Milstein and M. V. Tretyakov.  Numerical integration of stochastic differential equations with nonglobally Lipschitz coefficients.  SIAM J. Numer. Anal., 43(3):1139–1154, 2005.
  * [37] G.N. Milstein and M.V. Tretyakov.  Simulation of a space-time bounded diffusion.  Ann. Appl. Probab., 9(3):732–779, 1999.
  * [38] G.N. Milstein and M.V. Tretyakov.  Computing ergodic limits for Langevin equations.  Physica D, 229(1):81 – 95, 2007.
  * [39] G.N. Milstein and M.V. Tretyakov.  Stochastic numerics for mathematical physics.  Springer, 2nd edition, 2021.
  * [40] E. Nelson.  Dynamical theories of Brownian motion.  Princeton Univ. Press, 1967.
  * [41] F. Nier.  Boundary conditions and subelliptic estimates for geometric Kramers-Fokker-Planck operators on manifolds with boundaries.  Memoirs AMS, 252, 2013.
  * [42] J. W. Paisley, D. Blei, and M. I. Jordan.  Bayesian nonnegative matrix factorization with stochastic variational inference.  In Handbook of mixed membership models and their applications, 2014\.
  * [43] S. Sangalli, E. Erdil, A. Hoetker, O. F. Donati, and E. Konukoglu.  Constrained optimization to train neural networks on critical and under-represented classes.  In NeurIPS, 2021.
  * [44] K. Spiliopoulos.  A note on the Smoluchowski-Kramers approximation for the Langevin equation with reflection.  Stochastics and Dynamics, 07(02):141–152, 2007.
  * [45] D. Talay.  Stochastic Hamiltonian systems: exponential convergence to the invariant measure, and discretization by the implicit Euler scheme.  Markov Proc. Relat. Fields, 8:163–198, 2002.
  * [46] A.D. Wentzell.  A course in the theory of stochastic processes.  McGraw-Hill, 1981.
  * [47] T. Önskog.  Existence of pathwise unique Langevin processes on polytopes with perfect reflection at the boundary.  Stat. & Probab. Lett., 83(10):2211–2219, 2013.