Stability index of linear random dynamical systems

Given a homogeneous linear discrete or continuous dynamical system, its stability index is given by the dimension of the stable manifold of the zero solution. In particular, for the $n$ dimensional case, the zero solution is globally asymptotically stable if and only if this stability index is n. Fixed n, let p_k, k=0,1,...,n, denote the probabilities that the random variable that assigns to each linear random dynamical system its stability index takes the value k. In this paper we obtain either the exact values p_k, or their estimations by combining the Monte Carlo method with a least square approach that uses some affine relations among the values p_k,k=0,1,...,n. The particular case of n-order homogeneous linear random differential or difference equations is also studied in detail.


Introduction
Nowadays it is unnecessary to emphasize the importance of ordinary differential equations and discrete dynamical systems to model real world phenomena, from physics to biology, from economics to sociology. These dynamical systems, a concept that includes both continuous and discrete models (and even dynamic equations in time-scales), can have undetermined coefficients that in the case of real applications, must be adjusted to fixed values that serve to make good predictions: this is known as the identification process. Once these coefficients are fixed we obtain a deterministic model.
In recent years some authors have highlighted the utility of considering random rather than deterministic coefficients to incorporate effects due to errors in the identification process, natural variability in some of the physical parameters, or as a method to treat and to incorporate uncertainties in the model, see [4,5,16] for examples coming from biological modeling and [9] for examples coming from mechanical systems.
In the same aim that inspires some works like [1,11], in this paper we focus on giving a statistical measure of the stability for both discrete and continuous linear dynamical where both x, x k ∈ R n and A is an n × n real matrix.
More concretely, in the continuous (resp. discrete) case we define the stability index of the origin, s(A), as the number of eigenvalues, taking into account their multiplicities, of A with negative real part (resp. modulus smaller than 1). This index coincides with the dimension of the invariant stable manifold of the origin. Notice also that if s(A) = n (resp. s(A) = 0) the origin is a global stable attractor (resp. a global unstable repeller).
In this work we study the probabilities p k for a linear dynamic system (1) to have a given stability index k when the parameters of the matrix A are random variables with a given natural distribution. As we will see in Section 2, this distribution must be that all the elements of A are independent and identically distributed (i.i.d.) normal random variables with zero mean. We also will study the same question for linear n-th order differential equations and for linear difference equations.
We also remark that our results can be extrapolated to know a measure of the stability behaviour of critical or fixed points for general non-linear dynamical systems, because near them they can be written aṡ with f being a non-linear term vanishing at zero. Moreover, while the situation where the origin is non-hyperbolic is negligible, in the complementary one, the stability index of the 2 linear part coincides with the dimension of the local stable manifold at the point.
The key tool in the continuous case to know the stability index of a matrix is the Routh-Hurwitz criterion, see for instance [8, p. 1076]. This approach allows to know the number of roots of a polynomial with negative real part in terms of algebraic inequalities among its coefficients. Similarly, its counterpart for the discrete case is called Jury criterion. It is worth to observe that in fact both are equivalent and it is possible to get one from the other by using a Möbius transformation that sends the left hand part of the complex plane into the complex ball of radius 1.
In all the cases, when we do not know how to compute analytically the searched probabilities, we introduce a two steps approach to obtain estimations of them: • Step 1: We start using the celebrated Monte Carlo method. Recall that this computational algorithm relies on repeated random sampling and gives estimations of the searched probabilities based on the law of large numbers and the law of iterated logarithm, see [2,13]. It holds that using n samples this approach gives the searched value with an absolute error of order O ((log log n)/n) 1/2 , which practically behaves as O(n −1/2 ). Since in all our simulations we will work with n = 10 8 , our first approaches to the desired probabilities will have an absolute error of order 10 −4 .
• To have a flavour of the type or results that we will obtain we describe several consequences of some of our results for linear homogeneous differential or difference equations of order n with constant coefficients, see Sections 5 and 7. A first result is that in both cases the expected stability index is n/2. Moreover, let r n denote the probability of the 0 solution to be a global stable attractor (stability index equals n) for them. Then, for differential equations, r n ≤ 1/2 n . Furthermore, r 1 = 1/2, r 2 = 1/4, r 3 = 1/16 and our two steps approach gives that r 4 ≃ 0.00925, r 5 ≃ 0.00071, and that r k is smaller that 10 −4 for bigger k.

