Numerical integrators for confined Langevin dynamics
B. Leimkuhler A. Sharma M. V. Tretyakov


Understanding Privacy Risks of Embeddings Induced by Large Language Models
    Zhihao Zhu, Ninglu Shao, Defu Lian, Chenwang Wu, Zheng Liu, Yi Yang, Enhong Chen
    April 29, 2024
====================================================================================




 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.
     





§ 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 ⊂ℝ^d be a bounded domain with sufficiently smooth boundary ∂ G, n_G(q) be the outward normal unit vector at q ∈∂ G, and I_∂ G(q) be the indicator function of q ∈∂ G. Let σ be a positive constant, b: G̅×ℝ^d→ℝ^d satisfying certain conditions, and (·) denote scalar product. Consider the confined Langevin dynamics (CLD) governed by the  stochastic differential equations (SDEs):

    CLD
        Q(t) = Q(0) +    ∫_0^tP(s)ds,    
        P(t) = P(0) +∫_0^tb(P(s),Q(s))ds  + σ W(t) + ℛ(t),
    ℛ(t)    = -∑_0<s≤ t2(P(s^-)· n_G(Q(s)))n_G(Q(s))I_∂ G(Q(s)),     t ≥ 0,

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

The process Q(t) does not leave the set G̅ and it elastically reflects from the boundary ∂ G <cit.> which is achieved due to the cádlág process ℛ(t). At the collision of the position Q(t) with ∂ G,  ℛ(t) acts on the momentum P(t) to reverse the direction of its component orthogonal to ∂ G keeping the tangential component unchanged. 

If b(q,p) := -∇ U(q) - γ p  and  σ = √(2γ/β),   
i.e., consider the particular case of  confined Langevin dynamics (<ref>) defined in the set G̅×ℝ^d:

    ECLD
        Q(t)    = Q(0) + ∫_0^tP(s)ds, 
     
        P(t)     = P(0) -∫_0^t∇ U(Q(s))ds - γ∫_0^t P(s)ds + √(2γ/β)W(t) + ℛ(t),
    ℛ(t)    = -∑_0<s≤ t2(P(s^-)· n_G(Q(s)))n_G(Q(s))I_∂ G(Q(s)),

where β >0 is interpreted as the inverse temperature in molecular dynamics and γ >0 is a friction parameter.  Eq. (<ref>) is ergodic with the invariant density ρ(q,p) = 1/Ze^-β(U(q) + |p|^2/2), (q,p) ∈G̅×ℝ^d. 
The central theme of this paper is the efficient weak-sense numerical approximation of (<ref>)  including the computation of ergodic limits associated with (<ref>) 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 <cit.>, ordinal data models <cit.>, constrained Lasso and Ridge regression <cit.>, Latent Dirichlet Allocation <cit.>, non-negative matrix factorization <cit.>, and stochastic optimal control <cit.>. 
In <cit.> (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 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. <cit.>) which are also expected for the considered confined-space case. 

Molecular dynamics. To study molecular properties of a given substance under constant temperature 1/β, 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 <cit.> 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 - ∇ 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 <cit.>. Confined Langevin dynamics (<ref>) 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 <cit.>. 

Computational fluid dynamics. The SDEs (<ref>)  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 (<ref>) type allow to study, e.g., convection, polydispersity etc. related to mesoscopic particles suspended in turbulent flow (see e.g., <cit.> 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. <cit.>. One can exploit  (<ref>) to minimize a function U(q), where q satisfies constraints represented by domain G. Moreover, one can also consider replacing ∇ U with a noisy gradient ∇̃U to introduce batching of data for machine learning specific tasks. 

Solving specular boundary value problem. Via the Feynman-Kac formula <cit.>, the SDEs (<ref>) provide a probabilistic representation of the  specular boundary value problem for linear parabolic partial differential equations (PDEs) (see (<ref>)-(<ref>)). Solving such PDEs in high dimensions using deterministic methods is computationally infeasible. We can utilise the proposed here numerical methods for (<ref>) 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 (<ref>). Despite this,  very limited attention has been paid so far to numerical integration of confined Langevin dynamics. In this respect we note that in <cit.> the authors investigated a numerical scheme when G=ℝ^d-1×(0,∞) and proved its first order weak convergence. 

In this paper (Section <ref>), we derive first-order weak methods to numerically approximate (<ref>) as well as second-order schemes in the case of (<ref>). 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 <ref> 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 <ref> with necessary background material for (<ref>) and (<ref>) provided in Section <ref>. 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 (<ref>). In a deterministic setting, the discussed splitting approach would typically lead to first-order approximation, even in symmetric compositions <cit.>. Then stochastic numerics experience <cit.> suggests that finite-time convergence of a generalization of such deterministic methods to (<ref>) 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 <ref> we present results of numerical tests which confirm our theoretical predictions and compare the methods from  Section <ref>. In Section <ref> we summarize this work and suggest topics of interest for further research related to confined Langevin dynamics. 




§ PRELIMINARIES


Notation. 
For the boundary sets we write:
    Σ^+ = {  (q,p) ∈∂ G ×ℝ^d; (p· n_G(q)) >0} and Σ_t^+ = (0,t)×Σ^+. The space C^l_1,l_2,l_3([0,T]× A×ℝ^d) denotes the space of  functions whose continuous derivatives exist with respect to t ∈ [0,T], q ∈ A and p∈ℝ^d up to order l_1, l_2 and l_3, respectively; C^l_1,l_2(A×ℝ^d)  has  the same meaning but with respect to q and p;  C_b^l_1,l_2(A×ℝ^d)  is a space of functions in  C^l_1,l_2(A×ℝ^d) with bounded  derivatives up to order l_1 in q and l_2 in p; and  C_c^l_1,l_2(A×ℝ^d) are  functions in  C^l_1,l_2(A×ℝ^d) with compact support. 
As usual, |· | denotes the L^2-norm of a vector in ℝ^d.
For a differentiable vector-valued function  F(x) from G̅⊂ℝ^d or ℝ^d to ℝ^d, 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 <ref>, we need assumptions on the boundary ∂ G of G and the drift coefficient b(q,p).

We assume the boundary ∂ G of set G belongs to C^4.
  

The coefficient b(q,p) is a C^2,2(G̅×ℝ^d) function which grows at most linearly in p at infinity, i.e. there exists a constant K > 0 such that,

    |b(q,p)| ≤ K(1+|p|),

for all p ∈ℝ^d and q ∈G̅.


In <cit.>, existence of weak solutions and uniqueness in law of (<ref>) in G=(0,∞)×ℝ^d-1 with multiplicative noise and  b(q,p) = b(q) - γ p, γ > 0,  were established under boundedness assumptions on b(q) and  the diffusion coefficient σ(q). 
In <cit.>, existence of weak solution of (<ref>) in G=(0,∞)×ℝ^d-1  with bounded drift coefficient b was shown. 
This  well-posedness result was extended in <cit.> to the case of ∂ G ∈ C^3 being a compact submanifold of ℝ^d, by establishing the existence of a weak solution and its pathwise uniqueness for b(q,p) ≡ 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 <cit.>) as done in <cit.>.  Pathwise uniqueness of the solution of (<ref>) with b(q,p) ≡ 0 is proved in <cit.> when G is a convex polytope. 




 §.§ Backward Kolmogorov equation for (<ref>)


Consider the backward Kolmogorov equation

    ∂ u/∂ t  + σ^2/2Δ_pu + (p·∇_qu) + (b(q,p)·∇_pu)   = 0,     (t,q,p) ∈ (0,T)× G ×ℝ^d,

with terminal condition

    u(T,q,p) = φ(q,p),     (q,p) ∈ G×ℝ^d,

and specular boundary condition

    u(t,q,p) = u(t,q,p-2(p· n_G)n_G),     (t,q,p) ∈Σ_T^+,

where Δ_p is the Laplacian with respect to p, ∇_q and ∇_p are gradients with respect to q and p, respectively, and n_G:= n_G(q). 

The solution of (<ref>)-(<ref>) has the probabilistic representation  <cit.>:

    u(t_0,q_0,p_0) = 𝔼φ(Q_t_0,q_0,p_0(T),P_t_0,q_0,p_0(T)),

where (Q_t_0,q_0,p_0(t),P_t_0,q_0,p_0(t)), t ≥ t_0, solves (<ref>) with the starting point (t_0,q_0,p_0).

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

Take G =(0,∞) and ∂ G = { 0}, b(q,p) = (1/q+c + p), σ  =1, φ(q,p) = (q+c)^2e^-p^2, and boundary condition u(t,q,p) = u(t,q,-p) at Σ_T^+, where c ≥ 0. The solution is given by u(t,q,p) = (q+c)^2e^-p^2e^-(T-t).



Let G ⊂ℝ^d, U(q) ∈ C^1(G̅) and α, β  > 0. Take b(q,p) = -∇_qU(q) + α p, σ = √(2α/β), φ = e^-β(|p|^2/2 + U(q)) and boundary condition is u(t,q,p) = u(t,q,p-2(p· n_G)n_G), (t,q,p) ∈Σ_T^+. The solution is u(t,q,p)=exp(-β(|p|^2/2 + U(q))-α d(T-t)).




Let G ⊂ℝ^d, b(q,p) = -qcos|q|^2+2(d+2)p/2|p|^2 + 1, σ=√(2), and φ(q,p) = |p|^4 + |p|^2 + sin|q|^2. The boundary condition is u(t,q,p) = u(t,q,p-2(p· n_G)n_G),  (t,q,p) ∈Σ_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 φ and the solution of (<ref>)-(<ref>).

φ (q,p) ∈ C^4 ,4 (G̅×ℝ^d).


(i) The solution of (<ref>)-(<ref>) exists and  belongs to C^2,2,4((0,T)× G×ℝ^d) and C^0,1,0([0,T]×G̅×ℝ^d), and Hessian derivative of u(t,q,p) with respect to q exists on [0,T]×G̅×ℝ^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 (<ref>)-(<ref>) 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 <ref>. 
Also, it was shown in <cit.> that if φ∈ C_c^1,1(G ×ℝ^d) and b ∈ C_b^1,1(ℝ^d, ℝ^d), then u(t,q,p) ∈ C([0,T]×G̅×ℝ^d) and ∇_qu and ∇_pu exist and belong to C((0,T)×G̅×ℝ^d). 





 §.§ Ergodicity


Consider (<ref>) which is a particular case of (<ref>).


