Variational integrators for stochastic dissipative Hamiltonian systems

TL;DR

This paper develops variational integrators for stochastic dissipative Hamiltonian systems, enabling structure-preserving simulations that maintain key physical properties over long times, with applications to kinetic plasma models.

Contribution

It introduces a general methodology for deriving stochastic variational integrators based on a stochastic Hamiltonian framework, extending geometric numerical methods to stochastic dissipative systems.

Findings

Integrators preserve a discrete stochastic Lagrange-d'Alembert principle.

Numerical tests show superior stability and energy behavior.

Application to Vlasov-Fokker-Planck equation demonstrates effectiveness.

Abstract

Variational integrators are derived for structure-preserving simulation of stochastic forced Hamiltonian systems. The derivation is based on a stochastic discrete Hamiltonian which approximates a type-II stochastic generating function for the stochastic flow of the Hamiltonian system. The generating function is obtained by introducing an appropriate stochastic action functional and considering a stochastic generalization of the deterministic Lagrange-d'Alembert principle. Our approach presents a general methodology to derive new structure-preserving numerical schemes. The resulting integrators satisfy a discrete version of the stochastic Lagrange-d'Alembert principle, and in the presence of symmetries, they also satisfy a discrete counterpart of Noether's theorem. Furthermore, mean-square and weak Lagrange-d'Alembert Runge-Kutta methods are proposed and tested numerically to demonstrate their superior long-time numerical stability and energy behavior compared to non-geometric methods. The Vlasov-Fokker-Planck equation is considered as one of the numerical test cases, and a new geometric approach to collisional kinetic plasmas is presented.

Equations

$$d_{t}q$$

$$d_{t}p$$

$$d_{t}z=A(z)dt+B(z)dW(t),$$

$$\displaystyle A(z)=\begin{pmatrix}\phantom{-}\frac{\partial H}{\partial p}+\frac{1}{2}\sum_{i=1}^{m}\Big{[}\frac{\partial^{2}h_{i}}{\partial p\partial q}\frac{\partial h_{i}}{\partial p}+\frac{\partial^{2}h_{i}}{\partial p^{2}}\Big{(}f_{i}-\frac{\partial h_{i}}{\partial q}\Big{)}\Big{]}\\ -\frac{\partial H}{\partial q}+F+\frac{1}{2}\sum_{i=1}^{m}\Big{[}\Big{(}\frac{\partial^{2}h_{i}}{\partial q\partial p}-\frac{\partial f_{i}}{\partial p}\Big{)}\Big{(}\frac{\partial h_{i}}{\partial p}-f_{i}\Big{)}-\Big{(}\frac{\partial^{2}h_{i}}{\partial q^{2}}-\frac{\partial f_{i}}{\partial q}\Big{)}\frac{\partial h_{i}}{\partial p}\Big{]}\end{pmatrix},\qquad\quad B(z)=\begin{pmatrix}\phantom{-}\big{(}\frac{\partial h}{\partial p}\big{)}^{T}\\ -\big{(}\frac{\partial h}{\partial q}\big{)}^{T}+f\end{pmatrix},$$

$$C([t_{a},t_{b}])=\big{\{}(q,p):\Omega\times[t_{a},t_{b}]\longrightarrow T^{*}Q\,\big{|}\,\text{$q$, $p$ are almost surely continuous $\mathcal{F}_{t}$-adapted semimartingales}\big{\}}.$$

$$\mathcal{B}\big{[}q(\cdot),p(\cdot)\big{]}=p(t_{b})q(t_{b})-\int_{t_{a}}^{t_{b}}\Big{[}p\circ d_{t}q-H\big{(}q(t),p(t)\big{)}\,dt-\sum_{i=1}^{m}h_{i}\big{(}q(t),p(t)\big{)}\circ dW^{i}(t)\Big{]},$$

$$\delta\mathcal{B}\big{[}q(\cdot),p(\cdot)\big{]}\equiv\frac{d}{d\epsilon}\bigg{|}_{\epsilon=0}\mathcal{B}\big{[}q(\cdot)+\epsilon\delta q(\cdot),p(\cdot)+\epsilon\delta p(\cdot)\big{]}.$$

$$\delta\mathcal{B}\big{[}q(\cdot),p(\cdot)\big{]}-\int_{t_{a}}^{t_{b}}F\big{(}q(t),p(t)\big{)}\cdot\delta q(t)\,dt-\sum_{i=1}^{m}\int_{t_{a}}^{t_{b}}f_{i}\big{(}q(t),p(t)\big{)}\cdot\delta q(t)\circ dW^{i}(t)=0\,,$$

$$\displaystyle\delta\mathcal{B}\big{[}q(\cdot),p(\cdot)\big{]}$$

$$\displaystyle\phantom{=}+\int_{t_{a}}^{t_{b}}\bigg{[}\frac{\partial H}{\partial q}\big{(}q(t),p(t)\big{)}\,\delta q(t)+\frac{\partial H}{\partial p}\big{(}q(t),p(t)\big{)}\,\delta p(t)\bigg{]}\,dt$$

$$\displaystyle\phantom{=}+\sum_{i=1}^{m}\int_{t_{a}}^{t_{b}}\bigg{[}\frac{\partial h_{i}}{\partial q}\big{(}q(t),p(t)\big{)}\,\delta q(t)+\frac{\partial h_{i}}{\partial p}\big{(}q(t),p(t)\big{)}\,\delta p(t)\bigg{]}\circ dW^{i}(t),$$

$${\int }_{t_a}^{t_b}p(t)\circ d_t\delta q(t)=p(t_b)\delta q(t_b)-p(t_a)\delta q(t_a)-{\int }_{t_a}^{t_b}\delta q(t)\circ d_{t}p(t).$$

$$\displaystyle\delta\mathcal{B}\big{[}q(\cdot),p(\cdot)\big{]}$$

$$\displaystyle-\int_{t_{a}}^{t_{b}}\delta p(t)\bigg{[}\circ d_{t}q(t)-\frac{\partial H}{\partial p}\big{(}q(t),p(t)\big{)}\,dt-\sum_{i=1}^{m}\frac{\partial h_{i}}{\partial p}\big{(}q(t),p(t)\big{)}\circ dW^{i}(t)\bigg{]},$$

$$\displaystyle\delta\mathcal{B}\big{[}q(\cdot),p(\cdot)\big{]}-\int_{t_{a}}^{t_{b}}F\big{(}q(t),p(t)\big{)}\cdot\delta q(t)\,dt-\sum_{i=1}^{m}\int_{t_{a}}^{t_{b}}f_{i}\big{(}q(t),p(t)\big{)}\cdot\delta q(t)\circ dW^{i}(t)$$

$$\displaystyle=\underbrace{\int_{t_{a}}^{t_{b}}\delta q(t)\Bigg{[}\circ d_{t}p(t)+\bigg{(}\frac{\partial H}{\partial q}\big{(}q(t),p(t)\big{)}-F\big{(}q(t),p(t)\big{)}\bigg{)}\,dt+\sum_{i=1}^{m}\bigg{(}\frac{\partial h_{i}}{\partial q}\big{(}q(t),p(t)\big{)}-f_{i}\big{(}q(t),p(t)\big{)}\bigg{)}\circ dW^{i}(t)\Bigg{]}}_{A}$$

$$\displaystyle-\underbrace{\int_{t_{a}}^{t_{b}}\delta p(t)\bigg{[}\circ d_{t}q(t)-\frac{\partial H}{\partial p}\big{(}q(t),p(t)\big{)}\,dt-\sum_{i=1}^{m}\frac{\partial h_{i}}{\partial p}\big{(}q(t),p(t)\big{)}\circ dW^{i}(t)\bigg{]}}_{B}.$$

$$q(t)=q(t_a)+\underset{M_0(t)}{\underbrace{{\int }_{t_a}^t\frac{\partial H}{\partial p}(q(s),p(s))ds}}+\sum _{i=1}^m\underset{M_i(t)}{\underbrace{{\int }_{t_a}^t\frac{\partial h_i}{\partial p}(q(s),p(s))\circ d{W}^i(s)}},$$

$${\int }_{t_a}^{t_b}\delta p(t)\circ d_{t}q(t)$$

$$={\int }_{t_a}^{t_b}\delta p(t)\circ d_t{M}_0(t)+\sum _{i=1}^m{\int }_{t_a}^{t_b}\delta p(t)\circ d_t{M}_i(t)$$

$$={\int }_{t_a}^{t_b}\delta p(t)\frac{\partial H}{\partial p}(q(t),p(t))dt+\sum _{i=1}^m{\int }_{t_a}^{t_b}\delta p(t)\frac{\partial h_i}{\partial p}(q(t),p(t))\circ d{W}^i(t),$$

$$S(q_{a},p_{b})=\mathcal{B}\big{[}\bar{q}(\cdot;q_{a},p_{b}),\bar{p}(\cdot;q_{a},p_{b})\big{]},$$

$$F^{-}(q_a,p_b)$$

$$F^{+}(q_a,p_b)$$

$$q_b=D_{2}S(q_a,p_b)-F^{+}(q_a,p_b),p_a=D_{1}S(q_a,p_b)-F^{-}(q_a,p_b),$$

$$\frac{\partial S}{\partial q_a}(q_a,p_b)$$

$$\displaystyle\phantom{=}+\int_{t_{a}}^{t_{b}}\bigg{[}\bigg{(}\frac{\partial\bar{q}(t)}{\partial q_{a}}\bigg{)}^{T}\frac{\partial H}{\partial q}\big{(}\bar{q}(t),\bar{p}(t)\big{)}+\bigg{(}\frac{\partial\bar{p}(t)}{\partial q_{a}}\bigg{)}^{T}\frac{\partial H}{\partial p}\big{(}\bar{q}(t),\bar{p}(t)\big{)}\bigg{]}\,dt$$

$$\displaystyle\phantom{=}+\sum_{i=1}^{m}\int_{t_{a}}^{t_{b}}\bigg{[}\bigg{(}\frac{\partial\bar{q}(t)}{\partial q_{a}}\bigg{)}^{T}\frac{\partial h_{i}}{\partial q}\big{(}\bar{q}(t),\bar{p}(t)\big{)}+\bigg{(}\frac{\partial\bar{p}(t)}{\partial q_{a}}\bigg{)}^{T}\frac{\partial h_{i}}{\partial p}\big{(}\bar{q}(t),\bar{p}(t)\big{)}\bigg{]}\circ dW^{i}(t),$$

$$\int_{t_{a}}^{t_{b}}d_{t}\bigg{(}\frac{\partial\bar{q}(t)}{\partial q_{a}}\bigg{)}^{T}\circ\bar{p}(t)=\bigg{(}\frac{\partial\bar{q}(t_{b})}{\partial q_{a}}\bigg{)}^{T}p_{b}-\bar{p}(t_{a})-\int_{t_{a}}^{t_{b}}\bigg{(}\frac{\partial\bar{q}(t)}{\partial q_{a}}\bigg{)}^{T}\circ d_{t}\bar{p}(t),$$

$$\sum _{j=1}^N{\int }_{t_a}^{t_b}{\bar{p}}^j(t)\circ d_t\frac{\partial {\bar{q}}^j(t)}{\partial q_a^i},$$

$$\frac{\partial S}{\partial q_a}(q_a,p_b)=\bar{p}(t_a)$$

$$\displaystyle+\int_{t_{a}}^{t_{b}}\bigg{(}\frac{\partial\bar{p}(t)}{\partial q_{a}}\bigg{)}^{T}\bigg{[}\circ d_{t}\bar{q}-\frac{\partial H}{\partial p}\big{(}\bar{q}(t),\bar{p}(t)\big{)}\,dt-\sum_{i=1}^{m}\frac{\partial h_{i}}{\partial p}\big{(}\bar{q}(t),\bar{p}(t)\big{)}\circ dW^{i}(t)\bigg{]}$$

$$\displaystyle=\bar{p}(t_{a})+\int_{t_{a}}^{t_{b}}\bigg{(}\frac{\partial\bar{q}(t)}{\partial q_{a}}\bigg{)}^{T}\Big{[}F\big{(}\bar{q}(t),\bar{p}(t)\big{)}\,dt+\sum_{i=1}^{m}f_{i}\big{(}\bar{q}(t),\bar{p}(t)\big{)}\circ dW^{i}(t)\Big{]},$$

$$=\bar{p}(t_a)+F^{-}(q_a,p_b),$$

