# Symplectic Integrators to Stochastic Hamiltonian Dynamical Systems Derived from  Composition Methods

Authors: Tetsuya Misawa
Journal: Mathematical Problems in Engineering
 (2010)

---

## Abstract

“Symplectic” schemes for stochastic Hamiltonian dynamical systems are formulated through “composition methods (or operator splitting methods)” proposed by Misawa (2001). In the proposed methods, a symplectic map, which is given by the solution of a stochastic Hamiltonian system, is approximated by composition of the stochastic flows derived from simpler Hamiltonian vector fields. The global error orders of the numerical schemes derived from the stochastic composition methods are provided. To examine the superiority of the new schemes, some illustrative numerical simulations on the basis of the proposed schemes are carried out for a stochastic harmonic oscillator system.


---

## Body

## 1. Introduction

The theory of Hamiltonian dynamical systems is understood as one of the fundamental tools for the description of dynamical phenomena treated in physics, engineering, and economics. To investigate the behavior of such a system, numerical simulation is often carried out. Then, as one of the important subjects in the numerical approach, attention has been paid to how to construct numerical schemes preserving Hamiltonian dynamical structures, that is,  “symplectic integration schemes" (Hairer et al. [1] and references therein). The symplectic schemes are superior to numerical realization of the phase flows of Hamiltonian dynamical systems on long time intervals, and hence one can apply them to various dynamical systems described in the framework of classical Hamiltonian mechanics.As an extension of this topic, Milstein  et al. [2] have proposed symplectic integrators for  “stochastic Hamiltonian dynamical systems (e.g., Misawa [3])" which are governed by stochastic differential equations (SDEs) (Ikeda and Watanabe [4]). The stochastic Hamiltonian dynamical systems are viewed as “open" Hamiltonian systems perturbed by random fluctuations from an external world; from the viewpoint of finance, such systems may also be regarded as dynamical ones describing some risky assets. For the stochastic systems, they have formulated  “symplectic schemes of Runge-Kutta type" which are analogous to those in deterministic cases. On the other hand, symplectic integrators for deterministic Hamiltonian systems are often produced through the so-called  “composition methods (operator splitting methods)", which are formulated on the basis of Lie algebraic theory and Lie group theory (McLachlan and Quispel [5]). The method is outlined as follows: Let X denote vector fields on some space with coordinates x, with flows exp(tX), that is, the solutions of differential equations of the form ẋ(t)=X(x) are given in the form x(t)=exp(tX)(x(0)). Assume that one can write X=A+B in such a way that exp(tA) and exp(tB) can both be calculated explicitly. (More generally, this can be relaxed by approximations of the exponential maps.) In the most elementary case, the method gives the approximation for x(t) through(1.1)x̃(t)=exp(tA)exp(tB)(x)=x(t)+O(t2);
the last equality is obtained by using the Baker-Campbell-Hausdorff (BCH) formula, which is well known in Lie algebraic theory. The advantage of this method is that geometric properties of the true flow exp(tX) are often preserved. If X, A, and B are Hamiltonian vector fields, then both true flow and the approximating flow are symplectic. Therefore, it seems to be quite natural to derive symplectic integrators for stochastic Hamiltonian systems by using these methods.The purpose of the present paper is to propose some symplectic numerical integration schemes for stochastic Hamiltonian dynamical systems. Then, on the basis of the author's results on composition methods for SDEs (Misawa [6]), the numerical approximation errors of the new stochastic schemes are estimated.This paper is organized as follows. In Section2, we first review the composition methods for SDEs and the theorem on error estimation to the schemes. In Section 3, some symplectic schemes by composition methods are proposed for stochastic Hamiltonian dynamical systems described by SDE with a one-dimensional Wiener process. The approximation error of numerical solutions through the schemes is also estimated. In Section 4, in order to examine the superiority of the new schemes, we compare the numerical solutions from the proposed schemes with those from the stochastic Euler-Maruyama scheme through an illustrative example. Some concluding remarks are given in the end of the section.Finally, we would like to quote here another approach to composition methods for SDEs by K. Burrage and P. M. Burrage [7]. Furthermore, the other contributions related to our topic can also be found in the recent articles by Malham and Wiese [8] and Bou-Rabee and Owhadi [9].


