Stochastic discrete Hamiltonian variational integrators.

Variational integrators are derived for structure-preserving simulation of stochastic Hamiltonian systems with a certain type of multiplicative noise arising in geometric mechanics. 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 its corresponding variational principle. Our approach permits to recast in a unified framework a number of integrators previously studied in the literature, and presents a general methodology to derive new structure-preserving numerical schemes. The resulting integrators are symplectic; they preserve integrals of motion related to Lie group symmetries; and they include stochastic symplectic Runge-Kutta methods as a special case. Several new low-stage stochastic symplectic methods of mean-square order 1.0 derived using this approach are presented and tested numerically to demonstrate their superior long-time numerical stability and energy behavior compared to nonsymplectic methods.

Introduction

Stochastic differential equations (SDEs) play an important role in modeling dynamical systems subject to internal or external random fluctuations. Standard references include [5, 27–29, 42, 50]. Within this class of problems, we are interested in stochastic Hamiltonian systems, which take the form (see [6, 30, 43])where and are the Hamiltonian functions, W(t) is the standard one-dimensional Wiener process, and denotes Stratonovich integration. The system (1.1) can be formally regarded as a classical Hamiltonian system with the randomized Hamiltonian given by , where H(q, p) is the deterministic Hamiltonian and h(q, p) is another Hamiltonian, to be specified, which multiplies (in the Stratonovich sense, denoted as ) a one-dimensional Gaussian white noise, . 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. Particular examples include modeling synchrotron oscillations of particles in particle storage rings (see [17, 56]) and stochastic dynamics of the interactions of singular solutions of the EPDiff basic fluids equation (see [23]). More examples are discussed in Sect. 4. See also [31, 37, 46, 54, 57, 58, 61].

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 particular, their phase space flows (almost surely) preserve the symplectic structure. 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 [18, 41, 55] and the references therein) and symplectic integrators have been shown to demonstrate superior performance in long-time simulations of Hamiltonian systems, 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 were first proposed in [43, 44]. Stochastic generalizations of symplectic partitioned Runge–Kutta methods were analyzed in [13, 35, 36]. A stochastic generating function approach to constructing stochastic symplectic methods, based on approximately solving a corresponding stochastic Hamilton–Jacobi equation satisfied by the generating function, was proposed in [65, 66], and this idea was further pursued in [2, 4, 16]. Stochastic symplectic integrators constructed via composition methods were proposed and analyzed in [45]. A first order weak symplectic numerical scheme and an extrapolation method were proposed and their global error was analyzed in [3]. More recently, an approach based on Padé approximations has been used to construct stochastic symplectic methods for linear stochastic Hamiltonian systems (see [60]). Higher-order strong and weak symplectic partitioned Runge–Kutta methods have been proposed in [67, 68]. High-order conformal symplectic and ergodic schemes for the stochastic Langevin equation have been introduced in [25]. Other structure-preserving methods for stochastic Hamiltonian systems have also been investigated, see, e.g., [1, 15, 26], and the references therein.

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., [18] and the references therein). To the best of our knowledge, such rigorous analysis has not been attempted in the stochastic context as yet. However, the numerical evidence presented in the papers cited above is promising and suggests that stochastic symplectic integrators indeed possess the property of very accurately capturing the evolution of the Hamiltonian H over exponentially long time intervals (note that the Hamiltonian H in general does not stay constant for stochastic Hamiltonian systems).

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, and in the presence of a symmetry, satisfy a discrete version of Noether’s theorem. For an overview of variational integration for deterministic systems see [40]; see also [21, 32, 33, 47, 48, 53, 63, 64]. Variational integrators were introduced in the context of finite-dimensional mechanical systems, but were later generalized to Lagrangian field theories (see [39]) and applied in many computations, for example in elasticity, electrodynamics, or fluid dynamics; see [34, 49, 59, 62].

Stochastic variational integrators were first introduced in [8] and further studied in [7]. However, those integrators were restricted to the special case when the Hamiltonian function was independent of p, and only low-order Runge–Kutta types of discretization were considered. In the present work we extend the idea of stochastic variational integration to general stochastic Hamiltonian systems by generalizing the variational principle introduced in [33] and applying a Galerkin type of discretization (see [32, 33, 40, 47, 48]), which leads to a more general class of stochastic symplectic integrators than those presented in [7, 8, 35, 36]. Our approach consists in approximating a generating function for the stochastic flow of the Hamiltonian system, but unlike in [65, 66], we do make this discrete approximation by exploiting its variational characterization, rather than solving the corresponding Hamilton–Jacobi equation.

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

In Sect. 2 we introduce a stochastic variational principle and a stochastic generating function suitable for considering stochastic Hamiltonian systems, and we discuss their properties.

In Sect. 3 we present a general framework for constructing stochastic Galerkin variational integrators, prove the symplecticity and conservation properties of such integrators, show they contain the stochastic symplectic Runge–Kutta methods of [35, 36] as a special case, and finally present several particularly interesting examples of new low-stage stochastic symplectic integrators of mean-square order 1.0 derived with our general methodology.

In Sect. 4 we present the results of our numerical tests, which verify the theoretical convergence rates and the excellent long-time performance of our integrators compared to some popular non-symplectic methods.

Section 5 contains the summary of our work.

Variational principle for stochastic Hamiltonian systems

The stochastic variational integrators proposed in [7, 8] 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 are applicable to (1.1) only when the Hamiltonian function is independent of p and the Hamiltonian H is non-degenerate (i.e., the associated Legendre transform is invertible). However, in the case of general the paths q(t) of the system become almost surely nowhere differentiable, which poses a difficulty in interpreting the meaning of the corresponding Lagrangian. Therefore, we need a different sort of action functional and variational principle to construct stochastic symplectic integrators for (1.1). To this end, we will generalize the approach taken in [33]. 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.1).

