Drift-preserving numerical integrators for stochastic Poisson systems

We perform a numerical analysis of randomly perturbed Poisson systems. For the considered It\^o perturbation of Poisson differential equations, we show the longtime behavior of the energy and quadratic Casimirs for the exact solution. We then propose and analyze a drift-preserving splitting scheme for such problems with the following properties: exact drift preservation of energy and quadratic Casimirs, mean-square order of convergence one, weak order of convergence two. Finally, extensive numerical experiments illustrate the performance of the proposed numerical scheme.


Introduction
Deterministic Hamiltonian dynamics are universally applied as mathematical models to describe the dynamical evolution of various physical systems in science and engineering. If the Hamiltonian is not explicitly time dependent, then its value is constant along exact trajectories of the problem. This constant equals the total energy of the system. The recent years have seen plenty of research activities in the design and numerical analysis of energypreserving numerical integrators for deterministic Hamiltonian systems, see for instance [7,8,12,17,22,23,29,30,37,38,34,39,49] and references therein.
In the last few years, the above research has been extended to stochastic differential equations (SDEs) by, for instance, adding a random term to a deterministic Hamiltonian in an additive way, see [9,14,21,19,28,44,43,13]. Observe that extensions to stochastic partial differential equations have been studied in [3,4,15,18,41].
In the present publication, we propose and analyze a drift-preserving scheme for stochastic Poisson systems driven by an additive Brownian motion. Such problems are a direct generalization of the SDEs studied recently in [13], as well as in all the above references, but the proposed numerical integrator is not a trivial generalization of the one given in [13]. In particular, we present a novel numerical scheme that exactly satisfies a trace formula for the expected value of the Hamiltonian and of the Casimir for all times in Section 2. Such longtime behavior corresponds to the one of the exact solution of stochastic Poisson systems and can also be seen as a longtime weak convergence estimate. For the sake of completeness, under general assumptions, we show mean-square and weak convergence of the proposed numerical scheme in Section 3. The main properties of the proposed time integrators are then numerically illustrated in Section 4.

Drift-preserving scheme for stochastic Poisson problem
This section presents the problem, introduces the drift-preserving integrator and shows some of its main geometric properties.
2.1. Setting. For a fixed dimension d, let W ptq P R d denote a standard d-dimensional Wiener process defined for t ą 0 on a probability space equipped with a filtration and fulfilling the usual assumptions. For a fixed dimension m and a smooth potential V : R m Ñ R, let us consider the separable Hamiltonian function of the form We next set Xptq " ppptq, qptqq P R mˆRm and consider the following stochastic Poisson system with additive noise (2) dXptq " BpXptqq∇HpXptqq dt`ˆΣ 0˙d W ptq.
Here, BpXq P R 2mˆ2m is a smooth skew-symmetric matrix and Σ P R mˆd is a constant matrix. In addition, we assume that the initial value X 0 " pp 0 , q 0 q of the above SDE is either non-random or a random variable with bounded moments up to any order (and adapted to the filtration). For simplicity, we assume in the analysis of this paper that px, yq Þ Ñ Bpxq∇Hpyq is globally Lipschitz continuous on R 2mˆR2m and that H and B are C 7 , resp. C 6 -functions with all partial derivatives with at most polynomial growth. This is to ensure existence and uniqueness of solutions to (2) for all times t ą 0 as well as bounded moments at any orders of such solutions. These regularity assumptions on the coefficients B and H will also be used to show strong and weak convergence of the proposed numerical scheme for (2). We observe that one could weaken these assumptions, but this is not the aim of the present work. The present setting covers, for instance, the following examples: simplified versions of the stochastic rigid bodies studied in [45,47], the stochastic Hamiltonian systems considered in [13] by taking he constant canonical symplectic matrix, for which the SDE (2) yields dpptq "´∇V pqptqq dt`Σ dW ptq, dqptq " pptq dt, the Hamiltonian considered in [9] (where the matrix Σ is diagonal), the linear stochastic oscillator from [44], and various stochastic Hamiltonian systems studied in [36,Chap. 4], see also [35], or [42,50,27,26].
Remark 1. We emphasize that our analysis is not restricted to the above form of the Hamiltonian. Indeed, the results below as well as the proposed numerical scheme can be applied to the more general problem (no needed of partitioning the vector X neither to have the separable Hamiltonian (1)) dXptq " BpXptqq∇HpXptqq dt`ˆΣ 0˙d W ptq, as long as the Hessian of the Hamiltonian has a nice structure. One could for instance consider a (linear in p) term of the formṼ pqqp or most importantly the case when the Hamiltonian is quadratic as in the example of a stochastic rigid body. See below for further details.
Applying Itô's lemma to the function HpXq on the solution process Xptq of (2), one obtains Using the skew-symmetry of the matrix BpXq, we have ∇HpXq T BpXq∇HpXq " 0, and using that the partial Hessian ∇ 2 pp HpXq " Id m is a constant matrix, thanks to the separable form (1), rewriting the above equation in integral form and taking the expectation, one then gets the trace formula for the energy of the above problem (2): Hence the expected energy grows linearly with time for all t ą 0. We now would like to design a numerical scheme having the same longtime property.