$$\bar{q}(t_b)=D_{2}S(q_a,p_b)-F^{+}(q_a,p_b),\bar{p}(t_a)=D_{1}S(q_a,p_b)-F^{-}(q_a,p_b).$$

$${\Phi }_g^{TQ}(q,\dot{q})$$

$${\Phi }_g^{T^{\ast }Q}(q,p)$$

$$\displaystyle\xi_{Q}(q)=\frac{d}{d\lambda}\bigg{|}_{\lambda=0}\Phi_{\exp\lambda\xi}(q),\qquad\xi_{TQ}(q,\dot{q})=\frac{d}{d\lambda}\bigg{|}_{\lambda=0}\Phi^{TQ}_{\exp\lambda\xi}(q,\dot{q}),\qquad\xi_{T^{*}Q}(q,p)=\frac{d}{d\lambda}\bigg{|}_{\lambda=0}\Phi^{T^{*}Q}_{\exp\lambda\xi}(q,p).$$

Full text

1 Introduction

Stochastic differential equations (SDEs) play an important role in modeling dynamical systems subject to internal or external random fluctuations. Standard references include [10], [54], [62], [69], [89], [100]. Within this class of problems, we are interested in stochastic forced Hamiltonian systems, which take the form

[TABLE]

where $H=H(q,p)$ and $h_{i}=h_{i}(q,p)$ for $i=1,\ldots,m$ are the Hamiltonian functions, $F=F(q,p)$ and $f_{i}=f_{i}(q,p)$ are the forcing terms, $W(t)=(W^{1}(t),\ldots,W^{m}(t))$ is the standard $m$-dimensional Wiener process, and $\circ$ denotes Stratonovich integration. We use $d_{t}$ to denote the stochastic differential of stochastic processes (other than the Wiener process $W(t)$) to avoid confusion with the exterior derivative $d$ of differential forms. The system (1) can be formally regarded as a classical forced Hamiltonian system with the randomized Hamiltonian given by $\widehat{H}(q,p,t)=H(q,p)+\sum_{i=1}^{m}h_{i}(q,p)\circ\dot{W}^{i}(t)$, and the randomized forcing given by $\widehat{F}(q,p,t)=F(q,p)+\sum_{i=1}^{m}f_{i}(q,p)\circ\dot{W}^{i}(t)$, where $H(q,p)$ and $F(q,p)$ are the deterministic Hamiltonian and forcing, respectively, and $h_{i}(q,p)$, $f_{i}(q,p)$ represent the intensity of the noise. Equation (1) is a generalization of stochastic Hamiltonian systems considered in [13], [50], [72], and [93]. Such systems can be used to model, e.g., mechanical systems with uncertainty, or error, assumed to arise from random forcing, limited precision of experimental measurements, or unresolved physical processes on which the Hamiltonian of the deterministic system might otherwise depend. Applications arise in many models in physics, chemistry, and biology. Particular examples include molecular dynamics (see, e.g., [12], [55], [70], [112]), dissipative particle dynamics (see, e.g., [103]), investigations of the dispersion of passive tracers in turbulent flows (see, e.g., [108], [121]), energy localization in thermal equilibrium (see, e.g., [102]), lattice dynamics in strongly anharmonic crystals (see, e.g., [39]), description of noise induced transport in stochastic ratchets (see, e.g., [71]), and collisional kinetic plasmas ([61], [113]).

As occurs for other SDEs, most Hamiltonian SDEs cannot be solved analytically and one must resort to numerical simulations to obtain approximate solutions. In principle, general purpose stochastic numerical schemes for SDEs can be applied to stochastic Hamiltonian systems. However, as for their deterministic counterparts, stochastic Hamiltonian systems possess several important geometric features: in the case of systems without forcing, their phase space flows (almost surely) preserve the symplectic structure ([13], [92], [93]); when the forcing terms are present, then the solutions also satisfy the stochastic Lagrange-d’Alembert principle, as will be shown in Section 2, and in some special cases the phase space flow may be conformally symplectic (see [14], [51], [94]). When simulating these systems numerically, it is therefore advisable that the numerical scheme also preserves such geometric features. Geometric integration of deterministic Hamiltonian systems has been thoroughly studied (see [41], [88], [107] and the references therein) and symplectic integrators have been shown to demonstrate superior performance in long-time simulations of Hamiltonian systems without forcing, compared to non-symplectic methods; so it is natural to pursue a similar approach for stochastic Hamiltonian systems. This is a relatively recent pursuit. Stochastic symplectic integrators are discussed in [5], [7], [8], [9], [22], [25], [33], [52], [80], [81], [92], [93], [95], [118], [127], [128], [130], [132].

Long-time accuracy and near preservation of the Hamiltonian by symplectic integrators applied to deterministic Hamiltonian systems have been rigorously studied using the so-called backward error analysis (see, e.g., [41] and the references therein). To the best of our knowledge, such general rigorous results have not yet been proved for stochastic Hamiltonian systems, but backward error analysis for SDEs is currently an active area of research. Modified SDEs associated with some particular numerical schemes are considered in [1], [31], [32], [111], [129], and [133]. Backward error analysis for the Langevin equation with additive noise is studied for several integrators in [2], [63], and [64]. Recently, backward error analysis for a weak symplectic scheme applied to a stochastic Hamiltonian system has been presented in [6]. Asymptotic preservation of large deviation principles by stochastic symplectic methods is investigated in [29]. The numerical evidence and partial theoretical results to date are promising and suggest that stochastic geometric integrators indeed possess the property of very accurately capturing the evolution of the Hamiltonian $H$ over long time intervals.

An important class of geometric integrators are variational integrators. This type of numerical schemes is based on discrete variational principles and provides a natural framework for the discretization of Lagrangian systems, including forced, dissipative, or constrained ones. These methods have the advantage that they are symplectic when applied to systems without forcing, and in the presence of a symmetry, they satisfy a discrete version of Noether’s theorem. For an overview of variational integration for deterministic systems see [84]; see also [44], [56], [59], [74], [75], [97], [98], [106], [123], [126]. Variational integrators were introduced in the context of finite-dimensional mechanical systems, but were later generalized to Lagrangian field theories (see [83]) and applied in many computations, for example in elasticity, electrodynamics, or fluid dynamics; see [77], [99], [117], [122].

Stochastic variational integrators were first introduced in [16] and further studied in [15]. However, those integrators were restricted to the special case when the Hamiltonian functions $h_{i}=h_{i}(q)$ were independent of $p$, and only low-order Runge-Kutta types of discretization were considered. Stochastic discrete Hamiltonian variational integrators applicable to a general class of Hamiltonian systems were proposed in [50] by generalizing the variational principle for deterministic systems introduced in [75] and applying a Galerkin type of discretization; see also [48]. In the present work we extend the ideas put forth in [50] to forced systems of the form (1) and propose the corresponding Lagrange-d’Alembert variational integrators.

When the forcing terms in Eq. (1) are linear functions of the momentum variable $p$, then the stochastic flow of the system is conformally symplectic (see [94] and Section 2.4). Stochastic conformally symplectic integrators for such systems were proposed in [14], [17], and [51]. Quasi-symplectic integrators were introduced in [94] and further studied in [90]. These ideas are very interesting, but at present seem to be limited only to systems that exhibit a very special form, that is, systems with separable Hamiltonians, linear forcing terms, and additive noise. The stochastic Lagrange-d’Alembert variational integrators introduced in Section 3 are applicable to the general class of systems of the form (1) and preserve their underlying variational structure.

Main content

The main content of the remainder of this paper is, as follows.

In Section 2 we introduce a stochastic Lagrange-d’Alembert principle and a stochastic generating function suitable for considering stochastic forced Hamiltonian systems, and we discuss their properties.

In Section 3 we present a general framework for constructing stochastic Lagrange-d’Alembert variational integrators, prove the discrete stochastic Lagrange-d’Alembert principle, propose mean-square and weak stochastic Lagrange-d’Alembert Runge-Kutta methods, and present several particularly interesting examples of low-stage schemes. We also discuss connections with the idea of quasi-symplectic integrators.

In Section 4 we present the results of our numerical tests, which verify the excellent long-time performance of our integrators compared to some popular non-geometric methods. In particular, as one of the test cases we consider the Vlasov-Fokker-Planck equation, which is used as a model for collisional kinetic plasmas.

Section 5 contains the summary of our work.

2 Lagrange-d’Alembert principle for stochastic forced Hamiltonian systems

The stochastic variational integrators proposed in [16] and [15] were formulated for dynamical systems which are described by a Lagrangian and which are subject to noise whose magnitude depends only on the position $q$. Therefore, these integrators can be extended to (1) only if the Hamiltonian functions $h_{i}=h_{i}(q)$ are independent of $p$ and the Hamiltonian $H$ is non-degenerate (i.e., the associated Legendre transform is invertible). However, in the case of general $h_{i}=h_{i}(q,p)$ the paths $q(t)$ of the system become almost surely nowhere differentiable, which poses a difficulty in interpreting the meaning of the corresponding Lagrangian. To avoid these kind of issues, in [50] an action functional based on a phase space Lagrangian was introduced, and variational integrators for unforced Hamiltonian systems were constructed. In the present work we extend the approach taken in [50] to include forced Hamiltonian systems. To begin, in the next section, we will introduce an appropriate stochastic action functional and show that it can be used to define a type-II generating function for the stochastic flow of the system (1).

2.1 Stochastic Lagrange-d’Alembert principle

Let the Hamiltonian functions $H:T^{*}Q\longrightarrow\mathbb{R}$ and $h_{i}:T^{*}Q\longrightarrow\mathbb{R}$ for $i=1,\ldots,m$ be defined on the cotangent bundle $T^{*}Q$ of the configuration manifold $Q$, and let $(q,p)$ denote the canonical coordinates on $T^{*}Q$. The Hamiltonian forces $F:T^{*}Q\longrightarrow T^{*}Q$ and $f_{i}:T^{*}Q\longrightarrow T^{*}Q$ for $i=1,\ldots,m$ are fiber-preserving mappings with the coordinate representations $F(q,p)=(q,F(q,p))$ and $f_{i}(q,p)=(q,f_{i}(q,p))$, respectively, where by a slight abuse of notation we use the same symbol to denote the force and its local representation. For simplicity, in this work we assume that the configuration manifold has a vector space structure, $Q\cong\mathbb{R}^{N}$, so that $T^{*}Q=Q\times Q^{*}\cong\mathbb{R}^{N}\times\mathbb{R}^{N}$ and $TQ=Q\times Q\cong\mathbb{R}^{N}\times\mathbb{R}^{N}$. In this case, the natural pairing between one-forms and vectors can be identified with the scalar product on $\mathbb{R}^{N}$, that is, $\langle(q,p),(q,\dot{q})\rangle=p\cdot\dot{q}$, where $(q,\dot{q})$ denotes the coordinates on $TQ$. Let $(\Omega,\mathcal{F},\mathbb{P})$ be the probability space with the filtration $\{\mathcal{F}_{t}\}_{t\geq 0}$, and let $W(t)=(W^{1}(t),\ldots,W^{m}(t))$ denote a standard $m$-dimensional Wiener process on that probability space (such that $W(t)$ is $\mathcal{F}_{t}$-measurable). We will assume that the Hamiltonian functions and the forcing terms are sufficiently smooth and satisfy all the necessary conditions for the existence and uniqueness of solutions to (1), and their extendability to a given time interval $[t_{a},t_{b}]$ with $t_{b}>t_{a}\geq 0$. One possible set of such assumptions can be formulated by considering the Itô form of (1),

[TABLE]

with $z=(q,p)$ and

[TABLE]

where $\partial^{2}h_{i}/\partial q^{2}$, $\partial^{2}h_{i}/\partial p^{2}$, and $\partial^{2}h_{i}/\partial q\partial p$ denote the Hessian matrices of $h_{i}$, whereas $\partial h/\partial q$, $\partial h/\partial p$, $\partial f_{i}/\partial q$, and $\partial f_{i}/\partial p$ denote the Jacobian matrices of $h=(h_{1},\ldots,h_{m})$ and $f_{i}$, respectively, and the $n\times m$ forcing matrix $f$ is defined as $f=(f_{1},\ldots,f_{m})$. For simplicity and clarity of the exposition, throughout this paper we assume that (see [10], [54], [62], [69])

• (H1)

$H$ and $h_{i}$ for $i=1,\ldots,m$ are $C^{2}$ functions of their arguments,

• (H2)

