Stability Issues for Selected Stochastic Evolutionary Problems : A Review

We review some recent contributions of the authors regarding the numerical approximation of stochastic problems, mostly based on stochastic differential equations modeling random damped oscillators and stochastic Volterra integral equations. The paper focuses on the analysis of selected stability issues, i.e., the preservation of the long-term character of stochastic oscillators over discretized dynamics and the analysis of mean-square and asymptotic stability properties of θ-methods for Volterra integral equations.


Introduction
This paper provides a brief review of recent results obtained in the context of stability analysis for stochastic numerical methods.The treatise is essentially twofold: regarding stability properties as the preservation of qualitative features of the continuous problem as well as the numerical preservation of stable behaviors in the solution of the continuous problem.
The first highlighted aspect, i.e., stability as numerical preservation of qualitative properties, is here framed in the context of stochastic differential equations (SDEs), with special emphasis on problems describing stochastic oscillators [1][2][3][4].The perspective we follow consists of providing a long-term analysis of numerical methods for SDEs in terms of preserving invariance laws that characterize the dynamics provided by the exact solution.A first contribution in this sense was given in [5], where the author analyzed the invariance of asymptotic laws characterizing linear stochastic systems under given discretizations.The analysis of partitioned methods for linear oscillators in the presence of additive noise has been an object of [6], while the analysis of linear second order SDEs describing damped stochastic oscillators has been provided by Burrage and Lythe in [1,2] and inspired the paper [3] giving a more general two-step framework.More recent contributions have regarded stochastic Hamiltonian problems [7,8] and stochastic oscillators with multiplicative noise [9].
The second highlighted issue, i.e., numerical preservation of stable behaviors in the solution of the continuous problem, the attention is focused on stochastic Volterra integral equations (SVIEs).For such operators, whose numerical discretization has been considered by the recent literature (see, for instance, [10][11][12] and references therein), researchers have started to provide a parallel with the classical theory of the numerical approximation of SDEs [13,14].However, as far as the authors are aware, the first contribution in developing a stability analysis in the time-stepping numerics for SVIEs is the recent paper [15], where the analysis is given both in terms of mean-square stability properties as well as on asymptotic ones.
The treatise is divided into two parts: the first one, presented in Section 2, is regarding the approximation of stochastic differential equations modeling the dynamics of damped oscillators subject to both deterministic and random forcing terms; the second part, contained in Section 3, gives a glance on stochastic ϑ-methods for stochastic Volterra integral equations and, in particular, on the analysis of the stability properties with respect to suitable test equations.

The Problem
The motion of a particle constrained by a deterministic forcing and a stochastic one, in the presence of damping, can be modeled by the following scalar second order SDE [1,2] where f (x) is a deterministic forcing, η is the damping parameter, and εξ(t) is the stochastic forcing of amplitude ε.Clearly, assuming a linear forcing f (x) = −gX(t), Equation ( 1) admits the following first order It ô formulation where X(t) and V(t) are, respectively, position and velocity of the particle at time t whose dynamics is also governed by the occurences of the Wiener process W(t) [13,14].The long-term character of the problem described by Equation ( 2) is clearly highlighted by its stationary density [1,2] where the constant N 0 is the unknown of the equality In other terms, the stochastic dynamics described by Equation (2) has a Gaussian distributed velocity, uncorrelated with the position.An effective representation of such a long-term behavior is given by the following correlation matrix where