U(q) ∈ C^∞(G̅) and U(q) > 1 for all q ∈G̅.



Exploiting the approach used in proving Theorem 3.2 in <cit.> (see also <cit.>) one can establish under Assumptions <ref> and <ref> (note that ∇ U(q) is bounded thanks to G being bounded) that (<ref>) is exponentially ergodic, i.e.  
there exists a unique invariant measure μ(dq,dp) defined on G×ℝ^d and  the following inequality holds for any φ (q,p), (q,p) ∈G̅×^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:

    |𝔼(φ(Q(T),P(T))) - φ̅| ≤ Ke^-λ T,

where 
    φ̅ = ∫_ℝ^d∫_Gφ(q,p) μ(dq, dp).
  Moreover, the invariant measure μ has a probability density ρ_∞ with respect to Lebesgue measure. In <cit.>, geometric ergodicity is proved in the case of the semi-space G =(0,∞) ×ℝ^d-1.


The stationary Fokker-Plank equation for the density ρ_∞ has the form

    γ/βΔ_pρ_∞-(p·∇_qρ_∞) + (∇_qU(q)·∇_pρ_∞) + γ (∇_p·(pρ_∞))  = 0,     (q,p) ∈ G×ℝ^d,

with boundary condition

    ρ_∞(q,p) = ρ_∞
        (q,p-2(n_G· p)n_G),      (q,p) ∈Σ^+.

As one may check, the Gibbs density 

    ρ_∞(q,p) = 1/ℤexp(-β(|p|^2/2 + U(q)))

indeed satisfies the stationary Fokker-Plank equation (<ref>) and the boundary condition (<ref>). Here ℤ is the normalization constant so that ∫_ℝ^d∫_G̅ρ_∞(q,p)dq dp =1. 

In <cit.>, the stationary Fokker-Plank equation (<ref>)-(<ref>) in the case of G =(0,∞) ×ℝ^d-1 was considered.
In <cit.>, 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 (<ref>) introduced in Section <ref> equip us with the tools needed to approximately sample from a targeted density ρ(q), q ∈G̅, by choosing U(q) = -1/βln (ρ(q)), β >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 (<ref>)-(<ref>) related to (<ref>), i.e., of (<ref>)-(<ref>) with b(q,p) = -∇ U(q) - γ p.


The solution of (<ref>)-(<ref>) with b(q,p) = -∇ U(q) - γ p satisfies Assumption <ref> and the condition that 
there are positive constants C, λ and m independent of T such that the following bound holds for some m ≥ 1:

    ∑_l = 1^2∑_i + |j| = l    |D_t^iD_q^ju(t,q,p)| +  ∑_l = 1^4∑_|j| = l|D_p^ju(t,q,p)|≤ C(1 + |p|^2m)e^-λ (T-t),

where D_x^ju = ∂^|j|u/∂ x^j_1,…,∂ x^j_d, x = q ∈G̅ or  p ∈ℝ^d and j is a multi-index, and C>0 is independent of T and p.


We note that an estimate of the form (<ref>) in the case of Langevin equations in the whole space is proved in <cit.>.


Let us discuss the relationship between the confined Langevin dynamics (<ref>) and the reflected gradient SDE (Brownian dynamics with reflection or, in other words, the overdamped reflected Langevin dynamics) <cit.>:

    dX(t)    = -∇ U(X(t))dt + √(2/β) dW(t) + ν (X(t))I_∂ G(X(t)) dL(t),

 where L(t) is the  local time of the process X(t) on the boundary ∂ G.
 The solution X(t) of (<ref>) belongs to G̅ and it is instantaneously reflected along the internal normal ν to ∂ G when it hits the boundary ∂ G.
 Under the same conditions as we have imposed on the potential U and the boundary ∂ G above, the SDE (<ref>) is ergodic with the invariant measure 
 ρ_ RSDE(x) ∝exp(-β U(x)) <cit.>, which coincides with ρ_∞(q)=∫_ℝ^dρ_∞(q,p)dp, i.e., both (<ref>) and (<ref>) 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 ℝ^d space <cit.> (see also <cit.>), the solution Q(t) of (<ref>) tends in distribution to the solution X(t) of (<ref>) when γ→∞ <cit.>.
 In <cit.> we proposed and studied numerical algorithms for sampling from ρ_ RSDE(x), x ∈G̅, based on  (<ref>).
 Benefits of using Langevin equations vs  overdamped Langevin dynamics in the case of the whole space are well established (see e.g. <cit.>). Such a comparison in the case of confined spaces will be investigated in a future work. 






§ NUMERICAL METHODS


In Section <ref>, we present first-order weak methods for the SDEs (<ref>).
In Section <ref>, we propose splitting schemes for the SDEs (<ref>) which demonstrate (in numerical studies) second-order weak convergence.  In Section <ref>, we compare known numerical methods for (<ref>) with the schemes of Sections <ref> and <ref>, 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,…, N-1, and t_0=0. We take Q_0∈ G and P_0∈ℝ^d\{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. <cit.>).  
However, to construct numerical methods for (<ref>) and (<ref>), 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 ∂ G happen.

One step with possible collision.
Let the current state of a Markov chain approximating the solution of (<ref>) be (q, p) ∈G̅×ℝ^d. We construct the next state P_1, Q_1 of the chain as follows.

If q +ph ∈ G, 

    A_c(q,p;h)= (q + ph, p),

else (i.e., if q+ph ∈ G^c) A_c(q,p;h)= (Q_1,P_1), obtained from

    Q_c = q + pτ∈∂ G,
        P_1 = p - (n(Q_c)· p) n(Q_c),   
         Q_1 = Q_c + (h-τ)P_1.