A suitable probability space
In our approach, the starting point is to determine which is the "natural" election of the probability space and the distribution law of the coefficients of the linear dynamical system.

3
Only after this step is fixed we can ask for the probabilities of some dynamical features or some phase portraits.
For completeness, we start with some previous considerations and with an example, already considered in the literature, see [1,11,18]. Consider the planar linear differential system: ẋ where A, B, C, D are random variables, so that the sample space is Ω = R 4 . It is plausible to require that these real random variables are i.i.d. and continuous. Also, according to the principle of indifference (or principle of insufficient reason) [6], it would seem reasonable to impose that them were such that the random vector (A, B, C, D) had some kind of uniform distribution in R 4 . But there is no uniform distribution for unbounded probability spaces.
However still, there is a natural election for the distribution of the variables A, B, C and D.
Indeed, it is well known that the phase portrait of the above system does not vary if we multiply the right-hand side of both equations by a positive constant (which corresponds to a constant time scale). This means that in the space of parameters, R 4 , all the systems with parameters belonging to the same half-straight line passing through the origin are topologically equivalent and in particular have the same stability index. Hence, we can ask for a probability distribution density f of the coefficients such that the random vector has a uniform density on the sphere S 3 ⊂ R 4 , that is a compact set.
The question is, which are the probability densities f that give rise to a uniform distribution of the vector (3) on the sphere? The answer is that, just assuming that f is continuous and positive, f must be the density of a normal random variable with zero mean. Moreover, this result is true for arbitrary dimension, see next theorem. We remark that the reciprocal result is well known [12,14].
. . , X n be i.i.d. one-dimensional random variables with a continuous positive density function f . The random vector , has a uniform distribution in S n−1 ⊂ R n if and only if each X i is a normal random variable with zero mean.
Curiously, in the case that we cannot assign uniform distributions, there is an extension of the indifference principle which suggest to use those distributions that maximize the 4 entropy, i.e. the quantity h(f ) = − Ω f (x) ln(f (x))dx for any given density f . The onedimensional random variables with continuous probability density function f on Ω = R that maximize the entropy are again the Gaussian ones, [6,Thm 3.2].
Of course, if instead of properties concerning general dynamical systems one focuses on particular models in which the parameters have specific restrictions due to physical or biological reasons one must consider other type of distributions, see for instance [16].
Using Theorem 1, and going back to the initial motivating example, in order to study (2) we have to consider the probability space (Ω, F, P ) where Ω = R 4 , F is the σ-algebra da db dc dd, that is 1/2 by symmetry, as we will see.
Proof of Theorem 1. Let (X 1 , . . . , X n ) be the random vector associated with the with the random variables of the statement, with joint continuous density function g(x 1 , . . . , x n ). We claim that for some continuous function h.
In consequence, writing this fact in cartesian coordinates, we get that almost everywhere , for some continuous function h and the claim (4) follows. Now we complete the proof. Since X 1 , . . . , X n are i.i.d., with positive density f , we .
It is well known that all its continuous solutions are ϕ(x) = ax, for some a ∈ R. Hence all continuous solutions of (5) are H(x) = e ax .
As a consequence, f (x) = b e ax 2 for some (a, b) ∈ R 2 . Since f is a density function, The reciprocal part is straightforward and well known [12,14].

