Strong convergence of a halfexplicit Euler scheme for constrained stochastic mechanical systems
Abstract.
This paper is concerned with the numerical approximation of stochastic mechanical systems with nonlinear holonomic constraints. Such systems are described by second order stochastic differentialalgebraic equations involving an implicitly given Lagrange multiplier process. The explicit representation of the Lagrange multiplier leads to an underlying stochastic ordinary differential equation, the drift coefficient of which is typically not globally onesided Lipschitz continuous. We investigate a halfexplicit drifttruncated Euler scheme which fulfills the constraint exactly. Pathwise uniform convergence is established. The proof is based on a suitable decomposition of the discrete Lagrange multipliers and on norm estimates for the single components, enabling the verification of consistency, semistability and moment growth properties of the scheme. To the best of our knowledge, the presented result is the first strong convergence result for a constraintpreserving scheme in the considered setting.
AMSClassification. Primary 60H35, 74Hxx; Secondary 60H10, 58J65, 65C30.
Keywords. Stochastic differentialalgebraic equation; manifoldvalued stochastic differential equation; nonlinear constraint; numerical approximation; drifttruncated scheme; strong convergence.
1. Introduction
Both the numerical approximation of differentialalgebraic equations (DAEs) and of stochastic differential equations (SDEs) have been extensively studied in the literature. Convergence results for higher index DAEs can be found, e.g., in [3, 6, 8, 9, 25], and the convergence analysis of numerical schemes for SDEs with nonglobally Lipschitz continuous coefficients has been an active and rapidly evolving field of research within the last years, see, e.g., [1, 2, 5, 7, 10, 13, 14, 15, 16, 28, 29, 37, 38]. In contrast, the convergence analysis of numerical schemes for stochastic differentialalgebraic equations (SDAEs) is far less developed. In this article, we combine key concepts from both areas, numerics of DAEs and numerics of SDEs, and prove strong convergence of a constraintpreserving numerical scheme for a large class of second order SDEs with nonlinear algebraic constraints.
We are interested in the dynamics of constrained stochastic mechanical systems which are modelled by SDAEs of the type
(1.1) 
with dimensional position and velocity processes and , positive definite and symmetric mass matrix , dimensional Brownian motion , coefficient functions and , and a sufficiently smooth constraint function , . By we denote the transpose of the Jacobian matrix of at , which is assumed to be of full rank for all in a neighborhood of the constraint manifold The process is an dimensional semimartingale, implicitly given as the Lagrange multiplier to the holonomic constraint . It thus determines the constraint force, which, according to d’Alembert’s principle, acts orthogonal to the tangent space at the current position . Constrained SDEs of the form (1.1) occur in various applications, ranging from molecular dynamics [20, 19, 21, 22, 40, 41] to models for the dynamics of fibers in turbulent airflows in the context industrial production processes of nonwoven textiles [23, 24].
Our assumptions on the coefficients , and the constraint function in (1.1) are weak enough to cover a large variety of practically relevant examples. We assume that and are locally Lipschitz continuous, of polynomial and linear growth, respectively, and that the mapping satisfies a onesided linear growth condition. Note that the diffusion coefficient is allowed to depend on both position and velocity, which is motivated by [23, 24]. The derivatives of are assumed to satisfy suitable boundedness and nondegeneracy conditions in a neighborhood of the manifold. Our setting allows in particular for quadratic constraint functions, which are relevant in various applications, cf. Section 5. The precise assumptions are stated in Section 2.1. By a slight generalization of the existence result in [23], we know that for all initial conditions , there exists a unique global strong solution to the SDAE (1.1), see Section 2.2 for details.
There seem to be no strong convergence results available in the literature so far for constraintpreserving schemes for SDAEs of the type (1.1) with quadratic constraint functions . A particular class of such SDAEs are the constrained Langevintype equations considered in molecular dynamics. Various constraintpreserving numerical schemes have been proposed and applied in this context, see, e.g., [20, 21, 22, 40] and the references therein. Typically, in the corresponding sections of these works the focus lies mainly on the practical efficiency of the algorithms but not so much on a fully rigorous convergence analysis. Often the proposed partly implicit schemes are even known to be not always solvable, cf. [20, Section 2]. We are also not aware of fully completed proofs concerning weak convergence. Further works related to our problem concern SDEs which are given in an explicit form, without an implicit Lagrange multiplier process as in (1.1), and the analysis of structurepreserving algorithms for their numerical approximation, see, e.g., [11, 27, 32, 33, 34, 43]. The theoretical results in these works rely on classical global Lipschitz assumptions on the drift and diffusion coefficients. We note that the Lagrange multiplier process in (1.1) can be represented explicitly in terms of and , so that the SDAE can be equivalently reformulated as an inherent SDE which does not involve an implicit Lagrange multiplier, see Eq. (2.9) in Section 2.2. However, the drift coefficient appearing in this inherent SDE is typically neither globally Lipschitz continuous nor globally onesided Lipschitz continuous, even if the coefficients and in (1.1) are chosen to be constant, see Section 5.1 for a simple toy example. Thus, the results in the mentioned works are not applicable in our setting.
Our numerical scheme combines ideas from different areas: On the one hand, we consider a GearGuptaLeimkuhler (GGL) reformulation of (1.1) to simplify the approximation problem and use a halfexplicit method where only the algebraic variables are discretized in an implicit manner. These concepts are standard in DAE theory, cf. [3, 6, 8, 9, 26, 35]. On the other hand, we follow the taming or truncation approach developed in recent years in the context of numerical methods for SDEs with nonglobally Lipschitz continuous coefficients, cf. [14, 16, 28, 29, 37]. We refer to Section 3.1 for a short discussion of the concepts. The combination of these ideas motivates the following halfexplicit drifttruncated Euler scheme for the approximation of the SDAE (1.1). Given initial conditions , , a finite time interval , and a number of time steps , we set and approximate the solution processes and in (1.1) by the timediscrete processes and iteratively defined by
(1.2) 
Here is a scalar truncation term, chosen as
(1.3) 
with a suitable constant , see Section 3.1 for details and a discussion. The dimensional Lagrange multipliers and are implicitly determined by the constraints in the third and forth line of (1.2), and is the increment of the driving Brownian motion on the time interval . Since the Lagrange multiplier formally corresponds to the infinitesimal increment in (1.1), it is natural to use the timediscrete process defined by
(1.4) 
as an approximation of the Lagrange multiplier process in (1.1). We remark that the presence of the truncation factor in (1.2) fulfills two purposes: It ensures the solvability of the scheme and moreover allows for the derivation of moment bounds needed to obtain strong convergence.
In this article, we verify the strong convergence of the scheme (1.2) towards the system (1.1) in a pathwise uniform sense as described below. To the best of our knowledge, this is the first proof of strong convergence for a constraintpreserving scheme for SDAEs of the type (1.1) with possibly quadratic constraint functions . Moreover, it also seems to be the first proof of convergence at all for a constraintpreserving scheme in the considered general setting.
Our first main result, Theorem 3.4, states that the scheme (1.2) is uniquely solvable for all choices of and , , in the sense that there exists exactly one solution such that the Lagrange multiplier process satisfies a specific boundedness condition. This existence result is nontrivial, in particular in view of the fact mentioned above that alternative schemes considered in the literature are often known to be not always solvable. The proof of Theorem 3.4 is based on a homotopy argument and relies on the presence of the truncation term in (1.2). Note that we do not impose any condition on the stepsize and moreover no truncation of the noise term is needed. Besides the existence and uniqueness of a numerical solution, Theorem 3.4 also provides a suitable decomposition of the Lagrange multiplier as well as suitable norm estimates for the single components of and for , crucial for proving convergence of the scheme.
Our second main result, Theorem 4.1, concerns the strong convergence of the scheme (1.2). Let , and be defined by piecewise constant or piecewise linear interpolation of the corresponding timediscrete processes , and in (1.2) and (1.4). For instance, in the piecewise linear case we have for all , . Theorem 4.1 states that the timeinterpolated solution to (1.2), (1.4) converges strongly towards the solution to the SDAE (1.1) with initial conditions , in the sense that
(1.5) 
for all . The main idea of the proof is to use the existence, uniqueness and decomposition result from Theorem 3.4 to formally rewrite the scheme (1.2) as a fullyexplicit onestep approximation of the underlying inherent SDE, an explicit drifttruncated Euler scheme with an additional explicit perturbation term, and to use the norm estimates for the Lagrange multipliers and from Theorem 3.4 to verify suitable consistency, semistability and moment growth conditions for the discrete solution. This enables us to apply a general convergence result from [14] to obtain that . The pathwise uniform strong convergence of and the convergence of to the Lagrange multiplier process are then proven in separate steps, by further exploiting the assumptions on the coefficients and , the decomposition and norm estimates from Theorem 3.4 for the Lagrange multipliers , , and the equivalence of (1.1) and the underlying inherent SDE.
Let us shortly discuss the question of convergence rates. Most strong approximation results with rates for multidimensional SDEs are based on at least a global monotonicity assumption on the coefficients and thus in particular on a global onesided Lipschitz assumption on the drift coefficient, see, e.g., [10, 16, 37]. As already mentioned, the drift coefficient of the inherent SDE associated to the SDAE (1.1) typically fails to satisfy such a condition, cf. Section 2.2 and the example in Section 5.1. A recently developed strategy to obtain strong convergence rates in the case of SDEs with nonglobally monotone coefficients makes use of exponential integrability properties of suitably tamed/truncated schemes [13]. This approach might potentially also be useful in the context of SDAEs of the type (1.1) but lies beyond the scope of the present work.
To complete the picture, let us also note that the approximation of socalled index one stochastic differentialalgebraic equations has been analyzed in several papers, mainly in the context of electric circuits, cf. [17, 18, 36, 39, 42] and the references therein. These equations are of a different structure than (1.1). In particular, the standard assumptions used in the context of index one SDAEs, e.g., that the constraints are globally uniquely solvable for the algebraic variables [42], are not fulfilled in our setting. One can think of Theorem 3.4 as partly substituting such assumptions.
The article is organized as follows: In Section 2 we describe our setting in detail (Subection 2.1) and state an existence and uniqueness result which ensures the global strong solvability of the SDAE (1.1) as well as the equivalence of (1.1) and the above mentioned underlying inherent SDE (Subsection 2.2). In Section 3 we first specify and discuss our approximation scheme (Subsection 3.1) before we give a detailed analysis of its solvability in Theorem 3.4 and Corollary 3.6 (Subsection 3.2). The strong convergence result (1.5) is proven in Section 4; a complete formulation is given in Theorem 4.1. Its proof is based on the results from Section 3 as well as a suitable reformulation of the problem (Subsection 4.1) and a collection of auxiliary results from the literature (Subsection 4.2). We deduce the convergence of in a nonpathwise sense by verifying specific consistency, semistability and moment growth conditions (Subsection 4.3) and show the pathwise uniform convergence of in a separate step (Subsection 4.4). Concrete examples for SDAEs of the type (1.1) are given in Section 5. Finally, a proof of the global strong solvability of (1.1) is sketched in Appendix A and a globalized version of the implicit function theorem used in the proof of Theorem 3.4 is presented in Appendix B.
2. Preliminaries
General notation. is the set of natural numbers excluding zero. In the sequel, let . By and we denote the Euclidean norm and the corresponding inner product in finitedimensional real vector spaces. For instance, we have for , for , and for . If is a differentiable function, we write for its Jacobian matrix at a point and for the transpose thereof. By and , , we denote the spaces of linear operators from to and fold multilinear operators from to , respectively, with norms and . The Jacobian matrix of a sufficiently smooth function at a point is identified with the corresponding element in , and, accordingly, the higher derivatives , , are elements in . If is a random variable on a probability space and we write .
2.1. Setting and assumptions
Here we describe our setting concerning the SDAE (1.1) in detail. Let , , and be a symmetric and (strictly) positive definite matrix. The constraint function in (1.1) determines the constraint manifold
(2.6) 
and satisfies the regularity conditions stated in Assumption 2.2 below. We first introduce some notation associated with .
Notation 2.1 (Constraint manifold).
The following notation associated with the constraint manifold in (2.6) is used throughout the article.