The random collision time τ is found so that if q + ph ∉ G then τ ensures that the intermediate position Q_c∈∂ G (note that (p· n_G(Q_c)) > 0). Thereafter, the numerical method follows the step P_1 = p - (n(Q_c)· 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τ∈∂ G, we take the smallest positive value of τ.  Such a choice is possible due to Assumption <ref> (see <cit.>). For illustration, consider an annular domain in ℝ^2: G = { (q_1, q_2)  |  r_1^2 < q_1^2 + q_2^2 < r_2^2} and ∂ G = { (q_1, q_2)  |  q_1^2 + q_2^2 = r_1^2} ∪ { (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 +τ p ∈∂ G are τ_1 = (r_2 - r_1) /2, τ_2 = (3r_1 + r_2)/2, τ _3 = (r_1 + 3r_2)/2 and τ_4 = (r_1 - r_2)/2. Here, τ_4 is negative, so we take τ = min{τ_1, τ_2, τ_3} = τ_1 = (r_2 -r_1)/2. 



 §.§ First-order schemes


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

    δ(q,p  ; h,ξ) = hb(q,p) +  h^1/2σξ,

where ξ = (ξ^1,…,ξ^d)^⊺ is an ℝ^d-valued random variable such that ξ^i, i=1,…,d, are i.i.d. real-valued random variables with bounded all moments and satisfying 

    𝔼ξ^i = 𝔼(ξ^i)^3=0, 𝔼(ξ^i)^2 = 1.
    
The conditions (<ref>) are, for instance, satisfied if the {ξ^i} are drawn from the standard normal distribution or as i.i.d. discrete random variables with the distribution:

    P(ξ^i=± 1)=1/2.


In Fig. <ref> we present the proposed schemes. For the sake of clarity, we describe [PA_c] scheme in detail. As already mentioned, the aim is to construct a Markov chain (t_k,Q_k,P_k)_(k ≥ 0)  to approximate  (<ref>) 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 P̂_k+1 and  Q̂_k+1 respectively:

    P̂_k+1 = P_k + δ(Q_k, P_k  ; h, ξ_k+1),    Q̂_k+1 = Q_k + hP̂_k+1,

where ξ_k+1∼ξ and they are mutually independent. 
The next state of the Markov chain, Q_k+1 and P_k+1, is evaluated according to A_c(Q_k, P̂_k+1;h), which means that if Q̂_k+1∈ G, then P_k+1 = P̂_k+1 and Q_k+1 = Q̂_k+1, else P_k+1 = P̂_k+1 - (n_G(Q_c, k+1)·P̂_k+1)n_G(Q_c, k+1) and Q_k+1 = Q_c, k+1 + (h - τ)P_k+1, where Q_c, k+1 = Q_k + P̂_k+1τ∈∂ G.  



Schemes [PA_c] and [A_cP] have finite-time convergence with first weak order when they are applied to (<ref>). By way of illustration, we state below Theorem <ref> for scheme [PA_c] which is proved in Section <ref>. When these schemes are applied to the ergodic SDE (<ref>), they approximate the ergodic limit with error 𝒪 (h+exp(-λ T)).  We state below the corresponding Theorem <ref> for scheme [PA_c] which is proved in Section <ref>.


Under Assumptions <ref>-<ref>, the order of weak convergence of scheme [PA_c]  for (<ref>) is 𝒪(h), i.e. the following inequality holds, for sufficiently small h>0:

    |𝔼φ(Q_N,P_N) - 𝔼φ(Q(T),P(T))| ≤ Ch,

where 

(Q(t),P(t)) is from (<ref>) and C>0 is independent of h.



Let Assumptions <ref>, <ref>, <ref>-<ref> hold. For sufficiently small h>0 the following error estimate holds for scheme [PA_c] applied to (<ref>):

    |𝔼φ(Q_N,P_N) - φ̅| ≤ C(h+exp(-λ T)),

where C>0 and λ >0 are constants independent of T and h. 







 §.§ [A_c, 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. <cit.> and references therein). Hence, it is natural to extend those methods to the confined Langevin dynamics (<ref>). The [A_c, B, O] splitting schemes proposed in this subsection degenerate to the standard deterministic Störmer-Verlet schemes with collision step A_c when (<ref>) 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 A_c having local error of order O(h) <cit.>.
Because of this, stochastic numerics experience <cit.> suggests that finite-time weak convergence of the [A_c, B, O]  schemes for (<ref>) should be first order, too. However, the [A_c, B, O]  schemes for (<ref>) unexpectedly turned out to be of weak order 2, which we demonstrate in several carefully constructed and performed numerical experiments (Section <ref>) and theoretically justify when G is a half-space, i.e. G=(0,∞)×ℝ^d-1, in Section <ref>. 

Let ξ = (ξ^1,…,ξ^d)^⊺ be an ℝ^d valued random variable  such that ξ^i, i=1,…,d, are i.i.d. real-valued random variables with bounded all moments and satisfying 

    𝔼ξ^i = 𝔼(ξ^i)^3=𝔼(ξ^i)^5=0,   𝔼(ξ^i)^2 = 1,   𝔼(ξ^i)^4 = 3.
    
The conditions (<ref>) are, for instance, satisfied for ξ^i having the standard normal distribution or being i.i.d. discrete random variables with the distribution:

    P(ξ^i=±√(3))=1/6,    P(ξ^i=0)=2/3.

In Fig. <ref>, ξ_k∼ξ and ζ_k∼ξ  and they all are mutually independent. 
Introduce the notation 

    O(p;h,ξ,γ):=p e^- γ h +√(σ ^2/2γ(1-e^-2γ h))ξ .

The three proposed methods are presented in Fig. <ref>. It is straightforward to write down the other three possible symmetric concatenations, namely [BOA_cOB], [A_cBOBA_c], [A_cOBOA_c], which we test in Experiment <ref>.








 §.§ Comparison with other numerical integrators for (<ref>)


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

Let G = (0, ∞). 
In <cit.> (see also <cit.>), the authors studied a symmetrized scheme in which, given (q,p) at time t_0, the approximate position is always updated as 

    Q̅(t_0+h) = |q + ph|.

At every step,  θ = t_0 - q/p is calculated. A collision happens if θ∈ (t_0, t_0+h). In that case  

    Q̅(θ) = 0,    P̅(θ -) = p + b(q,p) (θ - t_0) + σ (W(θ) - W(t_0)),    P̅(θ) = - P̅(θ -), 
    P̅(t+h) = p + b(Q̅(θ),P̅(θ)) (t_0+h - θ) + σ (W(t_0+h) - W(θ)) .

If no collision occurs then P̅(t_0+h) = p + b(q,p) h + σ (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,∞)×ℝ^d-1, in <cit.>. 

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, θ) due to a realized large Wiener increment. In such a case the update P̅(θ) = - P̅(θ-) 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 [A_cP], [PA_c] and the splitting schemes of Section <ref> if implemented for G = (0,∞)×ℝ^d-1. We also observe that if one extends the symmetrized reflection of (<ref>) around the boundary in curved domains, then the symmetrized reflection around the boundary along the normal does not imitate specular reflection of  (<ref>), unlike the schemes of Sections <ref> and <ref>. 

 Further, some penalized schemes for (<ref>) are also proposed in <cit.> without analysis. It is known (see <cit.> and references therein) that symmetrized reflection and penalized schemes are good candidates to  approximate the reflected SDEs (<ref>),  in which the local time L(t) is responsible for keeping the trajectories X(t) in G̅. In contrast, the confinement of position Q(t) in G̅ is achieved in (<ref>) by the cádlág  process ℛ(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  ℛ(t), which, in particular, ensures  conservation of (conjugate) momenta in the event of collision, i.e.,  

    |P_k|^2 = |P̂_k|^2I_G̅(Q̂_k) + |P̂_k - 2(n^π·P̂_k)n^π|^2I_G̅^c(Q̂_k)= |P̂_k|^2,

where n^π = 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 <ref> and <ref>).  

Although we present our results in a bounded domain G with boundary ∂ G, the results can be extended to the half-space G = (0,∞)×ℝ^d-1 or any other unbounded domain with boundary satisfying Assumption <ref>. 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 ∇ U or replaced by using the concept of rejecting trajectories <cit.> 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) · p) (see A_c 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 A_c will not lie in G. We note that when ∂ 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 ∂ 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 A_c 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 [PA_c] scheme: depending on the update of the momentum P,  A_c 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 (q̃, p̃) the position and new velocity, respectively, after the P step in [PA_c].  Also, denote the continuous-time dynamics evolving according to the A_c  step by (Q̃(t),  P̃(t)). As already highlighted, the dynamics of Q̃(t) follows straight paths, and indeed is the exact solution of the following ODE system:
 
    dQ̃(t)    = P̃(t) dt,     Q̃(t_0) = q̃,
    
         d P̃(t)    = d ℛ(t),      P̃(t_0) = 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 <cit.>, that the Lebesgue measure of a set of configurations (q̃, 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 <ref>).     

 We observe that all the schemes of Sections <ref> and <ref> 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 ∇ U to deliberately push the dynamics towards the boundary to enhance the number of collisions (see Experiment <ref>). 




§ PROOFS


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




 §.§ Proof of Theorem <ref>


The proof is based on  two main ingredients: a one-step error estimate (Lemma <ref>) and an estimate related to the average number of collisions (Lemma <ref>). 
According to the lemmas, inside the domain, where we use the usual Euler-type walker, the one-step error is 𝒪(h^2) and the walker makes 𝒪(1/h) steps; the error in the collision steps is 𝒪(h) but the average number of collision steps is bounded and independent of the time step h. The finite-time error is  naturally 𝒪(h). 
We note that the position Q(t) in the continuous dynamics (<ref>) and its discrete approximation Q_k is bounded as we consider (<ref>) in a bounded domain G̅. At the same time, momentum P(t) lives in ℝ^d and hence, to prove convergence, we need that moments of its discrete approximation P_k are bounded (see Lemma <ref>). The proof relies on the fact that  the drift in  (<ref>) is growing at most linearly in P(t) and the noise is additive. 


Let Assumptions <ref>-<ref> hold, then the moments of the momentum P_k constructed according to the scheme [PA_c] applied to (<ref>) are bounded, i.e. for any m ≥ 1:

    𝔼(|P_k|^2m) ≤ C(1 + |P_0|^2m),

where  C>0 is a constant independent of h.


Using (<ref>), we get

    𝔼|P_k+1|^2m=𝔼|P̂_k+1|^2m= 𝔼|(P̂_k+1-P_k)+P_k|^2m
       ≤𝔼|P_k|^2m+m𝔼|P_k|^2m-2[2(P_k ·(P̂_k+1-P_k))+(2m-1)|P̂_k+1-P_k|^2] 
       + K ∑_l=3^2m𝔼 |P_k|^2m-l  |P̂_k+1-P_k|^l,

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

    |𝔼[(P̂_k+1-P_k)|P_k]| ≤ Kh|P_k|,    𝔼[|P̂_k+1-P_k|^l|P_k] ≤ Kh^l/2 (1+|P_k|^l),   l ≥ 2,

with K>0 dependent on m, b and the moments of ξ_k. 
We observe that

    𝔼|P_k+1|^2m≤𝔼|P_k|^2m+Kh 𝔼|P_k|^2m +Kh

and, using the Gronwall lemma, we get (<ref>).


Next we estimate the one-step weak error.  


Under Assumptions <ref>-<ref>, the one-step error of the scheme [PA_c] is estimated as

    |𝔼[u(t_k+1,Q_k+1,P_k+1) - u(t_k,Q_k,P_k)| Q_k,P_k]| 
       ≤ C(h^2(1+|P_k|^2m) + h𝔼[(n_G(Q_c,k+1)·P̂_k+1)(1+|P_k+1|^2m)I_G^c(Q̂_k+1)|Q_k,P_k]) ,

where u(t,q,p) is the solution of (<ref>)-(<ref>) and C>0 is a constant independent of h and P_k.


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 (<ref>) that 

    𝔼[|P_k+1|^2m| Q_k,P_k]≤ C(h+|P_k|^2m),

where C>0 is independent of h, P_k, Q_k.

We have

    u(   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,P̂_k+1))I_G(Q̂_k+1) 
    
        +(u(t_k+1,Q_k+1,P_k+1) - u(t_k+1,Q_k,P̂_k+1))I_G^c(Q̂_k+1) +
    (u(t_k+1,Q_k,P̂_k+1) - u(t_k,Q_k,P_k)) 
       =: A_1 + A_2 + A_3.

Consider the first term. Using Taylor's expansion and Assumption <ref>, we obtain 

    A_1   := (u(t_k+1,Q_k+1,P_k+1) - u(t_k+1,Q_k,P̂_k+1))I_G(Q̂_k+1) 
       
        = (u(t_k+1,Q_k+1,P̂_k+1) - u(t_k+1,Q_k,P̂_k+1))I_G(Q̂_k+1)
       
       = (h(P̂_k+1·∇_qu(t_k+1,Q_k,P̂_k+1)) + r_1)I_G(Q̂_k+1),