The Methodology: Indirect Stochastic Linear Two-Step Methods
In [1][2][3], the authors have analyzed the ability of some of the most common numerical methods for SDEs in preserving the correlation matrix over long times.The analysis given in [1,2] only involves one-step methods, while that provided in [3] is focused on the larger class of indirect two-step methods (ITS methods), i.e., the family of stochastic two-step methods applied to the first order representation given by Equation (2) of the second order SDE in Equation (1).ITS methods assume the form where the matrices R 1 (h), R 0 (h) ∈ R 2×2 and the vectors r 1 (h), r 0 (h) ∈ R 2 are the characteristic coefficients, collected in the Butcher tableau The corresponding discretized correlation matrix is then given by with where X n and V n are numerical solutions of Equation ( 2) generated by the ITS method given by Equation (5).As proved in [3], Σ satisfies the matrix equation which is the most relevant tool in analyzing the long-term dynamics along the solutions computed by the ITS method defined in Equation ( 5).Indeed, for a fixed ITS method, the matrix in Equation ( 7) can be a priori computed in order to appreciate how far it is from the exact correlation matrix defined in Equation ( 4).Alternatively, for those particularly interested in developing new methods, the matrix equality (8) can be used to derive the constraints that the coefficients of an ITS method have to fulfill in order to provide a measure-exact numerical scheme [1,2].
As an example, let us provide the analysis of the famous Euler-Maruyama method [13] that, in the form of the ITS method, assumes the form Solving Equation ( 8) with respect to the Σ gives We observe that Σ = Σ + O(h), i.e., the matrix Σ associated with the Euler-Maruyama method is a first order approximation of Σ.In order to appreciate the values of the errors over the damping coefficient η, we depict the corresponding graphs in Figure 1.One can recognize that, for increasing values of the damping parameter, position and velocity become less correlated, as happens for the continuous problem defined by Equation (2).
Let us now analyze a genuine two-step method, i.e., the following Adams-Moulton method [16] where Solving Equation ( 8) with respect to the Σ gives 2, show that the more the problem is damped, the more the long-term mean-square of the velocity and the position are approximated with slowly decreasing errors.In the long-term, numerical position and velocity of the particle computed by the Adams-Moulton method are only weakly correlated, as desired.A qualitative comparison arising from Figures 1 and 2 shows that the Euler-Maruyama method better approximates the long-term mean-square of the velocity and the position, while the Adams-Moulton method better captures the long term expectation of the product of the velocity and the position.
While the Euler-Maruyama method is able to reproduce the exact correlation matrix Σ as h goes to 0, this is not the case for the Adams-Moulton method.This aspect opens a relevant question on the difference between one-step and genuine multistep methods.The fact that multistep methods do not recover the exact invariants of the continuous problem in the limit as h tends to 0 is typical of the deterministic setting: for instance, multistep methods do not recover the Hamiltonian function of Hamiltonian problems [17].We conjecture that there are sources of parasitism to be properly addressed also in stochastic linear multistep methods, but this requires a very deep analysis based on the stochastic version of the backward error analysis [17], which will be an object of future investigations.

The Problem
We now focus our attention on stochastic Volterra integral Equations (SVIEs) in It ô form The next section introduces the general class of stochastic ϑ-methods [15], for which a stability analysis is provided according to the basic linear test equation, with λ, µ ∈ R, and the convolution test equation, with λ, µ, σ ∈ R.

The Methodology: ϑ-Methods for SVIEs
According to the classical theory on the numerical approximation of Volterra integral equations [18][19][20][21][22], we refer to the equidistant grid I h = {t n = nh, n = 0, ..., N, Nh = T} and, by evaluating Equation (10) in t n , we obtain The stochastic ϑ-method, introduced in [15], assumes the following form with Y 0 = X 0 , h = t i+1 − t i .G i is a standard Gaussian random variable, i.e., it is N (0, 1)-distributed.Convergence analysis of the stochastic ϑ-method, provided in [15] as well as in [10] for ϑ = 0, relies on the same hypothesis of the existence and uniqueness of the solution to Equation (10) [12,23].In particular, the authors proved that the stochastic ϑ-method is convergent with strong order 1/2, i.e., there exists a real constant C such that However, as is usual in the numerical approximation of SDEs, such an order of convergence can be improved by adding further terms in the expansion of the right-hand side [11,24,25], leading to the following improved stochastic ϑ-method where ∂ ∂x denotes partial differentiation with respect to the second argument, while G i,1 and G i,2 are mutually independent N (0, 1) random variables.Clearly, the numerical scheme can be made free from any derivative by suitable finite difference approximation.For instance, acting as in [11] leads to the following derivative free method One can prove (see [11,15]) that the strong order of convergence of the methods defined by Equations ( 15) and ( 16) is equal to 1, since there exists a real constant K such that