A preliminary result and methodology
We will investigate the probabilities of having a certain stability index for several linear dynamical systems with random coefficients. In particular we consider: (a) Differential systemsẋ = A x where x ∈ R n and A is a real constant n × n matrix, (b) Homogeneous linear differential equation of order n with constant coefficients: a n x (n) + a n−1 x (n−1) (d) Linear homogeneous difference equation of order n with constant coefficients a n x k+n + a n−1 x k+n−1 + · · · + a 1 x k+1 + a 0 x k = 0.
Notice that in the four situations the behaviour of the dynamical systems does not change if we multiply all the involved constants by the same positive real number. This fact situates the four problems in the same context that the motivating example (2). Hence, following the results of Section 2, in all the cases the coefficients will be i.i.d. normal variables with zero mean and standard deviation 1.
Hence in all cases we have a well-defined probability space (Ω, F, P ), where Ω = R m , with m = n 2 , n + 1, n 2 + 1 or n + 1 according we are in case (a), (b), (c) or (d), respectively, F is the σ-algebra generated by the open sets and for each A ∈ F, where a = (a 1 , a 2 , . . . , a m ), ||a|| 2 = m j=1 a 2 j and da = da 1 da 2 . . . da m . For instance the matrices A appearing in case (a) and (c) are the so called random matrices.
The use of Routh-Hurwitz algorithm is a very useful tool to count the number of roots of a polynomial with negative real parts and it is implemented in many computer algebra systems. These conditions are given in terms of algebraic inequalities among the coefficients of the polynomials. Let us recall how to use it to count the number of roots with modulus less than one of a polynomial and, hence, to obtain the so called Jury conditions. Given any polynomial Q(λ) = q n λ n + q n−1 λ n−1 + · · · + q 1 λ + q 0 with q j ∈ C, by using the conformal transformation λ = z+1 z−1 , we get the associated polynomial It is straightforward to observe that λ 0 ∈ C is a root of of Q(λ) such that |λ 0 | < 1 if and Hence, because Routh-Hurwitz and Jury conditions are semi-algebraic, in all cases the random variable that X that assigns to each dynamical system its stability index k, 0 ≤ k ≤ n, is measurable. Hence A k := {a ∈ R m : X(a) = k} ∈ F and its probability p k := P (A k ) is well defined. Observe also that the non-hyperbolic cases are totally negligible because for their characterization some algebraic equalities appear. In this paper we will either calculate or estimate in the four situations the values p k for k ≤ 10.

A preliminary result
In three of the above considered cases we will apply the following auxiliary result: Let (Ω, F, P ) be a probability space and let Y : Ω → R be a discrete random variable with image Im(Y ) = {0, 1, . . . , n}, and probability function (b) If n is even and n ≥ 2 then 2 If, additionally, n is even * and i odd p i = i even p i = 1 2 , then (c) If n 2 is even, then 2 Proof. We start proving that E(Y ) = n/2. Assume for instance that n is odd. Since When n is even the proof is similar.
The proof of all the four items is straightforward and we omit it.