2.2.
Definition of the numerical scheme. The numerical integrator studied in [13] cannot directly be applied to the stochastic Poisson system (2). Our idea is to combine a splitting scheme with one of the (deterministic) energy-preserving schemes from [17]. Observe that a similar strategy was independently presented in [20] in the particular context of the Langevin equation with other aims than here. We thus propose the following time integrator for problem (2), which is shown in Theorem 4 below to be a drift-preserving integrator for all times: where h ą 0 is the stepsize of the numerical scheme and t n " nh.

Remark 2. (Further extensions)
Let us observe that the (deterministic) energy-preserving scheme from [17] present in the term in the middle of (4) could be replaced by another (deterministic) energy-preserving scheme for (deterministic) Poisson systems, see for example: [8,6,48,10] or a straightforward adaptation of the energy-preserving Runge-Kutta schemes for polynomial Hamiltonians in [11]. Let us further remark that it is also possible to interchange the ordering in the splitting scheme by considering first half a step of the (deterministic) energy-preserving scheme, then a full step of the stochastic part, and finally again half a step of the (deterministic) energy-preserving scheme. Finally, let us add that one could add a damping term in the SDE (2) to compensate for the drift in the energy thus getting conservation of energy for such problems (either in average or a.s.). In this case, one would add the damping term in the first and last equations of the numerical scheme (4) in order to get a (stochastic) energy-preserving splitting scheme. An example of application is Langevin's equation, a widely studied model in the context of molecular dynamics. We do not pursue further this question.
We now show the boundedness along time of all moments of the numerical solution given by (4).
Consider the numerical discretization of the Poisson system with additive noise (2) by the drift-preserving numerical scheme (4) on the compact time interval r0, T s. One then has the following bounds for the numerical moments: for all stepsize h assumed small enough and all m P N, Proof. To show boundedness of the moments of the numerical solution given by (4), we use [36, Lemma 2.2, p. 102], which states that it is sufficient to show the following estimates: where C is independent of h and M n is a random variable with moments of all orders bounded uniformly with respect to all h small enough. The above estimates are a consequence of the definition of the drift-preserving integrator (4) and the linear growth property of the coefficients of the SDE (2) (as a consequence of their Lischitzness).

2.3.
Exact drift preservation of energy. We next show that the numerical scheme (4) satisfies, for all times, the same trace formula for the energy as the exact solution to the SDE (2), i. e., that it is a drift-preserving integrator.

Theorem 4. Consider the numerical discretization of the Poisson system with additive noise
(2) by the drift-preserving numerical scheme (4). Then, for all timestep h assumed small enough, the expected energy of the numerical solution satisfies the following trace formula for all discrete times t n " nh, where n P N.
Proof. The first step of the drift-preserving scheme can be rewritten as and an application of Itô's formula gives Since the second step of the drift-preserving scheme (4) is the deterministic energy-preserving scheme from [17], one then obtains Finally, as in the beginning of the proof, the last step of the numerical integrator provides A recursion now completes the proof.
The above result can also be seen as a longtime weak convergence result and provides a longtime stability result (in a certain sense) for the drift-preserving numerical scheme (4).