Stability Issues
The stability analysis provided in [15] relies on investigating the behavior of the above methods when applied to the linear test Equation ( 11) and the convolution one (12).As regards the linear test equation, mean-square and asymptotic stability properties of the exact solution respectively occur when The following result, proved in [15], occurs.
Moreover, the above methods are asymptotically stable if and only if .
There is a strict analogy between the above presented methods and some similar formulae for SDEs.Indeed, one can see that the stability properties of the ϑ-method in Equation ( 13) completely parallels those of the Euler-Maruyama ϑ-method for SDEs [26][27][28][29][30][31].If we remove from Equation ( 15) the sum involving the derivative ∂a ∂x , we obtain i.e., this scheme is analogous to the Milstein ϑ-method for SDEs [26], with Figures 3 and 4, respectively, show the regions of mean-square stability of the methods defined by Equations ( 13), ( 15), ( 16) and ( 21) with respect to the basic test Equation (11) for ϑ = 1/2 and ϑ = 3/4.Such methods, introduced in [15], can achieve unbounded stability regions in correspondence with the some values of the parameter ϑ ≥ 1/2.As visible from Figure 5, some methods have unbounded regions when ϑ > 1. Regions of asymptotic stability are depicted in Figures 6 and 7, respectively, for ϑ = 1/2 and ϑ = 3/4.We observe that the regions of asymptotic stability are depicted by computing, for each point in the rectangle [−4, 2] × [0, 8], the value of the expectation contained in the left-hand side of the inequality (20) and checking if it is a negative number.Such expectation has been computed over 500 paths.Let us now move to the analysis with respect to the convolution test Equation (12).The following result occurs.Theorem 2. Let x = hλ, y = hµ 2 and z = h 2 σ.The stochastic ϑ-methods given by Equations (13), ( 15), ( 16) and (21) are mean-square stable with respect to the convolution test Equation (12) if the spectral radius ρ(K) is less than 1, where and (i) ζ = η = ψ = 0 for the method given by Equation ( 13), (ii) ζ = 1 2 y, η = x √ y, ψ = z √ y for for the improved method given by Equation ( 15), for the derivative free method given by Equation( 16), (iv) ζ = 1 2 y, η = 0, ψ = z √ y for the method given by Equation (21), .
Mean-square and asymptotic stability regions of the above methods with respect to the convolution test Equation ( 12) are depicted in Figures 8 and 9, respectively.We observe that the region of asymptotic stability is depicted by computing, for each point in the rectangle [−2, 0] × [0, 4], the absolute value of the solution to Equation ( 12) and checking if it is smaller than a prescribed threshold, assumed equal to 1e − 4 in our implementations.We also observe that larger selection of stability regions, in correspondence with different values of z and ϑ, has been provided in [15].

Conclusions
We have presented a short review of some recent results regarding the analysis of stability issues for stochastic differential Equation (1) and stochastic Volterra integral Equation (10).
As regards SDEs, the attention has been focused on analyzing the long-term behavior of two-step methods when applied to the linear stochastic damped linear oscillator described by Equation (2).We observe that the evidence provided in Section 2 mostly focuses on the case η > 1, i.e., that of strongly damped oscillators.It is very relevant to focus also on the case η << g, characterizing weakly damped oscillators, whose relevance is high in many physical problems.This issue will be analyzed in detail in future contributions.The tool introduced in [1][2][3] gives the possibility to provide a priori analysis relying on simple computations on 2 by 2 matrices.The structure-preserving approach to SDEs will next be devoted to stochastic Hamiltonian problems [7,8] and the stochastic extension of existing deterministic approaches [32][33][34].A stochastic version of the non-polynomial fitting for oscillatory problems will also be addressed [35][36][37][38].
As regards SVIEs, the investigation has concerned the analysis of mean-square and asymptotic stability properties of ϑ-methods, which is going to be oriented on wider families of methods, such as stochastic Runge-Kutta methods for SVIEs, in future investigations.It is also worth observing that there is a connection between stochastic integral equations and stochastic differential equations.This connection yields, for instance, to the possibility of analyzing the properties of numerical methods for SVIEs through the corresponding SDE theory.Future contributions will focus on the analysis of stability properties of numerical methods for SVIEs that inherit the stability recursion of certain methods for SDEs; this analysis is helpful to assess a stability theory of numerical methods for SVIEs, which has been partially given in Section 3.
A further relevant issue regards the employment of other test equations for the stability analysis of numerical methods for SVIEs.Indeed, it is physically relevant and interesting to consider test equations depending on exponential kernels of the type λ exp(−σ(t − s)), characteristic of system with fast response, rather than the case of the convolution test Equation (12), typical of slowly responding systems.This issues will also be addressed for stochastic fractional differential equations [35].

Figure 3 .
Figure 3. Mean-square stability regions in the (x, y)-plane with respect to the basic test Equation (11), case ϑ = 1/2.The values of x and y are given in Theorem 1.