Stochastic variational principle

Let the Hamiltonian functions and be defined on the cotangent bundle of the configuration manifold Q, and let (q, p) denote the canonical coordinates on . For simplicity, in this work we assume that the configuration manifold has a vector space structure, , so that and . In this case, the natural pairing between one-forms and vectors can be identified with the scalar product on , that is, , where denotes the coordinates on TQ . Let be the probability space with the filtration , and let W(t) denote a standard one-dimensional Wiener process on that probability space (such that W(t) is -measurable). We will assume that the Hamiltonian functions H and h are sufficiently smooth and satisfy all the necessary conditions for the existence and uniqueness of solutions to (1.1), and their extendability to a given time interval with . One possible set of such assumptions can be formulated by considering the Itô form of (1.1),with andwhere , , and denote the Hessian matrices of h. For simplicity and clarity of the exposition, throughout this paper we assume that (see [5, 27–29])These assumptions are sufficient1 for our purposes, but could be relaxed if necessary. Define the spaceSince we assume , the space is a vector space (see [50]). Therefore, we can identify the tangent space . We can now define the following stochastic action functional, ,where denotes Stratonovich integration, and we have omitted writing the elementary events as arguments of functions, following the standard convention in stochastic analysis.

H and h are functions of their arguments

A and B are globally Lipschitz

Theorem 2.1

(Stochastic variational principle in phase space) Suppose that H(q, p) and h(q, p) satisfy conditions (H1)–(H2). If the curve in satisfies the stochastic Hamiltonian system (1.1) for , where , then the pair is a critical point of the stochastic action functional (2.4), that is,almost surely for all variations such that almost surely and .

Proof

Let the curve in satisfy (1.1) for . It then follows that the stochastic processes q(t) and p(t) are almost surely continuous, -adapted semimartingales, that is, (see [5, 50]). We calculate the variation (2.5) aswhere we have used the end point condition, . Since the Hamiltonians are 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 and integration with respect to t and W(t). Upon applying the integration by parts formula for semimartingales (see [50]), we findSubstituting and rearranging terms produces,where we have used . Since satisfy (1.1), then by definition we have that almost surely for all ,that is, q(t) can be represented as the sum of two semi-martingales and , where the sample paths of the process are almost surely continuously differentiable. Let us calculatewhere in the last equality we have used the standard property of the Riemann–Stieltjes integral for the first term, as is almost surely differentiable, and the associativity property of the Stratonovich integral for the second term (see [27, 50]). Substituting (2.10) in (2.8), we show that the second term is equal to zero. By a similar argument we also prove that the first term in (2.8) is zero. Therefore, , almost surely.

Remark It is natural to expect that the converse theorem, that is, if is a critical point of the stochastic action functional (2.4), then the curve satisfies (1.1), should also hold, although a larger class of variations may be necessary. A variant of such a theorem, although for a slightly different variational principle and in a different setting, was proved in Lázaro-Camí and Ortega [30]. Another variant for Lagrangian systems was proved by Bou-Rabee and Owhadi [8] in the special case when is independent of p. In that case, one can assume that q(t) is continuously differentiable. In the general case, however, q(t) is not differentiable, and the ideas of [8] cannot be applied directly. We leave this as an open question. Here, we will use the action functional (2.4) and the variational principle (2.5) to construct numerical schemes, and we will directly verify that these numerical schemes converge to solutions of (1.1).

Stochastic type-II generating function

When the Hamiltonian functions H(q, p) and h(q, p) satisfy standard measurability and regularity conditions [e.g., (H1)–(H2)], then the system (1.1) possesses a pathwise unique stochastic flow . It can be proved that for fixed this flow is mean-square differentiable with respect to the q, p arguments, and is also almost surely a diffeomorphism (see [5, 27–29]). Moreover, almost surely preserves the canonical symplectic form (see [6, 30, 44]), that is,where denotes the pull-back by the flow . We will show below that the action functional (2.4) can be used to construct a type II generating function for . Let be a particular solution of (1.1) on . Suppose that for almost all there is an open neighborhood of , an open neighborhood of , and an open neighborhood of the curve such that for all and there exists a pathwise unique solution of (1.1) which satisfies , , and for . (As in the deterministic case, for sufficiently close to one can argue that such neighborhoods exist; see [38].) Define the function aswhere the domain is given by . Below we prove that S generates2 the stochastic flow .

Theorem 2.2

The function is a type-II stochastic generating function for the stochastic mapping , that is, is implicitly given by the equationswhere the derivatives are understood in the mean-square sense.

Under appropriate regularity assumptions on the Hamiltonians [e.g., (H1)–(H2)], the solutions and are mean-square differentiable with respect to the parameters and , and the partial derivatives are semimartingales (see [5]). We calculate the derivative of S aswhere for notational convenience we have omitted writing and explicitly as arguments of and . Applying the integration by parts formula for semimartingales (see [50]), we findSubstituting and rearranging terms, we obtain the result,since is a solution of (1.1). Similarly we show . By definition of the flow, then .

We can consider as a function of time if we let vary. Let us denote this function as . Below we show that satisfies a certain stochastic partial differential equation, which is a stochastic generalization of the Hamilton–Jacobi equation considered in [33].

Proposition 2.1

(Type II stochastic Hamilton–Jacobi equation) Let the time-dependent type-II generating function be defined aswhere and as before. Then the function satisfies the following stochastic partial differential equationwhere denotes the stochastic differential of with respect to the t variable.