Experimental methodology
In all the cases considered in the paper, when we can not give an exact value of the probabilities p k we start estimating them by using the Monte Carlo method (further details * When n is odd the imposed equalities automatically hold. are given in each of the following subsections), [13]. From these obtained estimations, the observed relative frequencies, we improve them via the least squares method, by using the linear constraints given in Corollaries 4, 6 and 13.
A more detailed explanation of this second step is as follows: The probabilities p k satisfy some affine relations, like the ones in Lemma 2 or the ones in Proposition 12. Then, if we denote p = (p 0 , . . . , p n ) t ∈ R n+1 it is possible to write p = M q + b where q ∈ R k with k ≤ n is a vector whose components are different elements of p 0 , . . . , p n ; M ∈ M n×k (R), and b ∈ R k . Let p = ( p 0 , . . . , p n ) t be the vector with the estimated probabilities obtained by the observed relative frequencies using the Monte Carlo method. Then, we can find the least squares solution [15, Def. 6.1] of the the system, obtaining that see [7, p. 198] or [17, p. 200]. So we can find some improved estimations p, via: Some detailed examples are given in Sections 4, 5 and 7.
We write λ 1 λ 2 · · · λ n = (λ 1 λ 2 · · · λ k )(λ k+1 λ k+2 · · · λ n ) where λ 1 , λ 2 , . . . , λ k are all the real negative ones. Observe also that for complex eigenvalues λλ > 0. Hence λ k+1 λ k+2 · · · λ n > 0, From the above proposition it easily follows: N (0, 1) entries, let X be the random variable defined above and p k = P (X = k). Then the probabilities p k satisfy all the consequences of Lemma 2. In particular E(X) = n/2. Now we reproduce some experiments to estimate the probabilities p k for low dimensional cases. We apply the Monte Carlo method, that is, for each considered dimension n, we have For each considered dimension of the phase space n, and in order to take advantage of the relations stated in Corollary 4, we can refine the solutions using the least squares solutions of the inconsistent linear system associated to these relations when using the observed frequencies obtained by the Monte Carlo simulation.
We summarize the results of our experiments in the Table 1, where the observed relative frequencies and the estimates are presented only up to the fifth decimal (in the table, and in the whole paper, frequency stands for relative frequency) because as we already explained in the introduction, the expected absolute error will be of order 10   Table 1. Linear stability indices for linear random differential systems.
14 5 Linear random differential equations of order n In this section we consider linear random homogeneous differential equations of order n where x = x(t), the derivatives are taken in respect to t, and A j are again i.i.d. random variables with N (0, 1) distribution.
To get the stability index for them we only need to know the probability distributions of the number of roots with negative real part of its associated random characteristic polynomial: Let X be the random variable that counts the number of roots of Q(λ) with negative real parts and define p k = P (X = k) for k = 0, 1, . . . , n.
Proposition 5. Set p k = P (X = k), where X is the random variable defined above. Then (c) i odd p i = i even p i = 1 2 .
Proof. The proof of (a) is trivial. To prove (b) consider equation (10) with its characteristic polynomial Q(λ) and also the new differential equation with its characteristic polynomial Q * (λ) = Q(−λ) = (−1) n A n λ n + (−1) n−1 A n−1 λ n−1 + · · · − A 1 λ + A 0 . Since Q(λ) = 0 if and only if Q * (−λ) = 0 we get that p k = p * n−k where p * i the probability that Q * (λ) has i roots with negative real part. But also p k = p * k because the coefficients of the equations (10) and (11)  Similarly as in the proof of (c) of Proposition 3 we observe that the polynomial Q(λ) has an odd number of roots with negative real part if and only if A 0 · A n < 0, because we can neglect the case of having some roots with zero real part. Since the coefficients of (10) are symmetric independent random variables, the probability that Q(λ) has an odd number of roots with negative real part is Corollary 6. Consider the linear random homogeneous differential equation of order n (10), with all A i being i.i.d. N (0, 1) random variables, let X be defined above, and set p k = P (X = k). Then the probabilities p k satisfy all the consequences of Lemma 2. In particular E(X) = n/2.
For each n, let r n be the probability of the origin to be a global stable attractor (asymptotically stable equilibrium) for (10), that is r n = p n . By Proposition 5(b) this probability coincides with the one of being a repeller because p n = p 0 . Our results in Proposition 8 seem to indicate that r n decreases with n. Before proving this proposition we need a preliminary result. This is so, because for instance in the definition of A + , the last inequality can also be written as U > SV /T > 0 and from it we know that the condition U > 0 can be removed. Finally, interchanging U and V and S and T we get the same relations in the definitions of A + and A − . Since all variables are independent N (0, 1), both sets have the same probability and p + = p − , as we wanted to prove.
Proof of Proposition 8. Notice that r n is the probability that the characteristic polynomial Q(λ), associated to the random differential equation (10) is a Hurwitz stable polynomial, that is r n = P (Every root of Q(λ) belongs to R − ), where R − = {z ∈ C such that Re(z) < 0}. It is well-known that a necessary condition for a polynomial to have every root in R − is that all its coefficients have the same sign. This is so because it holds for polynomials of degree 1 and 2, and this property is preserved when we multiply two polynomials satisfying it. Hence, Since the variables A i are independent and symmetric As a consequence, and the first statement follows.
Let us prove that r 3 = p 3 = 1/16. By using the Routh-Hurwitz criterion [8, p. 1076], it can be seen that a 3 λ 3 + a 2 λ 2 + a 1 λ + a 0 has every root in R − if and only if all its coefficients have the same sign and moreover a 1 a 2 − a 0 a 3 > 0. Hence, , being all A i ∼ N (0, 1) and independent. Due to their symmetry, the random variables A i and −A i , for i = 0, . . . , 3 have the same distribution and hence p + 3 = p − 3 . Therefore p 3 = 2p + 3 . The result follows now by Lemma 7, which gives p + 3 = 1/32. Let us prove that r 4 < r 3 /2. To compare both probabilities, here it will be more convenient to write the coefficients of the polynomials with subscripts with increasing ordering, that is q n (x) = a 0 x n + a 1 x n−1 + · · · + a n−1 x + a n . With this notation, which also respects the traditional notation when writing the Hurwitz matrices, and when a 0 > 0, the Routh- Hence, these conditions when a 0 > 0 and for n = 3 are: a 1 > 0, a 1 a 2 − a 0 a 3 > 0 and a 3 > 0.
Consider now, for n = 3, 4, the random polynomials Q n (x) =Ã 0 x n +Ã 1 x n−1 + · · · + A n−1 x +Ã n , whereÃ i ∼ N (0, 1) and are independent (notice that with this notation each coefficientÃ i is the coefficient A n−i of the characteristic polynomial). For simplicity we denote with the same name the coefficients of Q 3 and Q 4 although they are different random variables. As above, r 3 = 2p + 3 , and r 4 = 2p + 4 , where p + k = P (A + k ), with Notice that if we define it is clear that P (B) = p + 3 /2 and, moreover A + 4 ⊂ B 3 , with a strict inclusion. Hence p + 4 = P (A + 4 ) < P (B) = p + 3 /2, and r 4 < r 3 /2, as we wanted to show.
Corollary 9. Consider a linear random homogeneous differential equation of order n = 3 and the random variable X defined above. Then p 0 = p 3 = 1/16 and p 1 = p 2 = 7/16.
Proof. By the above proposition, for n = 3, p 0 = p 3 = r 3 = 1/16. Hence, by Proposition 5, The computations in this case are similar to the ones of the previous section and the obtained results are summarized in Table 2. We only give some comments for the cases n = 8 and 10, where we have encountered that the vectors p have negative and very small , 13569 80000000 , 4631293 200000000 .

