Exponential mean-square stability properties of stochastic linear multistep methods

The aim of this paper is the analysis of exponential mean-square stability properties of nonlinear stochastic linear multistep methods. In particular it is known that, under certain hypothesis on the drift and diffusion terms of the equation, exponential mean-square contractivity is visible: the qualitative feature of the exact problem is here analysed under the numerical perspective, to understand whether a stochastic linear multistep method can provide an analogous behaviour and which restrictions on the employed stepsize should be imposed in order to reproduce the contractive behaviour. Numerical experiments confirming the theoretical analysis are also given.

the property of the continuous problem is detected and identified through sufficient conditions, we aim to provide their numerical counterparts, equivalent to conditional stability properties for the numerical method. Before entering into the details, we first recover here the deterministic case, with special emphasis on linear and nonlinear stability properties.

Linear and nonlinear stability issues for deterministic differential equations
For a well-posed deterministic initial value problem based on the ordinary differential equation y (t) = f (y(t)) linear stability properties are classically provided with respect to the scalar linear test equation where ξ ∈ C and Re(ξ ) ≤ 0, considered for the first time by Dahlquist [6]. The solution of this simple problem remains bounded when t goes to infinity and one needs to require that the numerical solution possesses an analogous stability property to that displayed by the exact solution: this fact is at the basis of the linear stability analysis of numerical methods (see [4] and references therein).
Based on the concept of one-sided Lipschitz continuity G. Dahlquist laid the foundation of the theory of nonlinear stability in [7] (this paper has also been celebrated by J. Butcher in the review paper [3]), in order to treat situations beyond the realm of linear stability analysis. Hence, the following nonlinear test problem y (t) = g t, y(t) , t ≥ 0, is considered, with g : R × R m → R m satisfying a one-sided Lipschitz condition of the form g(t, y 1 ) − g(t, y 2 ) T (y 1 − y 2 ) ≤ 0 ( 3 ) for all t ≥ 0 and y 1 , y 2 ∈ R m . Denote by y(t) andỹ(t) two solutions to (2) with initial conditions y 0 andỹ 0 , respectively. Then it is known (see [3,9,11] and references therein) that the condition (3) implies the contractivity of the trajectories, given by the inequality where · is any norm in R m , and the corresponding problem is said to be dissipative.
Extensive research has been carried out in order to provide numerical methods generating contractive numerical solutions for dissipative problems, giving rise to the notions of AN-stability, G-stability, algebraic stability and so on (see [4,15] and references therein).

Exponential mean-square stability of stochastic differential equations
We now move to the main issue of this paper, i.e. a nonlinear Itô stochastic differential equation (SDE) for t ∈ [0, T ], where f : R n → R n , g : R n → R n×m , and W (t) is an m-dimensional Wiener process. Following [12], we assume the following conditions on the functions f and g: (i) f, g ∈ C 1 ; (ii) the drift coefficient f satisfies a one-sided Lipschitz condition, i.e. there exists μ ∈ R such that (iii) the diffusion coefficient g is globally Lipschitz, i.e. there exists L > 0 such that Under the aforementioned assumptions, the following result holds [12,14].
As a consequence of this result, the following definition is provided [12].

Definition 1
The solution of any given SDE (5) satisfying (8) with α < 0 is said to be exponential mean-square contractive.
For other possible stability issues, compare for instance [1,10,18] and references therein.
We aim to investigate the numerical counterpart of (8) when trajectories are generated by stochastic linear multistep methods, in order to provide bounds on the stepsize whose values eventually reproduce such a contractive behaviour, according to the following definition.

Definition 2
Assuming that the integration interval [0, T ] of (5) is discretized in the following grid {t n = nΔt, n = 0, 1, . . . , N, Δt = T /N}, the numerical solution of any given SDE (5) satisfying (8) with α < 0 is said to be exponential mean-square contractive if there exists ν < 0 such that In other terms, (9) provides the discrete time counterpart of (8). The article proceeds as follows: Section 2 provides inequalities of type (9) when linear multistep methods (10) are used; Sections 3 and 4 are respectively devoted to the analysis of the exponential mean-square contractive behaviour of one and two-step methods; Section 5 contains the numerical results obtained by applying the analysed methods to scalar nonlinear equations as well as to systems of equations; Section 6 is devoted to giving some concluding remarks.