Choose an arbitrary pair and define the particular solution . Form the function and consider its total stochastic differential3 with respect to time. On one hand, the rules of Stratonovich calculus implywhere denotes the partial stochastic differential of with respect to the t argument. On the other hand, integration by parts in (2.17) impliesComparing (2.19) and (2.20), and using (2.13), we obtainThis equation is satisfied along a particular path , however, as in the discussion preceding Theorem 2.2, we can argue, slightly informally, that for almost all , and for t sufficiently close to , one can find open neighborhoods and which can be connected by the flow , i.e., given and , there exists a pathwise unique solution such that and . With this assumption, (2.21) can be reformulated as the full-blown stochastic PDE (2.18).

Remark Similar stochastic Hamilton–Jacobi equations were introduced in [65, 66], where they were used for constructing stochastic symplectic integrators by considering series expansions of generating functions in terms of multiple Stratonovich integrals. This was a direct generalization of a similar technique for deterministic Hamiltonian systems (see [18]). In this work we explore the generalized Galerkin framework for constructing approximations of the generating function in (2.13) by using its variational characterization (2.12). Our approach is a stochastic generalization of the techniques proposed in [33, 48] for deterministic Hamiltonian and Lagrangian systems.

Stochastic Noether’s theorem

Let a Lie group G act on Q by the left action . The Lie group G then acts on TQ and by the tangent and cotangent lift actions, respectively, given in coordinates by the formulas (see [22, 38])where and summation is implied over repeated indices. Let denote the Lie algebra of G and the exponential map (see [22, 38]). Each element defines the infinitesimal generators , , and , which are vector fields on Q, TQ, and , respectively, given byThe momentum map associated with the action is defined as the mapping such that for all the function is the Hamiltonian for the infinitesimal generator , i.e.,where . The momentum map J can be explicitly expressed as (see [22, 38])Noether’s theorem for deterministic Hamiltonian systems relates symmetries of the Hamiltonian to quantities preserved by the flow of the system. It turns out that this result carries over to the stochastic case, as well. A stochastic version of Noether’s theorem was proved in [6, 30]. For completeness, and for the benefit of the reader, below we state and provide a less involved proof of Noether’s theorem for stochastic Hamiltonian systems.

Theorem 2.3

(Stochastic Noether’s theorem) Suppose that the Hamiltonians and are invariant with respect to the cotangent lift action of the Lie group G, that is,for all . Then the cotangent lift momentum map associated with is almost surely preserved along the solutions of the stochastic Hamiltonian system (1.1).

Equation (2.26) implies that the Hamiltonians are infinitesimally invariant with respect to the action of G, that is, for all we havewhere dH and dh denote differentials with respect to the variables q and p. Let (q(t), p(t)) be a solution of (1.1) and consider the stochastic process , where is arbitrary. Using the rules of Stratonovich calculus we can calculate the stochastic differentialwhere we used (1.1), (2.24), and (2.27). Therefore, almost surely for all , which completes the proof.

Stochastic Galerkin Hamiltonian variational integrators

If the converse to Theorem 2.1 holds, then the generating function defined in (2.12) could be equivalently characterized bywhere the extremum is taken pointwise in the probability space . This characterization allows us to construct stochastic Galerkin variational integrators by choosing a finite dimensional subspace of and a quadrature rule for approximating the integrals in the action functional . Galerkin variational integrators for deterministic systems were first introduced in [40], and further developed in [21, 32, 33, 47, 48]. In the remainder of the paper, we will generalize these ideas to the stochastic case.

Construction of the integrator

Suppose we would like to solve (1.1) on the interval [0, T] with the initial conditions . Consider the discrete set of times for , where is the time step. In order to determine the discrete curve that approximates the exact solution of (1.1) at times we need to construct an approximation of the exact stochastic flow on each interval , so that . Let us consider the space defined asFor convenience, we will express q(t) in terms of Lagrange polynomials. Consider the control points and let the corresponding Lagrange polynomials of degree s be denoted by , that is, . A polynomial trajectory can then be expressed aswhere for are the control values, denotes the time derivative of , and denotes the derivative of the Lagrange polynomial with respect to its argument. The restriction of the action functional (2.4) to the space takes the formsince for differentiable functions , where for brevity . Next we approximate the integrals in (3.4) using numerical quadrature rules and , where are the quadrature points, and , are the corresponding weights. At this point we only make a general assumption that for each i we have or . More specific examples will be presented in Sect. 3.4. The approximate action functional takes the formwhere is the increment of the Wiener process over the considered time interval and is a Gaussian random variable with zero mean and variance . The way of approximating the Stratonovich integral above was inspired by the ideas presented in [8, 12, 36, 43, 44]. Note that since we only used in the above approximation, we can in general expect mean-square convergence of order 1.0 only. 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 (see [12, 43, 44]). We can now approximate the generating function with the discrete Hamiltonian function defined aswhere we denoted . The numerical scheme is now implicitly generated by as in (2.13):Equations (3.6) and (3.7) can be written together as the following system:

where in (3.8b), in (3.8d), and for brevity we have introduced the notationEquation (3.8a) corresponds to the second equation in (3.7), Eqs. (3.8b), (3.8c) and (3.8d) correspond to extremizing (3.6) with respect to , and , respectively, and finally (3.8e) is the first equation in (3.7). Knowing , the system (3.8) allows us to solve for : we first simultaneously solve (3.8a), (3.8b) and (3.8d) [ equations] for and [ unknowns]; then and (3.8c) is an explicit update rule for . When , then (3.8) reduces to the deterministic Galerkin variational integrator discussed in [48]. Note that depending on the choice of the Hamiltonians and quadrature rules, the system (3.8) may be explicit, but in the general case it is implicit (see Sect. 3.4). One can use the Implicit Function Theorem to show that for sufficiently small and it will have a solution. However, since the increments are unbounded, for some values of solutions might not exist. To avoid problems with numerical implementations, if necessary, one can replace in (3.8) with the truncated random variable defined aswhere is suitably chosen for the considered problem. See [14, 44] for more details regarding schemes with truncated random increments and their convergence. Alternatively, one could employ the techniques presented in, e.g., [51, 52, 67], where the unbounded increments have been replaced with discrete random variables.