2.4.
Splitting methods with deterministic symplectic integrators and backward error analysis: linear case. As symplectic integrators for deterministic Hamiltonian systems or Poisson integrators for deterministic Poisson systems have proven to be very successful [24, Chapters VI and VII], it may be tempting to use them in a splitting scheme for the SDE (2). One could for instance replace the (deterministic) energy-preserving scheme in the middle step of equation (4) by a symplectic or Poisson integrator, such as for instance the second order Störmer-Verlet method [25,Sect. 5] which turns out to be explicit in the context of a separable Hamiltonian (1). Using a backward error analysis, see [40,Chapter 10], [24,Chapter IX], [32,Chapter 5], or [5,Chapter 5], one arrives at the following result in the case of a linear Hamiltonian system with additive noise (2) (i. e. for a quadratic potential V ), where the proposed splitting scheme is drift-preserving for a modified Hamiltonian.
Proposition 5. For a quadratic potential V in (1), consider the numerical discretization of the Hamiltonian system with additive noise (2) (where Bpxq " J for ease of presentation) by the drift-preserving numerical scheme (4), where the energy-preserving scheme in the middle for all discrete times t n " nh, where n P N, and r Proof. By backward error analysis and the theory of modified equations, see for instance [24, Chapter IX], the symplectic Runge-Kutta method Y 1 Þ Ñ Y 2 solves exactly a modified Hamiltonian system with initial condition Y 1 and modified Hamiltonian r H h pxq " Hpxq`Oph p q given by a formal series which turns out to be convergent in the linear case for all h small enough (and with a quadratic modified Hamiltonian). Following the lines of the proof of Theorem 4 applied with the modified Hamiltonian r H h , and observing that the partial Hessian ∇ 2 pp r H h pxq is a constant matrix independent of x (as r H h is quadratic), we deduce the estimate (6) for the averaged modified energy.
Observe in (6) that the constant scalar 1 2 Tr`Σ J r σ h Σ˘" 1 2 Tr`Σ J Σ˘`Oph p q is independent of x and a perturbation of size Oph p q of the drift rate for the exact solution of the SDE in (5).
Finally, note that an analogous result in the nonlinear setting (with nonquadratic potential V in (1)) does not seem straightforward due in particular to the non-boundedness of the moments of the numerical solution over long times and the fact that the modified Hamiltonian r H h pp, qq is nonquadratic with respect to p in general for a nonquadratic potential V .
2.5. Exact drift preservation of quadratic Casimir's. We now consider the case where the ordinary differential equation (ODE) coming from (2), i. e. equation (2) with Σ " 0, has a quadratic Casimir of the form

Recall that a function
CpXq is called a Casimir if ∇CpXq J BpXq " 0 for all X. Along solutions to the ODE, we thus have CpXptqq " Const. This property is independent of the Hamiltonian HpXq.
In this situation, one can show a trace formula for the Casimir as well as a drift-preservation of this Casimir for the numerical integrator (4). Proposition 6. Consider the numerical discretization of the Poisson system with additive noise (2) with the Casimir CpXq by the drift-preserving numerical scheme (4). Then, (1) the exact solution to the SDE (2) has the following trace formula for the Casimir for all times t ą 0. (2) the numerical solution (4) has the same trace formula for the Casimir, for all timestep h assumed small enough, for all discrete times t n " nh, where n P N.
Proof. The above results can be obtain directly by applying Itô's formula and using the property of the Casimir function CpXq.
Stochastic models with such a quadratic Casimir naturally appear for a simplified version of a stochastic rigid body motion of a spacecraft from [45] which has the quadratic Casimir CpXq " X 2 2 or a reduced model for the rigid body in a solvent from [47]. See also the numerical experiments in Section 4.3.

Convergence analysis
In this section, we study the mean-square and weak convergence of the drift-preserving scheme (4) on compact time intervals under the classical setting of globally Lipschitz continuous coefficients. Theorem 7. Let T ą 0. Consider the Poisson problem with additive noise (2) and the drift-preserving integrator (4). Then, there exists h˚ą 0 such that for all 0 ă h ď h˚, the numerical scheme converges with order 1 in the mean-square sense: for all t n " nh ď T , where the constant C is independent of h and n.

Proof. A Taylor expansion of the exact solution to (2) gives
Xphq where the term REST is bounded by h 2 in the mean and by h 3{2 in the mean-square sense. Performing a Taylor expansion for the numerical solution (4) gives, after some lenghty computations, where the term REST is bounded by h 2 in the mean and mean-square sense. The above computations result in the following estimates for the local errors ErXphq´X 1 s " Oph 2 q and Er Xphq´X 1 2 s 1{2 " Oph 3{2 q.  (4). Then, there exists h˚ą 0 such that for all 0 ă h ď h˚, the numerical scheme converges with order 2 in the weak sense: for all Φ P C 6 P pR 2m , Rq, the space of C 6 functions with all derivatives up to order 6 with at most polynomial growth, one has |ErΦpXpt n qqs´ErΦpX n qs| ď Ch 2 , for all t n " nh ď T , where the constant C is independent of h and n.

Numerical experiments
This section presents various numerical experiments in order to illustrate the main properties of the drift-preserving scheme (4), denoted by DP below. In some numerical experiments, we will compare this numerical scheme with classical ones for SDEs such as the Euler-Maruyama scheme (denoted by EM) and the backward Euler-Maruyama scheme (denoted by BEM).
When needed, we use fixed-point iterations in the drift-preserving scheme (4), but one could use Newton iterations as well.

4.1.
The linear stochastic oscillator. The linear stochastic oscillator has extensively been used as a test model since the seminal work [44]. We thus first consider the SDE (2) with BpXq " J the constant 2ˆ2 Poisson matrix and the following Hamiltonian and with Σ " 1 and W scalars. We take the initial values pp 0 , q 0 q " p0, 1q. For this problem, the integral present in the drift-preserving scheme (4) can be computed exactly, resulting in an explicit time integrator. This numerical scheme is different from the one proposed in [13].
We compute the expected values of the energy Hpp, qq using the stepsize h " 5{2 4 , resp. h " 100{2 8 , and the time interval r0, 5s, resp. r0, 100s along various numerical solutions. For this problem, we also use the stochastic trigonometric method, denoted by STM, from [14]. This time integrator is know to preserve the trace formula for the energy for this problem.
The expected values are approximated by computing averages over M " 10 6 samples. The results are presented in Figure 1. One can observe the excellent longtime behavior of the drift-preserving scheme with exact averaged energy drift as stated in Theorem 4. One can also see that the expected value of the Hamiltonian along the Euler-Maruyama scheme drifts exponentially with time. Furthermore, the growth rate of this quantity along the backward Euler-Maruyama scheme is slower than the growth rate of the exact solution to the considered SDE, see [44] for details. These growth rates are qualitatively different from the linear growth rate in the expected value of the Hamiltonian of the exact solution (3), of the STM from [14], and of the drift-preserving scheme (5). In the next numerical experiment, we numerically illustrate the strong convergence rate of the drift-preserving scheme (4) stated in Theorem 7. To do this, we discretize the linear stochastic oscillator on the interval r0, 1s using step sizes ranging from 2´6 to 2´1 0 . The loglog plots of the errors are presented in Figure 2, where mean-square convergence of order 1 for the proposed integrator is observed. The reference solution is computed with the stochastic trigonometric method using h ref " 2´1 2 . The expected values are approximated by computing averages over M " 10 6 samples.
The following numerical experiment illustrates the weak convergence rate of the driftpreserving scheme (4) stated in Theorem 8. In order to avoid Monte Carlo approximations, we focus on weak errors in the first and second moments only, where all the expectations can be computed exactly. We use the same parameters as above except for Σ " 0.1 and step sizes ranging from 2´4 to 2´1 6 . The results are presented in Figure 3, where one can observe weak order 2 in the first and second moments for the drift-preserving scheme.
As symplectic integrators for deterministic Hamiltonian systems have proven to be very successful [24], it may be tempting to use them in a splitting scheme for the SDE (2). In this last numerical experiment, we compare the behavior, with respect to the trace formula, of the drift-preserving scheme and of the symplectic splitting strategies discussed in Section 2.4. We    use the classical Euler symplectic and Störmer-Verlet schemes for the deterministic Hamiltonian and integrate the noisy part exactly. These numerical integrators are denoted by SYMP, resp. ST below. As a comparison with non-geometric numerical integrators, we also use the classical Euler and Heun's schemes in place of a symplectic scheme. These numerical integrators are denoted by splitEULER and splitHEUN. We discretize the linear stochastic oscillator on the time interval r0, 100s with 2 7 stepsizes. The results are presented in Figure 4. The splitting with a non-symplectic schemes behaves as poorly as standard explicit schemes for SDEs. Although not having the exact growth rates, the two symplectic splitting schemes behave much better than the classical Euler-Maruyama scheme with a linear drift in the averaged energy with a perturbed rate, as predicted by Proposition 5.  Hpp, qq " 1 2 p 2´c ospqq and with Σ " 1 and W scalars. We take the initial values pp 0 , q 0 q " p1, ? 2q. We again compare the behavior, with respect to the trace formula, of the DP, SYMP, ST and splitEULER schemes. To do this, we discretize the stochastic mathematical pendulum on the time interval r0, 100s with 2 7 stepsizes. The results are presented in Figure 5. Again, we recover the fact that the drift-preserving scheme exhibits the exact averaged energy drift, as predicted in Theorem 4. Furthermore, one can still observe a good behavior of the symplectic strategies from Section 2.4 analogously to the linear case in Section 4.1, although the analysis in Proposition 5 is only valid for the linear case.

4.3.
Stochastic rigid body. We now consider an Itô version of the stochastic rigid body studied in [33,1,16] for instance. This system has the following Hamiltonian and the skew-symmetric matrix Here, we denote the angular momentum by X " pX 1 , X 2 , X 3 q J and take the moments of inertia to be I " pI 1 , I 2 , I 3 q " p0.345, 0.653, 1q. The initial value for the SDE (2) is given by Xp0q " p0.8, 0.6, 0q and we consider a scalar noise W ptq with Σ " 0.25 (acting on the first component X 1 only).
Observe that, even if the Hamiltonian has not the desired structure (1), one still has a linear drift in the energy since the Hamiltonian is quadratic and thus the Hessian matrix present in the derivation of the trace formula has the correct structure as noted in Remark 1.
We compute the expected values of the energy HpXq and the Casimir CpXq using N " 2 5 stepsizes on the time interval r0, 4s (in order to still see the behavior of the Euler-Maruyama scheme) along various numerical solutions. The expected values are approximated by computing averages over M " 10 6 samples. The results are presented in Figure 6. One can observe the excellent behavior of the drift-preserving scheme as stated in Theorem 4 and Proposition 6. As in the previous numerical experiment, one can also see that the growth rates in H and C of the Euler-Maruyama schemes are qualitatively different from the linear growth rates in the expected values of these quantities of the exact solution of the stochastic rigid body. As for the previous example, we now numerically illustrate the strong convergence rate of the drift-preservation scheme (4) for the stochastic rigid body. To do this, we discretize this SDE on the time interval r0, 0.75s using step sizes ranging from 2´6 to 2´1 0 . The loglog plots of the errors are presented in Figure 7, where mean-square convergence of order 1 for the proposed integrator is observed. The reference solution is computed with the drift-preserving scheme using h ref " 2´1 2 . The expected values are approximated by computing averages over M " 10 5 samples. The last numerical experiment illustrates the weak convergence rate of the drift-preserving scheme (4) stated in Theorem 8. To do this we numerically compute the weak errors in the first and second moments of the first component of the solutions using the parameters: Σ " 0.1, T " 1, M " 2500 samples, and step sizes ranging from 2´1 0 to 2´2 0 . The rest of the parameters are as in the previous numerical experiment. The results are presented in Figure 8. One can observe weak order 2 in the first and second moments for the drift-preserving scheme (up to Monte-Carlo errors).