where the remainder r_1 satisfies 

    |𝔼[r_1 I_G(Q̂_k+1)| Q_k,P_k]| ≤ 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 ∇_qu(t_k+1, Q_c,k+1, P_k+1) around Q_k,  we get for the second term in (<ref>):

    A_2:= (u(t_k+1,Q_k+1,P_k+1) - u(t_k+1,Q_k,P̂_k+1))I_G^c(Q̂_k+1)  
        = (u(t_k+1,Q_c,k+1,P_k+1) + (h - τ)(P_k+1·∇_qu(t_k+1,Q_c, k+1, P_k+1)) 
        - u(t_k+1,Q_k,P̂_k+1) +r_2)I_G^c(Q̂_k+1)
       
       =  (u(t_k+1,Q_c,k+1,P_k+1)  + (h - τ)(P_k+1·∇_qu(t_k+1,Q_k, P_k+1)) 
          - u(t_k+1,Q_k,P̂_k+1)+r_3)I_G^c(Q̂_k+1)
       
       =  (u(t_k+1,Q_c,k+1,P̂_k+1)  + (h - τ)(P̂_k+1·∇_qu(t_k+1,Q_k, P_k+1)) 
           -2(h-τ)(n_G(Q_c,k+1)  ·P̂_k+1) (n_G(Q_c,k+1) ·∇_qu(t_k+1,Q_k, P_k+1)) 
          - u(t_k+1,Q_k,P̂_k+1)+r_3)I_G^c(Q̂_k+1),

 where we have used P_k+1 = P̂_k+1 -2(n_G(Q_c, k+1)·P̂_k+1)n_G(Q_c, k+1)I_G^c(Q̂_k+1) and  u(t_k+1,Q_c, k+1,P_k+1) = u(t_k+1, Q_c, k+1, P̂_k+1) from (<ref>). Here r_i, i=2,3, satisfy 
 
    |𝔼[r_i I_G^c(Q̂_k+1)| Q_k,P_k]| ≤ Ch^2 (1 + |P_k|^2m)

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,P̂_k+1) around Q_k and ∇_qu(t_k+1, Q_k, P_k+1) around P̂_k+1, we obtain
 
    A_2    =   (τ (P̂_k+1·∇_qu(t_k+1, Q_k, P̂_k+1)) +   (h - τ)(P̂_k+1·∇_qu(t_k+1,Q_k, P_k+1)) 
           -2 (h-τ)(n_G(Q_c,k+1)  ·P̂_k+1) (n_G(Q_c,k+1) ·∇_qu(t_k+1,Q_k, P_k+1)) + r_4)I_G^c(Q̂_k+1)
        = (τ (P̂_k+1·∇_qu(t_k+1, Q_k, P̂_k+1)) +(h - τ)(P̂_k+1·∇_qu(t_k+1,Q_k, P̂_k+1)) 
           -2 (h-τ)(P̂_k+1·J_∇_qu(P̂_k, θ)n_G(Q_c,k+1))(P̂_k+1· n_G(Q_c,k+1)) 
           -2 (h-τ)(n_G(Q_c,k+1)  ·P̂_k+1) (n_G(Q_c,k+1) ·∇_qu(t_k+1,Q_k, P_k+1))  
          + r_4 )I_G^c(Q̂_k+1)
         
        =  (h (P̂_k+1·∇_qu(t_k+1, Q_k, P̂_k+1)) 
         -2 (h-τ)(P̂_k+1·J_∇_qu(P̂_k,θ)n_G(Q_c,k+1))(P̂_k+1· n_G(Q_c,k+1))
         
           -2 (h-τ)(n_G(Q_c,k+1)  ·P̂_k+1) (n_G(Q_c,k+1) ·∇_qu(t_k+1,Q_k, P_k+1)) 
         + r_4)I_G^c(Q̂_k+1),

 where J_∇_qu(P̂_k,θ) is the Jacobian of ∇_qu with respect to p evaluated at (t_k+1, Q_k, P̂_k,θ) with 
 
    P̂_k,θ = P̂_k+1 -2θ (n_G(Q_c,k+1)·P̂_k+1)n_G(Q_c,k+1),   θ∈ (0,1),
 
 and r_4 satisfies (<ref>). 
 
 Combining the obtained expressions for A_1 and A_2, we ascertain that
 
    A_1 + A_2   =  h (P̂_k+1·∇_qu(t_k+1, Q_k, P̂_k+1))  
           -2 ((h-τ)(P̂_k+1·J_∇_qu(P̂_k,θ)n_G(Q_c,k+1))(P̂_k+1· n_G(Q_c,k+1)) 
           + (h-τ)(n_G(Q_c,k+1)  ·P̂_k+1) (n_G(Q_c,k+1) ·∇_qu(t_k+1,Q_k, P_k+1)) 
    )I_G^c(Q̂_k+1) + r_5

 with r_5 satisfying 
  
    |𝔼[r_5 | Q_k,P_k]| ≤ Ch^2 (1 + |P_k|^2m)
 
 with C>0 being a constant independent of h, P_k, Q_k. 
Taking the expectation conditioned on (Q_k, P_k), we get

    𝔼(A_1 + A_2| Q_k,P_k) =  h 𝔼[(P̂_k+1·∇_qu(t_k+1, Q_k, P̂_k+1))|Q_k,P_k] 
          
           -2𝔼[( (h-τ)(P̂_k+1·J_∇_qu(P̂_k,θ)n_G(Q_c,k+1))(P̂_k+1· n_G(Q_c,k+1)) 
           + (h-τ)(n_G(Q_c,k+1)  ·P̂_k+1) (n_G(Q_c,k+1) ·∇_qu(t_k+1,Q_k, P_k+1)))I_G^c(Q̂_k+1)|Q_k, P_k] 
           + 𝔼[ r_5| Q_k,P_k]
        =  h 𝔼[(P_k·∇_qu(t_k, Q_k, P_k))|Q_k,P_k]  + r_6
         -2 𝔼[((h-τ)(P̂_k+1·J_∇_qu(P̂_k,θ)n_G(Q_c,k+1))(P̂_k+1· n_G(Q_c,k+1)) 
         + (h-τ)(n_G(Q_c,k+1)  ·P̂_k+1) (n_G(Q_c,k+1) ·∇_qu(t_k+1,Q_k, P_k+1)) )
         I_G^c(Q̂_k+1)|Q_k, P_k],

where we have expanded ∇_qu(t_k+1, Q_k, P̂_k+1) around P_k and then ∇_q u(t_k+1, Q_k, P_k) around t_k,
and r_6 satisfies  
  
    |r_6| ≤ Ch^2 (1 + |P_k|^2m)
 
 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. <cit.>), we can write

    𝔼(A_3| Q_k,P_k)     = h∂/∂ tu(t_k,Q_k, P_k) +  h(b_k·∇_pu(t_k,Q_k,P_k)) + hσ^2/2Δ_pu(t_k,Q_k,P_k)  +r_7,

where r_7 satisfies (<ref>).

Note that (see Assumption <ref>, (<ref>) and (<ref>)):

    |n_G(Q_c, k+1)·∇_qu(t_k+1, Q_k, P_k+1)|    ≤ C (1 + |P_k+1|^2m), 
    
    |P̂_k,θ| ≤ C|P̂_k+1| = C|P_k+1|, 
     
    |(P̂_k+1·J_∇_q u(P̂_k,θ)n_G(Q_c, k+1))| ≤ C(1 + |P_k+1|^2m).


Thus, using (<ref>), (<ref>) and (<ref>), we infer

    |𝔼(u(t_k+1, Q_k+1, P_k+1) - u(t_k,Q_k,P_k)| Q_k, P_k)|
    
        
           ≤ C(h^2(1+|P_k|^2m) + h𝔼[(n_G(Q_c,k+1)·P̂_k+1)(1+|P_k+1|^2m)I_G^c(Q̂_k+1)|Q_k,P_k]),

where C>0 is independent of h and P_k. 

 

The estimate in the above local error lemma contains the term of order 𝒪(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 Σ_+, 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):

    𝕋V(t,q,p) - V(t,q,p)    = -g(t,q,p),     (t,q,p) ∈ [T_0, T-h]×G̅×ℝ^d, 
    
        V(t,q,p)    = 0,          (t,q,p) ∈{T}×G̅×ℝ^d,

where 𝕋 is the one-step transition operator:

    𝕋V(t,q,p) = 𝔼V(t_1, Q_1, P_1),    (t_0,Q_0,P_0) = (t,q,p),

and g(t,q,p) ≥ 0. The solution to this problem is given by <cit.>:

    V(t_0,q,p) =  𝔼(∑_k=0^N-1g(t_k,Q_k,P_k)| Q_0 = q,P_0 = p).


 
Let  Assumptions <ref>-<ref> hold, then for any integer m ≥ 1:

    𝔼(∑_k=0^N-1(n_G(Q_c, k+1)·P̂_k+1)(1+|P_k+1|^2m)I_G^c(Q̂_k+1) )≤ C,

where Q̂_k+1 = Q_k + P̂_k+1h, k=0, …, N-1, are from [PA_c] scheme and C>0 is a constant independent of h.





 Let us take g(t_k,q,p) in (<ref>) as 

    g(t_k,q,p)= 𝔼[(n_G(Q_c, k+1)·P̂_k+1)(1+|P_k+1|^2m)I_G^c(Q̂_k+1)| Q_k=q,P_k=p]

then the solution of (<ref>)-(<ref>) is given by

    v   (t_0,q,p) =  
       𝔼(∑_k=0^N-1 g(t_k,Q_k, P_k) | Q_0 = q , P_0 = p)
       
        =𝔼(∑_k=0^N-1𝔼((n_G(Q_c,k+1)·P̂_k+1)(1+|P_k+1|^2m)I_G^c(Q̂_k+1)| Q_k,P_k) | Q_0 = q, P_0 = p )  
       
       =𝔼(∑_k=0^N-1(n_G(Q_c, k+1)·P̂_k+1)(1 + |P_k+1|^2m)I_G^c(Q̂_k+1)| Q_0 = q, P_0 = p ).


If we can find a function V(t,q,p) with g(t,q,p), satisfying (<ref>)-(<ref>), such that 

    g(t,q,p) ≥𝔼[(n_G(Q_c)·P̂_1)(1+|P_1|^2m)I_G^c(𝒬̂_1)],

then the inequality, v(t_0,Q_0,P_0) ≤ V(t_0,Q_0,P_0) will hold. Therefore, the aim is to construct 
a function g = -(𝕋V - V) satisfying (<ref>). 

Consider a boundary zone G_R_0⊂ G defined as G_R_0 := { q ∈G̅  ;  (q,∂ G) < R_0}. Due to <cit.>,  there is an R_0 such that the projection of q on ∂ G is unique whenever q ∈ G_R_0. Introduce the function

    V(t,q,p) =  
        0,    (t, q, p) ∈{T}×G̅×ℝ^d,
    
        (A + e^L(T-t))(1+| p|^2m +2)  +  (B(q)· p)(1+|p|^2m),     (t, q, p) ∈ [0, T-h]×G̅×ℝ^d,

where B(q) ∈ C^3(G̅) is  an ℝ^d valued function such that  B(q) : = n_G(q^π) when q ∈G̅_R_0; here q^π is the projection of q on ∂ G. Note that n_G(q) ∈ C^4(G̅) due to Assumption <ref> (see <cit.>). We extend the function B(q) from G_R_0  to G̅ so that B(q) ∈ C^3(G̅) and |B(q)|_G^2,α≤ C |n(q^π)|_G_R_0^2,α, where |·|^2,α is the Hölder norm and C>0 depends on G. Such an extension is possible due to <cit.>. 
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 <ref> implies for any ℓ≥ 1
 
    |P_1|^2ℓ = |P̂_1|^2ℓ  
         = |p|^2ℓ  +2ℓσ h^1/2(ξ· p)|p|^2ℓ-2 +  r_1,

 and
 
    𝔼[|P_1|^2ℓ |q,p] = |p|^2ℓ +  r_2,

 where r_1 and r_2 satisfy for any ℓ≥ 1:
 
    (𝔼(|r_1|^2ℓ))^1/(2ℓ)≤ Ch(1 + |p|^2ℓ)

 and
 
    |r_2| ≤ Ch(1+|p|^2ℓ)

 with C>0 being independent of h and p.
 