Linear random maps
In order to keep the approach of the preceding sections, we suggest to consider random linear discrete dynamical systems of the form where B and each of the n 2 entries of the random matrix A are i.i.d. N (0, 1) random variables. Then, given a linear discrete random system (13), its characteristic random polynomial associated to the matrix 1 B A is where each random variable Q j is a polynomial expression of 1/B, A 1,1 , . . . , A n,n with not an easy distribution law. We denote by X the random variables that assigns to each Q its number of roots with modulus smaller than 1, that is, the stability index of the matrix 1 B A. Also p k denotes the probabilities that X takes the value k.
As we will see in the examples, in this case the condition p k = p n−k is no more satisfied.
Among other reasons it happens that the entries of A −1 have not simple distributions. Since we do not know other relations on the probabilities p k than the trivial one n k=0 p k = 1, which is directly fulfilled by the observed relative frequencies, in this case we do not perform the least squares refinement.
The case n = 1 is the only one that we have been able to solve analytically. Notice that in this situation the only eigenvalue of Q(λ) is λ = A/B, with A and B independent and N (0, 1). Hence p 0 = P (|A/B| > 1) and p 1 = P (|A/B| < 1) = P (|B/A| > 1). Since A/B and B/A have the same distribution it holds that p 0 = p 1 = 1/2.
As in the other models, for each dimension n ≤ 10, we generate 10 8 discrete systems of the form (13). For each matrix 1 B A we have computed the characteristic polynomial Q and its associated polynomial Q ⋆ (see Equation (6)) and have counted the number of roots of this last polynomial by using the Routh-Hurwitz zero counter. For n ≥ 5 and in order to decrease the computation time we have directly numerically computed the eigenvalues of the matrix and counted the number of them with modulus less than one. The results obtained are shown in Table 3. 7 Linear random difference equations of order n Finally we consider difference equations of order n of type where all the coefficients are i.i.d. random variables with N (0, 1) distribution. In this situation, the stability index is given by the number of zeros with modulus smaller than 1 of the random characteristic polynomial Q(λ) = A n λ n + A n−1 λ n−1 + · · · + A 1 λ + A 0 . As in the preceding sections let X be the random variable that counts the number of roots of Q(λ) with modulus smaller than 1 and set p k = P (X = k) for k = 0, 1, . . . , n.
Before proving some relations among the probabilities p k , we give two preliminary lemmas. Let erf(x) = 2 √ π x 0 e −u 2 du be the error function. The following result is stated in [3]. We prove it for the sake of completeness.
Next result is a consequence of the previous lemma.
Proof. The joint density function of the random vector (U, V ) is f σ (u)f ρ (v), where f s (u) = e −u 2 /(2s 2 ) /( √ 2πs). Observe that the points (u, v) ∈ R 2 such that u 2 − v 2 > 0 is the region where −|u| < v < |u|, hence by symmetry, where in the last equality we have used Lemma 10.
(d) When n = 2k is even, i even p i = 2 π arctan k + 1 k and i odd p i = 2 π arctan k k + 1 .
The sum of all p i when i is odd can be obtained from the previous one, see also (14).
Corollary 13. (i) Consider the linear random homogeneous difference equation of order n, let X be the random variable defined above and p k = P (X = k). Then the probabilities p k satisfy all the consequences of Lemma 2. In particular E(X) = n/2.