NUMERICAL PRESERVATION OF LONG-TERM DYNAMICS BY STOCHASTIC TWO-STEP METHODS

. The paper aims to explore the long-term behaviour of stochastic two-step methods applied to a class of second order stochastic diﬀerential equa- tions. In particular, the treatment focuses on preserving long-term statistics related to the dynamics of a linear stochastic damped oscillator whose velocity, in the stationary regime, is distributed as a Gaussian variable and uncorrelated with the position. By computing the solution of a very simple matrix equal- ity, we a-priori determine the long-term statistics characterizing the numerical dynamics and analyze the behaviour of a selection of methods.

1. Long-term dynamics of a damped stochastic oscillator. Let us consider a damped stochastic oscillator related to the motion of a particle subject to both a deterministic and a stochastic forcing term described by the following second order differential equationẍ = f (x) − ηẋ + εξ(t), where f (x) is a deterministic forcing term related to a potential V(x) by the relation f (x) = −V (x), η is the damping parameter, ξ(t) is the stochastic forcing term of amplitude ε. As highlighted by [5,6], such a differential problem can be recast in terms of position and velocity of the particle by the following system of coupled differential equations where X(t) is a stochastic process describing the position of the particle at time t, whose time dynamics is described by a deterministic differential equation, V (t) is a stochastic process describing the velocity of the particle at time t and varies in time according to an Itô stochastic differential equation (SDE), being W (t) a Wiener process. More details on Itô SDEs and their numerical treatment can be found, for instance, in [9,10] and references therein.
We are now interested in underlining some long-term properties of this dynamical system. Therefore, we remind that the probability density at time t is given by and analyze the stationary density that, related to the solution of (2), has the following analytical form [5,6] where the constant N 0 can be computed by the condition Thus, the long-time behaviour of the stochastic dynamical system described by (2) has a Gaussian distributed velocity, which is totally uncorrelated with the position of the particle. If the forcing term is linear, problem (2) assumes the form with g > 0 constant value, and the stationary density (3) becomes Therefore, the following long-term statistics characterize the dynamical system (4) We collect these values in the following correlation matrix which highlights the main features of the long-term behaviour of (4). This paper, following the path of structure-preservation in stochastic dynamical systems along numerical solutions (see [1,3,4,11,12] and references therein), aims to give answer to the following question: how accurately is the correlation matrix (6) reproduced by time-stepping methods for SDEs? In other terms, the goal of the paper is to analyze the conservation of the structure of the stochastic damped oscillator (4) along the numerical dynamics carried out by selected stochastic methods, extending the analysis provided by [5,6] for one-step methods. The considered family of methods, i.e. the family of indirect stochastic two-step methods, is introduced in Section 2; the long-term behaviour of these methods when applied to (4) is then discussed in Section 3; some conclusions are provided in Section 4.

2.
Indirect stochastic linear two-step methods. The starting point of our treatise, for a general system of Itô SDEs [9,10] dX(t) = f (X(t))dt + g(X(t))dW (t) is given by the family of two-step stochastic methods investigated in [2] and references therein. The Wiener increments ∆W i and ∆W i−1 in (7) are normally distributed with expected value 0 and variance h. Two-step methods (7), specialized to (4), give rise to the family of methods which, following the notation introduced by Henrici [8], will be denoted as indirect two-step methods (ITS methods) in the remainder of the treatise. Their indirect nature is due to the fact that ITS methods (8) are applied to the first order version (4) of the second order problem (1). Let us compute X i+1 from the first equation in (8), obtaining We replace above expression for X i+1 in the second equation in (8), leading to where We finally replace (10) in (9), obtaining Finally, in compact notation we get

RAFFAELE D'AMBROSIO, MARTINA MOCCALDI AND BEATRICE PATERNOSTER
with The coefficients of the method (12), given by the matrices R 0 (h), R 1 (h) and the vectors r 0 (h), r 1 (h), are then collected into the Butcher tableau which will be used in the remainder of the treatise to denote any specific method (12).
3. Preservation of the correlation matrix. We now aim to analyze the properties of the numerical correlation matrix with being X n and V n numerical solutions of (4) generated by the ITS method (12). The following result holds.
Theorem 3.1. The correlation matrix Σ that results from ITS method (12) satisfies the matrix equation Proof. We know that problem (4) has a stationary density, given by (5), which is Gaussian with expected value 0 and correlation matrix Σ given by (6). Since the advancing rule given by the ITS method (12) is a linear transformation, certainly also the corresponding numerical solution has a stationary density which is Gaussian with expected value 0 and correlation matrix Σ given by (14). Then, in the stationary regime, is still Gaussian with expected value 0 and correlation matrix is Gaussian with expected value 0 and correlation matrix R 0 (h) ΣR 0 (h) T . Moreover, proceeding similarly for the Wiener increments in (12), in the stationary regime we have that r 1 (h)∆W i has correlation matrix r 1 (h)r 1 (h) T h and r 0 (h)∆W i−1 has correlation matrix r 0 (h)r 0 (h) T h. Then, regarding (12) in terms of the correlations in the stationary regime here computed gives the thesis.
We observe that the matrix equality (15) is the two-step extension of the analogous one-step version, derived in [5,6]. It can be treated in a twofold way: • as a matrix equality whose unknown is Σ. In this way, we can a priori compute the entries of the numerical correlation matrix Σ and analyze the long-term properties of a given method; • by imposing Σ = Σ in (15), we can employ it as a matrix equality whose unknowns are some of the coefficients of the method. This provides a constructive tool for methods (12) that would exactly preserve, within round-off error, the matrix Σ along the numerical dynamics. However, constructing new methods is out of the scope of this paper and, therefore, we will only use (15) to a priori analyze the long-term properties of a selection of ITS methods (12), as highlighted in the following subsection.
3.1. Analysis of selected methods. We now aim to analyze the long-term properties of some known methods within the family of ITS methods (12), by a priori computing the statistics contained in the numerical correlation matrix (14). Such analysis strongly relies on Theorem (3.1) and, in particular, on solving the matrix equation (15).