$F$ and $f_{i}$ for $i=1,\ldots,m$ are $C^{1}$ functions of their arguments,

• (H3)

$A$ and $B$ are globally Lipschitz.

These assumptions are sufficient for our purposes, but could be relaxed if necessary. Define the space

[TABLE]

Since we assume $T^{*}Q\cong\mathbb{R}^{N}\times\mathbb{R}^{N}$, the space $C([t_{a},t_{b}])$ is a vector space (see [100]). Therefore, we can identify the tangent space $TC([t_{a},t_{b}])\cong C([t_{a},t_{b}])\times C([t_{a},t_{b}])$. We can now define the following stochastic action functional, $\mathcal{B}:\Omega\times C([t_{a},t_{b}])\longrightarrow\mathbb{R}$,

[TABLE]

where $\circ$ denotes Stratonovich integration, and we have omitted writing the elementary events $\omega\in\Omega$ as arguments of functions, following the standard convention in stochastic analysis. For a given curve \big{(}q(t),p(t)\big{)} in $T^{*}Q$ and its arbitrary variation \big{(}\delta q(t),\delta p(t)\big{)}, we define the corresponding variation of the action functional as

[TABLE]

Theorem 2.1** **(Stochastic Lagrange-d’Alembert Principle in Phase Space).

Suppose that $H(q,p)$, $F(q,p)$, and $h_{i}(q,p)$, $f_{i}(q,p)$ for $i=1,\ldots,m$ satisfy conditions (H1)-(H3). If the curve \big{(}q(t),p(t)\big{)} in $T^{*}Q$ satisfies the stochastic forced Hamiltonian system (1) for $t\in[t_{a},t_{b}]$, where $t_{b}\geq t_{a}>0$, then it also satisfies the integral equation

[TABLE]

almost surely for all variations \big{(}\delta q(\cdot),\delta p(\cdot)\big{)}\in C([t_{a},t_{b}]) such that almost surely $\delta q(t_{a})=0$ and $\delta p(t_{b})=0$.

Proof.

Let the curve \big{(}q(t),p(t)\big{)} in $T^{*}Q$ satisfy (1) for $t\in[t_{a},t_{b}]$. It then follows that the stochastic processes $q(t)$ and $p(t)$ are almost surely continuous, $\mathcal{F}_{t}$-adapted semimartingales, that is, \big{(}q(\cdot),p(\cdot)\big{)}\in C([t_{a},t_{b}]) (see [10], [100]). We calculate the variation (2.5) as

[TABLE]

where we have used the end point condition, $\delta p(t_{b})=0$. Since the Hamiltonians are $C^{2}$ and the processes $q(t)$, $p(t)$ are almost surely continuous, in the last two lines we have used a dominated convergence argument to interchange differentiation with respect to $\epsilon$ and integration with respect to $t$ and $W(t)$. Upon applying the integration by parts formula for semimartingales (see [100]), we find

[TABLE]

Substituting and rearranging terms produces,

[TABLE]

where we have used $\delta q(t_{a})=0$. Therefore, we have

[TABLE]

Since \big{(}q(t),p(t)\big{)} satisfy (1), then by definition we have that almost surely for all $t\in[t_{a},t_{b}]$,

[TABLE]

that is, $q(t)$ can be represented as the sum of the semi-martingales $M_{i}(t)$ for $i=0,\ldots,m$, where the sample paths of the process $M_{0}(t)$ are almost surely continuously differentiable. Let us calculate

[TABLE]

where in the last equality we have used the standard property of the Riemann-Stieltjes integral for the first term, as $M_{0}(t)$ is almost surely differentiable, and the associativity property of the Stratonovich integral for the second term (see [100], [54]). Substituting (2.1) in the term $B$ of (2.1), we show that $B=0$. By a similar argument we also prove that $A=0$. Therefore, the left-hand side of (2.1) is equal to zero, almost surely.

∎

Remark: It is natural to expect that the converse theorem, that is, if \big{(}q(\cdot),p(\cdot)\big{)} satisfy the integral principle (2.6), then the curve \big{(}q(t),p(t)\big{)} is a solution to (1), should also hold, although a larger class of variations $(\delta q,\delta p)$ may be necessary. Variants of such a theorem for systems without forcing have been proved in Lázaro-Camí & Ortega [72] and Bou-Rabee & Owhadi [16]. We leave this as an open question. Here, we will use the action functional (2.4) and the Lagrange-d’Alembert principle (2.6) to construct numerical schemes, and we will directly verify that these numerical schemes converge to solutions of (1).

2.2 Stochastic type-II generating function and forcing

When the functions $H(q,p)$, $F(q,p)$, $h_{i}(q,p)$, and $f_{i}(q,p)$ satisfy standard measurability and regularity conditions (e.g., (H1)-(H3)), then the system (1) possesses a pathwise unique stochastic flow $F_{t,t_{0}}:\Omega\times T^{*}Q\longrightarrow T^{*}Q$. It can be proved that for fixed $t,t_{0}$ this flow is mean-square differentiable with respect to the $q$, $p$ arguments, and is also almost surely a diffeomorphism (see [10], [54], [62], [69]). We will show below that the action functional (2.4) can be used to construct a type II generating function for $F_{t,t_{0}}$. Let $(\bar{q}(t),\bar{p}(t))$ be a particular solution of (1) on $[t_{a},t_{b}]$. Suppose that for almost all $\omega\in\Omega$ there is an open neighborhood $\mathcal{U}(\omega)\subset Q$ of $\bar{q}(\omega,t_{a})$, an open neighborhood $\mathcal{V}(\omega)\subset Q^{*}$ of $\bar{p}(\omega,t_{b})$, and an open neighborhood $\mathcal{W}(\omega)\subset T^{*}Q$ of the curve $(\bar{q}(\omega,t),\bar{p}(\omega,t))$ such that for all $q_{a}\in\mathcal{U}(\omega)$ and $p_{b}\in\mathcal{V}(\omega)$ there exists a pathwise unique solution $(\bar{q}(\omega,t;q_{a},p_{b}),\bar{p}(\omega,t;q_{a},p_{b}))$ of (1) which satisfies $\bar{q}(\omega,t_{a};q_{a},p_{b})=q_{a}$, $\bar{p}(\omega,t_{b};q_{a},p_{b})=p_{b}$, and $(\bar{q}(\omega,t;q_{a},p_{b}),\bar{p}(\omega,t;q_{a},p_{b}))\in\mathcal{W}(\omega)$ for $t_{a}\leq t\leq t_{b}$. (As in the deterministic case, for $t_{b}$ sufficiently close to $t_{a}$ one can argue that such neighborhoods exist; see [82].) Define the function $S:\mathcal{Y}\longrightarrow\mathbb{R}$ as

[TABLE]

where the domain $\mathcal{Y}\subset\Omega\times Q\times Q^{*}$ is given by $\mathcal{Y}=\bigcup\limits_{\omega\in\Omega}\{\omega\}\times\mathcal{U}(\omega)\times\mathcal{V}(\omega)$. Define further the two functions $F^{\pm}:\mathcal{Y}\longrightarrow\mathbb{R}^{N}$ as

[TABLE]

Below we prove that the functions $S$ and $F^{\pm}$ generate111A generating function for the transformation $(q_{a},p_{a})\longrightarrow(q_{b},p_{b})$ is a function of one of the variables $(q_{a},p_{a})$ and one of the variables $(q_{b},p_{b})$. Therefore, there are four basic types of generating functions: $S=S_{1}(q_{a},q_{b})$, $S=S_{2}(q_{a},p_{b})$, $S=S_{3}(p_{a},q_{b})$, and $S=S_{4}(p_{a},p_{b})$. In this work we use the type-II generating function $S=S_{2}(q_{a},p_{b})$. the stochastic flow $F_{t_{b},t_{a}}$.

Theorem 2.2**.**

The function $S(q_{a},p_{b})$ is a type-II stochastic generating function and the functions $F^{\pm}(q_{a},p_{b})$ are type-II stochastic exact discrete forces for the stochastic mapping $F_{t_{b},t_{a}}$, that is, $F_{t_{b},t_{a}}:(q_{a},p_{a})\longrightarrow(q_{b},p_{b})$ is implicitly given by the equations

[TABLE]

where the derivatives are understood in the mean-square sense.

Proof.

Under appropriate regularity assumptions on the Hamiltonians and forces (e.g., (H1)-(H3)), the solutions $\bar{q}(t;q_{a},p_{b})$ and $\bar{p}(t;q_{a},p_{b})$ are mean-square differentiable with respect to the parameters $q_{a}$ and $p_{b}$, and the partial derivatives are semimartingales (see [10]). We calculate the derivative of $S$ as

[TABLE]

where for notational convenience we have omitted writing $q_{a}$ and $p_{b}$ explicitly as arguments of $\bar{q}(t)$ and $\bar{p}(t)$. Applying the integration by parts formula for semimartingales (see [100]), we find

[TABLE]

where the left-hand side integral is understood as a column vector with the components given by

[TABLE]

for each $i=1,\ldots,N$. Substituting and rearranging terms, we obtain

[TABLE]

since $(\bar{q}(t),\bar{p}(t))$ is a solution of (1). After performing similar manipulations for $\partial S/\partial p_{b}(q_{a},p_{b})$, together we obtain the result

[TABLE]

By definition of the flow, then $F_{t_{b},t_{a}}(q_{a},\bar{p}(t_{a}))=(\bar{q}(t_{b}),p_{b})$.

∎

2.3 Noether’s theorem for stochastic systems with forcing

Let a Lie group $G$ act on $Q$ by the left action $\Phi:G\times Q\longrightarrow Q$. The Lie group $G$ then acts on $TQ$ and $T^{*}Q$ by the tangent $\Phi^{TQ}:G\times TQ\longrightarrow TQ$ and cotangent $\Phi^{T^{*}Q}:G\times T^{*}Q\longrightarrow T^{*}Q$ lift actions, respectively, given in coordinates by the formulas (see [47], [82])

[TABLE]

where $i,j=1,\ldots,N$ and summation is implied over repeated indices. Let $\mathfrak{g}$ denote the Lie algebra of $G$ and $\exp:\mathfrak{g}\longrightarrow G$ the exponential map (see [47], [82]). Each element $\xi\in\mathfrak{g}$ defines the infinitesimal generators $\xi_{Q}$, $\xi_{TQ}$, and $\xi_{T^{*}Q}$, which are vector fields on $Q$, $TQ$, and $T^{*}Q$, respectively, given by

[TABLE]

The momentum map $J:T^{*}Q\longrightarrow\mathfrak{g}^{*}$ associated with the action $\Phi^{T^{*}Q}$ is defined as the mapping such that for all $\xi\in\mathfrak{g}$ the function $J_{\xi}:T^{*}Q\ni(q,p)\longrightarrow\langle J(q,p),\xi\rangle\in\mathbb{R}$ is the Hamiltonian for the infinitesimal generator $\xi_{T^{*}Q}$, i.e.,

[TABLE]

where \xi_{T^{*}Q}(q,p)=\big{(}q,p,\xi^{q}_{T^{*}Q}(q,p),\xi^{p}_{T^{*}Q}(q,p)\big{)}. The momentum map $J$ can be explicitly expressed as (see [47], [82])

[TABLE]

Noether’s theorem for deterministic Hamiltonian systems relates symmetries of the Hamiltonian to quantities preserved by the flow of the system (see [47], [82]). When the Hamiltonian system is subject to external forces that are orthogonal to the infinitesimal generators of the symmetry group, then the corresponding momentum maps are still conserved (see [84]). It turns out that this result carries over to the stochastic case, as well. A stochastic version of Noether’s theorem for systems without forcing was proved in [13], [50], and [72]. Below we state and provide a proof of Noether’s theorem for stochastic forced Hamiltonian systems.

Theorem 2.3** **(Noether’s theorem for stochastic systems with forcing).

Suppose that the Hamiltonians $H:T^{*}Q\longrightarrow\mathbb{R}$ and $h_{i}:T^{*}Q\longrightarrow\mathbb{R}$ for $i=1,\ldots,m$ are invariant with respect to the cotangent lift action $\Phi^{T^{*}Q}:G\times T^{*}Q\longrightarrow T^{*}Q$ of the Lie group $G$, that is,

[TABLE]

for all $g\in G$. If the forcing terms are orthogonal to the infinitesimal generators of $G$, that is,

[TABLE]

