Stochastic Variational Partitioned Runge-Kutta Integrators for Constrained Systems

Stochastic Variational Partitioned Runge-Kutta Integrators for Constrained Systems

Nawaf Bou-Rabee    Houman Owhadi Applied & Computational Mathematics (ACM), Caltech, Pasadena, CA 91125 (,

Stochastic variational integrators for constrained, stochastic mechanical systems are developed in this paper. The main results of the paper are twofold: an equivalence is established between a stochastic Hamilton-Pontryagin (HP) principle in generalized coordinates and constrained coordinates via Lagrange multipliers, and variational partitioned Runge-Kutta (VPRK) integrators are extended to this class of systems. Among these integrators are first and second-order strongly convergent RATTLE-type integrators. We prove strong order of accuracy of the methods provided. The paper also reviews the deterministic treatment of VPRK integrators from the HP viewpoint.

1 Introduction

Since the foundational work of Bismut [1981], the field of stochastic geometric mechanics is emerging in response to the demand for tools to analyze the structure of continuous and discrete mechanical systems with uncertainty [Ha2007, ; VaCi2006, ; CiLeVa2008, ; Bi1981, ; MiReTr2002, ; MiReTr2003, ; LaOr2007a, ; LaOr2007b, ; MaWi2007, ]. Within this context the goal of this paper is to develop efficient, structure-preserving integrators for long-time simulations of constrained, mechanical systems perturbed by white-noise forces and torques. Our strategy is to employ stochastic variational integrators (SVIs) [BoOw2007a, ].

Variational integration theory derives integrators for mechanical systems from discrete variational principles [Ve1988, ; Ma1992, ; WeMa1997, ; MaWe2001, ]. The theory includes discrete analogs of the Lagrangian, Noether’s theorem, the Euler-Lagrange equations, and the Legendre transform. Variational integrators can readily incorporate holonomic constraints (e.g., via Lagrange multipliers) and non-conservative effects (via their virtual work) [WeMa1997, ; MaWe2001, ]. Altogether, this description of mechanics stands as a self-contained theory of mechanics akin to Hamiltonian, Lagrangian or Newtonian mechanics.

Variational integrators are symplectic, i.e., the discrete flow map they define exactly preserves the continuous symplectic 2-form. Using backward error analysis one can show that symplectic integrators applied to Hamiltonian systems nearly preserve the energy of the continuous mechanical system for exponentially long periods of time and that the modified equations are also Hamiltonian [HaLuWa2006, ]. Variational integrators are also distinguished by their ability to compute statistical properties of mechanical systems, such as in computing Poincaré sections, the instantaneous temperature of a system, etc.

Stochastic variational integrators are an extension of variational integrators to so-called stochastic mechanical systems. These systems are simple mechanical systems subject to certain random perturbations, and have their origins in Bismut’s foundational work [Bi1981, ]. Bismut showed that the flow of these stochastic mechanical systems extremize an action integral whose domain is the space of semimartingales on configuration space. Bismut’s work was further enriched and generalized to manifolds by recent work [LaOr2007a, ; LaOr2007b, ]. Lazaro-Cami and Ortega show that this general class of stochastic Hamiltonian systems on manifolds extremizes a stochastic action defined on the space of manifold-valued semimartingales [LaOr2007a, ]. Moreover, it has been shown that for a subclass of these systems, one can prove a converse, namely, a.s. a curve satisfies so-called stochastic Hamilton’s equations if and only if it extremizes a stochastic action [BoOw2007a, ].

With this variational principle, one can design SVIs [BoOw2007a, ]. Like their deterministic counterparts, these methods have the advantage that they are symplectic, and in the presence of symmetry, satisfy a discrete Noether’s theorem. Moreover, symplectic methods for stochastic mechanical systems on vector spaces have been shown to capture the correct energy behavior even in the presence of dissipation [MiReTr2002, ; MiReTr2003, ]. In particular, the energy injected or dissipated by the symplectic integrator is not an artifical function of the timestep. Moreover, the energy behavior is accurately captured by the integrator.

These structure-preserving properties are quite crucial. Consider, for instance, simulation of a simple mechanical system subject to random forces with amplitude . Suppose further the unperturbed system preserves energy, momentum, and possesses a first integral. Consider simulating this system with a higher-order accurate method, a standard integrator with simultaneous projection onto energy, momentum, and first integral level sets, and a stochastic variational integrator. If , a stochastic (physical) perturbation in the energy, momentum, and first integrals will appear. However, it is not clear how to modify the projection-based method to correctly capture these stochastic perturbations when . Moreover, the higher-order accurate method requires a time-step smaller than the amplitude of the perturbation in order to accurately represent its effects. And even then, if the time-span of integration is long enough systematic, artificial drift in these quantities will appear.

On the other hand, the stochastic perturbation in energy and momentum captured by the stochastic variational integrator is mechanical, i.e., it is only due to the -random forces. This is because these schemes define flows of discrete mechanical systems. Moreover, even when the amplitude of perturbation is not small, due to symplecticity we conjecture that SVIs on manifolds not only possess good long-time energy behavior, but also perform well in computing statistical properties, such as autocorrelation functions and the empirical distribution. We will investigate these questions in future work.

As far as we can tell, the extension of these structure-preserving integrators to stochastic mechanical systems with holonomic constraints has not been completed. The main goal of this paper is to extend SVIs to holonomic constraints, and in particular, introduce constrained, stochastic variational partitioned Runge-Kutta (VPRK) methods for such systems. Within this family we exhibit a first-order strongly convergent, symplectic integrator for constrained mechanical systems. We use a technique due to Vanden-Eijnden and Cicotti [2006] to prove order of accuracy of the integrators that appear in this paper [VaCi2006, ].

Future work will consider extensions of these schemes to Langevin equations with holonomic constraints. Continuous Langevin processes have been generalized to submanifolds and shown to be ergodic [CiLeVa2008, ]. The generalization to constraint submanifolds was done by appending holonomic constraints to Langevin equations. One can then check that the infinitesimal generator of the constrained Langevin process commutes with the Gibbsian density restricted to to determine that the restricted Gibbsian measure is an invariant measure of the constrained Langevin process. To prove this measure is unique one uses standard arguments based on nondegeneracy of the diffusion and drift vector fields on the momentums. On this note it would be interesting to ascertain ergodicity of SVIs for constrained Langevin systems.

2 Constrained VPRK Integrators

This section is provided to fix notation and clarify some aspects of deterministic constrained VPRK integrators which will be relevant in the stochastic context. We also provide a novel proof showing that constrained VPRK integrators can be derived directly from a discrete variational principle without explicitly introducing a discrete Legendre transform. For more details on unconstrained VPRK integrators we refer the reader to [Bo2007, ].

The setting of this section is a real, -dimensional vector space and a mechanical system with smooth, holonomic constraint function, , , that has a regular value at . The mechanical system’s configuration space is given by the constraint submanifold: . We will introduce Lagrange multipliers to prove an equivalence between a Hamilton-Pontryagin (HP) variational principle on and a constrained HP variational principle on . We will then show how to discretize this system to obtain variational RATTLE integrators for constrained mechanical systems.

The variational and symplectic character of VPRK integrators is discussed in [Su1990, ; MaWe2001, ; HaLuWa2006, ]. In what follows we will explicitly use the HP perspective, and specifically, extend the results in Hairer et al. [2006] to constrained systems using the HP perspective. It is worth mentioning, the Hamiltonian counterparts of the constrained VPRK methods, so-called symplectic partitioned Runge-Kutta methods, are also well understood [Ja1996, ; Re1997, ; Ha2003, ].

Discretization of HP Action

We will adopt a HP viewpoint to describe the action integral of this mechanical system with constraints. The HP description unifies the Hamiltonian and Lagrangian descriptions of a mechanical system [YoMa2006a, ; YoMa2006b, ; Bo2007, ; BoMa2007, ].

Definition 2.1.

The Pontryagin bundle of a manifold is defined as . Fixing the interval and , define the HP path space on as

The HP path space is a smooth infinite-dimensional manifold. One can show that its tangent space at consist of maps such that .

Definition 2.2.

Fix . Define the unconstrained HP action as

To discretize we first discretize the kinematic constraint: . An s-stage RK method is employed to discretize the kinematic constraint. Let and be given and define the fixed step size and , . The reason for using an s-stage RK discretization of the kinematic constraint is that the theory on such methods (order conditions, stability, and implementation) is mature. See, for instance, [HaNoWa1993, ].

Definition 2.3.

Consider the first order differential equation


Let () and let . An s-stage RK approximation is given by


The vectors and are called external and internal stage vectors, respectively.

It follows that an s-stage RK method is fully determined by its -matrix and -vector which are typically displayed using the so-called Butcher tableau:

Suppose that , , is given. Then an s-stage RK approximant applied to yields:


In what follows for will be introduced as internal stage unknowns and will be determined as a critical point of a discrete action sum.

VPRK Integrator for Constrained Systems

The VPRK method will be derived from a discretization of the HP action integral in which the kinematic constraint over the kth-time step is replaced with its discrete approximant: (3), and the integral of the Lagrangian over the kth-time step is approximated by the following quadrature:

The constraint is enforced for all internal stage positions using Lagrange multipliers as follows.

Definition 2.4.

Fix two points and on and define the discrete VPRK path space as:

and the discrete constrained VPRK action sum by:

The following condition on the coefficients of the s-stage RK method will be important in obtaining a well-defined discrete update on phase space that also respects the constraints [Ja1996, ].

Condition 2.1.

Consider an s-stage RK method with given -vector and -matrix. Assume and set . The coefficients of the s-stage RK method satisfy:

Under this condition on the coefficients of the s-stage RK method, one can prove the following theorem.

Theorem 2.5.

Given an s-stage RK method with -vector and -matrix that satisfy condition 2.1, a Lagrangian system with smooth Lagrangian such that is invertible, and smooth holonomic constraint function . A discrete curve satisfies the following VPRK method:


for and , if and only if it is a critical point of the function , that is, . Moreover, there exist such that the discrete flow map defined by the above scheme, , preserves the canonical symplectic form on .


Under the assumptions of the theorem, the existence of a numerical solution of (4) and a discrete flow map is guaranteed. That is, given one can solve (4) for . In particular, the condition that is invertible, implies one can eliminate using the Legendre transform of . The condition 2.1 implies that one can determine so that [Ja1996, ].

The differential of in the direction is given by:

Collecting terms with the same variations and summation by parts using the boundary conditions gives,

Since if and only if for all , one arrives at the desired equations with the elimination of and the introduction of the internal stage variables for . Conversely, if satisfies (4) then .

We will employ the variational proof of symplecticity to prove that is symplectic. This proof is standard, however, we provide it here to emphasize that the symplectic form that is exactly preserved by the method is the canonical one on .

Consider the subset of given by solutions of (4). Let denote the restriction of to this space. Since each of these solutions is determined by an initial point on , one can identify this space with , and hence, . Since is restricted to solution space,

Observe that these boundary terms are the canonical one-forms on evaluated at and . Preservation of the symplectic form follows from . ∎

Variational RATTLE for Constrained Systems

The variational RATTLE integrator is the Lagrangian analog of the RATTLE algorithm originally proposed as a constrained version of Verlet in [RyCiBe1977, ]. It was shown to be symplectic in [LeSk1994, ], and was extended to general constrained Hamiltonian systems by [Ja1996, ]. It is defined by the following two-stage RK discretization of the kinematic constraint (implicit trapezoidal rule),

Given and , the method determines and two Lagrange multipliers by solving the following system of equations,


By theorem 2.5, variational RATTLE defines a symplectic scheme. Moreover, as is well-known it is second-order accurate.

Variational Euler for Constrained Systems

One can relax the conditions assumed in theorem 2.5 on the coefficients of the s-stage RK method [Ha2003, ]. For instance, if in the s-stage RK method one can still obtain a well-defined variational integrator. Consider as an example the following two-stage RK method,

The corresponding variational integrator is given by:


However, the corresponding discrete flow does not satisfy the “hidden” velocity constraint, and hence, does not define a map from to . To satisfy the hidden constraint a projection step is taken:


One can check that this projection step defines a symplectic map. Thus, the composite map is symplectic. It is the Lagrangian version of the constrained symplectic Euler method [HaLuWa2006, ].

Another example relaxing the assumptions in theorem 2.5 is given by regarding the 2-stage RK method as being single-stage implicit Euler,

The corresponding variational integrator is:


We conclude with a theorem summarizing the structure-preserving properties of these constrained variational Euler methods.

Theorem 2.6.

The composition of one step of (6) or (8) and the projection step (7), defines a symplectic integrator on . Moreover, these integrators are first-order accurate.

3 Constrained, Stochastic Mechanical Systems

Constrained, Stochastic HP Principle

The setting in this section is an n-manifold and a stochastic mechanical system with smooth, holonomic constraint function, , , that has a regular value . Its configuration space is given by the constraint submanifold: . In this section we will introduce Lagrange multipliers to prove an equivalence between a stochastic variational principle on and a constrained, stochastic variational principle on .

We will adopt a HP viewpoint to describe this mechanical system with random perturbations and will refer to this perturbed system as a stochastic mechanical system. The HP principle unifies the Hamiltonian and Lagrangian descriptions of a mechanical system [YoMa2006a, ; YoMa2006b, ; Bo2007, ; BoMa2007, ]. Roughly speaking, in the stochastic context it states the following critical point condition on (cf. definition 2.1),

where are varied arbitrarily and independently with endpoint conditions and fixed, and and for are stochastic potentials and Wiener processes respectively. This principle builds in a Legendre transform, stochastic Hamilton’s equations and stochastic Euler–Lagrange equations. The action integral in the above principle, consists of two Lebesgue integrals with respect to the Lebesgue measure and Stratonovich stochastic integrals. This action is random, i.e., for every sample point one will obtain a different, time-dependent Lagrangian system. However, each system possesses a variational structure as made precise in [BoOw2007a, ]. The following definitions will be useful to state the constrained, variational principle of HP for mechanical systems with holonomic constraints.

Definition 3.1.

Let be the inclusion mapping. Fixing the interval , define constrained HP path space as

and the unconstrained HP path space as

Let be a probability space, , , be a nondecreasing family of -subalgebras of , , and , , be independent, real-valued Wiener processes. In terms of these Wiener processes, we define the following.

Definition 3.2.

Set and for . Moreover define the unconstrained action as

and the constrained action as .

The unconstrained HP path space is a smooth infinite-dimensional manifold. One can show that its tangent space at consist of maps such that and are of class .

As opposed to using generalized coordinates on , we wish to describe the mechanical system using constrained coordinates on and introduce Lagrange multipliers to enforce the constraint. However, because of the stochastic component of the action, the standard Lagrange multiplier theorem will not apply directly and one cannot introduce Lagrange multipliers in the standard way. Instead, we will introduce the Lagrange multiplier using the following.

Definition 3.3.

Given and , define

Differentiability of the flow map on will be defined in the mean-squared sense. In the following we define mean-squared derivatives on a Banach space with the understanding that this notion can be extended to any manifold using a local representative of the flow map.

Definition 3.4 (Mean-squared Derivative).

The mean squared norm of is given by:

Using this norm one can define the derivative of in the standard way, i.e., is mean squared differentiable at if there is a bounded linear map that satisfies,

With the above definitions we prove the following.

Theorem 3.5 (Constrained, Stochastic HP Principle).

Given a stochastic mechanical system with Lagrangian such that is invertible, stochastic potentials for , and holonomic constraint with . Then the following are equivalent:


extremizes .


satisfies stochastic HP equations


There exists such that and extremize the augmented action where and .


There exists such that and satisfy the constrained, stochastic HP equations


Moreover, the flows of (9) and (10) are mean-square symplectic.

Remark 3.1.

The constrained, stochastic HP equations should be thought of in integral form. In particular, the Lagrange multipler term in (10) by definition 3.3 satisfies:

for any .


The stochastic HP principle states that (i) and (ii) are equivalent [BoOw2007a, ]. Assume that (iii) is true. Then for all ,


Since for and for ,

which implies (i) and (ii).

Moreover, expanding (11) and setting and yields,

Consider the terms involving , and . Since these variations are arbitrary, the following hold:111This follows from the basic lemma that if and is arbitrary then .

Since , implies that for all .

Collecting the variations with respect to in the differential gives,

Using definition 3.3, we introduce the following function ,

In the following it is shown that if for arbitrary of class then satisfies (10).

Let be a partition of unity on . Expand in terms of this partition of unity,

Since the curves and are compactly supported, only a finite number of the are nonzero. For each nonzero, the terms in the integral can be expressed in local coordinates. Observe that since , the Stratonovish-Ito conversion formula implies that,

for .

We will select to single out the th-component of the covector field in . Introduce the following function for this purpose:

Observe that , , and . Let be a basis for the model space of . Now fix , and define in local coordinates as:

Introduce the following label to simplify subsequent calculations,

In terms of , one can write

We will show in the mean squared norm,


Using this result and the Borel-Cantelli lemma, one can deduce there exists that converges to such that a.s. converges to . It follows that almost surely.

We proceed to prove (12). Since ,