Although the scheme (3.8) formally appears to be a straightforward generalization of its deterministic counterpart, it should be emphasized that the main difference lies in the fact that the increments are random variables such that , which makes the convergence analysis significantly more challenging than in the deterministic case. The main difficulty is in the choice of the parameters s, r, , , , so that the resulting numerical scheme converges in some sense to the solutions of (1.1). The number of parameters and order conditions grows rapidly, when terms approximating multiple Stratonovich integrals are added (see Sect. 3.6 and [10–12, 14]). In Sects. 3.2 and 3.3 we study the geometric properties of the family of schemes described by (3.8), whereas in Sects. 3.4 and 3.5 we provide concrete choices of the coefficients that lead to convergent methods.

Properties of stochastic Galerkin variational integrators

The key features of variational integrators are their symplecticity and exact preservation of the discrete counterparts of conserved quantities (momentum maps) related to the symmetries of the Lagrangian or Hamiltonian (see [40]). These properties carry over to the stochastic case, as was first demonstrated in [8] for Lagrangian systems. In what follows, we will show that the stochastic Galerkin Hamiltonian variational integrators constructed in Sect. 3.1 also possess these properties.

Theorem 3.1

(Symplecticity of the discrete flow) Let be the dicrete stochastic flow implicitly defined by the discrete Hamiltonian as in (3.7). Then is almost surely symplectic, that is,where is the canonical symplectic form on .

The proof follows immediately by observing that (see [33])where d in the above formula denotes the differential operator with respect to the variables q and p and is understood in the mean-square sense. The result holds almost surely, because Eq. (3.7) holds almost surely.

The discrete counterpart of stochastic Noether’s theorem readily generalizes from the corresponding theorem in the deterministic case.

Theorem 3.2

(Discrete stochastic Noether’s theorem) Let be the cotangent lift action of the action of the Lie group G on the configuration space Q. If the generalized discrete stochastic Lagrangian , where , is invariant under the action of G, that is,where is the projection onto , then the cotangent lift momentum map J associated with is almost surely preserved, i.e., a.s. .

See the proof of Theorem 4 in [33]. In our case the result holds almost surely, because Eq. (3.7) holds almost surely.

For applications, it is useful to know under what conditions the discrete Hamiltonian (3.6) inherits the symmetry properties of the Hamiltonians H and h. Not unexpectedly, this depends on the behavior of the interpolating polynomial (3.3) under the group action. We say that the polynomial is equivariant with respect to G if for all we have

Theorem 3.3

Suppose that the Hamiltonians and are invariant with respect to the cotangent lift action of the Lie group G, that is,for all , and suppose the interpolating polynomial is equivariant with respect to G. Then the generalized discrete stochastic Lagrangian corresponding to the discrete Hamiltonian (3.6), where , is invariant with respect to the action of G.

The proof is similar to the proofs of Lemma 3 in [33] and Theorem 3 in [48]. Let us, however, carefully examine the actions of G on Q, TQ, and . Let and , and let . First, note that for the stochastic discrete Hamiltonian (3.6), we havewhere we used (3.8e). Consider and for , and calculate (3.15) for the transformed values and :Let us perform a change of variables with respect to which we extremize. First, define , so that for . From (3.13) we have , which we use to define by for . Since and are bijective, extremization with respect to and is equivalent to extremization with respect to and , and implies . Moreover, from (3.13) and (2.22) we have that . Finally, the invariance of the Hamiltonians implieswhich completes the proof.

Remark One can easily verify that the interpolating polynomial (3.3) is in particular equivariant with respect to linear actions (translations, rotations, etc.), therefore the stochastic Galerkin variational integrator (3.8) preserves quadratic momentum maps (such as linear and angular momentum) related to linear symmetries of the Hamiltonians H and h.

Stochastic symplectic partitioned Runge–Kutta methods

A general class of stochastic Runge–Kutta methods for Stratonovich ordinary differential equations was introduced and analyzed in [10-12]. These ideas were later used by Ma et al. [35] and Ma and Ding [36] to construct symplectic Runge–Kutta methods for stochastic Hamiltonian systems. An s-stage stochastic symplectic partitioned Runge–Kutta method for (1.1) is defined in [36] by the following system:

where and for are the position and momentum internal stages, respectively, and the coefficients of the method , , , , , satisfy the symplectic conditions

for . We now prove that in the special case when , the stochastic Galerkin variational integrator (3.8) is equivalent to the stochastic symplectic partitioned Runge–Kutta method (3.18).

Theorem 3.4

Let and let for denote the Lagrange polynomials of degree associated with the quadrature points . Moreover, let the weights be given byand assume for . Then the stochastic Galerkin Hamiltonian variational integrator (3.8) is equivalent to the stochastic partitioned Runge–Kutta method (3.18) with the coefficients

for .

The proof follows the main steps of the proof of Theorem 2.6.2 in [40]. The time derivative is a polynomial of degree . Therefore, it can be uniquely expressed in terms of the Lagrange polynomials asUpon integrating with respect to time, we findwhere we have used . For this giveswhere we have used and (3.20). Define the internal stages . Then, upon using (3.8d), Eq. (3.24) becomes (3.18c). For Eq. (3.23) giveswhere is defined by (3.21a). Upon substituting (3.8d), Eq. (3.25) becomes (3.18a), where is defined by (3.21c). Next, sum Eqs. (3.8a), (3.8b), and (3.8c). Noting that , this gives Eq. (3.18d). Finally, we note that for each we have a unique decompositionsince the left-hand side is a polynomial of degree s, and therefore it can be uniquely expressed as a linear combination of the Lagrange polynomials with the coefficients . Evaluating this identity at , , and differentiating it with respect to yield the following three equations, respectively,We form a linear combination of Eqs. (3.8a), (3.8b) and (3.8c) with the coefficients , , and , respectively. By using the identities (3.27) and rearranging the terms, we obtain (3.18d), where and are defined by (3.21b) and (3.21d), respectively. One can easily verify that the coefficients (3.21) satisfy the conditions (3.19).