## 2. Composition Methods for Numerical Integration of SDEs

Now, we start with a review of composition methods for numerical integration of SDEs proposed by Misawa [6]. Let (Ω,ℱ,P) be a probability space. We consider an n-dimensional stochastic differential equation of the Stratonovich type associated with a one-dimensional Wiener process W as follows:(2.1)dS(t)=b(S(t))dt+g(S(t))∘dW(t),
where b and g are n-dimensional C∞ vector functions, respectively. By using these functions, we define the differential operators X0 and X1 as(2.2)X0=∑i=1nbi∂i,X1=∑i=1ngi∂i,
where ∂i=∂/∂Si. According to Kunita [10], a solution of the SDE under an initial value S(0)=s can be formally represented in terms of an exponential map derived from a  “stochastic" vector field Yt as follows:(2.3)S(t)=(expYt)(s),
where Yt(ω) is the vector field for each t and almost surely (a.s.) ω∈Ω given by(2.4)Yt=tX0+W(t)X1+12(∫0tudW(u)-∫0tW(u)du)[X0,X1]+∑J;3≤|J|{∑ΔJ*cΔJWΔJ(t)}XJ,
and XJ=[⋯[Xj1,Xj2]⋯]Xjm](J=(j1,…,jm)), where [X0,X1] is a Lie bracket defined as X0X1-X1X0. Here, ∑ΔJ* is the sum for all single and double divided indices of a multi-index J, and WΔJ(t) and cΔJ are multiple Wiener-Stratonovich integrals on time interval [0,t] and the coefficients, respectively, which are determined through a rule for the divided indices of J. Moreover, |J| denotes the length of J. We note that (2.2) with (2.4) means that the solution S(t,ω) equals ϕ(1,s,ω) a.s., where ϕ(τ,s,ω) is the solution of the ordinary differential equation(2.5)dϕ(τ)dτ=Yt(ω)(ϕ(τ)),ϕ(0)=s,
regarding t and ω as parameters.Now, we proceed to a numerical integration of the stochastic equation (2.1) on the discretized time series by using the representation of solutions to SDEs mentioned above. It adopts a uniform discretization of the time interval [0,T] with step size (2.6)Δt=TN
for a fixed natural number N. Let tn=nΔt(n=0,1,2,…,N) be the nth step point. Then, for all n∈{0,…,N}, we abbreviate Sn=S(tn) and set S0=S(0)=s0. Moreover, we use ΔWn for n=0,1,…,N to denote the increments W(tn+1)-W(tn); they are independent Gaussian random variables with mean 0 and variance Δt, that is, N(0,Δt)-distributed random variables.On account of (2.3), we find a series of numerical solutions Sn(n=0,1,…,N) to SDE (2.1) by using(2.7)Sn+1=exp(YnΔt)(Sn)(n=0,1,2,…,N-1),
formally, where YnΔt is a vector field derived by replacing all the multiple Wiener integrals on time [0,t] by those on time interval Δt in (2.4). In the theory of ordinary differential equations, exp(YnΔt)(·) is often called the time-Δt map or exponential map. However, it is usually difficult to find the explicit form of the exponential map, and hence, we need to build an approximation for (2.7). Hence, we formulate a new stochastic numerical scheme as the following two procedures, which are composed of the truncation of the vector field (2.4) and a composition method (or operator-splitting method) applied to the exponential map derived from the truncated vector field.Procedure.
We replace the vector fieldYnΔt in (2.7) by a  “truncated” vector field ŶnΔt with a truncation order γ which is defined as
(2.8)ŶnΔt=∑J;1≤|J|≤γHJ(Δt)XJ,
where HJ(Δt) is a polynomial function of multiple Wiener integrals for the vector field XJ with a multi-index J. Then, we define a numerical sequence (Ŝn)n=0N through
(2.9)Ŝn+1=exp(ŶnΔt)(Ŝn)(n=0,1,…,N-1),
where Ŝ0=S(0)=s0.Procedure.
ForŜn+1=exp(ŶnΔt)(Ŝn), we apply a  “composition method" in a way analogous to that in the theory of ordinary differential equations. Suppose that the vector field ŶnΔt is of the form
(2.10)ŶnΔt=AnΔt+BnΔt,
where exp(AnΔt) and exp(BnΔt) can both be explicitly calculated through (2.15). Then an approximation to the exponential map ŶnΔt is given by exp(AnΔt)exp(BnΔt). Hence, the sequence of (Ŝn)n=0N in Procedure 1 is approximated by
(2.11)S̃n+1=exp(AnΔt)exp(BnΔt)(S̃n)(n=0,1,…,N-1),
where S̃0=S(0)=s0.
We regard(S̃n)n=0N as a numerical approximation to the exact discretized solutions (Sn)n=0N.We turn to estimate the approximation error of the numerical scheme described previously. In stochastic numerics, it is usually measured through  “global errors in the mean-square sense" for the scheme, and hence we first review the following definition (Gard [11], Kloeden and Platen [12]).Definition 2.1.
Suppose thatS(t) and (S¯n)n=0N are an exact solution and the numerical approximation solutions to SDE (2.1), respectively. Moreover, let E0,ξ be the expectation conditioned on starting at ξ at  “initial time”τ. Then the global error orderλ is defined by
(2.12)E0,s[|S¯N-SN|2]=O((Δt)2λ)(Δt↓0),
where SN=S(NΔt), NΔt=T, and |·| denotes the Euclidean norm on the space Rd. Here, the condition in the expectation (2.12) means S0=S¯0=s.It is obvious that the accuracy of the numerical scheme improves with increasing global order.We remark that in the framework of  “global error order" defined by Saito and Mitsui [13], the global order for the numerical solutions S¯n satisfying (2.12) is given by 2λ.It is evident that such an error estimation order depends on a truncated order ofYΔt and the way of decomposition of the truncated vector field ŶΔt. In [6], Misawa gives a way to calculate the global error order of stochastic numerical schemes determined through the two procedures mentioned above.Theorem 2.2 (see [6, Theorem 3.2]).
Letγ, α, and β be a truncation order of ŶnΔt described in Procedure 1, the least values of “mean square orders" of multiple Wiener integrals which are included in coefficients of the operators AnΔt and BnΔt, respectively. The mean square order s of a multiple Wiener integral I(Δt) on time interval Δt mentioned above is defined as E[|I(Δt)|2]=O((Δt)s)(Δt↓0). Then
(2.13)E[|S¯N-SN|2∣X0=X¯0=x]=O((Δt)δ-1)(Δt↓0)
holds, where δ=min(γ+2,α+β). Therefore, the global error order of stochastic numerical schemes provided by Procedures 1 and 2 is equal to (δ-1)/2.As a preparation for later discussions, we review some examples of our schemes and the global error orders given in [6].Example 2.3.
Suppose that a truncated vector fieldŶnΔt in Procedure 1 is given by
(2.14)ŶnΔt=ΔtX0+ΔWnX1.
We further set AnΔt=ΔtX0 and BnΔt=ΔWnX1 in the decomposition in Procedure 2 and assume that the explicit forms of both exponential maps for them are obtained. Then, the numerical scheme (2.11) can be put into the following form.Scheme 1.
One has(2.15)S̃n+1=exp(ΔtX0)exp(ΔWnX1)(S̃n).
In this case, γ, α and β are equal to 1, 2 and 1, respectively, and therefore the above theorem provides the global order of Scheme 1 as 1.Example 2.4.
ForŶn(Δt) in Example 2.3, we set AnΔt =ΔtX0A+ΔWnX1A and BnΔt=ΔtX0B+ΔWnX1B, where X0=X0A+X0B and X1=X1A+X1B. We assume that [X1A,X1B]≠0 and that the explicit forms of both exponential maps for them are obtained. Then, the numerical scheme (2.11) can be put into the following form.Scheme 2.
One has(2.16)S̃n+1=exp(ΔtX0A+ΔWnX1A)exp(ΔtX0B+ΔWnX1B)(S̃n).
In this case, α and β are both equal to 1; hence the global order of the scheme (2.10) for the above AnΔt and BnΔt becomes 0.5. Note that this order is equal to that of the Euler-Maruyama scheme for SDEs which is the most famous scheme in stochastic numerics.Remark 2.5.
By further manipulating the BCH formula to eliminate higher-order terms, we can obtain various schemes which give higher-order approximations to the exponential map. For example, the scheme corresponding to  “leapfrog", which is well known in deterministic numerical analysis, is given by(2.17)exp((Δt)(X+Y))=exp(ΔtY2)exp(ΔtX)exp(ΔtY2)+O((Δt)3).
In a way analogous to that in (2.11), we define a stochastic leapfrog scheme as follows:
(2.18)S̃n+1=exp(BnΔt2)exp(AnΔt)exp(BnΔt2)(S̃n)(n=0,1,…,N-1).
Then, using the BCH formula, we can produce numerical schemes having a better error order than that of scheme (2.11). In detail, see Remark 3.3 and Example 3.2 in [6].