for all $\xi\in\mathfrak{g}$ and $(q,p)\in T^{*}Q$, then the cotangent lift momentum map $J:T^{*}Q\longrightarrow\mathfrak{g}^{*}$ associated with $\Phi^{T^{*}Q}$ is almost surely preserved along the solutions of the stochastic forced Hamiltonian system (1).

Proof.

Equation (2.25) implies that the Hamiltonians are infinitesimally invariant with respect to the action of $G$, that is, for all $\xi\in\mathfrak{g}$ we have

[TABLE]

where $dH$ and $dh$ denote differentials with respect to the variables $q$ and $p$. Let $(q(t),p(t))$ be a solution of (1) and consider the stochastic process $J_{\xi}(q(t),p(t))$, where $\xi\in\mathfrak{g}$ is arbitrary. Using the rules of Stratonovich calculus we can calculate the stochastic differential

[TABLE]

where we used (1), (2.23), (2.24), and (2.27). Therefore, if (2.26) holds, then J_{\xi}\big{(}q(t),p(t)\big{)}=\text{const} almost surely for all $\xi\in\mathfrak{g}$, which completes the proof.

∎

Remark.

When the external forces are not all orthogonal to the infinitesimal generators of the symmetry group, formula (2.3) provides the rate of change of the momentum map.

2.4 Conformal symplecticity and phase space volume

The flow $F_{t,t_{0}}$ for stochastic Hamiltonian systems without forcing almost surely preserves the canonical symplectic two-form

[TABLE]

that is, $F^{*}_{t,t_{0}}\Omega_{T^{*}Q}=\Omega_{T^{*}Q}$, where $F^{*}_{t,t_{0}}$ denotes the pull-back by the flow $F_{t,t_{0}}$ (see [93], [13], [72]). This property does not hold for the general stochastic forced Hamiltonian system (1). However, for certain choices of the forcing terms, the flow may be conformally symplectic, which means that for all $t\geq t_{0}$ there exists a constant (possibly random) $c_{t,t_{0}}\in\mathbb{R}$ such that

[TABLE]

Deterministic conformally symplectic systems are considered in [87]. Conformal symplecticity for the special case of (1) with a separable Hamiltonian, an additive noise, and the forcing terms equal to $F(q,p)=-\nu p$ with a real parameter $\nu$, and $f_{i}(q,p)=0$ for $i=1,\ldots,m$, was considered in [14] and [51]. Below we demonstrate that the property of conformal symplecticity persists for more general cases.

Theorem 2.4** **(Conformal symplecticity).

Suppose that $H(q,p)$, $F(q,p)$, and $h_{i}(q,p)$, $f_{i}(q,p)$ for $i=1,\ldots,m$ satisfy conditions (H1)-(H3). If the forcing terms have the form

[TABLE]

for real parameters $\nu_{i}$, then the stochastic flow $F_{t,t_{0}}$ for (1) is almost surely conformally symplectic with the parameter $c_{t,t_{0}}$ in (2.30) given by

[TABLE]

for all $t\geq t_{0}$.

Proof.

For fixed $(q,p)\in T^{*}Q$, the stochastic process $F_{t,t_{0}}(q,p)$ satisfies the system (1), which can be written as

[TABLE]

where $X$ and $Y_{i}$ are vector fields on $T^{*}Q$, and are given by, respectively,

[TABLE]

Let us calculate the stochastic differential of $F^{*}_{t,t_{0}}\Omega_{T^{*}Q}$. Using the stochastic generalization of the dynamic definition of the Lie derivative (see Theorem 1.2 in [48]), we can write

[TABLE]

where $\pounds_{X}$ and $\pounds_{Y_{i}}$ denote the Lie derivatives with respect to the vector fields $X$ and $Y_{i}$, respectively. Using Cartan’s magic formula (see, e.g., [3]) we have that

[TABLE]

since $d\Omega_{T^{*}Q}=0$, where $i_{X}$ denotes the interior product with the vector field $X$. Substituting (2.34), (2.31), and (2.29), we obtain

[TABLE]

since the Hamiltonian function $H$ is $C^{2}$. In a similar fashion we show that $\pounds_{Y_{i}}\Omega_{T^{*}Q}=-\nu_{i}\Omega_{T^{*}Q}$. Plugging this in (2.35), we obtain a stochastic differential equation of the form

[TABLE]

It is straightforward to verify that the solution of (2.38) that satisfies the initial condition $F^{*}_{t_{0},t_{0}}\Omega_{T^{*}Q}=\Omega_{T^{*}Q}$ has the form

[TABLE]

with $c_{t,t_{0}}$ given by (2.32), which proves the conformal symplecticity of the flow $F^{*}_{t,t_{0}}$. It holds almost surely, since the solution of the SDE (2.38) is pathwise unique (see [10], [54], [62], [69]).

∎

The evolution of stochastic Hamiltonian systems without forcing preserves volumes in phase space, that is, for the standard volume form on $T^{*}Q$ defined as

[TABLE]

we have that $F^{*}_{t,t_{0}}\mu=\mu$. This is a direct consequence of the symplecticity of the flow. Phase space volume preservation does not hold for the general forced system (1), although for certain choices of the forcing terms the flow $F^{*}_{t,t_{0}}$ may possess a property similar to (2.30). Such a property was proved for the special case of (1) with a separable Hamiltonian, an additive noise, and the forcing terms equal to $F(q,p)=-\Gamma p$ with a constant $N\times N$ matrix $\Gamma$, and $f_{i}(q,p)=0$ for $i=1,\ldots,m$ (see [13], [51], [90], [91], [92], [93], [94]). Below we demonstrate that this property holds also for more general cases.

Theorem 2.5** **(Phase space volume evolution).

Suppose that $H(q,p)$, $F(q,p)$, and $h_{i}(q,p)$, $f_{i}(q,p)$ for $i=1,\ldots,m$ satisfy conditions (H1)-(H3). If the forcing terms have the form

[TABLE]

for constant $N\times N$ matrices $\Gamma_{i}$, then the phase space volume form $\mu$ for $t\geq t_{0}$ almost surely evolves according to the formula

[TABLE]

where

[TABLE]

and $F_{t,t_{0}}$ is the stochastic flow for (1).

Proof.

This theorem is a special case of, e.g., Lemma 4.3.1 in [69]. We briefly outline an alternative geometric proof, analogous to the proof of Theorem 2.4. Similar to (2.35), we can write

[TABLE]

Using the property of the divergence operator (see, e.g., [3]), we calculate

[TABLE]

where we have used (2.34) and (2.41), and the fact that the Hamiltonian function $H$ is $C^{2}$. In a similar way we show that $\pounds_{Y_{i}}\mu=-(\operatorname{tr}\Gamma_{i})\cdot\mu$. Therefore, we obtain the SDE of the form

[TABLE]

It is straightforward to verify that the solution that satisfies the initial condition $F^{*}_{t_{0},t_{0}}\mu=\mu$ is given by (2.42) with $b_{t,t_{0}}$ as in (2.43). The formula (2.42) holds almost surely, because the solution of the SDE is pathwise unique (see [10], [54], [62], [69]).

∎

3 Stochastic Lagrange-d’Alembert variational integrators

Suppose we would like to solve (1) on the interval $[0,T]$ with the initial conditions $(q_{0},p_{0})\in T^{*}Q$. Consider the discrete set of times $t_{k}=k\cdot\Delta t$ for $k=0,1,\ldots,K$, where $\Delta t=T/K$ is the time step. In order to determine the discrete curve $\{(q_{k},p_{k})\}_{k=0,\ldots,K}$ that approximates the exact solution of (1) at times $t_{k}$ we need to construct an approximation of the exact stochastic flow $F_{t_{k+1},t_{k}}$ on each interval $[t_{k},t_{k+1}]$, so that $(q_{k+1},p_{k+1})\approx F_{t_{k+1},t_{k}}(q_{k},p_{k})$. A numerical method respecting the underlying Lagrange-d’Alembert principle (2.6) can be constructed by approximating the generating function and forcing terms in (2.15). Let the discrete Hamiltonian function $H^{+}_{d}(q_{a},p_{b};t_{a},t_{b})$ be an approximation of the generating function (2.13), and let the discrete forces $F^{\pm}_{d}(q_{a},p_{b};t_{a},t_{b})$ be approximations of the forcing terms (2.2). The approximate numerical flow $F^{+}_{t_{k+1},t_{k}}:(q_{k},p_{k})\longrightarrow(q_{k+1},p_{k+1})$ is now generated as in (2.20):

[TABLE]

If there is no risk of confusion, we will omit writing the time arguments of $H^{+}_{d}$ and $F^{\pm}_{d}$. We will refer to the scheme (3) as a stochastic Lagrange-d’Alembert variational integrator.

3.1 Discrete stochastic Lagrange-d’Alembert principle

The advantage of the integrator (3) is that it follows from a discrete version of the stochastic Lagrange-d’Alembert principle (2.6). The discrete Lagrange-d’Alembert principle for deterministic Lagrangian systems was proposed in [59]; see also [84]. Below we generalize it to the stochastic case in the setting of Hamiltonian systems defined on the phase space $T^{*}Q$. Define the discrete random curve space $C_{d}$ as

[TABLE]

On that space define the discrete action functional, $\mathcal{B}_{d}:\Omega\times C_{d}\longrightarrow\mathbb{R}$,

[TABLE]

Note that $\mathcal{B}_{d}$ is an approximation of the stochastic action functional (2.4) on the interval $[0,T]$.

Theorem 3.1** **(Discrete stochastic Lagrange-d’Alembert Principle in Phase Space).

Suppose the discrete Hamiltonian $H^{+}_{d}$ is almost surely continuously differentiable, and the discrete forces $F^{\pm}_{d}$ are almost surely continuous with respect to their arguments. The discrete random curve $\{(q_{k},p_{k})\}_{k=0,\ldots,K}$ satisfies the set of equations

[TABLE]

almost surely for $k=1,\ldots,K-1$, if and only if it almost surely satisfies the variational equation

[TABLE]

for all variations $\{(\delta q_{k},\delta p_{k})\}_{k=0,\ldots,K}$ such that $\delta q_{0}=0$ and $\delta p_{K}=0$ almost surely.

Proof.

Consider an arbitrary random curve $\{(q_{k},p_{k})\}_{k=0,\ldots,K}$. Let us calculate the variation $\delta\mathcal{B}_{d}$ corresponding to the arbitrary variation $\{(\delta q_{k},\delta p_{k})\}_{k=0,\ldots,K}$ with $\delta q_{0}=0$ and $\delta p_{K}=0$ (almost surely). We have

[TABLE]

where in the second equality we have shifted the summation index in the $\delta q_{k+1}$ term and used the fact that $\delta q_{0}=0$. It is now straightforward to see that if the set of equations (3.1) is satisfied, then the variational equation (3.5) holds almost surely. Conversely, if the variational equation (3.5) holds for all variations $\{(\delta q_{k},\delta p_{k})\}_{k=0,\ldots,K}$ with $\delta q_{0}=0$ and $\delta p_{K}=0$, then the set of equations (3.1) has to be satisfied almost surely.

∎

3.2 Discrete Noether’s theorem for stochastic systems with forcing

Another advantage of the integrator (3) is that one can prove a discrete counterpart of Theorem 2.3. If the discrete system inherits the symmetries of the continuous problem, then the evolution of the momentum maps will be accurately captured by the numerical solution. Discrete Noether’s theorem for systems described by a type-II generating function was first proved for deterministic systems in [75], and later generalized to the stochastic case in [50]. Discrete Noether’s theorem for deterministic Lagrangian systems with forcing was first proposed in [84]. Below we combine these ideas and formulate a version of discrete Noether’s theorem applicable to discrete systems described by (3). Let $R_{d}:\Omega\times Q\times T^{*}Q\longrightarrow\mathbb{R}$ be the generalized discrete stochastic Lagrangian defined as

[TABLE]

Consider the action of the Lie group $G$ on $Q\times T^{*}Q$ given by

[TABLE]

For any $\xi\in\mathfrak{g}$ the corresponding infinitesimal generator on $Q\times T^{*}Q$ is then given by

[TABLE]

Theorem 3.2** **(Discrete Noether’s theorem for stochastic systems with forcing).

Suppose the generalized discrete stochastic Lagrangian $R_{d}:\Omega\times Q\times T^{*}Q\longrightarrow\mathbb{R}$ is invariant under the action of the Lie group $G$, that is,

[TABLE]

If the discrete forces $F^{\pm}_{d}$ satisfy the condition