Examples

In the construction of the integrator (3.8) we may choose the degree s of the approximating polynomials and the quadrature rules and . In the deterministic case, the higher the degree of the polynomials and the higher the order of the quadrature rule, then the higher the order of convergence of the resulting integrator (see [48]). In our case, however, as explained in Sect. 3.1, we cannot in general achieve mean-square order of convergence higher than 1.0, because we only used in (3.5). Since the system (3.8) requires solving equations for variables, from the computational point of view it makes sense to only consider methods with low values of s and r. In this work we focus on the following classical numerical integration formulas (see [18-20]):In [48] notation PsNrQu was proposed to denote a Galerkin variational integrator based on polynomials of degree s and a quadrature rule of order u with r quadrature points. We adopt a similar notation, keeping in mind that u denotes the classical order of the used quadrature rule—when the rule is applied to a stochastic integral, as in (3.5), its classical order is not attained in general. We also use a three-letter code to identify which integration formula is used. For example, P2N2Q4Gau denotes the integrator defined by (3.8) using polynomials of degree 2 and the Gauss–Legendre quadrature formula of classical order 4 with 2 quadrature points for both the Lebesgue and Stratonovich integrals in (3.5). If two different quadrature rules are used, we first write the rule applied to the Lebesgue integral, followed by the rule applied to the Stratonovich integral, e.g., P1N1Q2GauN2Q2Lob. Below we give several examples of integrators obtained by using polynomials of degree and the quadrature rules listed above.

Gauss–Legendre quadratures (Gau): midpoint rule (), etc.

Lobatto quadratures (Lob): trapezoidal rule (), Simpson’s rule (), etc.

Open trapezoidal rule (Otr; )

Milne’s rule (Mil; )

Rectangle rule (Rec; )—only in the case when .

General Hamiltonian function h(q, p)

For a general Hamiltonian , Eq. (3.8d), which represents the discretization of the Legendre transform, needs to contain both and terms to correctly approximate the continuous system. Therefore, we only consider methods with for all . A few examples of interest are listed below.

P1N1Q2Gau (Stochastic midpoint method)

Using the midpoint rule (, , ) together with polynomials of degree gives a stochastic Runge–Kutta method (3.18) with . Noting that and , this method can be written as

The stochastic midpoint method was considered in [36, 44]. 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, and , then from the second equation can be substituted into the first one. In that case only N nonlinear equations have to be solved for .

P2N2Q2Lob (Stochastic Störmer–Verlet method)

If the trapezoidal rule (, , , , ) is used with polynomials of degree , we obtain another partitioned Runge–Kutta method (3.18) with , , , , , . Noting that , , and , this method can be more efficiently written as

This method is a stochastic generalization of the Störmer–Verlet method (see [18]) and was considered in [36]. 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, this method becomes fully explicit.

P1N2Q2Lob (Stochastic trapezoidal method)

This integrator is based on polynomials of degree with control points , , and the trapezoidal rule. Equations (3.8) take the form

This integrator is a stochastic generalization of the trapezoidal method for deterministic systems (see [40]). One may easily verify that if the Hamiltonians are separable, that is, and , then and (3.30) is equivalent to the Störmer–Verlet method (3.29) and is fully explicit.

P1N3Q4Lob

If we use Simpson’s rule (, , , , , , , ), the resulting integrator (3.8) requires solving simultaneously 4N nonlinear equations, so it is computationally expensive in general. However, if the Hamiltonians H and h are separable, then (3.8d) implies , and the integrator can be rewritten as

where

and and . In this case only the first equation needs to be solved for , and then the second equation is an explicit update.

P1N2Q2Otr

Like the method (3.30), this integrator is based on polynomials of degree with control points , , but uses the open trapezoidal rule (, , , , , ). Equations (3.8) take the form

In general one has to solve the first, third, and fourth equation simultaneously (3N equations for 3N variables). In case of separable Hamiltonians we have and one only needs to solve N nonlinear equations: can be explicitly calculated from the first equation and substituted into the third one, and the resulting nonlinear equation then has to be solved for .

P2N2Q2Otr

If the open trapezoidal rule is used with polynomials of degree , we obtain yet another partitioned Runge–Kutta method (3.18) with , , , , , . Inspecting Eq. (3.18) we see that, for example, is explicitly given in terms of and , therefore one only needs to solve 3N equations for the 3N variables , , , and the remaining equations are explicit updates. This method further simplifies for separable Hamiltonians H and h: and are now explicitly given in terms of and , and the nonlinear equation for can be solved separately from the nonlinear equation for .

P1N3Q4Mil

A method similar to (3.31) is obtained if we use Milne’s rule (, , , , , , , ) instead of Simpson’s rule. The resulting integrator is also computationally expensive in general, but if the Hamiltonians H and h are separable, then (3.8d) implies , and the integrator can be rewritten as

where

and and . In this case only the first equation needs to be solved for , and then the second equation is an explicit update.

Hamiltonian function h(q) independent of momentum

In case the Hamiltonian function is independent of the momentum variable p, the term does not enter Eq. (3.8d), and therefore we can allow a choice of quadrature rules such that or for some i. If , however, the system (3.8) becomes underdetermined, but at the same time the corresponding does not enter any of the remaining equations, therefore we can simply ignore it. To simplify the notation, let be the set of indices such that , and denote , for . Similarly, let be the set of indices such that , and denote , for . In (3.8) leave out the terms and equations corresponding to or , and replace , , and r by , , , , and , accordingly. In other words, this is equivalent to using the quadrature rules and in (3.6). We then simultaneously solve (3.8a), (3.8b) and (3.8d) [ equations] for and [ unknowns]. A few examples of such integrators are listed below.