## 3. Symplectic Integrators of Stochastic Hamiltonian DynamicalSystems

Now, we proceed to derive symplectic integrators for stochastic Hamiltonian dynamical systems in terms of composition methods mentioned in the previous section. We start with the definition of stochastic Hamiltonian dynamical systems defined in [3]. Consider the following 2ℓ-dimensional stochastic dynamical system with one Wiener process:(3.1)d(xi(t)xℓ+i(t))=(∂ℓ+iH0(x(t))-∂iH0(x(t)))dt+(∂ℓ+iH1(x(t))-∂iH1(x(t)))∘dW(t)(i=1,…,ℓ),
where x=(xk)k=12ℓ and ∂j=∂/∂xj(j=1,2,…,2ℓ), respectively. In (3.1), Hα(x)(α=0,1) are smooth scalar functions on R2ℓ. Formally, one may regard this as a Hamiltonian dynamical system (3.2)ddt(xixℓ+i)=(∂ℓ+iĤ(x)-∂iĤ(x))(i=1,…,ℓ)
with a  “randomized" Hamiltonian Ĥ given by (3.3)Ĥ=H0+H1γt,
where γt is a one-dimensional Gaussian white noise. Hence, we call (3.1) and Hα(α=0,1) an (ℓ-dimensional) stochastic Hamiltonian dynamical system and the Hamiltonian, respectively. In the following, for simplicity, we set ℓ=1 and denote x1(t) and x2(t) by q(t) and p(t), respectively, that is, we treat the following stochastic Hamiltonian dynamical system:(3.4)d(q(t)p(t))=(∂∂pH0(p(t),q(t))-∂∂qH0(p(t),q(t)))dt+(∂∂pH1(p(t),q(t))-∂∂qH1(p(t),q(t)))∘dW(t).Lettn=nΔt(n=0,1,2,…,N) be the nth step-point in a uniform discretization of the time interval [0,T] with step size Δt=T/N for fixed natural number N. Suppose that {(pn,qn)}n=0N is a series of the numerical solutions for the above Hamiltonian system derived by some numerical scheme on the above time discretization. In the same manner as that in the case of the numerical scheme on deterministic Hamiltonian systems, we call the numerical scheme for stochastic Hamiltonian dynamical systems  “symplectic", which produces numerical solutions satisfying (3.5)det(∂(pn+1,qn+1)∂(pn,qn))=1
on any step point n. In consideration of the advantages of using symplectic schemes for usual deterministic systems, we also may expect to obtain more stable numerical orbits by using these schemes than by using ordinary numerical schemes, for example, the Euler scheme.Now, we produce the symplectic schemes by using the composition methods mentioned in Section2. First, as an example of symplectic scheme on the basis of Scheme 1, we can provide the following scheme.Scheme 3.
One has(3.6)(q̃n+1p̃n+1)=exp(AnΔt)exp(BnΔt)(qnpn),
where a truncated vector field ŶnΔt shown by (2.14) and the splitting vecter fields AnΔt and BnΔt are given by
(3.7)ŶnΔt=ΔtX0+ΔWnX1=(∂pH0∂q-∂qH0∂p)Δt+(∂pH1∂q-∂qH1∂p)ΔWn,(3.8)AnΔt=ΔtX0=(∂pH0∂q-∂qH0∂p)Δt,(3.9)BnΔt=ΔWnX1=(∂pH1∂q-∂qH1∂p)ΔWn,
respectively.It is easy to see that the above composition scheme is symplectic, since for each fixed stochastic parameterω, AnΔt, and BnΔt are Hamiltonian vector fields with the Hamiltonians H0Δt and H1ΔWn on time interval [tn,tn+1], respectively, and hence the composition of their exponential maps is symplectic (McLachlan and Quispel [5]). Moreover, as mentioned in Section 2, the global error order of this scheme is 1.Next, we give an example of the symplectic scheme on the basis of Scheme2.Scheme 4.
One has(3.10)(q̃n+1p̃n+1)=exp(AnΔt)exp(BnΔt)(qnpn),
where for the same truncated vector field ŶnΔt as Scheme 3, we choose AnΔt and BnΔt as follows:
(3.11)AnΔt=ΔtX0A+ΔWnX1A=Δt(∂pH0A∂q-∂qH0A∂p)+ΔWn(∂pH1A∂q-∂qH1a∂p),(3.12)BnΔt=ΔtX0B+ΔWnX1B=Δt(∂pH0B∂q-∂qH0B∂p)+ΔWn(∂pH1B∂q-∂qH1B∂p),
under the conditions H0A+H0B=H0 and H1A+H1B=H1.
We can also prove that this scheme is symplectic, since for each fixed stochastic parameterω, AnΔt, and BnΔt are also Hamiltonian vector fields with the Hamiltonians ΔtH0A+ΔWnH1A and ΔtH0B+ΔWnH1B on time interval [tn,tn+1]. We note that the global error order of this scheme is 0.5 as mentioned in Section 2.Remark 3.1.
If concrete HamiltoniansH0 and H1 are given, one can also check the condition (3.5) for each scheme mentioned above by direct calculation.