Euler-Maruyama method
We first analyze long-term properties of the Euler-Maruyama method [9] that, specialized to problem (4), assumes the form Recast in the compact notation (12), the corresponding Butcher tableau (13) assumes the form We solve the matrix equation (15) with respect to the matrix Σ defined by (14), whose entries are given by .

RAFFAELE D'AMBROSIO, MARTINA MOCCALDI AND BEATRICE PATERNOSTER
We observe that above values are coherent with those derived in [5,6]. Moreover, thought as functions of the stepsize h, they satisfy the following properties . Figure 1 shows the patterns of the errors | σ 2 x,EM − σ 2 x |, | σ 2 v,EM − σ 2 v | and | µ EM − µ| over the damping coefficient η. It is visible from the figure that, the more the problem is damped, the more the error in the long term mean-square of the velocity remains constant. The error associated to µ EM decreases when the problem becomes more damped. As a consequence, in the long-term, numerical position and velocity of the particle computed by the Euler-Maruyama method tend to become less correlated, in accordance with the corresponding property of the continuous problem. Moreover, also the error on σ 2 x,EM decreases for increasing values of η.

Trapezoidal method
Let us analyze the properties of the trapezoidal method applied to (4) Regarded in the compact notation (12), the corresponding Butcher tableau (13) assumes the form Through the matrix equation (15), solved with respect to Σ defined by (14), we obtain the following long-term statistics . Figure 2 shows the patterns of the errors | σ 2 x,T RAP − σ 2 x |, | σ 2 v,T RAP − σ 2 v | and | µ T RAP − µ| over the damping coefficient η. One can recognize that, the more the problem is damped, the more the long term mean-square of the velocity and the position are approximated with low accuracy. The error associated to µ T RAP remains essentially constant when the problem becomes more damped. Hence, in the long-term, numerical position and velocity of the particle computed by the trapezoidal method mostly remain uncorrelated, as desired.

A BDF method
We now consider the following BDF method [2] applied to (4) corresponding to the following Butcher tableau Through the matrix equation (15), solved with respect to Σ defined by (14), we obtain the following long-term statistics σ 2 x,BDF = 9ε 2 h 2g(9 + 2gh 2 + 6hη) , . Figure 4 shows the patterns of the errors | σ 2 x,BDF − σ 2 x |, | σ 2 v,BDF − σ 2 v | and | µ BDF − µ| over the damping coefficient η. One can recognize that, the more the problem is damped, the more the long term mean-square of the velocity and the position are approximated with low accuracy, with slowly decreasing errors. The error associated to µ BDF remains essentially constant when the problem becomes more damped. Hence, in the long-term, numerical position and velocity of the particle computed by the BDF method mostly remain uncorrelated, as desired.

4.
Conclusions. The treatise presents a long-term analysis of two-step methods for second order stochastic differential equations. In particular, we have focused on the problem of preserving the long-term properties of a linear stochastic damped linear oscillator, whose position and velocity are uncorrelated in the long-term and the velocity has Gaussian distribution. The tool we have introduced, following [5,6], is given by the simple matrix equation (15), which enables to perform an a-priori analysis of the long-term properties of two-step stochastic methods (12). As shown on selected examples of methods, it is not automatic that a method preserves the long-term properties of the given problem and methods have to be accurately chosen to achieve the prescribed behaviour. Following this path, future investigations will be oriented to employing Equation (15) to derive measure exact methods within the class (12), exploiting the larger number of degrees of freedom in the method with respect to the one-step case described in [5,6] in order to achieve a better balance among accuracy, stability and conservation of qualitative properties of the continuous problem. Moreover, it is relevant to analyze the ability of stochastic multistep methods in preserving other types of properties of stochastic differential problems: this is the case, for instance, of stochastic Hamiltonian problems [3,4], where a long-term behaviour in the mean-square of the Hamiltonian function is a-priori known and it looks relevant to merge this information in the numerical scheme, in order to accurately reproduce the same property along the numerical dynamics over long times.