[TABLE]

for all $(q_{k},q_{k+1},p_{k+1})\in Q\times T^{*}Q$, then the cotangent lift momentum map $J$ associated with $\Phi^{T^{*}Q}$ is almost surely preserved along the solutions of the discrete equations (3), i.e., a.s. $J(q_{k+1},p_{k+1})=J(q_{k},p_{k})$.

Proof.

Since the generalized discrete Lagrangian $R_{d}$ is invariant with respect to the actions of $G$, for an arbitrary $\xi\in\mathfrak{g}$ we have

[TABLE]

where we have used the fact that $\xi^{q}_{T^{*}Q}(q_{k+1},p_{k+1})=\xi_{Q}(q_{k+1})$. Assume that $q_{k}$, $q_{k+1}$, and $p_{k+1}$ satisfy the discrete evolution equation (3). By substituting (3) in (3.2), we obtain

[TABLE]

This can be rewritten as

[TABLE]

where we have used the definition of the cotangent lift momentum map (2.24). If the condition (3.11) holds, then we have $J_{\xi}(q_{k+1},p_{k+1})=J_{\xi}(q_{k},p_{k})$. The result holds almost surely, because equation (3) is satisfied almost surely. ∎

Remark.

When the discrete forces do not satisfy the condition (3.11), equation (3.14) provides the rate of change of the momentum map, which mimicks formula (2.3) in the continuous case.

3.3 Mean-square Lagrange-d’Alembert partitioned Runge-Kutta methods

3.3.1 Construction

Partitioned Runge-Kutta methods for deterministic forced Hamiltonian systems have been proposed in [57] and [84]. A general class of stochastic mean-square Runge-Kutta methods for Stratonovich ordinary differential equations was introduced and analyzed in [19], [20], and [21]. These ideas were later used by Ma & Ding & Ding [80] and Ma & Ding [81] to construct symplectic Runge-Kutta methods for stochastic Hamiltonian systems without forcing; see also [50]. Below we combine these ideas and introduce mean-square Lagrange-d’Alembert partitioned Runge-Kutta methods for stochastic forced Hamiltonian systems of the form (1).

Definition 3.3**.**

An $s$-stage mean-square Lagrange-d’Alembert partitioned Runge-Kutta method for the system (1) is given by

[TABLE]

where $\Delta t$ is the time step, $\Delta W=(\Delta W^{1},\ldots,\Delta W^{m})$ are the increments of the Wiener process, $Q_{i}$ and $P_{i}$ for $i=1,\ldots,s$ are the position and momentum internal stages, respectively, and the coefficients of the method $a_{ij}$, $\bar{a}_{ij}$, $\hat{a}_{ij}$, $b_{ij}$, $\bar{b}_{ij}$, $\hat{b}_{ij}$, $\alpha_{i}$, $\hat{\alpha}_{i}$, $\beta_{i}$, and $\hat{\beta}_{i}$ satisfy the conditions

[TABLE]

for $i,j=1,\ldots,s$.

The partitioned Runge-Kutta method (3.15) can be represented by the tableau

[TABLE]

where $a=(a_{ij})_{i,j=1\ldots s}$, $\alpha=(\alpha_{i})_{i=1\ldots s}$, etc. The set of equations (3.15) forms a one-step numerical scheme. Knowing $q_{k}$ and $p_{k}$ at time $t_{k}$, one can solve Equations (3.15a)-(3.15) for the internal stages $Q_{i}$ and $P_{i}$, and then use (3.15c)-(3.15) to determine $q_{k+1}$ and $p_{k+1}$ at time $t_{k+1}$. If given $q_{k}$ and $p_{k+1}$ instead, one can also solve (3.15) for the remaining variables $Q_{i}$, $P_{i}$, $q_{k+1}$ and $p_{k}$. Note that since we have only used $\Delta W^{r}=\int_{t_{k}}^{t_{k+1}}dW^{r}(t)$ in (3.15), we can in general expect mean-square convergence of order 1.0 at most. To obtain mean-square convergence of higher order we would also need to include higher-order multiple Stratonovich integrals, e.g., to achieve convergence of order 1.5 we would need to include terms involving $\Delta Z^{r}=\int_{t_{k}}^{t_{k+1}}\int_{t_{k}}^{t}dW^{r}(\xi)\,dt$ (see [21], [92], [93]). Below we prove that the Runge-Kutta method (3.15) with the conditions (3.16) is indeed a stochastic Lagrange-d’Alembert method of the form (3).

Theorem 3.4**.**

The $s$-stage mean-square partitioned Runge-Kutta method (3.15) with the conditions (3.16) is a stochastic Lagrange-d’Alembert variational integrator of the form (3) with the discrete Hamiltonian

[TABLE]

and the discrete forces

[TABLE]

where $q_{k+1}$, $p_{k}$, $Q_{i}$, and $P_{i}$ satisfy the system of equations (3.15) and are understood as functions of $q_{k}$ and $p_{k+1}$.

Proof.

The proof involves straightforward, although rather lengthy and tedious algebraic manipulations. Therefore, for the clarity and brevity of the exposition, we only consider the one-dimensional noise case $m=1$ and point out the key steps of the derivations. Let us introduce the following shorthand notation:

[TABLE]

Differentiate each of the equations (3.15) with respect to $q_{k}$ and $p_{k+1}$ to express the Jacobians $\partial Q_{i}/\partial q_{k}$, $\partial P_{i}/\partial q_{k}$, $\partial q_{k+1}/\partial q_{k}$, $\partial p_{k}/\partial q_{k}$, and analogous Jacobians with respect to $p_{k+1}$, in terms of the derivatives of the terms (3.3.1). For instance, we have

[TABLE]

where $I$ denotes the $N\times N$ identity matrix. Let us now calculate the derivative of the discrete Hamiltonian (3.18) with respect to $p_{k+1}$. After substituting the Jacobians (3.21) and using (3.15) to replace $p_{k+1}$, we obtain the expression

[TABLE]

After using (3.16a)-(3.16d) in the last four terms (e.g., $\alpha_{i}\alpha_{j}-\alpha_{j}a_{ji}=\alpha_{i}\bar{a}_{ij}$), and substituting (3.15) for $P_{i}$, we get

[TABLE]

By using the conditions (3.16e)-(3.16h) and collecting terms, we finally arrive at

[TABLE]

In a similar fashion we derive

[TABLE]

which completes the proof.

∎

3.3.2 Convergence

Mean-square convergence concentrates on pathwise approximations of the exact solutions (see [62], [89]). Let $\bar{z}(t)=(\bar{q}(t),\bar{p}(t))$ be the exact solution to (1) with the initial conditions $q_{0}$ and $p_{0}$, and let $z_{k}=(q_{k},p_{k})$ denote the numerical solution at time $t_{k}$ obtained by applying (3.15) iteratively $k$ times with the constant time step $\Delta t$. The numerical solution is said to converge in the mean-square sense with global order $r$ if there exist $\delta>0$ and a constant $C>0$ such that for all $\Delta t\in(0,\delta)$ we have

[TABLE]

where $T=K\Delta t$, as defined before, and $E$ denotes the expected value. In principle, in order to determine the mean-square order of convergence of the Lagrange-d’Alembert partitioned Runge-Kutta method (3.15) we need to calculate the power series expansions of $q_{k+1}$ and $p_{k+1}$ in terms of the powers of $\Delta t$ and $\Delta W^{i}$, and compare them to the Stratonovich-Taylor expansions for the exact solution $\bar{q}(t_{k}+\Delta t)$ and $\bar{p}(t_{k}+\Delta t)$ (see [21], [62], [89]). As mentioned in Section 3.3.1, the mean-square order of the method (3.15) cannot exceed 1.0. Below we provide the conditions that have to be satisfied by the coefficients of the method (3.15) in order for it to be convergent.

Theorem 3.5**.**

Suppose that, in addition to conditions (H1)-(H3), the functions $H(q,p)$, $F(q,p)$, and $h_{i}(q,p)$, $f_{i}(q,p)$ for $i=1,\ldots,m$ have all the necessary partial derivatives. Let the coefficients of the method (3.15) satisfy the conditions

[TABLE]

If the noise is commutative, that is, if the following conditions are satisfied

[TABLE]

where the vectors $\Gamma_{ij}$ and $\Lambda_{ij}$ for each $i,j=1,\ldots,m$ are defined as

[TABLE]

then the method (3.15) is convergent with mean-square order 1.0. If the noise is noncommutative, then the method (3.15) is convergent with mean-square order 0.5.

Proof.

General order conditions for stochastic non-partitioned Runge-Kutta methods have been analyzed in [20] and [21]. Conditions for mean-square convergence of order 1.0 for stochastic partitioned Runge-Kutta methods with a one-dimensional noise have been derived in [81]. However, the method (3.15) is more general, as we allow a multidimensional noise, and different coefficients are applied to the Hamiltonian and forcing terms, but the method of proof is similar to the proof of Theorem 2.1 in [81], therefore we only present the main steps. To simplify the notation, denote $\alpha=(\alpha_{1},\ldots,\alpha_{s})^{T}$, $b=(b_{ij})_{i,j=1,\ldots,s}$, and similarly for the remaining coefficients of the method. Let also $e=(1,1,\ldots,1)^{T}$ be an $s$-dimensional vector. Then the conditions (3.5) can be written more compactly, e.g., $\alpha^{T}e=1$ or $\beta^{T}be=1/2$. We first determine power expansions of the internal stages $Q_{i}$ and $P_{i}$ in terms of the powers of $\Delta t$ and $\Delta W^{i}$. We plug in series expansions for $Q_{i}$ and $P_{i}$ in Equations (3.15a)-(3.15), and determine their coefficients by expanding the derivatives of the Hamiltonians and forcing terms into Taylor series around $(q_{k},p_{k})$. Then we plug in thus found series expansions into Equations (3.15c)-(3.15), and again expand the derivatives of the Hamiltonians and forcing terms into Taylor series around $(q_{k},p_{k})$. This way we obtain the series expansions of $q_{k+1}$ and $p_{k+1}$ as

[TABLE]

where the vectors $\bar{\Gamma}_{ij}$ and $\bar{\Lambda}_{ij}$ for each $i,j=1,\ldots,m$ are defined as

[TABLE]

and the forcing terms and the derivatives of the Hamiltonians are evaluated at $(q_{k},p_{k})$. Let $\bar{q}(t;q_{k},p_{k})$ and $\bar{p}(t;q_{k},p_{k})$ denote the exact solution of (1) such that $\bar{q}(t_{k};q_{k},p_{k})=q_{k}$ and $\bar{p}(t_{k};q_{k},p_{k})=p_{k}$. Using (1) we calculate the Stratonovich-Taylor expansions for $\bar{q}(t_{k+1};q_{k},p_{k})$ and $\bar{p}(t_{k+1};q_{k},p_{k})$ as (see [62])

[TABLE]

where $J_{ij}=\int_{t_{k}}^{t_{k+1}}\int_{t_{k}}^{t}dW^{i}(\tau)\circ dW^{j}(t)$ denotes a double Stratonovich integral, $\Gamma_{ij}$ and $\Lambda_{ij}$ have been defined in (3.5), and the forcing terms and the derivatives of the Hamiltonians are again evaluated at $(q_{k},p_{k})$. Assuming the conditions (3.5) are satisfied, we have that $\bar{\Gamma}_{ij}=\Gamma_{ij}$ and $\bar{\Lambda}_{ij}=\Lambda_{ij}$, but comparing (3.3.2) and (3.3.2), we find that in the general case of noncommutative noise not all first order terms agree, and therefore we only have the local error estimates

[TABLE]

Theorem 1.1 from [89] then implies that the method (3.15) has mean-square order 0.5. However, if the noise is commutative, then using the property $J_{ij}+J_{ji}=\Delta W^{i}\Delta W^{j}$ (see [62], [89]), one can easily show

[TABLE]

In that case all first-order terms in the expansions (3.3.2) and (3.3.2) agree, and we have the local error estimates

[TABLE]

Theorem 1.1 from [89] then implies that the method (3.15) has mean-square order $1.0$.

∎

In the case of a one-dimensional noise the commutation condition (3.28) is trivially satisfied, therefore we have the following corollary.

Corollary 3.6**.**

Under the assumptions of Theorem 3.5, the method (3.15) is convergent with mean-square order 1.0 for systems driven by a one-dimensional noise.

3.3.3 Examples