Exponential mean-square stability analysis
We develop a general analysis of the exponential mean-square stability properties of the general family of stochastic two-step methods [2] In our analysis, the following lemma, reproduced from [12,13], will be particularly useful. (6) and (7),

Lemma 1 Under conditions
Then, a 1 and a 2 exist, are unique and satisfy the inequality For two-step methods (10), the following three-terms inequality holds.

Theorem 2 Under conditions (i)-(iii), any two numerical solutions {X n } n and
{Y n } n with initial values such that E|X 0 | 2 < ∞ and E|Y 0 | 2 < ∞, generated by applying a two-step method (10) with stepsize Δt satisfy where β T S and γ T S are the expressions defined for values of Δt such that 1 − 2β 2 μΔt = 0 and being Proof We consider two numerical solutions {X n } n and {Y n } n with E|X 0 | 2 < ∞ and E|Y 0 | 2 < ∞, computed by (10) with stepsize Δt. Lemma 1 leads to Squaring the right-hand side of this inequality, passing to expectations and applying assumptions (6) and (7) leads to Let us now provide estimates for the last six summands in (14): Applying the mean value theorem to the function f leads to Proceeding as above leads to -estimate of 2E| X n − Y n , X n−1 − Y n−1 |. We use the well-known inequality 2ab ≤ a 2 + b 2 , for any a, b ∈ R, leading to Hence, We proceed as above, obtaining We proceed as above, obtaining Matching (15)-(20) with (14) leads to the thesis.
This result first allows us to provide a complete analysis of the exponential meansquare properties of one-step methods. Indeed, the following result holds. (6) and (7), any one-step stochastic method

Theorem 3 Under conditions
satisfies the following exponential mean-square stability inequality with defined for values of Δt such that 1 − 2β 2 μΔt = 0.
Proof We apply the thesis of Theorem 2 to the family of one-step methods (21), with β OS defined as in (23) and, by recursion, The thesis immediately follows by assuming ν = 1 Δt log β OS .
We now aim to provide an exponential mean-square stability inequality for the complete two-step case (10). For this purpose, we need to specify the one-step method we use in order to recover the missing starting values X 1 and Y 1 . Once this aspect is clarified, the following result holds. (6) and (7), any two-step stochastic method (10) satisfies the following exponential mean-square stability inequality over a single step

Theorem 4 Under conditions
where ζ n+1 is recursively defined by the three-term formula with ζ 1 = β OS and ζ 2 = β TS β OS +γ TS . We recall β TS and γ TS are given in (12), while β OS is given in (23) and characterizes the one-step method employed to compute the missing starting values X 1 and Y 1 .
Proof The thesis immediately follows by recursively applying the thesis of Theorem 2, and applying the inequality for the first step.

Contractivity of one-step methods
We now employ the results in Section 2, especially Theorem 3, to study the exponential mean-square contractive behaviour of one-step methods (21). According to Definition 1, a one-step method is exponential mean-square contractive if ν < 0 in in (23). Let us analyse this condition on the case by case basis.

Stochastic trapezoidal rule
We now aim to study the behaviour of the stochastic trapezoidal rule when applied to an exponentially mean-square stable SDE (5) generating contractive solutions. For this method, the value of β OS in (23) is given by Let us now illustrate the above result with a numerical example, given by a scalar SDE (5) with In this case, (6) holds with μ = −4 and (7) holds with L = 1. Since α = 2μ + L = −7 < 0, the problem exhibits an exponentially mean-square contractive behaviour, according to Definition 1. In order to compute the stepsize restriction corresponding to (26), we estimate the value of M (13) as follows where X EM n is the numerical solution computed by the Euler-Maruyama method with initial value X 0 , stepsize δt, over N + 1 grid points, while Y EM n is the numerical solution computed by the Euler-Maruyama method with initial value Y 0 . The result of the computation of (29) over 1000 paths with X 0 = 1, Y 0 = 2, N = 2 15  In order to provide a numerical evidence of the sharpness of this bound, we apply the trapezoidal method with X 0 = 1 and Y 0 = 2 and obtain the patterns reported in Fig. 1, where the pointwise mean-square of the error over 1000 paths is computed and depicted, in correspondence of Δt = 75/128 and Δt = 75/64. The method applied with Δt = 75/128 shows the exponential decay of the mean-square error, while such a decay is totally lost for Δt = 75/64, confirming the proved theoretical result.

