SYMPLECTIC INTEGRATION OF HAMILTONIAN SYSTEMS WITH ADDITIVE NOISE

Hamiltonian systems with additive noise possess the property of preserving symplectic structure. Numerical methods with the same property are constructed for such systems. Special attention is paid to systems with separable Hamiltonians and to second-order differential equations with additive noise. Some numerical tests are presented.


Introduction.
The problem of preserving integral invariants in numerical integration of deterministic differential equations is of great current interest (see, e.g., [13,3,5,12,15,16] and references therein).The phase flows of some classes of stochastic systems possess the property of phase-volume preservation or possess the stronger property of preserving symplectic structure (symplecticness) [2,8].For instance, Hamiltonian equations with additive noise make up a rather wide and important class of systems having these properties.In the present paper, we construct special numerical methods which preserve symplectic structure in such stochastic systems.
Here P, Q, f, g, σ r , γ r are n-dimensional column-vectors.Consider the differential 2-form We are interested in systems (1.2) such that the transformation (p, q) → (P, Q) preserves symplectic structure [1] as follows: dP ∧ dQ = dp ∧ dq, (1.4) i.e., when the sum of the oriented areas of projections of a two-dimensional surface onto the coordinate planes (p 1 , q 1 ), . . ., (p n , q n ) is an integral invariant.As a consequence, all external powers of the 2-form are invariant for such systems as well.The case of the nth external power gives preservation of phase volume.To avoid confusion, we note that the differentials in (1.2) and (1.4) have different meanings.In (1.2), P, Q are treated as functions of time, and p, q are fixed parameters, while differentiation in (1.4) is made with respect to the initial data p, q.
A one-step mean-square approximation X(t + h; t, x), t 0 ≤ t < t + h ≤ t 0 + T of (1.1) is constructed, depending on t, x, h, and {w 1 (ϑ)−w 1 (t), . . ., w m (ϑ)−w m (t); t ≤ ϑ ≤ t+h}.Using the one-step approximation, we recurrently obtain the approximation For simplicity, we take t k+1 − t k = h = T/N.Note that X 0 may be random.A meansquare method for (1.2) based on the one-step approximation P = P (t+h; t, p, q), Q = Q(t + h; t, p, q) preserves symplectic structure (said to be symplectic or Hamiltonian) if In section 3, we construct some symplectic methods of mean-square order 1 and order 3/2 for general Hamiltonian systems with additive noise.We propose more effective methods for systems with Hamiltonians of a special form.We consider the case of separable Hamiltonians H(t, p, q) = V (p) + U (t, q) in section 4 and the case of Hamiltonians H(t, p, q) = 1 2 p M −1 p + U (t, q) with M a constant, symmetric, invertible matrix in section 5.In addition, symplectic methods for Hamiltonian systems with small additive noise can be found in the preprint [10].These methods are constructed using the results from [11].Let us emphasize that all the derived methods are efficient with respect to simulation of the used random variables.Section 6 is devoted to numerical experiments which demonstrate the superiority of the proposed symplectic methods over long periods of time in comparison with nonsymplectic methods.
Symplectic methods for Hamiltonian systems with multiplicative noise are in progress.Hamiltonian weak methods will be considered in later publications.
2. Preservation of symplectic structure.Consider system (1.2).Our aim is to indicate a class of stochastic systems, which preserve symplectic structure, i.e., satisfy condition (1.4).
Using the formula of change of variables in differential forms, we obtain Hence, the phase flow of (1.2) preserves symplectic structure if and only if Introduce the notation For a fixed k, we obtain that P ik p , Q ik p , i = 1, . . ., n, obey the following system of SDEs: Analogously, for a fixed k, P ik q , Q ik q , i = 1, . . ., n, satisfy the system • dw r , (2.5) The coefficients in (2.4) and (2.5) are calculated at (t, P, Q) with P = P (t) = [P 1 (t; t 0 , p, q), . . ., P n (t; t 0 , p, q)] , where Consider condition (2.1).Clearly, It is not difficult to check that if the functions f i (t, p, q), g i (t, p, q) are such that Then condition (2.3) is fulfilled if and only if Using the same arguments again, we prove that the relations (2.9)-(2.10)ensure this condition as well.
Corollary 2.2.The phase flow of a Hamiltonian system with additive noise preserves symplectic structure.