Given (t_k,Q_k, P_k) = (t,q,p), we now aim at calculating -(𝕋V(t,q,p) - V(t,q,p)). 
Using Taylor's theorem, we obtain 

    B(Q_1)    = (B(q)  + hJ_B(q +ϵ_1 hP̂_1)P̂_1)I_G(Q̂_1)
         
           +  (B(Q_c) + (h-τ)J_B(Q_c+ ϵ_2(h- τ)P_1)P_1)I_G^c(Q̂_1)

 and
 
    B(Q_c)    = B(q) +  τJ_B(q +ϵ_3τP̂_1)P̂_1,
 
where ϵ_i∈ (0,1), i = 1,2, 3, and J_B(x) denotes the Jacobian matrix of B evaluated at x. 
Therefore, we get

    (B(Q_1)   · P_1) = (B(Q_1)·P̂_1) I_G(Q̂_1)  +  (B(Q_1) · P_1)I_G^c(Q̂_1) 
         = ((B(q)·P̂_1) + h(J_B(q +ϵ_1 hP̂_1)P̂_1·P̂_1) )I_G(Q̂_1) 
            +   ((B(Q_c)·  P_1) + (h-τ)(J_B(Q_c+ ϵ_2(h- τ)P_1)P_1· P_1))I_G^c(Q̂_1) )
        
         = ((B(q)·P̂_1) + h(J_B(q +ϵ_1 hP̂_1)P̂_1·P̂_1) )I_G(Q̂_1) 
            +   ((B(Q_c)·P̂_1) + (h-τ)(J_B(Q_c+ ϵ_2(h- τ)P_1)P_1· P_1))I_G^c(Q̂_1) )
           - 2 (P̂_1· n_G(Q_c))(B(Q_c)· n_G(Q_c))I_G^c(Q̂_1).

Since (B(Q_c)· n_G(Q_c)) = 1 and using (<ref>), we obtain

    (B(Q_1)   · P_1)   = ((B(q)·P̂_1) + h(J_B(q +ϵ_1 hP̂_1)P̂_1·P̂_1) )I_G(Q̂_1) 
            +   ((B(Q_c)·P̂_1) + (h-τ)(J_B(𝒬_c+ ϵ_2(h- τ)P_1)P_1· P_1))I_G^c(Q̂_1) )
           - 2 (P̂_1· n_G(Q_c))I_G^c(Q̂_1) 
       
    = (B(q)·P̂_1) + h(J_B(q +ϵ_1 hP̂_1)P̂_1·P̂_1) )I_G(Q̂_1)  
            +   ( τ (J_B(q +ϵ_3τP̂_1)P̂_1·P̂_1) + (h-τ)(J_B(𝒬_c+ ϵ_2(h- τ)P_1)P_1· P_1))I_G^c(Q̂_1) 
           - 2 (P̂_1· n_G(Q_c))I_G^c(Q̂_1).

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

    (J_B(q +ϵ_1 hP̂_1)P̂_1·P̂_1)   ≤ C|P̂_1|^2, 
    
        (J_B(𝒬_c+ ϵ_2(h- τ)P_1)P_1· P_1)    ≤  C|P_1|^2 =C|P̂_1|^2, 
    
        (J_B(q +ϵ_3τP̂_1)P̂_1·P̂_1)    ≤ C|P̂_1|^2,

where C>0 is a generic constant independent of P̂_1. Using these estimates, (<ref>) becomes 

    (B(Q_1)   · P_1) = (B(q)·P̂_1) - 2 (P̂_1· n_G(Q_c))I_G^c(Q̂_1) + r_3,

where |r_3| ≤ Ch|P̂_1|^2 with C >0 being independent of h and P̂_1. 
Using (<ref>) and linear growth of b in p, we obtain

    𝔼[(B(Q_1)· P_1)(1 + |P_1|^2m)] = 𝔼[ ((B(q) ·P̂_1) - 2 (P̂_1· n_G(Q_c))I_G^c(Q̂_1)   + r_3)(1 + |P_1|^2m)] 
         = 𝔼[(B(q) ·P̂_1) (1 + |P_1|^2m)] - 2𝔼[ (P̂_1· n_G(Q_c))(1 + |P_1|^2m)I_G^c(Q̂_1) ]  + r_4,

where  |r_4| ≤ Ch[1 + |p|^2m+2] with C >0 being independent of h and p. Using (<ref>) and (<ref>), we get

    𝔼[(B(q)    ·P̂_1) (1 + |P_1|^2m)] = (B(q)· p)𝔼[(1 + |P_1|^2m)] + h(B(q)· b)𝔼[(1 + |P_1|^2m)] 
           + h^1/2𝔼[(B(q)·ξ)(1 + |P_1|^2m)]
        = (B(q)· p)(1 + |p|^2m) + r_5,

where |r_5|≤ Ch(1 + |p|^2m+1) with C>0 being independent of h and p. From (<ref>) and (<ref>), we have

    𝔼[(B(Q_1)· P_1)(1    + |P_1|^2m)]
        = (B(q)· p)(1 + |p|^2m)- 2𝔼[ (P̂_1· n_G(Q_c))(1 + |P_1|^2m)I_G^c(Q̂_1) ]  + r_6,

where |r_6|≤ Ch(1 + |p|^2m+2) with C>0 being independent of h and p.

Using (<ref>) with ℓ = m+1, we have

    (A + e^L(T-t- h))(1 + 𝔼[|P_1|^2m +2])  =  (A + e^L(T-t)(1 - Lh + L^2r_7)(1 + |p|^2m +2 +  r_2) 
        =  (A + e^L(T-t))(1 + |p|^2m +2 + r_2) 
        -L h e^L(T-t)(1 + |p|^2m +2) +r_8 e^L(T-t),

where r_7 ≤ Ch^2 and r_8 ≤ Ch^2(1+L^2)(1 + |p|^2m +2) with C>0 independent of h, L and p.
Recall that r_2≤ Ch(1+|p|^2m +2), where C>0 is independent of h, L and p.
 Therefore, using (<ref>) and (<ref>), we get

    𝔼[   V(t_1,Q_1, P_1)] - V(t,q,p) =  (A + e^L(T-t- h))(1 + 𝔼[|P_1|^2m +2]) 
           + 𝔼[(B(Q_1)· P_1)(1 + |P_1|^2m)]  - (A + e^L(T-t))(1 + |p|^2m + 2) 
    - (B(q)· p)(1 + |p|^2m) 
        = (A + e^L(T-t))r_2-L h e^L(T-t)(1 + |p|^2m +2) +r_8 e^L(T-t)+r_6
            - 2 𝔼[(P̂_1· n_G(Q_c))( 1+ |P_1|^2m) I_G^c(Q̂_1)| q, p]

Hence, there are sufficiently large L>0 and sufficiently small h>0 so that g(t,q,p) := - (𝕋V-V) ≥𝔼[(P̂_1· n_G(Q_c))(1+|P_1|^2m)I_G^c(Q̂_1)|q,p]. This proves that (<ref>) holds with a constant C>0  independent of h (it depends on P_0 and T). 



 Using Lemma <ref>, we have 

    |𝔼(φ(Q_N,P_N)) - 𝔼(φ(Q(T),P(T)))| = |𝔼(∑_k=0^N-1𝔼(u(t_k+1,Q_k+1,P_k+1) -u(t_k,Q_k,P_k) | Q_k,P_k))| 
    ≤ Ch^2𝔼(∑_k=0^N-1|P_k|^2m) + Ch𝔼(∑_k=0^N-1(n_G(Q_c,k+1)·P̂_k+1)I_G^c(Q̂_k+1) (1+|P_k+1|^2m)),

which on applying Lemma <ref> and Lemma <ref> gives the required result.




 §.§ Proof of Theorem <ref>


The proof of Theorem <ref> has the same ingredients as the proof of Theorem <ref> in the previous section with the major difference that the global error bound in Theorem <ref> should not depend on the integration time T since here we are interested in ergodic limits associated with (<ref>), while the error bound in Theorem <ref> for the approximation of (<ref>) on the finite time interval [0,T] can depend on T. We start (Lemma <ref>) by estimating moments of the momentum P(t) uniformly in T. Next, in Lemma <ref>,  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 <ref>. 


Let Assumptions <ref> and <ref> hold and γ h < 1, then the moments of the momentum P_k constructed according to the scheme [PA_c] applied to (<ref>) are bounded, i.e. for any m ≥ 1:

    𝔼(|P_k|^2m) ≤ C,

where  C>0 is independent of h and T.

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

    𝔼|P_k+1|^2m=𝔼|P̂_k+1|^2m≤𝔼
    |P_k|^2m+m𝔼|P_k|^2m-2[2(P_k· (P̂
    _k+1-P_k))+(2m-1)|P̂_k+1-P_k|^2] 
        +K∑_l=3^2m𝔼|P_k|^2m-l|P̂_k+1-P_k|^l,

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

    𝔼[P_k· (P̂_k+1-P_k)|P_k]    ≤   -γ
    h|P_k|^2+Kh|P_k|,   𝔼[|P̂_k+1-P_k|^2|P_k] ≤ Kh^2|P_k|^2+
    2γ/βh, 
    𝔼[|P̂_k+1-P_k|^l|P_k]    ≤   Kh^l|P_k|^l+Kh^l/2, l≥ 3,

with K>0 independent of h (it depends on max_q∈G̅|∇
U(q)|, γ , β , l and moments of ξ ). Using the
generalized Young's inequality, i.e. ab  ≤ϵ^ia^i/i + ϵ^-jb^j/j; a, b ≥ 0, ϵ>0; i,j >1; 1/i + 1/j =1, with a = |P_k|^2m-1, b = K, i = 2m/(2m -1), j = 2m and ϵ = γ^2m-1/2m, we get  

    2m𝔼[ |P_k|^2m-2(P_k· (P̂_k+1-P_k))] 
    ≤ -2mγ h 𝔼|P_k| ^2m + 2mh K𝔼(|P_k|^2m-1) 
          ≤ -2mγ h 𝔼|P_k| ^2m + (2m -1)γ h 𝔼|P_k|^2m + h K^2m(1/γ)^2m -1
          ≤ -γ h𝔼|P_k|^2m  + h K^2m(1/γ)^2m -1.

In the similar manner, we have

    (2m-1)𝔼|P_k|^2m-2[|P̂_k+1-P_k|^2] ≤ Kh^2
    𝔼|P_k|^2m+γ h/2𝔼|P_k|^2m + Kh,

and similarly with some new K>0

    K𝔼|P_k|^2m-l|P̂_k+1-P_k|^l≤γ h^l/2/
    2m𝔼|P_k|^2m+Kh^l/2, l≥ 3.

Next note that

    𝔼|P_k+1|^2m≤𝔼|P_k|^2m-γ h/2(1-
    √(h))𝔼|P_k|^2m+Kh

and, using the Gronwall lemma, we get 

    𝔼|P_k|^2m≤( 1-γ h/2(1-√(h)))
    ^k|P_0|^2m+2K/γ (1-√(h))[ 1-( 1-γ h
    /2(1-√(h))) ^k] .






Under Assumptions <ref>, <ref>, <ref>-<ref>, the one-step error of numerical scheme [PA_c] applied to (<ref>) is estimated as 

    |𝔼[u(t_k+1,Q_k+1,P_k+1) - u(t_k,Q_k,P_k)| Q_k,P_k]| ≤ C(h^2(1+|P_k|^2m)e^-λ(T - t_k)
            + he^-λ(T - t_k)𝔼[(n_G(Q_c,k+1)·P̂_k+1)(1+|P_k+1|^2m)I_G^c(Q̂_k+1)|Q_k,P_k]) ,

where u(t,q,p) is the solution of (<ref>)-(<ref>) with b(q,p) = - ∇ U(q) - γ p and C,λ>0  are constants independent of h, T and P_k.

The proof follows exactly the same steps as in the proof of Lemma <ref> with the modification that Assumption <ref> is used for derivatives' decay to bound the terms in the Taylor expansion.

 
Let  Assumptions <ref> and <ref> hold. For sufficiently small h>0 the following inequality holds for any integer m ≥ 1:

    𝔼(∑_k=0^N-1e^λ t_k(n_G(Q_c, k+1)·P̂_k+1)(1+|P_k+1|^2m)I_G^c(Q̂_k+1) )≤ C e^λ T,

where Q̂_k+1 = Q_k + P̂_k+1h and Q_c, k+1, Q_k+1, P̂_k+1, k=0, …, N-1, are according to the  Markov chain constructed via scheme [PA_c]  and C>0 is a constant independent of T and h.  



If in (<ref>)  g(t,q,p)=𝔼[e^λ t_k(n_G(Q_c)·P̂_1)(1+|P_1|^2m)I_G^c(Q̂_1)| q,p] then the solution of (<ref>)-(<ref>) is 

    v(t_0,q,p)    =  
       𝔼(∑_k=0^N-1 g(t_k,Q_k, P_k) | Q_0 = q , P_0 = p)
       
        =𝔼(∑_k=0^N-1𝔼(e^λ t_k(n_G(Q_c,k+1)·P̂_k+1)(1+|P_k+1|^2m)I_G^c(Q̂_k+1)| Q_k,P_k) | Q_0 = q, P_0 = p )  
       
       =𝔼(∑_k=0^N-1e^λ t_k(n_G(Q_c, k+1)·P̂_k+1)(1 + |P_k+1|^2m)I_G^c(Q̂_k+1)| Q_0 = q, P_0 = p ).

If we can find a function V(t,q,p) with g(t,q,p), satisfying (<ref>)-(<ref>), such that 

    g(t,q,p) ≥𝔼[e^λ t(n_G(Q_c)·P̂_1)(1+|P_1|^2m)I_G^c(Q̂_1)|q,p],

then the inequality v(t_0,Q_0,P_0) ≤ V(t_0,Q_0,P_0) holds. 


As in the proof of Lemma <ref>,  we consider the boundary zone G_R_0. Introduce the function 

    V(t,q,p)=
    
    0,     (t,q,p)∈{T}×G̅×ℝ^d, 
     
    Ae^λ t|p|^2ℓ+L(e^λ T-e^λ t)    
     
    +(B(q)· p)(1+|p|^2m)e^λ t,     (t,q,p)∈ 0,T-h]×G̅×ℝ^d,