P1N1Q1Rec (Stochastic symplectic Euler method)

The rectangle rule (, , ) does not yield a convergent numerical scheme in the general case, but when , the Itô and Stratonovich interpretations of (1.1) are equivalent, and the rectangle rule can be used to construct efficient integrators. In fact, applying the rectangle rule to both the Lebesgue and Stratonovich integrals and using polynomials of degree yield a method which can be written as

This method is a straightforward generalization of the symplectic Euler scheme (see [18, 40]) and is particularly computationally efficient, as only the first equation needs to be solved for , and then the second equation is an explicit update. Moreover, in case the Hamiltonian H is separable, the method becomes explicit.

P1N1Q1RecN2Q2Lob

The accuracy of the stochastic symplectic Euler scheme above can be improved by applying the trapezoidal rule to the Stratonovich integral instead of the rectangle rule. The resulting integrator takes the form

where

While having a similar computational cost, this method yields a more accurate solution than (3.36) (see Sect. 4 for numerical tests). Moreover, in case the Hamiltonian H is separable, the method becomes explicit.

P1N1Q1RecN1Q2Gau

Similarly, if we apply the midpoint rule instead of the trapezoidal rule, we obtain the following modification of the stochastic symplectic Euler method:

where

This method demonstrates a similar performance as (3.37) (see Sect. 4 for numerical tests). It becomes explicit if the Hamiltonian H is separable and the noise is additive, i.e., .

P2N2Q2LobN1Q1Rec

A modification of the stochastic Störmer–Verlet method (3.29) is obtained if we use the rectangle rule to approximate the Stratonovich integral:

This integrator has a similar computational cost as the stochastic Störmer–Verlet method (see Sect. 4), but it yields a slightly less accurate solution (see Sect. 4). Moreover, in case the Hamiltonian H is separable, the method becomes explicit.

P1N1Q2GauN2Q2Lob

This integrator is a modification of the stochastic midpoint method (3.28). We apply the midpoint rule (, , ) to the Lebesgue integral in (3.4), and the trapezoidal rule (, , , , ) to the Stratonovich integral. The resulting numerical scheme can be written as

where

This method is fully implicit, but similar to (3.28), simplifies when the Hamiltonian H is separable.

P1N2Q2LobN1Q2Gau

If instead we apply the trapezoidal rule to the Lebesgue integral and the midpoint rule to the Stratonovich integral in (3.4), we obtain a modification of the stochastic trapezoidal rule (3.30):

This method becomes explicit when the Hamiltonian H is separable and the noise is additive, i.e., .

Convergence

Various criteria for convergence of stochastic schemes have been suggested in the literature (see [28, 42]). Some criteria concentrate on pathwise approximations of the exact solutions (mean-square convergence, strong convergence), while others focus on approximations of some functionals instead (weak convergence). We are here primarily interested in mean-square convergence. Let be the exact solution to (1.1) with the initial conditions and , and let denote the numerical solution at time obtained by applying (3.8) iteratively k times with the constant time step . The numerical solution is said to converge in the mean-square sense with global order r if there exist and a constant such that for all we havewhere , as defined before, and E denotes the expected value. In principle, in order to determine the mean-square order of convergence of the Galerkin variational integrator (3.8) we need to calculate the power series expansions of and in terms of the powers of and , and compare them to the Stratonovich–Taylor expansions for the exact solution and (see [12, 28, 42]). It is quite a tedious task to do in the general case, therefore we only discuss the examples presented in Sect. 3.4. For instance, in case of the stochastic trapezoidal method (3.30) we plug in series expansions for , , and , and determine their coefficients by expanding the derivatives of the Hamiltonians into Taylor series around and comparing the terms corresponding to the like powers of and . We find thatwhere the derivatives of the Hamiltonians are evaluated at . Let and denote the exact solution of (1.1) such that and . Using (1.1) we calculate the Stratonovich–Taylor expansions for and , and comparing them to (3.46) we find thatUsing Theorem 1.1 from [42], we conclude that the stochastic trapezoidal method (3.30) has mean-square order of convergence . In a similar fashion we prove that all methods presented in Sect. 3.4 are convergent with mean-square order 1. We further verify these results numerically in Sect. 4.1.

Remark

For simplicity and clarity of the exposition, in this work we are mainly concerned with a one-dimensional noise in (1.1). However, all of the constructions and results presented in Sects. 2 and 3 generalize in a straightforward manner, when a multidimensional noise , together with the corresponding Hamiltonian functions , is considered in (1.1), except that the integrators derived in Sect. 3.4 in general do not attain mean-square order 1.0 of convergence, unless the noise is commutative. Indeed, if we repeat the procedure described above for the stochastic trapezoidal method, we will obtain the following power series expansions in terms of the powers of and :where the vectors and for each are defined asand the derivatives of the Hamiltonians are evaluated at . On the other hand, the Stratonovich–Taylor expansions for and read, respectively,where denotes a double Stratonovich integral. Comparing (3.49) and (3.50), we find that in the general case not all first order terms agree, and therefore we only have the local error estimates

Theorem 1.1 from [42] then implies that the stochastic trapezoidal method has mean-square order 1 / 2. However, if the noise is commutative, that is, if the following conditions are satisfiedthen using the property (see [28, 42]), one can easily showIn that case all first-order terms in the expansions (3.48) and (3.50) agree, and we again have the local error estimates (3.47), meaning that the scheme has mean-square order 1.0. Similar analysis holds for all the methods presented in Sect. 3.4. It should be noted that the commutation conditions (3.52) hold for two important special cases:The latter in particular means that the methods presented in Sect. 3.4.2 retain their mean-square order of convergence for multidimensional noises.