Symplectic mean-square methods for general Hamiltonian systems with additive noise.
In this section we consider the general Hamiltonian system with additive noise where P, Q, f, g, σ r , γ r are n-dimensional column-vectors, w r (t), r = 1, . . ., m, are independent standard Wiener processes, and H(t, p, q) is a Hamiltonian.
In this paper we suppose that H is sufficiently smooth and that first derivatives of f and g are bounded.Moreover, we also require boundedness of some higher order derivatives of f and g.At the same time, we emphasize that these conditions are not necessary and the methods obtained are more widely applicable.
The following lemma guarantees the unique solvability of (3.3) with respect to P, Q for any P k , Q k , and sufficiently small h.Lemma 3.1.Let F (x; c, s) be a continuous d-dimensional vector-function depending on x ∈ R d , c ∈ R d , and s ∈ S, where S is a set from an R l .Suppose F has the first partial derivatives ∂F i /∂x j , i, j = 1, . . ., d, which are uniformly bounded in Then there is an h 0 > 0 such that the equation The solution of (3.5) can be found by the method of simple iteration with an arbitrary initial approximation.
The proof of this lemma is not difficult and is therefore omitted.The next lemma is true for system (3.1) with arbitrary f and g (i.e., f and g may not obey condition (3.2)).
Lemma 3.2.The mean-square order of the methods The proof is based on comparison of the one-step approximation of methods (3.3)-(3.4)with the one-step approximation of the Euler method (see details in [10]).
As remarked in our introduction, the method based on a one-step approximation P = P (t + h; t, p, q), Q = Q(t + h; t, p, q) preserves symplectic structure if its one-step approximation satisfies condition (1.6).The one-step approximation P , Q of methods The relations for P, Q coincide with those for the one-step approximation corresponding to the deterministic symplectic method [16].Therefore, methods (3.3)-(3.4)are symplectic as well.From here and Lemma 3.2, we get the theorem.
Theorem 3.3.Methods (3.3)-(3.4)for system (3.1)-(3.2) preserve symplectic structure and have the first mean-square order of convergence.Now consider another generalization of the same family of deterministic symplectic methods to system (3.1), For sufficiently small h, equations (3.6) are uniquely solvable with respect to P k+1 , Q k+1 according to Lemma 3.1.
Proof.Comparing the one-step approximation of method (3.6) with the one-step approximation of the Euler method, one can establish that the mean-square order of method (3.6) is equal to 1. Now we check symplecticness of the method.Let P , Q be the one-step approximation corresponding to method (3.6).Introduce We have The relations for P , Q coincide with the one-step approximation corresponding to the symplectic deterministic method.Therefore,

Methods of order 3/2.
For i = 1, . . ., s, consider the relations where ϕ i , ψ i , η, ζ do not depend on p and q, the parameters α ij and β i satisfy the conditions and c i are arbitrary parameters.
If Proof.We generalize the proof of Theorem 6.1 in [13].Introduce the temporary notation ) and form the exterior products Now using (3.11)-(3.12),find the expressions for df i ∧ dq and dp ∧ dg j and substitute them in (3.10) as follows: The last term in the right-hand side vanishes owing to (3.9).Consider the second term in the right-hand side of (3.13).We have Taking into account that the exterior product is skew-symmetric and f and g satisfy condition (1.5), it is not difficult to see that this expression vanishes.Returning to (3.13), we obtain The next lemma is used in Theorem 3.7 for the Hamiltonian system (3.1)-(3.2).However, this lemma is of interest for arbitrary systems with additive noise as well (see Remark 3.2 below).So, consider the system with additive noise and introduce the following parametric family of one-step approximations for (3.14): (3.15) where ∆w r := w r (t + h) − w r (t), ) and the parameters α, λ 1 , λ 2 , µ 1 , µ 2 are such that For example, the following set of parameters satisfies (3.17)-(3.18): Note that random variables ∆w r and J r0 are of the same mean-square order O(h 1/2 ).
Proof.Due to properties of the Wiener process and Itô integrals, we get (3.20) Let ∆X i := X i − x, i = 1, 2. We have (3.21) Expand (3.15) as follows: Using (3.20)-(3.21),one can obtain where the coefficients and their derivatives are calculated at (t, x) and where ρ satisfies the same relations as ρ (see (3.26)).
Remark Remark 3.2.In a way similar to how method (3.15) was obtained, it is not difficult to construct new explicit RK methods of mean-square order 3/2 for an arbitrary system of differential equations with additive noise (3.14).For instance, we obtain the following explicit RK method of order 3/2 for (3.14): Note that if we apply this method, as well as any other explicit RK method, to (3.1)-(3.2), it will not preserve symplectic structure.Now consider the parametric family of methods for the Hamiltonian system with additive noise (3.1), where the parameters α, λ 1 , λ 2 , µ 1 , µ 2 satisfy (3.17)-(3.18).Under σ r ≡ 0, γ r ≡ 0, r = 1, . . ., m, method (3.30) is reduced to the wellknown second-order symplectic RK method for deterministic Hamiltonian systems (see, e.g., [13, p. 101]).Let us note that, using this deterministic method with α = 0 (the midpoint rule), another implicit 3/2-order method for Hamiltonian systems with noise was proposed in [17]-however, without preserving symplectic structure.
The one-step approximation corresponding to method (3.30) is of the form (3.15).Therefore, due to Lemma 3.6, method (3.30) is of mean-square order 3/2.Moreover, this one-step approximation is of the form (3.7) with s = 2 and This set of parameters α ij , β i , i, j = 1, 2, satisfies conditions (3.9).Then due to Lemma 3.5, the method (3.30) is symplectic.Thus we obtain the following theorem.
Remark 3.3.Formula (3.30) contains the random variables ∆ k w r (h), (J r0 ) k , (I 0r ) k , the joint distribution of which is Gaussian.They can be simulated at each step by 2m independent N (0, 1)-distributed random variables ξ rk and η rk , r = 0, . . ., m as follows: So, method (3.30) can be rewritten in the constructive form.