where A=max_q∈G̅|B(q)|; integer ℓ≥ 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 <ref>
(i.e., B(q)=n(q^π) for q∈ G_R_0 and q^π is projection
of q on ∂ G). Note that V(t,q,p)≥ 0.

Using Taylor's expansion, we get for any l≥ 1

    𝔼|P_1|^2l   =𝔼|P̂_1|^2l=|p|^2l-2lhγ
    (p· p)|p|^2l-2 - 2lh(∇ U(q) · p)|p|^2l-2
        +2γ/βlh(2(l-1)+d)|p|^2l-2+r_1,

where r_1 satisfies 

    |r_1|≤ Cl^4h^2(1+|p|^2l)

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

    𝔼e^λ (t+h)|P_1|^2ℓ-e^λ t|p|^2ℓ
       ≤ e^λ t[ -ℓ hγ |p|^2ℓ+ℓ hγ[ 
    2ℓ/β (ℓ -1)(2(ℓ -1)+d)] ^ℓ+(λ+1) h|p|^2ℓ
    +h (2 ℓ |∇ U(q)|)^2l] +r_2,

where r_2 satisfies (<ref>) with l=ℓ. 

Applying similar arguments used to obtain (<ref>), we get 

    𝔼[(B(Q_1)· P_1)(1+|P_1|^2m)]=(B(q)· p)(1+|p|^2m)-2
    𝔼[(P̂_1· n_G(Q_c))(1+|P_1|^2m)I_G^c(Q̂_1)]+r_3,

where r_3 satisfies 

    |r_3|≤ Ch(1+|p|^2m)

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

From (<ref>) and (<ref>) and using e^λ
(t+h)≥ e^λ t+λ he^λ t, it follows that 

    𝕋V-V=𝔼Ae^λ (t+h)|P_1|^2ℓ+L(e^λ T-e^λ
    (t+h))+𝔼[(B(Q_1)· P_1)(1+|P_1|^2m)]e^λ (t+h)
          -Ae^λ t|p|^2ℓ-L(e^λ T-e^λ t)-(B(q)·
    p)(1+|p|^2m) 
       ≤   Ae^λ t[ -ℓ hγ |p|^2ℓ+ℓ hγ[ 
    2ℓ/β (ℓ -1)(2(ℓ -1)+d)] ^ℓ+(λ +1) h|p|^2ℓ
    +h (2 ℓ |∇ U(q)|)^2l] 
          -Le^λ tλ h -2𝔼[(P̂_1· n_G(Q_c))(1+|P_1|^2m)I_G^c(
    Q̂_1)]e^λ t+r_2+r_4,

where r_4 satisfies (<ref>). Choosing a sufficiently 
large ℓ and, depending on the choice of ℓ , choosing a sufficiently small h>0 and large L, we get g(t,q,p)=-(𝕋V-V)≥𝔼[
e^λ t(n_G(Q_c)·P̂_1)(1+|P_1|^2m)I_G^c(Q̂
_1)]. 
Hence, the lemma is proved.




 Using Assumption <ref> and Lemma <ref>, we have 

    |𝔼φ(Q_N,P_N) - φ̅| ≤ |𝔼(φ(Q(T),P(T))) - φ̅| +  |𝔼(φ(Q_N,P_N)) - 𝔼(φ(Q(T),P(T)))| 
       ≤ Cexp(-λ T)+ |𝔼(∑_k=0^N-1𝔼(u(t_k+1,Q_k+1,P_k+1) -u(t_k,Q_k,P_k) | Q_k,P_k))|   
       ≤ Cexp(-λ T)+ Ch^2𝔼(|P_k|^2m)∑_k=0^N-1e^-λ (T - t_k)
        + Che^-λ T𝔼(∑_k=0^N-1e^λ t_k(n_G(Q_c,k+1)·P̂_k+1)(1+|P_k+1|^2m)I_G^c(Q̂_k+1) ),

which together with Lemmas <ref> and <ref> gives the desired result.





 §.§ One-step lemma for scheme [OBA_cBO]


Consider G=(0,∞ )× R^d-1. In this case the one step of scheme
[OBA_cBO] from Section <ref> takes the form:

    𝒫_1 =P_0e^-γ h/2+√(σ ^2/2γ(1-
    e^-γ h))ξ,   𝒫_2 =𝒫_1+h/2f(Q_0),   𝒬 =Q_0+h𝒫_2; 
          If 𝒬^1 <0,  τ =-Q_0^1/𝒫_2^1
    , 𝒫_3^1=-𝒫_2^1, 𝒬_c^1=0=Q_0^1-
    τ𝒫_3^1=Q_0^1+τ𝒫_2^1, 
          𝒬_c^i =Q_0^i+τ𝒫_3^i=Q_0^i+τ𝒫
    _2^i , i≠ 1;   
    Q^1 =(h-τ )𝒫_3^1=( h-Q_0^1/𝒫
    _3^1) 𝒫_3^1, 
          Q^i =𝒬_c^i+(h-τ )𝒫
    _3^i=Q_0^i+h𝒫_2^i, i≠ 1; 
          else Q =𝒬, 𝒫_3=𝒫_2;
          𝒫_4 =𝒫_3+h/2f(Q),   
    P =𝒫_4e^-γ h/2+√(σ ^2/2γ(1-
    e^-γ h))ζ.

Here the components ξ^i of ξ and ζ^i of ζ are i.i.d. with the standard normal distribution and 𝒬^i, i= 1,…, d, denotes i-th component of 𝒬 and likewise for q and 𝒫_j, j= 2, 3.