Hamiltonian functions linear in q and p for all , i.e. additive noise

Hamiltonian functions simultaneously independent of one of the variables q or p for all

Methods of order 3 / 2

In order to construct stochastic Galerkin variational integrators of higher order one needs to include higher order terms in the discretization of the Stratonovich integral in (3.5). For example, a method of mean-square order 3/2 must include terms involving (see [12, 43, 44]). Inspired by the theory presented in [12], we can add extra terms to the discrete Hamiltonian (3.6) and write it asThe random variables and have a Gaussian joint distribution (see [28, 44]), and at each time step they can be simulated by two independent -distributed random variables and asIn order to achieve mean-square convergence of order 3/2 one needs to determine appropriate values for the parameters s, r, , , , and . However, we will not attempt to achieve this in the present work. Instead, we will show that some known stochastic symplectic integrators can be derived as stochastic Galerkin variational integrators.

Suppose the Hamiltonian is separable, i.e., , and the Hamiltonian function does not depend on momentum. Consider the discrete Hamiltonianwhere different weights and were applied to the potential U(q) and kinetic T(p) terms, respectively. Similar to (3.8), the corresponding stochastic variational integrator takes the formwhere in the second equation, and in the fourth equation. In the special case when andwe can show, similar to Theorem 3.4, that the stochastic Galerkin variational integrator (3.57) is equivalent to the stochastic partitioned Runge–Kutta methodwith the coefficientswhere we assume and for all i. Partitioned Runge–Kutta methods of type (3.59) were considered in [44]. In particular, it was shown that for the choice of the coefficientsgives a method of mean-square order 3/2 (see Theorem 4.3 in [44]).

Numerical experiments

In this section we present the results of our numerical experiments. We verify numerically the convergence results from Sect. 3.5 and investigate the conservation properties of our integrators. In particular, we show that our stochastic variational integrators demonstrate superior behavior in long-time simulations compared to some popular general purpose non-symplectic stochastic algorithms.

Numerical convergence analysis

Kubo oscillator

In order to test the convergence of the numerical algorithms from Sect. 3.4.1 we performed computations for the Kubo oscillator, which is defined by and , where is the noise intensity (see [44]). The Kubo oscillator is used in the theory of magnetic resonance and laser physics. The exact solution is given bywhere and are the initial conditions. Simulations with the initial conditions , and the noise intensity were carried out until the time for a number of time steps . In each case 2000 sample paths were generated. Let denote the numerical solution. We used the exact solution (4.1) as a reference for computing the mean-square error , where . The dependence of this error on the time step is depicted in Fig. 1. We verified that our algorithms have mean-square order of convergence 1.0. The integrators P1N3Q4Lob, P1N3Q4Mil, P1N2Q2Lob (stochastic trapezoidal method), and P2N2Q2Lob (stochastic Störmer–Verlet method) turned out to be the most accurate, with the latter two having least computational cost.

Fig. 1

The mean-square error at the time as a function of the time step for the simulations of the Kubo oscillator with the initial conditions , and the noise intensity . In each case 2000 sample paths were generated. The tested integrators proved to be convergent of order 1.0 in the mean-square sense. Note that the plots for P2N2Q2Lob and P1N2Q2Lob, as well as for P1N3Q4Lob and P1N3Q4Mil, overlap very closely

Synchrotron oscillations of particles in storage rings

We carried out a similar test for the numerical schemes from Sect. 3.4.2. We performed computations for the stochastic Hamiltonian system defined by and , where is the noise intensity. Systems of this type are used for modeling synchrotron oscillations of a particle in a storage ring. Due to fluctuating electromagnetic fields, a particle performs stochastically perturbed oscillations with respect to a reference particle which travels with fixed energy along the design orbit of the accelerator; in this description p corresponds to the energy deviation of the particle from the reference particle, and q measures the longitudinal phase difference of both particles (see [17, 56] for more details). Simulations with the initial conditions , and the noise intensity were carried out until the time for a number of time steps . In each case 2000 sample paths were generated. The mean-square error was calculated with respect to a high-precision reference solution generated using the order 3/2 strong Taylor scheme (see [28], Chapter 10.4) with a very fine time step . The dependence of this error on the time step is depicted in Fig. 2. We verified that our algorithms have mean-square order of convergence 1.0.

Fig. 2

The mean-square error at the time as a function of the time step for the simulations of the synchrotron oscillations of a particle in a storage ring with the initial conditions , and the noise intensity . In each case 2000 sample paths were generated. The tested integrators proved to be convergent of order 1.0 in the mean-square sense. Note that the plots for P1N1Q1Rec, P1N1Q1RecN2Q2Lob, and P1N1Q1RecN1Q2Gau, as well as for P2N2Q2Lob and P1N2Q2LobN1Q2Gau, overlap very closely

Long-time energy behavior

One can easily check that in the case of the Kubo oscillator the Hamiltonian H(q, p) stays constant for almost all sample paths, i.e., almost surely. We used this example to test the performance of the integrators from Sect. 3.4.1. Simulations with the initial conditions , , the noise intensity , and the relatively large time step were carried out until the time (approximately 160 periods of the oscillator in the absence of noise) for a single realization of the Wiener process. For comparison, similar simulations were carried out using non-symplectic explicit methods like Milstein’s scheme and the order 3/2 strong Taylor scheme (see [28]). The numerical value of the Hamiltonian H(q, p) as a function of time for each of the integrators is depicted in Fig. 3. We find that the non-symplectic schemes do not preserve the Hamiltonian well, even if small time steps are used. For example, we find that Milstein’s scheme does not give a satisfactory solution even with , and though the Taylor scheme with yields a result comparable to the variational integrators, the growing trend of the numerical Hamiltonian is evident. On the other hand, the variational integrators give numerical solutions for which the Hamiltonian oscillates around the true value (one can check via a direct calculation that the stochastic midpoint method (3.28) in this case preserves the Hamiltonian exactly; of course this does not necessarily hold in the general case).