Symplectic mean-square methods in the case of a separable Hamiltonian.
In this section we consider the Hamiltonian system with additive noise (3.1), a Hamiltonian of which has the special structure H(t, p, q) = V (p) + U (t, q).(4.1)Such Hamiltonians are called separable.We note that it is not difficult to consider a slightly more general separable Hamiltonian H(t, p, q) = V (t, p) + U (t, q), but we restrict ourselves here to (4.1).In the case of separable Hamiltonian (4.1), system (3.1)takes the partitioned form where Obviously, the implicit symplectic methods from the previous section can be applied to the partitioned system (4.2), and these methods take a more simple form in this case (we do not write them here).We recall that there are no explicit symplectic RK methods for the general system (3.1)-(3.2).However, for the partitioned system (4.2) it is possible to construct explicit symplectic methods just as in the deterministic case [16,12,13].

Explicit first-order methods.
On the basis of the known family of deterministic partitioned RK (PRK) methods [16,12,13], we construct the family of explicit partitioned methods for stochastic system (4.2) as follows: Since the expressions for dP k+1 , dQ k+1 coincide with those corresponding to the deterministic symplectic method, then method (4.3) is symplectic.Further, it is not difficult to show that method (4.3) has the first mean-square order of accuracy.As a result, we obtain the following theorem.Theorem 4.1.The explicit partitioned method (4.3) for system (4.2) preserves symplectic structure and has the first mean-square order of convergence.
Remark 4.1.In the special cases of α = 0 and α = 1, method (4.3) takes a more simple form.In these cases it requires evaluation of each of the coefficients f, g only once per step.
Remark 4.2.It is possible to propose other symplectic first-order methods for (4.2) on the basis of the same deterministic PRK methods as above.For instance, the method is of first mean-square order and symplectic.
Introduce the relations (cf.(3.7)-(3.8)) where ϕ i , ψ i , η, ζ do not depend on p and q; the parameters α ij , αij , β i , and βi satisfy the conditions and c i are arbitrary parameters.
It is not difficult to see that method (4.8)-(4.9)has the form of (4.5)-(4.6)and that its parameters satisfy conditions (4.7).Then, Lemma 4.2 implies that this method preserves symplectic structure.Using ideas of the proof of Lemma 3.6, we establish that method (4.8)-(4.9)with (4.10)-(4.11) is of mean-square order 3/2.We have thus proved the following theorem.
Remark 4.3.Attracting other explicit deterministic second-order PRK methods from [13,16], it is possible to construct other explicit symplectic methods of order 3/2 for system (4.2).For instance, by swapping the roles of p and q in method (4.8)-(4.9),we can obtain another 3/2-order symplectic PRK method.

Symplectic methods in the case of Hamiltonian H(t, p, q) =
1 2 p M −1 p + U (t, q).Here we propose symplectic methods for the Hamiltonian system (4.2), when γ r (t) = 0 and the separable Hamiltonian has the special form with M a constant, symmetric, invertible matrix (i.e., the kinetic energy V (p) in (4.1) is equal to 1  2 p M −1 p).In this case, system (4.2) reads This system can be written as a second-order differential equation with additive noise Clearly, the symplectic methods from sections 3 and 4 can be applied to (5.2)-(5.3).Due to specific features of this system, these methods have a more simple form here.Moreover, one can prove that method (4.8)-(4.9),when applied to (5.2)-(5.3), is of mean-square order 2. In this section we restrict ourselves to new explicit methods of orders 2 and 3.
Other methods of order 2 are given in [10].In [14] a symplectic method of mean-square order 1 for (5.2)-(5.3) is proposed on the basis of the Störmer-Verlet method.There, some physical applications of stochastic symplectic integrators are also discussed.
Proof.It is not difficult to check that dP k+1 ∧dQ k+1 = dP 3 ∧dQ 3 .The expression for dP 3 ∧ dQ 3 coincides with that corresponding to the deterministic third-order symplectic method from [15,16,13].This implies that method (5.10)-(5.11) is symplectic.By Lemma 5.2 we get that the method has mean-square order 3.

