Stochastic differential equations (SDE) of the second and higher order with or without time-varying delays, naturally appear in multiple applications, where deterministic models are perturbed by the white noise or its generalizations. A classical example is the Langevin equation (see e.g.  ). Liénard-type second-order stochastic equations were examined in multiple papers (see, for example,    and  ). In  explicit conditions for asymptotic stability of the second-order equation under additive white noise were obtained. In  boundedness and exponential stability conditions for second-order SDEs with a constant delay were examined. Other examples deal with the study of sensibility on stochastic perturbations of regenerative vibrations in milling process described by the second order linear differential equation with delays  ; a planar inverted pendulum on a cart, operating under modelling uncertainties and stochastic perturbations, modeled by the two-dimensional subsystem in  .
Stochastic high-order models of processes related to abrasive waterjet milling or fluid energy milling (batch grinding) are well-known as well (see e.g.  ). Stochastic high-order feedforward nonlinear system with time-varying delay was suggested to model many physical devices, such as the ball-beam with a friction term and the cart-pendulum system  . Large-scale stochastic high-order systems with time-varying delay are described by a series of interconnected subsystems in  . A hyperjerk system is a dynamical system governed by an n-th order ordinary differential equation with describing the time evolution of a single scalar variable (see e.g.  ), which can also be perturbed by a stochastic noise.
Several definitions of stochastic Lyapunov stability are used in the literature, e.g. stability in probability, stability in the mean and almost sure stability, stability of the p-th mean (p-stability), and even more. For applications to real systems, stability properties that are close to deterministic stability (almost sure sample stability) are the most desired, while conditions for p-stability are technically easier to obtain.
In this paper we study the global p-stability of the linear n-th order Itô delay equation
and its non-autonomous generalization
where are independent scalar Brownian motions defined on a probability space.
Stability of stochastic first-order differential equations with delays, as well as systems of equations, has been extensively studied (see       ,    and the references therein). The main tool for studying the global stability is the Lyapunov functional method and its stochastic modifications (see e.g.  and the references therein). While this method has been and remains the leading technique, numerous difficulties with the theory and applications to specific problems, even in the deterministic case, persist. It is, in particular, the case if one examines stability of high-order differential equations. Of course, one can always apply the Lyapunov method after reducing a high-order equation to a first-order system in the well-known way, and this technique does work in specific situations. Yet, this approach may also fail in many cases. That is why, very few papers in recent literature have examined stability of high-order stochastic differential equations with time-varying delays as such.
In the recent paper  , a new, more flexible algorithm of reducing a high-order deterministic differential equation with delay to a first-order system was suggested. The approach goes back to the theory of M-matrices. This idea is not new: absolute stability via M-matrices was studied in the monograph  , while in the more recent paper  this technique was applied to stability of neural networks. The efficiency of this method in connection with high-order deterministic equations was demonstrated in  . In this paper we claim that this approach is applicable to the SDE with delays as well, provided that the analysis based on M-matrices is combined with the regularization method from  . The latter method differs from the classical Lyapunov technique, which presupposes the existence of suitable Lyapunov functionals. Rather, the method from  requires the existence of a suitable auxiliary equation which is used to regularize the original equation and subsequently to check solvability of a regularized equation in a carefully chosen space of stochastic processes.
In conclusion, we stress that even if this paper studies stochastic linear equations, the various linearization criteria for a nonlinear stochastic differential equation (see e.g.  and the references therein) in combination with the tests obtained in this paper, can be used to examine the local stability of a nonlinear stochastic differential equations.
Let be a stochastic basis, where Ω is set of elementary probability events, F is a σ-algebra of all events on Ω, is a right continuous family of s-subalgebras of F, P is a probability measure on F; all the above s-algebras are assumed to be complete w.r.t. P, i.e. containing all subsets of zero measure; the symbol E stands below for the expectation related to the probability measure P. The expectation (the integral w.r.t. the measure P) is denoted by E.
We will use the following notations:
- is an arbitrary yet fixed norm in Rn, being the associated matrix norm.
- is the Lebesgue measure on .
- is the norm in a normed space X.
- p is an arbitrary real number satisfying .
- is the standard m-dimensional Brownian motion (i.e. the scalar Brownian motions are all independent).
Recall that the classic Marcinkiewicz-Zygmund inequality
where are independent random variables with the zero mean, can be extended to the integral form
for any predictable stochastic process ( ), any and any component ( ) of the Brownian motion B. The inequality (4) is often used in this paper. In 1988 D.L. Burkholder proved (see for example  ,  ) that in the Marcinkiewicz-Zygmund inequality (3) the constant is best possible for all for .
Equation (2) is assumed to be equipped with the initial conditions
1) , are Lebesgue measurable functions defined on ; in addition, we assume that μ―almost everywhere for some positive constants , μ―almost everywhere for some positive constants , and μ―almost everywhere for some positive constants .
2) are Lebesgue measurable functions defined on and satisfying the estimates μ―almost everywhere for some positive constants , , .
3) φ is an ―measurable, scalar stochastic process defined on , where .
4) is an ―measurable random variable for .
We define a solution of the initial value problem (2), (5), (6) to be a predictable stochastic process , which is ―times differentiable on and which satisfies the initial conditions (5), (6) and the integral equation
where the integrals are understood in the Lebesgue and the Itô sense, respectively, and if .
The initial value problem (2), (5), (6) has a unique (up to the natural P-equivalency) solution (see e.g.  ). In other words, the stochastic process satisfies Equation (2) and the initial conditions ( ), ( ).
We will write , where denotes the linear space of all n―dimensional, ―measurable random values. In addition, we define the following normed space:
Definition 1  We say that Equation (2) is exponentially p-stable ( ) if there are positive numbers such that all solutions of the initial value problem (2), (5), (6) satisfy
The analysis of the exponential p-stability of Equation (2) will be performed via an equivalent first order system of Itô equations. The technique of reduction of a high-order linear differential equation to a system by the substitution is quite common, and for system (2) it yields
where the first component of the solution of initial value problem (10), (11), (12) coincides with the solution x of the initial value problem (2), (5), (6), so that the exponential p-stability of Equation (2) follows from the exponential p-stability of system (10); and the latter can be, at least in the theory, studied by the Lyapunov-Razumikhin method of the stability analysis of stochastic delay equations. This method is based on finding a suitable Lyapunov function satisfying special conditions (see e.g.  ), which guarantee the stability properties in question. However, practical implementation of this technique seems to be difficult.
Below we use the generalized reduction technique based on a set of positive parameters, which can be chosen arbitrarily. Adapting this set to the coefficients of the given stochastic equation considerably increases, and this will be shown in the paper, flexibility of the reduction method. In addition, we will combine this technique with the method of stability analysis based on positive invertible matrices  being, to our opinion, a more efficient alternative to the Lyapunov-Razumikhin algorithm, at least in the case of stochastic linear equations with delay.
Let ( ) be some positive numbers. Consider the following generalization of system (10):
where ( ), ( , ) and
for , and the other entries are obtained from Equation (2). System (13) is supposed to be equipped with the initial conditions (11), (12).
Let us make some comments on this reduction technique. According to the paper  , the solution of the deterministic counterpart of the initial value problems (13), (11), (12) (i.e. in the absence of all Brownian motions) gives the solution of the (deterministic) problem (2), (5), (6) if one puts . Replacing the chain rule by Itô’s formula leads to the same conclusion for the stochastic initial value problem (2), (5), (6). In particular, the exponential p-stability of Equation (2) follows from the exponential stability of system (13) for any .
Lemma 1 Let be a scalar function which is square integrable on , be a predictable stochastic process satisfying . Then
Proof. Once we prove the inequality (15), the inequality (16) can be justified similarly.
3. Main Result
An -matrix is called nonnegative if , , and positive if , .
Definition 2 A matrix is called an M-matrix if for , and one of the following conditions is satisfied:
- G has a positive inverse matrix G−1;
- the principal minors of the matrix G are positive.
Now we define an -matrix G which plays a crucial role in the theorem below. Let
Here and for all , . These numbers can be expressed via the constants , and from assumption 1 in Section 2. Thus, the matrix G becomes
Theorem 1 Assume that and there exist positive numbers such that and
Then system (13), and hence Equation (2), is exponentially 2p-stable.
Proof. First of all, we observe that the determinant of the matrix G is equal to the left-hand side of the equality (19), while the other principal minors are all equal to 1. Hence G is an M-matrix.
Now, system (13) with the conditions (11) can be rewritten as follows:
where is an unknown scalar predictable stochastic process on such that for , and a known scalar predictable stochastic process on such that for and outside the interval .
Let be the solution of (20) satisfying the initial conditions (12). A straightforward calculation shows that coincides with the solution of the initial value problem (11), (12), (13) for (but not necessarily for , of course).
We choose a positive number for all and make the following substitution into system (20): , where is an unknown predictable stochastic process defined on . By the definition, for and , thus
Let . Then by using (12), we rewrite system (21)
Denote , and
From the first equations in (22) we obtain
The estimate (4) and the last equation in (22) yield
the inequality (24) yields
Denote and define the --matrix by putting
Then from (27) we obtain the componentwise vector inequality
where is an n-dimentional column vector.
Since is an M-matrix, is also an M-matrix for small l. Therefore there exists a number such that the matrix is positive invertible. The inequality (29) justifies
Based on the inequality (30), we conclude that the solution of the initial value problem (13), (11), (12) satisfies
where , . Therefore system (13) is exponentially 2p-stable. Theorem 1 is proven.
4. Some Corollaries
In this section we consider a second order equation (a particular case of Equation (10) if n = 2)
that we transform into
where q1 is some positive number and , for .
The matrix G is now defined as
and , , for all .
Corollary 1 Assume that and there exists a positive number such that . Then system (33) is exponentially 2p-stable.
Proof. The statement follows from Theorem 1 and the observation that the determinant of the matrix G is equal to .
To demonstrate the efficiency of Corollary 1, let us consider the following particular case of Equation (33):
where B is the standard scalar Brownian motion, are positive real numbers. Then a straightforward application of Corollory 1 yields.
Corollary 2 Assume that and there exists a positive number such that
Then Equation (36) is exponentially 2p-stable.
5. Conclusions and Outlook
We studied Lyapunov stability of high-order linear stochastic Itô equations with delay using a non-Lyapunov approach, which combines the method described in the review paper  with the reduction technique based on the so-called M-matrices (see e.g.  ). This gave us an opportunity to find several efficient stability conditions in terms of the coefficients of the equations.
Solution of the following problems will complement the results of the present paper:
-Find explicit stability conditions for the 3rd and higher order stochastic delay equations in terms of their coefficients.
-Find sufficient conditions for the stability of the linear hybrid SDE with delay
where is a Markov chain with its state space S, which is independent of the Brownian motions and which represents random switchings between different delay equations.
-Generalize the suggested framework to the case of high-order SDEs driven by an arbitrary semimartingale, rather than by the Brownian motion.
-For the van der Pol-type SDE delayed equation under perturbations of white noise
examine stability by using the linearization criteria introduced in   in combination with the tests obtained in this paper.
The work of the first and the third author was partially supported by the grant DDG-2015-00046 of NSERC. The work of the second and the third author was partially supported by the grant \#239070 of the Norwegian Research Council.