## 4. Illustrative Examples

In this section, as an illustrative example of Schemes3 and 4, we will show the numerical solutions to a stochastic harmonic oscillator system by the symplectic schemes and compare them with the numerical results by the Euler scheme.Let us consider the following stochastic Hamiltonians dynamical system (3.4) with the Hamiltonian H0=(1/2)p2+(1/2)q2 and H1=σH0:(4.1)d(q(t)p(t))=(p(t)-q(t))dt+σ(p(t)-q(t))∘dW(t),
where σ is a constant. It is evident that this is a stochastic system with the  “conserved quantity" H0 [3, 6]; that is, the solution (q(t),p(t)) randomly moves on a circle determined by H0(q(t),p(t))=H0(q(0),p(0))=(1/2)(q(0)2+p(0)2) on p-q plane. We call such a system a stochastic harmonic oscillator system.Now, in order to investigate the advantage of our symplectic schemes, we will produce numerical solutions to this system by the new schemes and the  “Euler scheme" [11–13]. Then, by straightforward calculation of the definitions of the schemes together with the exponential maps (2.5), the Euler scheme, Schemes 3 and 4 under setting H0A, H1A, H0B, and H1B as H0A=(1/σ)H1A=(1/2)p2 and H0B=(1/σ)H1B=(1/2)q2 for this system are given as follows.Euler Scheme:
(4.2)(qn+1pn+1)=(qnpn)+(pn-12qn-qn-12pn)Δt+(pn-qn)σΔWn.Scheme 5 (Scheme3 for (4.1)).
One has(4.3)(qn+1pn+1)=(cos(Δt)-sin(Δt)sin(Δt)cos(Δt))(cos(σΔWn)-sin(σΔWn)sin(σΔWn)cos(σΔWn))(qnpn).Scheme 6 (Scheme4 for (4.1)).
One has(4.4)(qn+1pn+1)=(p̂n(Δt+σΔWn)+q̂np̂n),
where
(4.5)(q̂np̂n)=(qn-qn(Δt+σΔWn)+pn).
Here we note that we can directly examine the symplectic condition (3.5) which holds for Schemes 5 and 6mentioned above.In these schemes,ΔWn are numerically realized in terms of the independent N(0,1) random numbers γn as [12](4.6)ΔWn=γnΔt.Moreover, we set σ, Δt, and N as σ=0.2, Δt=0.05, and N=1000, respectively, and we choose the initial values of p0 and q0 as 1 and 0, respectively.Figure1 denotes the numerical time series of “p" by these schemes, and Figures 2(a)–2(c) display sample paths of numerical solutions from Euler scheme andSchemes 6 and 5, respectively. Here, as mentioned, it should be remembered that the exact solution to the original stochastic system should run on the circle H0(q(t),p(t))=(1/2)(q(0)2+p(0)2)=1, that is, a unit circle on the p-q plane. Figures 1 and 2(a), however, indicate that a series of the numerical solutions by the Euler scheme gradually goes far from the circle. In contrast, Figures 1, 2(b), and 2(c) show that the numerical solutions by our symplecticSchemes 5 and 6move around the circle stably, and especially, a series of the solution byScheme 5runs completely on the circle. Thus, as in the deterministic case, the new stochastic symplectic schemes also numerically preserve the Hamiltonian dynamical structure; this is the superiority of the new stochastic schemes over the ordinary stochastic schemes.Figure 1
Numerical sample paths of “p” of a stochastic harmonic oscillator system (SHOS) via Euler scheme, Composition 1 (Scheme 6) and Composition 2 (Scheme 5) (dt = 0.05, sigma = 0.2).(a) A numerical sample path of SHOS in phase space via Euler scheme. (b) A numerical sample path of SHOS in phase space via Composition 1 (Scheme6). (c) A numerical sample path of SHOS in phase space via Composition 2 (Scheme 5).
(a)(b)(c)Finally, we give some concluding remarks on our method.(1)
In this paper, we have treated the stochastic Hamiltonian dynamical systems with a one-dimensional Wiener process. However, it often happens that the error order of a numerical method collapses, if there is more than one Wiener process. Hence, it would be important to investigate the error orders of our schemes in the case of stochastic Hamiltonian dynamical systems with a  “multidimensional Wiener process". Therefore, we should be able to improve our composition methods and offer new schemes with high order for SDEs with a multi-dimensional Wiener process.(2)
Moreover, we have referred only to a single composition method to SDEs so far. It seems that one would readily be able to extend the method to a multi-composition one, since such compositions are obtained by iterations of single compositions. Using the multicompositions, we may easily produce a symplectic scheme for the more complicated stochastic Hamiltonian systems.In future papers, the topics mentioned above will be reported.



---

*Source: 384937-2010-08-09.xml*