In this section we need the following assumptions. The potential U(q)
satisfies Assumption <ref>, and ∇ U(q) is globally Lipschitz and
derivatives of ∇ U(q) are bounded. To ensure ergodicity (see <cit.> and <cit.>[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 α_1 >0 and 0<α_2 <1 such that 

    1/2(∇ U(q),q)≥α_2 U(q)+γ ^2α_2 (2-α_2 )/8(1-α_2)|q|^2+ α_1.

Further, we assume that the function φ (q,p)∈ C^6,6(G̅×ℝ^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 (<ref>)-(<ref>) with b(q,p)=-∇ U(q)-γ p is sufficiently smooth and, in particular, satisfies the
condition that there are positive constants C and λ 
independent of T such that the following bound holds for some m_1,m_2≥ 1 (cf. (<ref>) in Assumption <ref>): 

    ∑_l=1^3∑_i+|j|=l|D_t^iD_q^ju(t,q,p)|+∑
    _l=1^6∑_|j|=l|D_p^ju(t,q,p)|≤
    C(1+|q|^2m_1+|p|^2m_2)e^-λ (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 [OBA_cBO] are bounded uniformly in time (cf. Lemma <ref>). 

 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 [OBA_cBO] has
the following estimate 

    |𝔼[ u(t+h,Q,P)-u(t,Q_0,P_0)] |≤ Ce^-λ (T-t)(
    h^3+h^2𝔼[ I_G^c(𝒬)] ),

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


Discussion The number of steps is T/h while we expect (cf. Lemmas <ref> and <ref>) that the average number of steps with collisions is independent of h. This together with (<ref>) leads to the claim that [OBA_cBO] (and analogously other schemes of Section <ref>) has second weak convergence order in h, both at finite time and for the ergodic limit. We recall again that the deterministic counterpart of [OBA_cBO] 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 τ being equal to h/2 up to terms of order 𝒪 (h^2).  



We write the one step error as

    R :=[ u(t+h,Q,P)-u(t,Q_0,P_0)]  
          =[ u(t+h,Q,P)-u(t+h/2,Q,𝒫_3)]  
    +[ u(t+h/2,Q,𝒫_3)-u(t+h/2,𝒬_c,𝒫
    _3)] I_G^c(𝒬) 
          +[ u(t+h/2,𝒬_c,𝒫_2)-u(t+h/2,Q_0,𝒫
    _2)] I_G^c(𝒬) 
          +[ u(t+h/2,Q,𝒫_2)-u(t+h/2,Q_0,𝒫_2)]
    I_G(𝒬) 
    +[ u(t+h/2,Q_0,𝒫_2)-u(t+h/2,Q_0,𝒫_1)]  
          +[ u(t+h/2,Q_0,𝒫_1)-u(t+h/2,Q_0,P_0)] +[
    u(t+h/2,Q_0,P_0)-u(t,Q_0,P_0)] .

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

    R    =   [ (h/2-τ )𝒫_3·∇ _qu(t+h/2,𝒬_c,
    𝒫_3)I_G^c(𝒬)] +[ (τ
    -h/2)𝒫_2·∇ _qu(t+h/2,Q_0,𝒫
    _2)I_G^c(𝒬)]  
          +r_1+r_2,

where |r_1|≤ Ce^-λ (T-t)h^3
and |r_2|≤ Ce^-λ (T-t)h^2[ I_G^c(𝒬)
].
Introduce π so that π𝒫_3=𝒫_2, i.e., π𝒫_3=(-
𝒫_3^1,𝒫_3^2,… ,𝒫_3^d).
We have

    R    =   [ (h/2-τ )( 𝒫_3·∇
    _qu(t+h/2,𝒬_c,𝒫_3)-π𝒫_3·∇
    _qu(t+h/2,𝒬_c,π𝒫_3)) I_G^c(𝒬
    )]  
          +r_1+r_3,

where r_3 satisfies the same estimate as r_2 above. 

Let us first consider the following component of (<ref>).
Using the expression for τ , we get 

    R_1   =   𝔼τ𝒫_3^1[ ∂/∂ q^1
    u(t+h/2,0,Q_0^2,… ,Q_0^d,𝒫_3)+∂/∂
    q^1u(t+h/2,0,Q_0^2,… ,Q_0^d,π𝒫_3)] I_G^c(𝒬) 
          =𝔼[Q_0^1v_1(Q_0^2,… ,Q_0^d,𝒫_3) I_G^c(𝒬)],

where 

    v_1(q^2,… ,q^d,p):=∂/∂ q^1
    u(t+h/2,0,q^2,… ,q^d,p)+∂/∂ q^1
    u(t+h/2,0,q^2,… ,q^d,π p).

Let ρ (q,p) be the joint density for  Q_0, 𝒫_3. Assume
that all the required moments with respect to this density exist.
We have 

    R_1=∫_ℝ^2d-2∫_0^∞∫_0^p^1hq^1v_1(q^2,… ,q^d,p)ρ (q,
    p)dq^1dp^1dq^2⋯ dq^ddp^2⋯ dp^d/∫_R^2d-2
    ∫_0^∞∫_0^p^1hρ (q, p)dq^1dp^1dq^2⋯
    dq^ddp^2⋯ dp^d.

Note that

    ∫_0^p^1hq^1ρ (q, p)dq^1=1/2( p^1)
    ^2h^2ρ (0,q^2,… ,q^d, p)+ O(h^3),


    ∫_0^p^1hρ (q, p)dq^1=p^1hρ (0,q^2,… ,q^d,
    p)+ O(h^2).

Hence 

    R_1=h/2∫_ℝ^2d-2∫_0^∞( p^1)
    ^2ρ (p,0,q^2,… ,q^d)U_1(p,q^2,…
    ,q^d)dp^1dq^2⋯ dq^ddp^2⋯ dp^d+ O(h)/
    ∫_R^2d-2∫_0^∞p^1ρ (p,0,q^2,…
    ,q^d)dp^1dq^2⋯ dq^ddp^2⋯ dp^d+ O(h).

Next consider the following component of (<ref>):

    R_2   =   h/2𝔼(𝒫_3^1[ ∂/∂
    q^1u(t+h/2,0,Q_0^2,… ,Q_0^d,𝒫_3)+∂/
    ∂ q^1u(t+h/2,0,Q_0^2,… ,Q_0^d,π𝒫_3)]I_G^c(𝒬))  
       =   h/2∫_ℝ^2d-2∫_0^∞∫_0^p^1hp^1v_1(q^2,… ,q^d,p)ρ
    (q,p)dq^1dp^1dq^2⋯ dq^ddp^2⋯ dp^d/∫_R^2d-2
    ∫_0^∞∫_0^p^1hρ (q,p)dq^1dp^1dq^2⋯
    dq^ddp^2⋯ dp^d
       =   h/2∫_ℝ^2d-2∫_0^∞( p^1) ^2ρ
    (0,q^2,… ,q^d,p)v_1(q^2,… ,q^d,p)dp^1dq^2⋯
    dq^ddp^2⋯ dp^d+ O(h)/∫_R^2d-2∫_0^∞p^1ρ
    (0,q^2,… ,q^d,p)dp^1dq^2⋯ dq^ddp^2⋯ dp^d+ O(h).

Thus |R_1-R_2| ≤ Ch^2 e^-λ (T-t)𝔼[ I_G^c(𝒬)].
The terms for i≠ 1

    A_i   =𝔼( τ -h/2) 𝒫_3^i
       ×[ ∂/
    ∂ q^iu(t+h/2,0,Q_0^2,… ,Q_0^d,𝒫_3)-∂/∂ q^iu(t+h/2,0,Q_0^2,… ,Q_0^d,π𝒫
    _3)]  I_G^c(𝒬)

are analyzed analogously. 





§ NUMERICAL TESTS

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

The following notation is used in this section. The Monte Carlo estimator for 𝔼(φ(Q_N, P_N)) is φ̌_M = 1/M∑_k=1^Mφ(Q_N^(k), P_N^(k)), where (Q_N^(k), P_N^(k)) are independent realizations of (Q_N, P_N); χ denotes the number of collisions along a trajectory and  χ̌_M is the Monte-Carlo estimator for the average number of collisions; and D_M denotes the usual sample variance of Monte Carlo estimators.




In this experiment we consider finite-time weak convergence.
Consider Example <ref> with G := { q^2_1 + q^2_2 < 4 }] and U(q_1, q_2) = -q_1^2 - q_2^2, β = 1, α = 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 = φ̌_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 [PA_c] and two second-order schemes [OBA_cBO] and [BA_cOA_cB]. Tables <ref> and <ref> and Fig. <ref> confirm second-order accuracy of [OBA_cBO] and [BA_cOA_cB]. The first-order method [PA_c] had slightly higher accuracy than first order but clearly lower than second order. 
In this experiment O step in [OBA_cBO] and [BA_cOA_cB] schemes is

    O(p;h,ξ,α):=p e^α h +√(σ ^2/2α(e^2α h - 1))ξ .




tableExperiment <ref>. The global error e and estimated average number of collisions χ̌_M in interval [0,4] for  [PA_c] scheme. 
 

 h     M     φ̌_M± 2 √(D_M)     e      χ̌_M± 2√(D_M)     
 [0.5ex] 

  0.08     4 × 10^8     1.13004 ±  3.5 × 10^-4     0.139994     3.739 ± 1.3 × 10^-4      
 
 0.04    	4 × 10^8    	1.03974 ± 3.2 × 10^-4    	0.049692	    3.679 ± 	1.2 × 10^-4 

0.02    	4 × 10^8     1.00826 ±  3.1 × 10^-4     	0.018207		   	3.660	± 1.2 × 10^-4 

0.01    	4 × 10^8     	0.99762 ± 3.1 × 10^-4      0.007568	     3.654 ± 1.2 × 10^-4 

[1ex]




tableExperiment <ref>. The global error  e and the estimated average number of collisions χ̌_M in interval [0,4] for [OBA_cBO] scheme.
 

 h     M     φ̌_M± 2 √(D_M)     e      χ̌_M± 2√(D_M)     
 [0.5ex] 


  0.08     4 × 10^8     1.08582 ±  3.4 × 10^-4     0.095775      3.714 ± 1.2 × 10^-4      
 
 0.04    	4 × 10^8    	1.01649 ± 3.2 × 10^-4    	0.026441	    3.665 ± 	1.2 × 10^-4 

0.02    	4 × 10^8     0.99683 ±  3.1 × 10^-4    	0.006776	   	3.653	± 1.2 × 10^-4 

0.01    	4 × 10^8     	0.99164 ± 3.1 × 10^-4    	0.001591	     3.650 ±  1.2 × 10^-4   

[1ex]




tableExperiment <ref>. The global error e and the estimated average number of collisions χ̌_M in interval [0,4] for [BA_cOA_cB] scheme. 
 

 h     M     φ̌_M± 2 √(D_M)     e      χ̌_M± 2√(D_M)     
 [0.5ex] 



  0.08     4 × 10^8     1.07393 ±  3.4 × 10^-4     0.083880     3.710 ± 1.2 × 10^-4      
 
 0.04    	4 × 10^8    	1.01308 ± 3.2 × 10^-4    	0.023029	    3.664 ± 	1.2 × 10^-4 

0.02    	4 × 10^8     0.99586 ±   3.1 × 10^-4     	0.005807		   	3.653	± 1.2 × 10^-4 

0.01    	4 × 10^8     	0.99148 ± 3.1 × 10^-4     0.001431	     3.650 ± 1.2 × 10^-4 

[1ex]








 








tableExperiment <ref>. Dependence of φ̌_M on h for the [A_cBOBA_c] and [A_cOBOA_c] schemes. 
 

          [A_cBOBA_c]     [A_cOBOA_c]

 h     M     φ̌_M± 2√(D_M)     φ̌_M± 2√(D_M)    
 [0.5ex] 


  0.5     5 × 10^6     1.32450 ± 7.8852 × 10^-4     1.34389 ± 8.0906× 10^-4  

0.4      5 ×  10^6     1.30362	± 7.7085 × 10^-4     1.31691 ± 7.8433 × 10^-4 

0.2     5 × 10^6     1.27357 ±  7.4647 × 10^-4     1.27655 ± 7.4981 × 10^-4 

0.1     5 × 10^7     1.26542 ± 2.3408 × 10^-4     1.26605 ± 2.3426 × 10^-4

0.08     5 × 10^7     1.26437 ±  2.3383 × 10^-4     1.26497 ± 2.3400 × 10^-4

[1ex]








 


tableExperiment <ref>. Dependence of φ̌_M on h for [BA_cOA_cB] and [BOA_cOB] schemes. 
 

      [BA_cOA_cB]          [BOA_cOB]

 h     M     φ̌_M± 2√(D_M)      M    φ̌_M± 2√(D_M)    
 [0.5ex] 


  0.5     5 × 10^6     1.30502 ± 7.5645 × 10^-4     5 × 10^6    1.42944 ± 8.6166 × 10^-4  

0.4      5 ×  10^6     1.28807	± 7.4837 × 10^-4     5 × 10^6    1.36950 ± 8.1676 × 10^-4 

0.2     5 × 10^6     1.26871 ± 7.4072 × 10^-4     5 × 10^6    1.28988 ± 7.5830 × 10^-4 

0.1     5 × 10^8     1.26403 ± 7.3866 × 10^-5     5× 10^7     1.26922 ± 2.3492 × 10^-4

0.08     5 × 10^8     1.26352 ±  7.3854 × 10^-5     5× 10^7    1.26685 ± 2.3440 × 10^-4

[1ex]










tableExperiment <ref>. Dependence of φ̌_M on h for  [OA_cBA_cO] and [OBA_cBO] schemes. 
  

      [OA_cBA_cO]           [OBA_cBO]    

 h     M     φ̌_M± 2√(D_M)       M    φ̌_M± 2√(D_M)    
 [0.5ex] 


  0.5     5 × 10^6     1.22339 ± 7.0554 × 10^-4     5 × 10^6    1.42944 ± 8.4456 × 10^-4  

0.4      5 ×  10^6     1.23709	± 7.1800 × 10^-4     5 × 10^6    1.3695 ± 8.4456 × 10^-4 

0.2     5 × 10^6     1.25579 ± 7.3281 × 10^-4     5 × 10^6    1.28988 ± 7.5295 × 10^-4 

0.1     5 × 10^8     1.26094 ± 7.3689 × 10^-5     5× 10^7     1.26922 ± 2.3457 × 10^-4

0.08     5 × 10^8     1.26151 ±  7.3736 × 10^-5     5× 10^7    1.26685 ± 2.3417 × 10^-4

[1ex]










In this experiment we test convergence of second-order integators [A_cBOBA_c], [A_cOBOA_c], [BA_cOA_cB], [BOA_cOB], [OA_cBA_cO] and [OBA_cBO] to ergodic limits.

Take G := (1,∞) and ∂ G := {1}.  Consider (<ref>) with the potential

    U(q) = q^2/2

and β = 1, γ = 1.  For numerical simulation, we choose T = 20 as a sufficiently large time to achieve accurate values of the ergodic limit of φ(q) = q^2/2 with respect to this choice of (<ref>). The exact value of φ̅ is 1.2626 (4 d.p.). The starting point of Markov chain is Q_0 = 2 and P_0 = -0.1.  We denote e:=|φ̌_M - φ̅|.
The results are presented in Tables <ref>-<ref> and Fig. <ref>, which clearly show second order convergence for all the integrators. As expected <cit.>,   [BA_cOA_cB] and [OA_cBA_cO] are the most accurate schemes. 











tableExperiment <ref>. The estimate φ̌_M  for [OBA_cBO] scheme.
 

 h     M     φ̌_M± 2 √(D_M)      
 [0.5ex] 


  0.1     5 × 10^6     -4.09714 ± 1.2002× 10^-3  
 
  0.075     5 × 10^7     -4.13173 ± 3.7841 × 10^-4 

  0.05     5 × 10^7     -4.15843 ± 3.7753 × 10^-4  

  0.04     5 × 10^7     -4.16630 ± 3.7727 × 10^-4 

[1ex]







tableExperiment <ref>. The estimate φ̌_M  for [BA_cOA_cB] scheme. 
 

 h     M     φ̌_M± 2 √(D_M)    
 [0.5ex] 

  0.1     5 × 10^6     -4.14494 ± 1.1933 × 10^-3  
 
  0.075     5 × 10^7      -4.16074 ± 3.7709 × 10^-4  

  0.05     5 × 10^7     -4.17151 ± 3.7709  × 10^-4

  0.04     5 × 10^8     -4.17511 ± 1.1920 × 10^-4 

[1ex]









In this experiment we test convergence of second-order integators  [BA_cOA_cB] and [OBA_cBO] to ergodic limits.
Take G:= { q_1^2 + q_2^2 < 4}. Consider (<ref>) with the potential

    U(q_1,q_2) = 1/2(q_1 - q_2)^2 + q_1^2(q_1^2-12)/12 + q_2^2(q_2^2-24)/12

and β = 1, γ = 4. We choose 

    φ(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 φ̅. We take T = 12. We calculate the integral φ̅ 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:=|φ̌_M - φ̅|. The results are presented in Tables <ref>-<ref> and Fig. <ref>, which clearly show second-order convergence for both integrators with [BA_cOA_cB] being more accurate.










tableExperiment <ref>. Dependence of global error e on h for [OBA_cBO] scheme with T = 20.
   
    
    h      M    φ̌_M±  2 √(D_M)     e     χ̌_M±  2 √(D_M)

0.25      10^4     1.697 ± 0.024     2.102     32.52 ± 0.08      

     0.2         10^4     1.912 ± 0.023     1.888     34.71 ± 0.10 

     0.1        10^4     2.860 ± 0.017       0.940      53.01 ± 0.19

     0.08        10^4     3.103 ± 0.014     0.697      61.16 ± 0.21 

     0.05        10^4     3.4514 ± 0.0092     0.3486      78.58 ± 0.26

     0.04        10^4     3.5569 ± 0.0077      0.2431     86.92 ± 0.30 

     0.025        10^4     3.6859 ± 0.0058     0.1141      102.86 ±  0.37

     0.02        10^4     3.7238 ± 0.0053     0.0762      108.89 ± 0.40 
 
     0.01       10^5     3.7785 ± 0.0014      0.0215      120.97 ± 0.15

     0.008       10^7     3.78550 ± 0.00013     0.01450     123.12 ± 0.016 

     0.004       10^7     3.79630 ± 0.00013     0.003670      126.33 ± 0.017

     0.002       10^7     3.79907 ±  0.00013     0.000928     127.24 ±  0.017 

     
    






tableExperiment <ref>. Dependence of global error e on h for [BA_cOA_cB] scheme with T = 20.  
   
    
    h      M    φ̌_M±  2 √(D_M)     e     χ̌_M±  2 √(D_M)

0.25         10^4     1.712 ± 0.024     2.088     33.52 ± 0.08 

    0.2         10^4     1.945 ± 0.023     1.855     34.97 ± 0.10 

    0.1         10^4     2.868 ± 0.017     0.932     53.68 ± 0.19

    0.08         10^4     3.124 ± 0.014     0.676     61.86 ± 0.21 

    0.05         10^4     3.4689 ± 0.0090     0.331     79.92 ± 0.27

    0.04         10^4     3.5690 ± 0.0076     0.2310     88.07 ± 0.31 

    0.025     10^4     3.6859 ± 0.0058      0.1141     103.70 ± 0.38

    0.02         10^4     3.7292 ± 0.0051     0.0708     109.94 ± 0.41 

    0.01         10^5     3.7798 ± 0.0014      0.0202     121.56 ± 0.16

    0.008         10^7     3.78663 ± 0.00013      0.01337     123.538 ± 0.016 

    0.004         10^7     3.79644 ± 0.00013     0.00356     129.494 ± 0.017

    0.002         10^7     3.79910 ± 0.00013     0.00090     127.279 ± 0.017

    

    





tableExperiment <ref>. Stability of  [OBA_cBO] and [BA_cOA_cB]  schemes with T = 20 for large values of h.
    
    
             [BA_cOA_cB]      [OBA_cBO]      

           h      M    φ̌_M±  2 √(D_M)     φ̌_M±  2 √(D_M)

   0.25     10^4      1.712 ± 0.024     1.697 ± 0.024 

   0.3         10^4     1.588 ± 0.024       – 

       0.35         10^4     1.496 ± 0.024     – 

       0.42         10^4     1.469 ± 0.024      –

       0.5         10^4     –      –

    

    



We take G = { q_1^2 + q_2^2 < 4} and ∂ G = { q_1^2 + q_2^2 = 4}. We sample from distribution given by 

    ρ(q,p) = 1/ℤe^K |q|^2 - |p|^2,     q = (q_1,q_2) ∈G̅,  p = (p_1, p_2) ∈ℝ^d,

and compute 

    φ̅ = ∫_G (q_1^2 + q_2^2)e^K (q_1^2 + q_2^2)dq_1 dq_2

via its Monte Carlo estimator based on two schemes [OBA_cBO] and [BA_cOA_cB] which are used to discretize (<ref>). 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  [OBA_cBO] and [BA_cOA_cB] 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 [OBA_cBO] and [BA_cOA_cB] 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 <ref> - <ref> and Fig. <ref> that both numerical schemes [OBA_cBO] and [BA_cOA_cB] have 𝒪(h) convergence then 𝒪(h^3/2) convergence and, as h decreases, we clearly observe 𝒪(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 <ref> that [BA_cOA_cB] is more stable than [OBA_cBO]; it allows the use of step sizes of around 0.42, whereas [OBA_cBO] becomes unstable for a step size greater than about 0.30.






§ 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 [A_c, 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 A_c having local error of order 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.

plain