In the construction of the integrator (3.15) we may choose the number of stages $s$. In the deterministic case, the higher the number of stages, the higher order of convergence can be achieved (see [41], [42], [43]). In our case, however, as explained earlier, we cannot in general achieve mean-square order of convergence higher than 1.0, because we only used $\Delta W^{r}$ in (3.15). Since the system (3.15a)-(3.15) requires solving $2sN$ equations for $2sN$ variables, from the computational point of view it makes sense to only consider methods with low values of $s$. In this work we focus on the following classical numerical integration formulas (one can easily verify that the conditions (3.16) and (3.5) are satisfied for the discussed methods).

Stochastic midpoint method

Using the midpoint rule we obtain a one-stage non-partitioned Runge-Kutta method represented by the tableau

[TABLE]

Noting that $Q_{1}=(q_{k}+q_{k+1})/2$ and $P_{1}=(p_{k}+p_{k+1})/2$, this method can be written as

[TABLE]

The stochastic midpoint method was considered in [93] and [81] in the context of symplectic integrators for stochastic Hamiltonian systems without forcing; see also [50]. This example demonstrates that the stochastic midpoint method retains its geometric properties also for forced systems. It is an implicit method and in general one has to solve $2N$ equations for $2N$ unknowns. However, if the Hamiltonians are separable, that is, $H(q,p)=T_{0}(p)+U_{0}(q)$ and $h_{i}(q,p)=T_{i}(p)+U_{i}(q)$, then $q_{k+1}$ from the first equation can be substituted into the second one. In that case only $N$ nonlinear equations have to be solved for $p_{k+1}$. 2. 2.

Stochastic Störmer-Verlet method

A generalization of the classical Störmer-Verlet method can be obtained by choosing the tableau

[TABLE]

Noting that $Q_{1}=q_{k}$, $Q_{2}=q_{k+1}$, and $P_{1}=P_{2}$, this method can be more efficiently written as

[TABLE]

This method was considered in [81] in the context of symplectic integrators for stochastic Hamiltonian systems without forcing; see also [50]. It is particularly efficient, because the first equation can be solved separately from the second one, and the last equation is an explicit update. Moreover, if the Hamiltonians are separable, the second equation becomes explicit. If in addition the forcing terms $F$ and $f_{i}$ have special forms, then further improvements in efficiency are possible. For instance, if the forcing terms depend linearly on $p$, as is often the case in practical applications, then the first equation is a linear equation for $P_{1}$, and can be solved using linear solvers. In case the forcing terms are independent of $p$ altogether, then the whole method becomes fully explicit. 3. 3.

2-stage stochastic DIRK method

In order to reduce the computational cost of solving nonlinear equations, diagonally implicit Runge-Kutta (DIRK) methods use lower-triangular tableaus (see [41], [42], [43]). One can easily verify that the most general family of 2-stage stochastic DIRK methods that satisfy the conditions (3.16) and (3.5) has a tableau of the form

[TABLE]

where $\lambda\in\mathbb{R}$ is an arbitrary parameter. One can check that for $\lambda=0$ and $\lambda=1$, this method reduces to the stochastic midpoint method (1). For other choices of $\lambda$, one needs to solve equations (3.15a) and (3.15), first for $i=1$ ($2N$ equations) in order to calculate the internal stages $Q_{1}$ and $P_{1}$ ($2N$ variables), and then for $i=2$ ($2N$ equations) to find the internal stages $Q_{2}$ and $P_{2}$ ($2N$ variables). If the Hamiltonians are separable, then equations (3.15a) can be substituted into equations (3.15), and the problem is reduced to solving two systems of $N$ equations each.

Note that the methods (1), (2), and (3.40) are in general implicit. One can use the Implicit Function Theorem to show that for sufficiently small $\Delta t$ and $|\Delta W^{i}|$, the relevant nonlinear equations will have a solution. However, since the increments $\Delta W^{i}$ are unbounded, for some values of $\Delta W^{i}$ solutions might not exist. To avoid problems with numerical implementations, if necessary, one can replace $\Delta W^{i}$ in equations (1) and (2) with the truncated random variables $\overline{\Delta W^{i}}$ defined as

[TABLE]

where $A>0$ is suitably chosen for the considered problem. See [23] and [93] for more details regarding schemes with truncated random increments and their convergence.

3.4 Weak Lagrange-d’Alembert Runge-Kutta methods

3.4.1 Construction

A general class of weak stochastic Runge-Kutta methods for Stratonovich ordinary differential equations was introduced and analyzed in [104] and [105]. These ideas were later used by Wang & Hong & Xu [130] to construct weak symplectic Runge-Kutta methods for stochastic Hamiltonian systems without forcing. Below we combine these ideas and introduce weak Lagrange-d’Alembert Runge-Kutta methods for stochastic forced Hamiltonian systems of the form (1).

Definition 3.7**.**

An $s$-stage weak Lagrange-d’Alembert Runge-Kutta method for the system (1) is given by

[TABLE]

where $\Delta t$ is the time step, $\hat{I}_{1},\ldots,\hat{I}_{m}$ are independent three-point distributed random variables with $P(\hat{I}_{r}=\pm\sqrt{3\Delta t})=1/6$ and $P(\hat{I}_{r}=0)=2/3$ , $Q^{(0)}_{i}$, $Q^{(l)}_{i}$, $P^{(0)}_{i}$, and $P^{(l)}_{i}$ for $i=1,\ldots,s$ and $l=1,\ldots,m$ are the position and momentum internal stages, respectively, and the coefficients of the method $a^{(0)}_{ij}$, $a^{(1)}_{ij}$, $b^{(0)}_{ij}$, $b^{(1)}_{ij}$, $b^{(3)}_{ij}$, $\alpha_{i}$, $\beta_{i}$ satisfy the conditions

[TABLE]

for $i,j=1,\ldots,s$.

The Runge-Kutta method (3.42) can be represented by the tableau

[TABLE]

where $a^{(0)}=(a^{(0)}_{ij})_{i,j=1\ldots s}$, $\alpha=(\alpha_{i})_{i=1\ldots s}$, etc. The set of equations (3.42) forms a one-step numerical scheme. Knowing $q_{k}$ and $p_{k}$ at time $t_{k}$, one can solve Equations (3.42a)-(3.42) for the internal stages $Q^{(0)}_{i}$, $Q^{(l)}_{i}$, $P^{(0)}_{i}$ and $P^{(l)}_{i}$, and then use (3.42e)-(3.42) to determine $q_{k+1}$ and $p_{k+1}$ at time $t_{k+1}$. Depending on the choice of the coefficients, the method (3.42) is in general implicit. However, since the random variables $\hat{I}_{l}$ are bounded, one can show that for sufficiently small $\Delta t$, the relevant nonlinear equations will have a solution. Below we prove that the Runge-Kutta method (3.42) with the conditions (3.43) is indeed a stochastic Lagrange-d’Alembert method of the form (3).

Theorem 3.8**.**

The $s$-stage weak Runge-Kutta method (3.42) with the conditions (3.43) is a stochastic Lagrange-d’Alembert variational integrator of the form (3) with the discrete Hamiltonian

[TABLE]

and the discrete forces

[TABLE]

where $q_{k+1}$, $p_{k}$, $Q^{(0)}_{i}$, $Q^{(r)}_{i}$, $P^{(0)}_{i}$, and $P^{(r)}_{i}$, satisfy the system of equations (3.42) and are understood as functions of $q_{k}$ and $p_{k+1}$.

Proof.

The proof is analogous to the proof of Theorem 3.4. ∎

Remark.

For stochastic Hamiltonian systems without forcing, i.e. $F\equiv 0$, $f_{r}\equiv 0$, the method (3.42) reduces to a weak symplectic Runge-Kutta method of the type introduced in [130]. Therefore, in that case Theorem 3.8 also provides a type-II generating function for such a family of methods, and consequently an alternative proof of their symplecticity.

3.4.2 Convergence

Rather than precisely approximating each sample path, weak convergence concentrates on approximating the probability distribution and functionals of the exact solution (see [62], [89]). Let $\bar{z}(t)=(\bar{q}(t),\bar{p}(t))$ be the exact solution to (1) with the initial conditions $q_{0}$ and $p_{0}$, and let $z_{k}=(q_{k},p_{k})$ denote the numerical solution at time $t_{k}$ obtained by applying (3.42) iteratively $k$ times with the constant time step $\Delta t$. The numerical solution is said to converge weakly with weak global order $r$ if for each $\varphi\in C^{2(r+1)}_{P}(T^{*}Q,\mathbb{R})$ there exists $\delta>0$ and a constant $C>0$ such that for all $\Delta t\in(0,\delta)$ we have

[TABLE]

where $T=K\Delta t$, and $C^{\alpha}_{P}(T^{*}Q,\mathbb{R})$ denotes the space of all $\varphi\in C^{\alpha}(T^{*}Q,\mathbb{R})$ with polynomial growth, i.e., there exists a constant $A>0$ and $\gamma\in\mathbb{N}$ such that $|\partial^{\beta}_{z}\varphi(z)|\leq A(1+\|z\|^{2\gamma})$ for all $z\in T^{*}Q$ and any partial derivative of order $\beta\leq\alpha$. Weak convergence of the Runge-Kutta methods of type (3.42) has been analyzed, and the relevant order conditions for the coefficients have been derived in [105].

3.4.3 Examples

In [130] a number of weak symplectic Runge-Kutta methods for stochastic Hamiltonian systems without forcing have been proposed. Since the symplecticity conditions derived in [130] are equivalent to the conditions (3.43), these methods become Lagrange-d’Alembert integrators when applied to systems with forcing. In this work, we particularly focus on two methods, namely $SRKw1$ and $SRKw2$, as dubbed in [130].

SRKw1

The family of 1-stage $SRKw1$ methods is defined by the tableau

[TABLE]

where $\lambda\in\mathbb{R}$ is an arbitrary parameter. This method is weakly convergent with order 1.0 (see [105], [130]). Since $b^{(1)}=b^{(3)}$, equations (3.42) and (3.42) imply that $Q^{(1)}_{1}=\ldots=Q^{(m)}_{1}$ and $P^{(1)}_{1}=\ldots=P^{(m)}_{1}$. Therefore, in general one has to solve the system (3.42a)-(3.42) for the $4N$ variables $Q^{(0)}_{1}$, $P^{(0)}_{1}$, $Q^{(1)}_{1}$, and $P^{(1)}_{1}$. However, for several choices of the parameter $\lambda$ the computational cost can be reduced. If $\lambda=0$, then one can first solve the $2N$ equations (3.42a)-(3.42) for the $2N$ variables $Q^{(0)}_{1}$, $P^{(0)}_{1}$, and then the $2N$ equations (3.42)-(3.42) for the remaining $2N$ variables $Q^{(1)}_{1}$, $P^{(1)}_{1}$. Moreover, if the Hamiltonians are separable, that is, $H(q,p)=T_{0}(p)+U_{0}(q)$ and $h_{i}(q,p)=T_{i}(p)+U_{i}(q)$, then equation (3.42a) can be substituted into equation (3.42), and equation (3.42) can be substitted into equation (3.42), thus reducing the complexity to solving two systems of $N$ equations each. A similar situation occurs for $\lambda=1$. For $\lambda=\frac{1}{2}$ we further have $Q^{(0)}_{1}=Q^{(1)}_{1}=(q_{k}+q_{k+1})/2$ and $P^{(0)}_{1}=P^{(1)}_{1}=(p_{k}+p_{k+1})/2$, and the $SRKw1$ method takes the form of the stochastic midpoint method (1) with $\Delta W^{i}$ replaced by $\hat{I}_{i}$. 2. 2.

SRKw2

For systems driven by a single noise ($m=1$) we can consider methods with $b^{(3)}\equiv 0$. The family of 4-stage $SRKw2$ methods is defined by the tableau

[TABLE]