The tangent space at a point and the tangent bundle are given by

For we denote by
the Gram matrix associated with the constraint, where is the positive definite and symmetric mass matrix appearing in (1.1).

For all such that is invertible we set
where . Note that for the matrix represents the orthogonal projection of onto , corresponding to the inner product .

For we introduce the environments
of the constraint manifold and the tangent bundle . Here is the open ball in with radius and center .
With this notation at hand, we are able to state the regularity assumptions on the constraint function in detail.
Assumption 2.2 (Constraint function ).
The constraint function in (1.1) is three times continuously differentiable. There exists such for all the Jacobian matrix has full rank, and
Moreover, the higher order derivative mappings and are bounded.
We assume that is a complete probability space, a filtration of subalgebras of satisfying the usual conditions, and the process in (1.1) is an dimensional standard Wiener process on , for some . For the coefficient functions and we assume the following.
Assumption 2.3 (Coefficient functions and ).
The mappings and in (1.1) fulfill the following conditions.

Local Lipschitz continuity: For all there exist a constant such that, for all ,
Here is the open ball in with radius and center .

Growth conditions for : The function satisfies the following onesided linear growth and polynomial growth conditions. There exist constants and such that, for all ,

Linear growth of : There exists a constant such that, for all ,
Finally, we specify the the intitial conditions and for the position and velocity processes and in (1.1).
Assumption 2.4 (Initial conditions).
The initial conditions , are measurable random variables, integrable for all , such that .
2.2. Solvability of the SDAE
Under the assumptions in Subsection 2.1 there exists a unique global strong solution to the SDAE (1.1). We specify the notion of a solution as follows.
Definition 2.5.
The following result concerning the global strong solvability of SDAEs of the type (1.1) is a slight generalization of [23, Theorem 3.1]. Compared to [23] we consider weaker assumptions on the drift coefficient , a more general class of constraint functions , is not assumed to be the identity matrix, and finite moments of all orders are established.
Theorem 2.6 (Existence and uniqueness).
Let the assumptions in Subsection 2.1 be fulfilled. Then there exists a unique (up to indistinguishability) global strong solution to the SDAE (1.1) with initial conditions in the sense of Definition 2.5. With probability one, the following equality holds for all :
(2.8) 
Moreover, for all the th moment is finite.
The proof of Theorem 2.6 is mostly analogous to that of [23, Theorem 3.1]. A sketch of the key steps of the proof and the differences to [23] can be found in Appendix A.
As a direct consequence of Theorem 2.6 we know that, under the assumptions in Subsection 2.1, the following holds: If is the unique global strong solution to the SDAE (1.1) with initial conditions , then solves the socalled inherent or underlying SDE
(2.9) 
Conversely, consider arbitrary locally Lipschitz continuous extensions of the coefficients in the inherent SDE (2.9) to the whole space , and let be a global strong solution to (2.9) with initial conditions . Then it is not difficult to check that and, if further is the continuous valued semimartingale defined by (2.8), then is the unique global strong solution to the SDAE (1.1).
3. Analysis of the numerical scheme
We shortly discuss our numerical scheme in Subsection 3.1 before we present a detailed analysis of its solvability in Subection 3.2. The results in this section are essential for the convergence analysis in Section 4.
3.1. A halfexplicit drifttruncated Euler scheme
Suppose that the assumptions in Section 2.1 are fulfilled and recall from Section 1 the halfexplicit drifttruncated Euler scheme, which we specify and slightly reformulate as follows: Given we are searching for valued processes , and valued Lagrange multiplier processes , defined recursively by
(3.10a)  
(3.10b) 
Here is the step size, is a Brownian increment, and the truncation function is given by
(3.11) 
with depending on the constant introduced in (2.7). Let us remark that the choice of considering both and in the maximum in (3.11), and not solely , is mainly for convenience reasons as it guarantees a precise control of the norm of also for small values of , independently of the step size ; alternatively one could impose a suitable step size restriction.
In the following remarks we sketch the main concepts the scheme (3.10) is based on. As these are wellknown in the respective scientific communities, we refer to the mentioned literature and the references therein for more detailed expositions. The concepts presented in the first two remarks are standard in the context of DAEs.
Remark 3.1 (GearGuptaLeimkuhler formulation).
Compared to the SDAE (1.1) the scheme (3.10) involves the additional Lagrange multiplier in (3.10a) as well as the additional constraint in (3.10b). It is clear that the solution to (1.1) also solves the SDAE
(3.12)  
if we choose the integrator process to be identically zero. This is the socalled GearGuptaLeimkuhler (GGL) reformulation of (1.1). Observe that the scheme (3.10) is a discrete version of (3.12), with and corresponding to the infinitesimal increments and . The GGL stabilization is a standard index reduction technique for deterministic mechanical systems, compare [6, 8] and [9, Chapter VII]. Loosely speaking, it reduces the influence that a perturbation of the constraint has on the Lagrange multiplier, see [8, Chapter 1].
Remark 3.2 (Halfexplicit schemes).
The main idea of halfexplicit methods for deterministic DAEs is to discretize the differential variables in an explicit manner and only the algebraic variables in an implicit manner, see [26, 35] and [9, Section VII.6]. Thereby, the dimension of the system of equations which has to be solved implicitly in each time step is kept minimal. In our reformulated SDAE (3.12) one can interpret , as the “differential variables” and the formal time derivatives of , as the “algebraic variables”. As a particularly useful consequence of our halfexplicit approach, the position in (3.10) is stochastically independent of the Brownian increment .
Next, we explain the concept of truncating or taming, which has been used and studied extensively in the last years in the context of SDEs with nonglobally Lipschitz continuous coefficients, cf. [1, 5, 10, 13, 14, 16, 28, 29, 37].
Remark 3.3 (Truncated and tamed schemes).
The concept of truncation or taming is used to obtain strongly convergent explicit methods for SDEs whose coefficients do not fulfill global Lipschitz conditions. In the onedimensional setting it has been shown that the classical EulerMaruyama scheme may diverge if the coefficient functions are of superlinear polynomial growth, see [15] for details. The divergence follows from the existence of a sequence of events of exponentially small probability on which the numerical approximations grow at least doubleexponentially fast. The idea of tamed or truncated schemes is to adjust the coefficient functions in such a way that this growth behaviour is avoided while the adjustment is negligible with high probability. In our setting, the presence of the truncation term in (3.10) enables the derivation of suitable moment bounds for the numerical solution, but it also ensures the existence of a numerical solution at all, compare the example in Section 5.1. Let us further note that the standard drifttruncated explicit Euler scheme for the inherent SDE (2.9) involves the truncation function
compare [14, Section 3.6]. The choice (3.11) for our halfimplicit scheme is a simplification of this truncation function, making use of the boundedness properties of the constraint function from Assumption 2.2. This simplification reduces the computational effort significantly since the calculations of , and are avoided in each time step.
3.2. Solvability of the scheme and Lagrange multiplier estimates
Here we verify the unique solvability of the halfexplicit drifttruncated Euler scheme (3.10) and derive a suitable decomposition as well as norm estimates for the Lagrange multipliers , . Let us stress that, although the main result in this subsection, Theorem 3.4, is formulated in a deterministic setting, it is tailormade for the analysis of our stochastic problem, cf. Remark 3.7. Its application to the scheme (3.10) is described in Corollary 3.6.
For , and , consider the system of equations
(3.13a)  
(3.13b) 
whose solution consists of the points and the Lagrange multipliers .
Theorem 3.4 (Solvability, Lagrange multiplier estimates).
Let Assumptions 2.2 and 2.3 be fulfilled, let be the constant given by (2.7) and set , . Let be a continuous function satisfying
(3.14) 
Then there exists an open neighborhood of such that, for all , the system (3.13) has a unique solution
(3.15) 
in , where is the open ball in with radius and center zero. The solution (3.15) depends continuously on , and the Lagrange multiplier can be represented as
(3.16) 
with remainder terms , depending continuously on . Moreover, there exist constants such that
(3.17) 
for all .
Proof.
For the sake of clarity we divide the proof into several steps: After defining a suitable neighborhood of in Step 1, we show in Step 2 the existence of a solution to the system (3.13a) as well as the first estimate in (3.17). The uniqueness of this solution is verified in Step 3. In Step 4 we show the existence of a unique solution to the system (3.13b). The decomposition (3.16) of and the corresponding estimates in (3.17) are derived in Step 5. In the final Step 6 we prove that and depend continuously on . For the sake of readability, we use the notation and for the minimum and maximum of numbers .
Step 1: We begin by choosing an open neighborhood of is such a way that
(3.18) 
Note that the term on the left hand side of (3.18) is continuous as a function of and equals zero for . Thus we can define as the preimage of the open interval with respect to this function. The inequality (3.18) is used in Step 2 and in Step 5. We remark that similar assumptions are used in the analysis of numerical schemes for deterministic DAEs, see, e.g., [9, Theorem VII.4.1].
Unless stated otherwise, we assume throughout this proof that is given and fixed.
Step 2: In order to find a solution to (3.13a), we use a homotopy ansatz as commonly used in the DAE context, see for instance [8, Theorem 4.1]. It is particularly fruitful in our setting due to the presence of the truncation function . Note that the system (3.13a) is not influenced by . Our goal is to apply the globalized implicit function theorem (gIFT), see Theorem B.2, to the function defined by
(3.19) 
where is an arbitrary small positive number. To this end, we are going to show that, for all with , the Jacobian matrix is invertible and
(3.20) 
Since , Theorem B.2 (gIFT) then implies that there exists a continuously differentiable function such that and for all . Thus, we obtain a solution to (3.13a) by setting . In view of (3.20) and the identity we also have the inequality