Numerical tests.
In this section we consider the following Hamiltonian system with additive noise: We have, for the solution X = (X 1 , X 2 ) of (6.1), where When applied to (6.1), the explicit symplectic method (4.4) with α = 1 takes the form Method (6.3) can be written as Our aim is to analyze propagation of the error r k := X k − X(t k ).We get Proposition 6.1.Suppose T and h are such that T h 2 is sufficiently small.Then for k = 0, 1, . . ., N, T = Nh, the following inequality holds: , and the columns of the matrix G are eigenvectors of H corresponding to the eigenvalues λ 1 , λ 2 .We write the matrices Λ and G in the form diag(e kiϕ − e kih , e −kiϕ − e −kih ).It is not difficult to show that (the norms of matrices are Euclidean) Taking into account that ϕ = arcsin h 1 − h 2
Consequently, the Euler method can be used on the interval [0, T E ] if T

3/2
E h is sufficiently small.Due to Proposition 6.2, the error of the symplectic method (6.3) on [0, T S ] with T S = T 2 E is equal to O(T E h + T 3 E h 2 ); i.e., the symplectic method is applicable on longer time intervals than the Euler method.Of course, the Euler method possesses properties worse than the symplectic method since the absolute values of the eigenvalues of H are greater than 1.
Finally, consider the following optimal method from [9, p. 62] (this method also uses only the increments ∆ k w as the information regarding w(t) but it uses this information optimally): Xk+1 = Ĥ Xk + vk (6.14) = cos h sin h − sin h cos h Xk + Evidently, this method is symplectic.Also, as Ĥ = F, it has no error in the absence of noise.We get for its error Consequently, the error of the optimal method is estimated as O(T 1/2 h).This implies that method (6.14) is more applicable on the longer time interval [0, T O ] = [0, T 3 E ] than the symplectic method (6.3).
In the numerical tests we simulate system (6.1) by (i) the exact formulae (6.2), (ii) the symplectic method (6.3), and (iii) the Euler method (6.13).Figure 1 corresponds to the time interval [0, 128] which contains approximately 20 oscillations of (6.1) (note that the period of free oscillations of (6.1) is equal to 2π).The results clearly demonstrate that the Euler method is unacceptable for simulation of the Hamiltonian system (6.1) on a long time interval.After 10 oscillations (Figure 1) the norm of its error is already half the norm of the solution, and after 200 oscillations (see Figure 2) the amplitude of oscillations simulated by the Euler method is 50,000 times greater than the exact amplitude.
In contrast to the Euler method, the symplectic method reproduces oscillations of the system (6.1)quite accurately.After 10 oscillations (Figure 1) the norm of its error is approximately 2% of the norm of the solution.But it is more astonishing that the relations (3.7)-(3.8)coincide with a general form of s-stage Runge-Kutta (RK) methods for deterministic differential equations.It is known (see, e.g., Theorem 6.1 in [13]) that the symplectic condition d P ∧d Q = dp∧dq holds for P , Q from (3.7)-(3.8)with (3.9) and ϕ i = ψ i = η = ζ = 0. Let us check the case of arbitrary ϕ i , ψ i , η, ζ.Lemma 3.5.Relations (3.7)-(3.8)with condition (3.9) preserve symplectic structure, i.e., d P ∧ d Q = dp ∧ dq.
.25) |Eρ| = O(h 3 ), Eρ 2 = O(h 5 ).(3.26) Substituting (3.22)-(3.23) in (3.24) and using (3.17), we get X = x + m r=1 b r ∆w r + m r=1 b r 3.1.It is very unusual that the direct expansion of (3.15) does not contain the habitual term h 2 A similar term in the expansion contains some combinations of ∆w r and J r0 instead of h.This is necessary for a method conserving symplectic structure.At the same time, this new reception allows us to construct new RK methods for general (not only Hamiltonian) stochastic systems with additive noise (see the next remark).
2, we get Let T and h be such that T h 2 is sufficiently small.SupposeE |X 0 | 2 ≤ C.Then the mean-square error is estimated as k+1 ,