where $\lambda_{1},\lambda_{2},\lambda_{3}\in\mathbb{R}$ are arbitrary parameters. This method is weakly convergent with order 2.0 (see [105], [130]). Note that $\beta_{3}=\beta_{4}=0$, so the values of the internal stages $Q^{(1)}_{3}$, $Q^{(1)}_{4}$, $P^{(1)}_{3}$, and $P^{(1)}_{4}$ are not needed in (3.42e) and (3.42) to calculate $q_{k+1}$ and $p_{k+1}$, respectively. Moreover, equations (3.42) and (3.42) for $i=3,4$ are explicit updates, therefore there is no need to solve for or calculate the values of these internal stages. In fact, the choice of the parameters $\lambda_{1}$, $\lambda_{2}$, and $\lambda_{3}$ has no effect on the values of $q_{k+1}$ and $p_{k+1}$, therefore we can set them to zero for convenience. Consequently, the system of equations (3.42a) and (3.42) for $i=1,2,3,4$, and equations (3.42) and (3.42) for $i=1,2$ ($12N$ equations) has to be solved for the internal stages $Q^{(0)}_{1},\ldots,Q^{(0)}_{4}$, $P^{(0)}_{1},\ldots,P^{(0)}_{4}$, $Q^{(1)}_{1}$, $Q^{(1)}_{2}$, $P^{(1)}_{1}$, and $P^{(1)}_{2}$ ($12N$ variables). If the Hamiltonians are separable, then equations (3.42a) and (3.42) can be substituted into equations (3.42) and (3.42), and the resulting system of $6N$ equations can be solved for $P^{(0)}_{1},\ldots,P^{(0)}_{4}$, $P^{(1)}_{1}$, and $P^{(1)}_{2}$ ($6N$ variables).

3.5 Quasi-symplecticity

The idea of quasi-symplectic integrators has been proposed in [94] as an attempt to construct numerical methods that at least to some extent emulate the special time evolution of the symplectic and volume forms, as pointed out in Theorem 2.4 and Theorem 2.5, respectively. The authors considered a special form of the stochastic forced Hamiltonian system, namely

[TABLE]

where $M$ is an $N\times N$ constant positive definite matrix, $\Gamma$ is an $N\times N$ constant matrix, and $\sigma_{i}$ are constant vectors. The authors call a numerical integrator $F^{+}_{t_{k+1},t_{k}}:(q_{k},p_{k})\longrightarrow(q_{k+1},p_{k+1})$ quasi-symplectic if it satisfies the following two conditions when applied to the system (3.5):

• (QS1)

it degenerates to a symplectic method when the forcing term vanishes, i.e., $\Gamma=0$,

• (QS2)

the Jacobian

[TABLE]

does not depend on $q_{k}$ and $p_{k}$.

The condition (QS2) is natural, since the exact Jacobian (2.43) does not depend on the phase space variables. Several quasi-symplectic numerical methods have been proposed and tested in [94]; see also [90]. Below we demonstrate that the idea of quasi-symplecticity can be extended to more general systems than (3.5).

The methods presented in Section 3.3.3 and Section 3.4.3 preserve the underlying variational structure of the general system (1), as has been shown in Theorem 3.1. These methods also naturally reduce to symplectic methods, when the forcing terms $F$ and $f_{i}$ vanish (see [50], [80], [81], [93], [130]). Below we show that the Störmer-Verlet method satisfies the condition (QS2) for a much broader class of systems than (3.5).

Theorem 3.9**.**

Suppose that $H(q,p)$, $F(q,p)$, and $h_{i}(q,p)$, $f_{i}(q,p)$ for $i=1,\ldots,m$ satisfy conditions (H1)-(H3). If the Hamiltonians are separable, that is,

[TABLE]

and the forcing terms have the form

[TABLE]

for constant $N\times N$ matrices $\Gamma_{i}$, then the Jacobian $J$ of the discrete flow $F^{+}_{t_{k+1},t_{k}}:(q_{k},p_{k})\longrightarrow(q_{k+1},p_{k+1})$ defined by the Störmer-Verlet method (2) does not depend on $q_{k}$ and $p_{k}$, and is almost surely equal to

[TABLE]

where $I$ is the $N\times N$ identity matrix, $\gamma=\Delta t\Gamma_{0}+\sum_{i=1}^{m}\Delta W^{i}\Gamma_{i}$, and we assume that the matrix $I-\frac{1}{2}\gamma$ is almost surely invertible.

Proof.

With the separable Hamiltonians (3.52) and the linear forcing terms (3.53), the first equation in (2) is linear, and $P_{1}$ can be expressed as

[TABLE]

We then plug in $P_{1}$ into the second and third equations in (2) to obtain expressions for $q_{k+1}$ and $p_{k+1}$ as functions of $q_{k}$ and $p_{k}$. Let us introduce the notation

[TABLE]

Using this notation, the Jacobian $J$ of the mapping $(q_{k},p_{k})\longrightarrow(q_{k+1},p_{k+1})$ can be expressed as

[TABLE]

Let us transform this determinant into a block upper triangular form by performing basic linear manipulations on its columns and rows. First, multiply the upper and lower right blocks by $\frac{1}{2}B$ on the right, and add the results to the upper and lower left blocks, respectively. Then, multiply the upper left and right blocks by $\frac{1}{2}C$ on the left, and add the results to the lower left and right blocks, thus obtaining a block upper triangular form. Writing out these steps explicitly, we have

[TABLE]

which completes the proof. ∎

Remark.

In case the matrix $\eta=I-\frac{1}{2}\gamma$ is not almost surely invertible, one can replace $\Delta W^{i}$ with the suitably chosen truncated increments (3.41).

4 Numerical experiments

In this section we present the results of our numerical experiments. We have tested the performance of the stochastic Lagrange-d’Alembert integrators presented in Section 3, namely the midpoint method (1), the Störmer-Verlet method (2), the DIRK method (3.40) with $\lambda=1/2$, the $SRKw1$ method (3.48) with $\lambda=0$, and the $SRKw2$ method (3.49), and compared it to the performance of some popular general purpose non-geometric explicit stochastic integrators, namely the mean-square Heun method ([24], [62]), the mean-square $R2$ and $E1$ methods (see [19], [20], [21], [24]), and the weak $RS1$ and $RS2$ methods ([105]). The Lagrange-d’Alembert integrators have demonstrated superior behavior in long-time simulations in all of the examples described below. In the case of the midpoint, Störmer-Verlet, and DIRK methods, we used unbounded increments $\Delta W^{i}$, but observed no numerical issues. In principle, one should use truncated increments of the form (3.41), but for the chosen parameters in the examples below, the probability of encountering a singularity was negligible. All computations have been performed in the Julia programming language with the help of the GeometricIntegrators.jl library (see [67]).

4.1 Long-time energy behavior

The Kubo oscillator is a stochastic Hamiltonian system with the Hamiltonians given by $H(q,p)=p^{2}/2+q^{2}/2$ and $h(q,p)=\beta(p^{2}/2+q^{2}/2)$, where $\beta$ is the noise intensity (see [93]). It is an example of an oscillator with a fluctuating frequency and it was first introduced in the context of the line-shape theory (see [4], [68]), but later also found many other applications in connection with mechanical systems, turbulence, laser theory, wave propagation (see [124] and the references therein), magnetic resonance spectroscopy, nonlinear spectroscopy (see [96] and the references therein), single molecule spectroscopy ([58]), and stochastic resonance ([26], [27], [28], [38]). The Kubo oscillator serves as a prototype for multiplicative stochastic processes, and since its solutions can be calculated analytically, it is often used for validation of numerical algorithms (see, e.g., [36], [81], [93], [118]). Here we consider the damped Kubo oscillator with the forcing terms given by $F(q,p)=-\nu p$ and $f(q,p)=-\beta\nu p$, where $\nu$ is the damping coefficient. It is straightforward to verify that the exact solution is given by

[TABLE]

where $q_{0}$ and $p_{0}$ are the initial conditions, the angular frequency is $\omega=\frac{1}{2}\sqrt{4-\nu^{2}}$, and we have assumed the underdamped case $0\leq\nu<2$. Note that (4.1) is the solution of the deterministic damped harmonic oscillator with the time argument shifted by $\beta W(t)$. Given that $W(t)\sim N(0,t)$ is normally distributed, one can explicitly calculate the expected value of the Hamiltonian $H$ as a function of time as

[TABLE]

where

[TABLE]

Simulations with the initial conditions $q_{0}=2$, $p_{0}=0$, and the parameters $\beta=0.5$ and $\nu=0.001$ were carried out until the time $T=5000$ (approximately 800 periods of the oscillator in the absence of noise). In each case 50000 sample paths were generated. The numerical value of the mean Hamiltonian $E(H)$ as a function of time is depicted in Figure 4.1 and Figure 4.2 for the mean-square and weak integrators, respectively. We see that the Lagrange-d’Alembert integrators capture the exponential decay of $E(H)$ very accurately even when relatively large time steps $\Delta t$ are used. The explicit Heun and $R2$ methods fail to reproduce that behavior even for the significantly smaller time step. While the explicit $E1$, $RS1$, and $RS2$ methods capture the qualitative decay of $E(H)$, still much smaller time steps would be needed to reach the level of accuracy of the Lagrange-d’Alembert integrators, thus rendering them inefficient. The accuracy of the Monte Carlo approximation of $E(H)$ at each time step was controlled by estimating the relative error $\sigma(E(H))/E(H)$, where $\sigma(E(H))$ denotes the standard deviation of the mean. The maximum relative error for the Störmer-Verlet method was $2.87\cdot 10^{-3}$, and for all other methods it did not exceed $5.26\cdot 10^{-4}$.

4.2 Ergodic limits

In many cases of practical interest the system (1) is ergodic, which means that

• (1)

it possesses a unique invariant measure represented by the probability density function $\rho_{\infty}(\xi,\zeta)$ with $(\xi,\zeta)\in T^{*}Q$, i.e. a stationary solution of the corresponding Fokker-Planck equation (see [37])

• (2)

for any function $\varphi:T^{*}Q\longrightarrow\mathbb{R}$ with polynomial growth at infinity, its ergodic limit, i.e. the expected value with respect to the invariant measure, can be calculated as the limit

[TABLE]

where $(\bar{q}(t),\bar{p}(t))$ is an arbitrary solution of (1) with arbitrary initial conditions.

For more information about ergodic systems and ergodic numerical schemes see, e.g., [14], [17], [51], [85], [86], [90], [119]. For many applications, it is interesting to compute the mean of a given function with respect to the invariant law of the diffusion, but the explicit form of the invariant measure is often not known. If the considered system is ergodic, then the ergodic limit can be approximated as

[TABLE]

by choosing a sufficiently large time $T$. One can then use numerical integrators to approximate $\bar{q}(T)$ and $\bar{p}(T)$. However, formula (4.5) requires integration of the system over comparatively long time intervals, which poses a significant computational difficulty. Below we compare the performance of the geometric integrators introduced in Section 3 with the performance of explicit schemes. Note that we do not make any claims about the ergodicity of the used schemes and defer this issue to future work.

In recent years the analysis of nonlinear oscillators subjected to random excitations has been of significant interest, for instance in the context of stochastic resonance and stochastic bifurcation theory. The van der Pol oscillator is one of the most extensively studied systems in nonlinear dynamics and has a long history of being used in physical and biological sciences (see, e.g., [40]). It possesses a trivial fixed point and a limit cycle attractor. Various stochastic extensions of the van der Pol oscillator have been considered to test the effect of external noises on its self-sustaining mechanism, the period and lifetime of its oscillations, and the attraction basins of its fixed point and limit cycle (see, e.g., [34], [53], [76], [78], [79], [109], [114], [120]). A numerical study of such stochastic extensions requires long integration times and serves as an interesting testbed for numerical algorithms. Consider van der Pol’s equation with additive noise (see [90]), which is a stochastic forced Hamiltonian system of the form (1) with

[TABLE]

where $\nu\geq 0$ and $\sigma\geq 0$ are parameters. The explicit form of the invariant measure for this system is unknown, however, it is interesting to compute the ergodic value of the energy. Note that the forcing term $F(q,p)$ is not globally Lipschitz, therefore this example also tests the Lagrange-d’Alembert integrators in the situation when the assumption (H3) from Section 2.1 is not satisfied. Simulations with the initial conditions $q_{0}=1$, $p_{0}=1$, and the parameters $\sigma=0.05$ and $\nu=0.001$ were carried out until the time $T=5000$. In each case $10^{6}$ sample paths were generated. The numerical value of the mean Hamiltonian $E(H)$ as a function of time is depicted in Figure 4.3 for the DIRK, Heun, and $E1$ methods. As the reference value we take $H^{\text{erg}}=2.3165$, which was calculated in [90] using a second-order weak quasi-symplectic method at the time $T_{\text{ref}}=10000$ with the time step $\Delta t=0.05$ and $4\times 10^{6}$ sample paths. We see that the DIRK method accurately reproduces the reference value even with the relatively large time step $\Delta t=0.2$, while the Heun and $E1$ methods require the much smaller time step $\Delta t=0.02$ to reach that level of accuracy. The situation is similar for the other Lagrange-d’Alembert and explicit integrators. Figure 4.4 depicts the behavior of $E(H)$ near the reference value on the time interval $[4000,5000]$ for each of the tested integrators. The maximum relative Monte Carlo error $\sigma(E(H))/E(H)$ did not exceed $7.17\cdot 10^{-4}$ in any of the simulations.