Contractivity of two-step methods
We now move to analysing the contractivity of two-step methods (10), specializing the result in Theorem 4 to a selection of methods. According to Definition 1, a twostep method is exponential mean-square contractive if η < 0 in (24), i.e. if, after performing n steps, 0 < ζ n < 1, with ζ n defined in (25). Let us analyse this condition on the case by case basis.

Adams-Moulton method
Adams-Moulton method is given by (10) with For this method, we have In order to check condition (30), we need to compute the values of ζ n for large enough values of n, assuming to employ the trapezoidal method as a starting method. As a numerical illustration, we still consider problem (5) with drift and diffusion coefficients given by (28). For n = 30, (30) leads to the following restriction on the stepsize Δt ≤ 0.87, in order to have exponential mean-square contractivity. The numerical evidence originated by applying the Adams-Moulton method with X 0 = 1 and Y 0 = 2 is reported in Fig. 2, where the pointwise mean-square of the error over 1000 paths is /computed and depicted, in correspondence of Δt = 25/64 and Δt = 125/128. The method applied with Δt = 25/64 shows the exponential decay of the mean-square error, while such a decay is totally lost for Δt = 125/128, confirming what we proved.

Milne-Simpson
Milne-Simpson method is given by (10) with For this method, we have In order to check condition (30), we need to compute the values of ζ n for large enough values of n. By symbolic manipulation we have realized that, for large enough values of n, condition (30) is never fulfilled by positive values of Δt; thus, Milne-Simpson method does not show an exponential mean-square stable character.

BDF2
BDF2 method is given by (10) with For this method, we have .
In order to check condition (30), we need to compute the values of ζ n for large enough values of n, assuming to employ the trapezoidal method as a starting method. As a numerical illustration, we still consider problem (5) with drift and diffusion coefficients given by (28). For n = 30, (30) does not lead to any stepsize restriction; in other terms, the method is nearly unconditionally stable for exponential mean-square contractivity. The numerical evidence originated by applying the Adams-Moulton method with X 0 = 1 and Y 0 = 2 is reported in Fig. 3 the exponential decay of the mean-square error is visible and, clearly, more accurate for smaller values of Δt.

Numerical experiments
As a test for the pertinence of the analysed methods in solving exponential meansquare stable SDEs (5), we provide the numerical evidence associated with the following nonlinear problems. We first consider the following nonlinear SDE [16] with 2a − 1 ≥ b 2 and b = 0 to ensure an exponential mean-square contractive behaviour, as proved in [16]. In our experiments, we choose a = 2 and b = 1. The evidence, reported in Fig. 4, shows the capability of trapezoidal, Adams-Moulton and BDF2 methods to preserve the expected mean-square contractivity. We next consider the following nonlinear system [17] dX whose solution is exponentially mean-square stable (see [17]). Also in this case, the results reported in Fig. 5 show the capability of trapezoidal and BDF2 methods to nicely retain the mean-square contractivity behaviour, for Δt = 25/64. Adams-Moulton method, subject to a more severe stepsize restriction, is not able to provide a mean-square contractivity behaviour for Δt = 25/64.

Conclusions
We have analysed the behaviour of stochastic linear multistep methods when applied to nonlinear stochastic differential (5), in presence of exponential mean-square contractivity, under the conditions highlighted in Theorem 1. The theoretical analysis has provided sharp bounds on the stepsize to be employed in order to retain numerically the property of contractivity, as described in Sections 2-4. If a starting procedure is needed (and this is the case of two-step methods), the obtained bound also depends on the property of the employed starter method. Future issues of this research, which can be assumed as a starting point in assessing a theory of nonlinear stability for stochastic numerical methods, regard the possibility of extending the notion of algebraic stability for stochastic Runge-Kutta methods [8], in order to understand its possible connection with geometric numerical integration of stochastic Hamiltonian problems [5].