Fig. 3

The numerical Hamiltonian for the simulations of the Kubo oscillator with the initial conditions , and the noise intensity . Top The results obtained with Milstein’s scheme and the order 3/2 strong Taylor scheme. We see that the Hamiltonian tends to blow up despite using small time steps. Bottom The results obtained with the integrators derived in Sect. 3.4.1. For comparison, the solution obtained with the Taylor scheme for is also included. Note that for clarity the same color code is applied when the plots for some integrators overlap very closely (color figure online)

Anharmonic oscillator

In general the Hamiltonian H(q, p) does not stay constant for stochastic Hamilton equations. To determine how well our integrators perform in such cases we considered the anharmonic oscillator defined by and , where is the noise intensity and is a parameter. One can calculate the expected value of the Hamiltonian analytically asthat is, the mean value of the Hamiltonian grows linearly in time (see [56]). Simulations with the initial conditions , , the parameter , and the noise intensity were carried out until the time (approximately 100 periods of the oscillator in the absence of noise). In each case 10,000 sample paths were generated. The numerical value of the mean Hamiltonian E(H) as a function of time for each of the integrators is depicted in Fig. 4. We see that the variational integrators accurately capture the linear growth of E(H), whereas the Taylor scheme fails to reproduce that behavior even when a smaller time step is used. It is worth noting that the integrators P1N1Q1RecN2Q2Lob and P1N1Q1RecN1Q2Gau yield a very accurate solution, while being computationally efficient, as discussed in Sect. 3.4.2.

Fig. 4

Top The numerical value of the mean Hamiltonian E(H) for the simulations of the anharmonic oscillator with the initial conditions , , the parameter , and the noise intensity is shown for the solutions computed with the order 3/2 strong Taylor scheme using the time step and the variational integrators derived in Sect. 3.4.1 using the time step or . The variational integrators accurately capture the linear growth of E(H), whereas the Taylor scheme fails to reproduce that behavior. Middle The difference between the numerical value of the mean Hamiltonian E(H) and the exact value (4.2) is shown for the integrators derived in Sect. 3.4.1. Bottom Same for the integrators derived in Sect. 3.4.2. The integrators P1N1Q1RecN2Q2Lob and P1N1Q1RecN1Q2Gau prove to be particularly accurate, while having a low computational cost

Remark

One can verify by a direct calculation that when the P2N2Q2Otr integrator (example 6 in Sect. 3.4.1) is applied to the Kubo oscillator, then the corresponding system of Eqs. (3.18) does not have a solution when . To avoid numerical difficulties, one could in principle use the truncated increments (3.9) with, e.g., (for ). However, given the negligible probability that for the parameters used in Sects. 4.1.1 and 4.2.1, we did not observe any numerical issues, even though we did not use truncated increments. In the case of all the other numerical experiments presented in Sect. 4, the applied algorithms either turned out to be explicit, or the corresponding nonlinear systems of equations had solutions for all values of . Nonlinear equations were solved using Newton’s method and the previous time step values of the position and momentum were used as initial guesses.

Summary

In this paper we have presented a general framework for constructing a new class of stochastic symplectic integrators for stochastic Hamiltonian systems. We generalized the approach of Galerkin variational integrators introduced in [33, 40, 48] to the stochastic case, following the ideas underlying the stochastic variational integrators introduced in [8]. The solution of the stochastic Hamiltonian system was approximated by a polynomial of degree s, and the action functional was approximated by a quadrature formula based on r quadrature points. We showed that the resulting integrators are symplectic, preserve integrals of motion related to Lie group symmetries, and include stochastic symplectic Runge–Kutta methods introduced in [35, 36, 44] as a special case when . We pointed out several new low-stage stochastic symplectic methods of mean-square order 1.0 for systems driven by a one-dimensional noise, both for the case of a general Hamiltonian function and a Hamiltonian function independent of p, and demonstrated their superior long-time numerical stability and energy behavior via numerical experiments. We also stated the conditions under which these integrators retain their first order of convergence when applied to systems driven by a multidimensional noise.

Our work can be extended in several ways. In Sect. 3.6 we indicated how higher-order stochastic variational integrators can be constructed and showed that a type of stochastic symplectic partitioned Runge–Kutta methods of mean-square order 3/2 considered in [44] can be recast in that formalism. It would be interesting to derive new stochastic integrators of order 3/2 by choosing appropriate values for the parameters in (3.54) or (3.56). It would also be interesting to apply the Galerkin approach to construct stochastic variational integrators for constrained (see [7]) and dissipative (see [9]) stochastic Hamiltonian systems, and systems defined on Lie groups (see [32]). Another important problem would be stochastic variational error analysis. That is, rather than considering how closely the numerical solution follows the exact trajectory of the system, one could investigate how closely the discrete Hamiltonian matches the exact generating function. In the deterministic setting these two notions of the order of convergence are equivalent (see [40]). It would be instructive to know if a similar result holds in the stochastic case. A further vital task would be to develop higher-order weakly convergent stochastic variational integrators. As mentioned in Sects. 3.1 and 3.6, higher-order methods require inclusion of higher-order multiple Stratonovich integrals, which are cumbersome to simulate in practice. In many cases, though, one is only interested in calculating the probability distribution of the solution rather than precisely approximating each sample path. In such cases weakly convergent methods are much easier to use (see [28, 42]). 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 [24].