4.3 Vlasov equation

In recent years there has been a growing interest in applying geometric integration to particle-in-cell (PIC) simulations of the Vlasov equation in plasma physics. The results to date concern almost entirely collisionless cases (see [18], [35], [45], [101], [110], [115], [116], [131]). The first step towards a geometric description of collision operators, using the so-called metriplectic formulation, has been recently made in [46]. Below we demonstrate that stochastic forced Hamiltonian systems provide an alternative structure-preserving description, and further consider two examples, namely the Lenard-Bernstein and the Lorentz collision operators, to test the long-time behavior of the stochastic Lagrange-d’Alembert integrators.

4.3.1 Lenard-Bernstein collision operator

The following two-dimensional Vlasov-Fokker-Planck equation

[TABLE]

has been studied in [61] and [113] as a model for collisional kinetic plasmas, where $\rho=\rho(x,v,t)$ denotes the particle distribution function in the position-velocity phase space, $E(x)=-\phi^{\prime}(x)$ is the external electric field with the electrostatic potential $\phi(x)$, and $\nu>0$, $\mu>0$, $D>0$ are real parameters. The right-hand side of (4.7) is the so-called Lenard-Bernstein collision operator, which models small-angle collisions and was originally used to study longitudinal plasma oscillations (see [73]). A stochastic split particle-in-cell (PIC) method for the numerical simulation of (4.7) has been proposed in [113], whereby the advection part is solved using the standard PIC method, and the diffusion part is modeled by a stochastic differential equation. Below we demonstrate a structure-preserving approach to solving (4.7). When $\rho$ is interpreted as a probability density function, then (4.7) is the Fokker-Planck equation for the two-dimensional stochastic process $(X(t),V(t))$ whose evolution is governed by the stochastic differential equation (see [37], [62])

[TABLE]

driven by the one-dimensional Wiener process $W(t)$. This equation is a stochastic forced Hamiltonian system (1) with

[TABLE]

It can be easily verified that the stationary solution of (4.7) is given by the Gibbs measure

[TABLE]

where $Z$ is the normalizing constant such that $\int\int\rho_{\infty}(x,v)\,dv\,dx=1$. Let us consider (4.7) on the domain $(x,v)\in[0,1]\times\mathbb{R}$ with periodic boundary conditions in $x$, and with the electrostatic potential

[TABLE]

where $E_{0}>0$ is the maximum magnitude of the electric field $E(x)=-\phi^{\prime}(x)$. One can check that the system (4.3.1) with the potential (4.11) is ergodic (see Theorem 3.2 in [85]). As the initial condition, we take the probability distribution of the form

[TABLE]

where $\rho_{X}(x)$ for $\epsilon>0$ describes a perturbation of the uniform distribution along the spatial direction $x$, and $\rho_{V}(v)$ for $a>0$ is the so called bump-on-tail distribution in the velocity space, where the bump is centered at $v_{0}$ with the standard deviation $\sigma>0$. Simulations with the parameters $\nu=0.01$, $\mu=1$, $D=\sqrt{2}$, $E_{0}=3$, $\epsilon=0.25$, $a=0.5$, $v_{0}=4$, and $\sigma=0.5$ were carried out until the time $T=1000$. In each case $10^{7}$ sample paths were generated. The initial conditions $X_{0}$ and $V_{0}$ were randomly drawn from the probability distribution (4.12) using rejection sampling (see Figure 4.5). The exact ergodic value $H^{\text{erg}}$ of the Hamiltonian can be calculated using the invariant probability density (4.10) as

[TABLE]

The numerical value of the mean Hamiltonian $E(H)$ as a function of time is depicted in Figure 4.6 for the DIRK, Heun, and $E1$ methods. We see that the DIRK method accurately reproduces the ergodic limit even with the relatively large time step $\Delta t=0.15$, while the $E1$ method requires the much smaller time step $\Delta t=0.02$ to reach a comparable level of accuracy. The Heun method yields a less accurate result even for $\Delta t=0.02$. The situation is similar for the other Lagrange-d’Alembert and explicit integrators. Figure 4.7 depicts the behavior of $E(H)$ near the exact ergodic limit on the time interval $[500,1000]$ for each of the tested integrators. The maximum relative Monte Carlo error $\sigma(E(H))/E(H)$ did not exceed $8.24\cdot 10^{-4}$ in any of the simulations. The numerical probability density at the final time $T=1000$ calculated with each of the mean-square and weak methods is depicted in comparison to the exact invariant measure (4.10) in Figure 4.8 and Figure 4.9, respectively.

4.3.2 Lorentz collision operator

The following four-dimensional Vlasov-Fokker-Planck equation

[TABLE]

has been used in [11] to study the electron-ion collision effects on the damping of electron plasma waves, where $\rho=\rho(x,y,v_{x},v_{y},t)$ denotes the particle distribution function in the position-velocity phase space, $E_{x}(x,y)=-\frac{\partial\phi}{\partial x}(x,y)$ and $E_{y}(x,y)=-\frac{\partial\phi}{\partial y}(x,y)$ are the components of the external electric field with the electrostatic potential $\phi(x,y)$, and $\nu>0$ is a real parameter. The right-hand side of (4.14) is the so-called Lorentz collision operator, which models electron-ion interactions via pitch-angle scattering. The primary effect of this type of scattering is a change of the direction of the electron’s velocity with negligible energy loss. More information about the Lorentz collision operator can be found in [60]. Below we demonstrate a structure-preserving approach to solving (4.14). When $\rho$ is interpreted as a probability density function, then (4.14) is the Fokker-Planck equation for the four-dimensional stochastic process $(X(t),Y(t),V_{x}(t),V_{y}(t))$ whose evolution is governed by the stochastic differential equation (see [37], [62])

[TABLE]

driven by the one-dimensional Wiener process $W(t)$. This equation is a stochastic forced Hamiltonian system (1) with

[TABLE]

Let us consider (4.14) on the domain $(x,y,v_{x},v_{y})\in[0,1]^{2}\times\mathbb{R}^{2}$ with periodic boundary conditions in $x$ and $y$, and with the electrostatic potential

[TABLE]

where $E_{0}>0$ is the maximum magnitude of the electric field $E(x,y)=-\nabla\phi(x,y)$. As the initial condition, we take the probability distribution of the form

[TABLE]

where the parameters $\epsilon_{1},\epsilon_{2}>0$ describe a perturbation of the uniform distribution along the spatial directions $x$ and $y$, and the velocity part is Maxwellian. The Lorentz collision operator by construction preserves the total energy of the plasma, that is,

[TABLE]

where E(H)\equiv E\Big{(}H\big{(}X(t),Y(t),V_{x}(t),V_{y}(t)\big{)}\Big{)} for short (see [60]). Moreover, in the stochastic description (4.3.2) the Hamiltonian is almost surely preserved for each sample path, which can be easily verified by calculating the stochastic differential

[TABLE]

where the last equality follows from (4.3.2) and (4.3.2).

Mean-square integrators aim to approximate each sample path of the exact solution, and therefore they should also approximate the stronger energy preservation condition (4.20). In order to test the long-time performance of the mean-square integrators discussed in Section 3.3.3, simulations with the parameters $\nu=0.005$, $E_{0}=3$, and $\epsilon_{1}=\epsilon_{2}=0.25$ were carried out for a single sample path until the time $T=100000$. For each integrator the same random initial condition, drawn from the probability distribution (4.18), and the same realization of the Brownian motion were used. The numerical value of the Hamiltonian $H$ as a function of time is depicted in Figure 4.10. Even with relatively large time steps the mean-square Lagrange-d’Alembert integrators preserve energy much more accurately than the non-geometric explicit methods.

On the other hand, weak integrators aim to approximate the probability distribution and functionals of the exact solutions rather than each sample path, therefore they may not preserve energy on each sample path, but nevertheless they should approximate the mean energy (4.19). In order to test the long-time performance of the weak integrators discussed in Section 3.4.3, simulations with the same parameters as above were carried out until the time $T=10000$. In each case $10^{6}$ sample paths were generated. The initial conditions were randomly drawn from the probability distribution (4.18). The exact mean energy can be calculated by substituting (4.3.2), (4.17), and (4.18) into (4.19). For the chosen parameters, we have $E(H)_{\text{exact}}=1$.The numerical value of the mean Hamiltonian $E(H)$ as a function of time is depicted in Figure 4.11 for the $SRKw1$, $SRKw2$, $RS1$, and $RS2$ methods. We see that the weak Lagrange-d’Alembert methods accurately reproduce the mean energy conservation even with the relatively large time steps, while the non-geometric methods require much smaller time steps to reach a comparable level of accuracy. The maximum relative Monte Carlo error $\sigma(E(H))/E(H)$ did not exceed $0.001$ in any of the simulations.

5 Summary and future work

In this paper we have presented a general framework for constructing a new class of stochastic variational integrators for stochastic forced Hamiltonian systems. We have extended the approach taken in [50] by considering the stochastic Lagrange-d’Alembert principle and constructing the corresponding structure-preserving schemes, which we have dubbed stochastic Lagrange-d’Alembert variational integrators. We have shown that in the presence of a symmetry such integrators satisfy a discrete version of Noether’s theorem. We have further considered certain classes of mean-square and weak Runge-Kutta methods previously known in the literature, and determined the conditions under which such methods become Lagrange-d’Alembert integrators. We have finally pointed out several examples of low-stage Runge-Kutta methods of that type, and demonstrated their superior long-time numerical performance via numerical experiments. In particular, as one of the test cases we have considered the Vlasov-Fokker-Planck equation and proposed a new geometric approach to the simulation of collisional kinetic plasmas.

Our work can be extended in several ways. The mean-square partitioned Runge-Kutta methods introduced in Section 3.3 only use the increments $\Delta W^{r}=\int_{t_{k}}^{t_{k+1}}dW^{r}(t)$, therefore their mean-square order of convergence cannot exceed 1.0 (see [21], [92], [93]). To obtain mean-square convergence of higher order one can extend the definitions of the discrete Hamiltonian (3.18) and the discrete forces (3.4) to include higher-order multiple Stratonovich integrals, e.g., to achieve convergence of order 1.5 we would need to include terms involving $\Delta Z^{r}=\int_{t_{k}}^{t_{k+1}}\int_{t_{k}}^{t}dW^{r}(\xi)\,dt$; see [50] for an example how this can be done for unforced Hamiltonian systems. Another aspect worth a more detailed investigation is the issue of ergodicity of the Lagrange-d’Alembert methods. In Section 4.2 and Section 4.3 we have experimentally demonstrated the usefulness of our integrators in calculating the ergodic limits, but have not formally proved their ergodicity. It would be beneficial to determine under what conditions Lagrange-d’Alembert integrators can be ergodic in the sense discussed in, e.g., [85], [86], or [119], when applied to ergodic Hamiltonian systems. It would also be interesting to extend the idea of Lagrange-d’Alembert integrators to stochastic Hamiltonian systems that are both forced and constrained. Structure-preserving numerical methods for such systems would be of great interest in molecular dynamics (see [15], [30], [125]). Yet another direction of great practical significance would be a further study of the geometric approach to collisional kinetic plasmas presented in Section 4.3 and application of more realistic collision operators that preserve the total energy and momentum, as well as an extension to the self-consistent Maxwell-Vlasov equations (see [65], [66]). Finally, one may extend the idea of variational integration to stochastic multisymplectic partial differential equations such as the stochastic Korteweg-de Vries, Camassa-Holm or Hunter-Saxton equations. Theoretical groundwork for such numerical schemes has been recently presented in [49].

Acknowledgements

We would like to thank Christopher Albert, Darryl Holm, Katharina Kormann, Omar Maj, Philip Morrison, Bruce Scott, Cesare Tronci, and Udo von Toussaint for useful comments and references. We are particularly endebted to Eric Sonnendrücker for pointing out the connections between stochastic systems and the Vlasov equation. The study is a contribution to the Reduced Complexity Models grant number ZT-I-0010 funded by the Helmholtz Association of German Research Centers.