This is a preprint of: “Stability index of linear random dynamical systems”, Anna Cima, Armengol Gasull, Vı́ctor Mañosa, Electron. J. Qual. Theo., vol. 2021(15), 1–27, 2021. DOI: [10.14232/ejqtde.2021.1.15] Stability index of linear random dynamical systems∗ Anna Cima(1) , Armengol Gasull(1,2) and Vı́ctor Mañosa(3) (1) Departament de Matemàtiques, Facultat de Ciències, Universitat Autònoma de Barcelona, 08193 Cerdanyola del Vallès, Barcelona, Spain cima@mat.uab.cat, gasull@mat.uab.cat (2) Centre de Recerca Matemàtica, Edifici Cc, Campus de Bellaterra, 08193 Cerdanyola del Vallès, Barcelona, Spain (3) Departament de Matemàtiques, Institut de Matemàtiques de la UPC-BarcelonaTech (IMTech), Universitat Politècnica de Catalunya Colom 11, 08222 Terrassa, Spain victor.manosa@upc.edu February 19, 2021 Abstract 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 X be the random variable that assigns to each linear random dynamical system its stability index, and let pk with k = 0, 1, . . . , n, denote the probabilities that P (X = k). In this paper we obtain either the exact values pk , or their estimations by combining the Monte Carlo method with a least square approach that uses some affine relations among the values pk , k = 0, 1, . . . , n. The particular case of n-order homogeneous linear random differential or difference equations is also studied in detail. Mathematics Subject Classification 2010: 37H10, 34F05, 39A25, 37C75. Keywords: Stability index; random differential equations; random difference equations; random dynamical systems. ∗ The authors are supported by Ministry of Science and Innovation–State Research Agency of the Spanish Government through grants PID2019-104658GB-I00 (MICINN/AEI, first and second authors) and DPI2016- 77407-P (MICINN/AEI/FEDER, UE, third author). The first and second authors are also supported by the grant 2017-SGR-1617 from AGAUR, Generalitat de Catalunya. The third author acknowledges the group’s research recognition 2017-SGR-388 from AGAUR, Generalitat de Catalunya. 1 1 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 con- tinuous and discrete models (and even dynamic equations in time-scales), can have unde- termined 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 pro- cess, natural variability in some of the physical parameters, or as a method to treat and to incorporate uncertainties in the model, see [5, 6, 21] for examples coming from biological modeling and [11] for examples coming from mechanical systems. In the same aim that inspires some works like [1, 7, 14], in this paper we focus on giving a statistical measure of the stability for both discrete and continuous linear dynamical systems, ẋ = A x or xk+1 = A xk , (1) where both x, xk ∈ Rn 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 pk 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 as ẋ = A x + f (x), or xk+1 = A xk + f (xk ), 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 situation, the stability index of 2 the linear part coincides with the dimension of the local stable manifold at the point. In the continuous case, the key tool to know the stability index of a matrix is the Routh- Hurwitz criterion, see for instance [10, 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 the Jury criterion. It is worth observing 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 true probabilities, we introduce a two step approach to obtain estimations of them: • Step 1: We start using the celebrated Monte Carlo method. Recall that this compu- tational algorithm relies on repeated random sampling and gives estimations of the true probabilities based on the law of large numbers and the law of iterated logarithm, see [2, 3, 13, 18]. It is the case that using M samples this approach gives the true value  with an absolute error of order O ((log log M )/M )1/2 , which practically behaves as O(M −1/2 ), where O stands for the usual Landau notation. In all our simulations we will work with M = 108 , so our first approaches to the desired probabilities will have an asymptotic absolute error of order 10−4 . More detailed explanations of the sharpness of our estimations for this value of M are given in Section 3.2 by using the Chebyshev inequality and the Central limit theorem. We have used the default in-built pseudo-random number generator in the Statistics package of Maple in our simulations∗ . This procedure use the Mersenne Twister method with period 219937 − 1 to generate uniformly-distributed pseudo-random num- bers, and then the Ziggurat method, which is a kind of rejection sampling algorithm, to obtain the normally-distributed pseudo-random numbers, see [16] and [17]. Observe that our sample size M = 108 is much smaller than the period of the pseudo-random number generator, which is greater that 106001 . • Step 2: Since the results of the plain Monte Carlo simulations do not satisfy certain linear constrains concerning the true probabilities, we propose to correct them by using a least squares approach. We take as final estimates of the true probabilities the least squares solution ([20, Def. 6.1]) of the inconsistent overdetermined system obtained when relative frequencies of the Monte Carlo simulation are forced to satisfy these linear constrains. See Section 3.2 for more details. We would like to remark that there are other options to improve plain Monte Carlo simulations like variance reduction ∗ Concretely, we use the commands RandomVariable(Normal(0,1)) and Sample. 3 and quasi-Monte Carlo methods [2, 13]. To have a flavour of the type of results that we will obtain we describe several conse- quences 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 rn denote the probability of the 0 solution to be a global stable attractor (stability index equals n) for them. Then, for differential equations, rn ≤ 1/2n . Furthermore, r1 = 1/2, r2 = 1/4, r3 = 1/16 and our two step ap- proach gives that r4 ' 0.00925, r5 ' 0.00071, and that rk is smaller that 10−4 for bigger k. √ In the case of difference equations we prove that r1 = 1/2 and r2 = π1 arctan( 2) ' 0.304. 2 A suitable probability space In our approach, the starting point is to determine which is the natural choice of the probability space and the distribution law of the coefficients of the linear dynamical system. 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, 14, 23]. Consider the planar linear differential system: ! ! ! ẋ A B x = (2) ẏ C D y where A, B, C, D are random variables, so we can set the sample space to be Ω = R4 . It is plausible to require that these real random variables are independent and identically distributed (i.i.d.) and continuous. Also, according to the principle of indifference (or principle of insufficient reason) [8], it would seem reasonable to impose that these variables were such that the random vector (A, B, C, D) had some kind of uniform distribution in R4 . But there is no uniform distribution for unbounded probability spaces. Nevertheless, there is a natural choice 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 proportional change in the time scale). This means that in the space of parameters, R4 , 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   A B C D p , , , , with S = A2 + B 2 + C 2 + D 2 , (3) S S S S 4 has a uniform distribution on the sphere S3 ⊂ R4 . This achieves our objective, since S3 is a compact set. The question is: which are the probability densities f that give rise to a uniform dis- tribution 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 the next theorem. We remark that the converse result is well–known [15, 19]. Theorem 1. Let X1 , X2 , . . . , Xn be i.i.d. one-dimensional random variables with a contin- uous positive density function f . The random vector   X n 1/2 X1 X2 Xn , ,..., , with S= Xi2 , S S S i=1 has a uniform distribution in Sn−1 ⊂ Rn if and only if each Xi 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 suggests to use those distributions that maximize the R entropy, i.e. the quantity h(f ) = − Ω f (x) ln(f (x))dx for any given density f . The one- dimensional random variables with continuous probability density function f on Ω = R that maximize the entropy are again the Gaussian ones, [8, 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 [21]. 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 Ω = R4 , F is the σ-algebra generated by the open sets of R4 and P : F → [0, 1] is the probability function with 1 −(a2 +b2 +c2 +d2 )/2 density 4π 2 e , where for simplicity we take variance 1 in each marginal density function. For instance, assume that we want to compute the probability α of system (2) to have exactly one eigenvalue with negative real part. Next, we observe that the probability of having one null eigenvalue is zero. This is because the event which characterizes this pos- sibility is a subset of an event which is itself described by an algebraic equality between the random variables A, B, C, D. This subset has Lebesgue measure zero and therefore, by virtue of the fact that the joint distribution is continuous, the probability of this event, and therefore the event characterizing the null eigenvalue, must also be zero. Thus we have that α coincides with the probability of having a saddle (stability index 1) at the origin, i.e. 5 AD − BC < 0. Then, the open set U := {(a, b, c, d) ∈ R4 : ad − bc < 0} belongs to F and Z 1 2 2 − a +b +c 2 +d2 α = P (AD − BC < 0) = e 2 da db dc dd, 4π 2 U which is 1/2 by symmetry, as we will see. Proof of Theorem 1. Let (X1 , . . . , Xn ) be the random vector associated with the random variables of the statement, with joint continuous density function g(x1 , . . . , xn ). We claim that g(x1 , . . . , xn ) = h(x21 + · · · + x2n ), (4) for some continuous function h. Taking spherical coordinates, we consider the new random vector (R, Θ) ∈ Rn where R = (X12 +· · ·+Xn2 )1/2 and Θ = (Θ1 , . . . , Θn−1 ). We have X1 = R cos Θ1 , X2 = R sin Θ1 cos Θ2 , . . . Xn−1 = R sin Θ1 sin Θ2 · · · sin Θn−2 cos Θn−1 and Xn = R sin Θ1 sin Θ2 · · · sin Θn−2 sin Θn−1 . By the change of variables theorem, the joint density function of (R, Θ) is gR,Θ (r, θ) = g(r cos(θ1 ), . . . , r sin(θ1 ) · · · sin(θn−1 )) rn−1 sinn−2 (θ1 ) sinn−1 (θ2 ) · · · sin(θn−2 ) · χ where θ = (θ1 , . . . , θn−1 ), and n−2 Y χ := χ[0,∞) (r) · χ[0,2π) (θn−1 ) · χ[0,π) (θi ), i=1 where χA stands for the characteristic function of the set A. The density function of (R, Θ) conditioned to R, gΘ|R , is gR,Θ (r, θ) gΘ|R (r, θ) := , gR (r) where gR (r) is the marginal density of R: Z π Z π Z 2π gR (r) := ··· g(r cos(θ1 ), . . . , r sin(θ1 ) · · · sin(θn−1 )) dS, 0 0 0 where dS = rn−1 sinn−2 (θ1 ) sinn−1 (θ2 ) · · · sin(θn−2 )dθn−1 · · · dθ1 is the n-dimensional surface element in spherical coordinates. To prove the statement, we need to characterize which are the joint density functions g(x1 , . . . , xn ) such that when we fix R = r, the probability on the (n − 1)-dimensional sphere of radius r, denoted by Sn−1 (r), is uniformly distributed. In such a case the partial spherical segment Σr = {R = r, θi ∈ [αi , βi ] for i = 1, . . . , n − 1} must have probability P (Σr ) = S(Σr )/S(Sn−1 (r)) where S denotes the surface area. Set α = (α1 , . . . , αn−1 ) and β = (β1 , . . . , βn−1 ). Notice that Z β1 Z βn−1 S(Σr ) = ··· dS = rn−1 A(α, β) α1 αn−1 6 where Z β1 Z βn−1 A(α, β) = ··· sinn−2 (θ1 ) sinn−1 (θ2 ) · · · sin(θn−2 )dθn−1 · · · dθ1 α1 αn−1 n 2π 2 and S(Sn−1 (r)) = Γ( n rn−1 . Hence, on the one hand, 2)  n  A(α, β) P (Σr ) = Γ n , 2 2π 2 which does not depend on r. On the other hand, Z β1 Z βn−1 P (Σr ) = ··· gΘ|R dθ α1 αn−1 where dθ = dθn−1 · · · dθ2 dθ1 . This implies that Z Z β1 βn−1 g(r cos(θ1 ), . . . , r sin(θ1 ) · · · sin(θn−1 )) rn−1 sinn−2 (θ1 ) sinn−1 (θ2 ) · · · sin(θn−2 ) · χ ··· dθ α1 αn−1 gR (r)  Z Z Γ n2 β1 βn−1 = n ··· sinn−2 (θ1 ) sinn−1 (θ2 ) · · · sin(θn−2 )dθ, 2π 2 α1 αn−1 for all αi , βi ∈ [0, π) for i = 1, . . . , n − 2 with αi < βi and αn−2 , βn−2 ∈ [0, 2π) with αn−2 < βn−2 . This last equality implies that almost everywhere  Γ n2 g(r cos(θ1 ), . . . , r sin(θ1 ) · · · sin(θn−1 )) rn−1 n = , 2π 2 gR (r) and therefore g(r cos(θ1 ), . . . , r sin(θ1 ) · · · sin(θn−1 )) is a function that only depends on r. In consequence, writing this fact in Cartesian coordinates, we get that almost everywhere g(x1 , . . . , xn ) = h(x21 + · · · + x2n ), for some continuous function h and the claim (4) follows. Now we complete the proof. Since X1 , . . . , Xn are i.i.d. with positive density f , we know that g(x1 , · · · , xn ) = f (x1 ) . . . f (xn ). So equation (4) can be expressed as f (x1 ) · · · f (xn ) = h(x21 + · · · + x2n ) for all (x1 , . . . , xn ) ∈ Rn where h is a positive function. Taking x2 = · · · = xn = 0 we have that f (x1 ) f (0)n−1 = h(x21 ) and h(0) = (f (0))n > 0. Thus, h(x21 ) h(x2n ) h(x21 ) h(x2n ) f (x1 ) . . . f (xn ) = · · · = · · · (f (0))n = h(x21 + · · · + x2n ). (f (0))n−1 (f (0))n−1 (f (0))n (f (0))n Hence, using that h(0) = (f (0))n > 0, h(x21 ) h(x2n ) h(x21 + · · · + x2n ) ··· = . h(0) h(0) h(0) 7 Taking H(ξ) := h(ξ)/h(0), and ui = x2i , it holds that H(u1 ) · · · H(un ) = H(u1 + · · · + un ) with H(0) = 1. (5) Hence, ϕ(u) = log(H(u)) is a continuous function that satisfies Cauchy’s functional equation ϕ(u1 ) + · · · + ϕ(un ) = ϕ(u1 + · · · + un ) with ϕ(0) = 0. 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) = eax . 2 As a consequence, f (x) = b eax for some (a, b) ∈ R2 . Since f is a density function, R∞ 2 p a < 0. Moreover, using −∞ beax dx = b −π/a = 1, and setting a = −1/(2σ 2 ), we get that 1 x2 f (x) = √ e− 2σ2 , 2πσ 2 so each variable Xi is a normal random variable N(0, σ 2 ). The converse part is straightforward and well–known [15, 19]. Remark 1. The continuity condition for f in Theorem 1 is relevant since Equation (5) also admits non-continuous solutions that can be constructed, for instance, from non-continuous solutions of Cauchy’s functional equation known for n = 2, see [12]. 3 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 ∈ Rn and A is a real constant n × n matrix, (b) Homogeneous linear differential equation of order n with constant coefficients: an x(n) + an−1 x(n−1) + · · · + a1 x0 + a0 x = 0, (c) Linear discrete systems b xk+1 = A xk where xk ∈ Rn , b ∈ R; and A is a real constant n × n matrix, (d) Linear homogeneous difference equation of order n with constant coefficients an xk+n + an−1 xk+n−1 + · · · + a1 xk+1 + a0 xk = 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 8 the results of Section 2, in all the cases, we may take the coefficients to be i.i.d. random normal variables with zero mean and variance 1. Hence in all cases we have a well–defined probability space (Ω, F, P ), where Ω = Rm , with m = n2 , n + 1, n2 + 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, Z 1 2 P (A) = √ e−||a|| /2 da, (6) ( 2π) m A P where a = (a1 , a2 , . . . , am ), ||a||2 = m 2 j=1 aj and da = da1 da2 . . . dam . 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(λ) = qn λn + qn−1 λn−1 + · · · + q1 λ + q0 with qj ∈ C, by using z+1 the conformal transformation λ = z−1 , we get the associated polynomial Q? (z) = qn (z + 1)n + qn−1 (z + 1)n−1 (z − 1) + . . . + q0 (z − 1)n . (7) It is straightforward to observe that λ0 ∈ C is a root of of Q(λ) such that |λ0 | < 1 if and only if z0 = (λ0 + 1)/(λ0 − 1) is a root of Q? (z) such that Re(z0 ) < 0. 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 Ak := {a ∈ Rm : X(a) = k} ∈ F and its probability pk := P (Ak ) is well–defined. Observe also that the non-hyperbolic cases are totally negligible because in their characterization some algebraic equalities appear. In this paper we will either calculate or estimate in the four situations the values pk for k ≤ 10. 3.1 A preliminary result In three of the above considered cases we will apply the following auxiliary result: Lemma 2. 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 mass function pk = P (Y = k) P such that pk = pn−k for all k = 0, . . . , n. Then E(Y ) = nk=0 kpk = n/2. Moreover P n−1 (a) If n is odd then 2 2 k=0 pk = 1. In particular, when n = 1, p0 = p1 = 12 . P n2 −1 (b) If n is even and n ≥ 2 then 2 k=0 pk + p n2 = 1. 9 P P If, additionally, n is even† and i odd pi = i even pi = 12 , then P n2 −2 P n2 −1 (c) If n2 is even, then 2 k=0, k even p k + p n = 2 1 2 , and 2 1 k=1, k odd pk = 2 . In particular, when n = 4, p1 = p3 = 14 , p2 = 1 2 − 2p0 and p4 = p0 . n P n2 −1 1 P n2 −2 1 (d) If 2 is odd, then 2 k=0, k even pk = 2 , and k=1, k odd pk + p n2 = 2. In particular, when n = 2, p0 = p2 = 14 and p1 = 12 . Proof. We start proving that E(Y ) = n/2. Assume for instance that n is odd. Since pk = pn−k , its holds that kpk + (n − k)pn−k = npk , for each k ≤ (n − 1)/2. Hence, n  E(Y ) = np0 + np1 + · · · + np n−1 = 2p0 + 2p1 + · · · + 2p n−1 2 2 2 n  n = (p0 + pn ) + (p1 + pn−1 ) + · · · + (p n−1 + p n+1 ) = . 2 2 2 2 When n is even the proof is similar. The proof of all the four items is straightforward and we omit it. 3.2 Experimental methodology In all the cases considered in the paper, when we can not give an exact value of the probabil- ities pk we start estimating them by using the Monte Carlo method, see [18]. The estimates obtained (namely, the observed relative frequencies) are then improved via the least squares method, by using the linear constraints given in Corollaries 4, 6 and 13. In all the cases we will use Monte Carlo method with M = 108 to obtain an estimation, say pe, for a probability p := P (A) like the one given in equality (6) for different measur- able sets A. Further details for each concrete situation are given in each of the following subsections. In brief, recall that pe is given by the proportion of samples that are in A. For studying, for a given M, how close are p and pe, let Bj , j = 1, . . . , M be i.i.d. Bernoulli random variables, where each one of them that takes the value 1 with probability p and the value 0 with probability 1 − p. 1 PM Define PM = M e is j=1 Bj . Then, the value obtained for the random variable PM , p the approximation of p given by the Monte Carlo method. Let us see, by using Cheby- shev inequality or the Central limit theorem, that with very high probability, pe is a good approximation of p. Notice first that E(PM ) = p and due to the independence of the Bj ,   M X 1 1 p(1 − p) 1 Var(PM ) = Var  Bj  = 2 M Var(B1 ) = ≤ , M M M 4M j=1 † When n is odd the imposed equalities automatically hold. 10 because p(1 − p) ≤ 1/4. Recall also that for each ε > 0 and any random variable X, with E(X 2 ) < ∞, the Chebyshev inequality reads as Var(X) P (|X − E(X)| < ε) ≥ 1 − . ε2 Hence, applying the Chebyshev inequality to X = PM we get that p(1 − p) 1 P (|PM − p| < ε) ≥ 1 − ≥1− . M ε2 4M ε2 Taking M = 108 , as in our computations, denoting pe = P108 , and considering ε = 10−3 we get that the above probability gives the following conservative estimate of the reliability of the method  1 399 P |ep − p| < 10−3 ≥ 1 − = = 0.9975. 400 400 Let us see, by using the Central limit theorem, that the above probability seems to be much bigger. By this theorem we know that for M big enough, and p(1 − p)M also big enough, the distribution of the random variable PM − E(PM ) PM − p p =q Var(PM ) p(1−p) M can be practically considered to be a random variable Z with distribution N(0, 1). In fact, in Statistics it is usually imposed that p(1 − p)M > 18. Hence, unless p is very close to 0 or 1, the value M = 108 is big enough. Hence √ √ ! √ ! M |PM − p| ε M ε M P (|PM − p| < ε) =P p

2Φ 2ε M − 1, p(1 − p) where Φ is the distribution function of a N(0, 1) random variable. Taking again M = 108 and ε = 10−3 or ε = 2 × 10−4 we get  p − p| < 10−3 & 2Φ (20) − 1 > 1 − 10−88 , P |e  p − p| < 2 × 10−4 & 2Φ (4) − 1 > 0.99993. P |e In fact, for instance looking at the values pk of Table 2 for n = 2 in Section 4, that can pk − pk | ≤ 6 × 10−5 , for k = 0, 1, 2. So, the actual also be obtained analytically, we get that |e bound is smaller that the bounds obtained above. Finally, to illustrate how the error decays when the sample size increases, we show the evolution of the errors in one case where the true probabilities are known. We consider the 11 second order difference equation A2 xk+2 + A1 xk+1 + A0 xk = 0 where Ai are i.i.d. random variables with N(0, 1) distribution. The stability index is given by the number of zeroes with modulus smaller than 1 of the characteristic polynomial Q(λ) = A2 λ2 + A1 λ + A0 . Let X be the random variable that counts the number of roots with modulus smaller than 1 of Q(λ), and pk = P (X = k) for k = 0, 1, 2. The true value of the probabilities pk is obtained in Corollary 13. Performing Monte Carlo simulations with M = 10m with m = 2, . . . , 10 we obtain the observed frequencies pe2 (m) shown in Table 1. These frequencies are the estimated probabilities for the origin to be asymptotically stable. Notice that √ in Proposition 12 and in Corollary 13 we prove that p0 = p2 = arctan( 2)/π and, of √ course, p1 = 1 − p0 − p2 = 2 arctan(1/ 2)/π). For M = 10m we denote the absolute error em = |e p2 (m) − p2 |: M = 102 M = 103 M = 104 pe2 (2) = 0.37 pe2 (3) = 0.319 pe2 (4) = 0.3102 e2 ≈ 0.065913276015 e3 ≈ 0.014913276015 e4 ≈ 0.006113276015 M = 105 M = 106 M = 107 pe2 (5) = 0.30416 pe2 (6) = 0.303892 pe2 (7) = 0.3041241 e5 ≈ 0.000073276015 e6 ≈ 0.000194723985 e7 ≈ 0.000037376015 M= 108 M= 109 M = 1010 pe2 (8) = 0.30406079 pe2 (9) = 0.304076699 pe2 (10) = 0.304079098 e8 ≈ 0.000025933985 e9 ≈ 0.000010024985 e10 ≈ 0.000007625985 Table 1. Observed frequency and absolute error of p2 for second order difference √ equations, using that p2 = arctan( 2)/π ≈ 0.304086723985. With the above results, the regression line of Y = log(em ) versus X = log(M ) = m log(10) is Y = −0.505 X − 1.260 with R2 = 0.893. The slope is therefore −0.505 ≈ −1/2 as was expected a priori since, in practice, the absolute error behaves as O(M −1/2 ) as M → ∞ (see the Step 2 in the Introduction). A more detailed explanation of the second step, about the improvement of the Monte Carlo estimations using the least squares method, is as follows: the probabilities pk satisfy some affine relations, like the ones in Lemma 2 or the ones in Proposition 12 below. Then, if we denote p = (p0 , . . . , pn )t ∈ Rn+1 it is possible to write p = Mq + b where q ∈ Rk with k ≤ n is a vector whose components are different elements of p0 , . . . , pn ; M ∈ Mn×k (R); and b ∈ Rk . Let p p0 , . . . , pen )t be the vector with the estimated probabilities obtained e = (e by the observed relative frequencies using the Monte Carlo method. Then, we can find the least squares solution [20, Def. 6.1] of the the system, e = Mb p q + b, (8) 12 which is b = (Mt · M)−1 · Mt · (e q p − b), (9) b , via see [9, p. 198] or [22, p. 200]. So we can find some improved estimations p b = Mb p q + b. (10) Some detailed examples are given in Sections 4, 5 and 7. 4 Linear random differential systems Consider linear differential systems ẋ = A x where x ∈ Rn , A ∈ Mn×n (R), where A is a random matrix whose entries are i.i.d. random variables with N(0, 1) distribution. Let X be the random variable that counts the number of eigenvalues of A with negative real part, s(A). Proposition 3. With the above notations, set pk = P (X = k). The following holds: Pn (a) k=0 pk = 1. (b) For all k ∈ {0, 1, . . . , n}, pk = pn−k . P P (c) i odd pi = i even pi = 21 . Proof. The assertion (a) is trivial. To prove (b) we observe that if a matrix A has k eigenvalues with negative real part, then B = −A has n − k eigenvalues with negative real part. Calling qm the probability that B has m eigenvalues with negative real part, we get that pm = qm . This is so, because if X ∼ N(0, 1) then −X ∼ N(0, 1) and as a consequence the entries of A and B have the same distribution. Then, qk = pn−k and the result follows. To see (c) we claim that s(A) is even if and only if the determinant of A is positive P and, moreover, P (det(A) > 0) = 1/2. From this claim we get the result because i even pi is the probability of s(A) being even. To prove the first part of the claim notice first that we can assume that 0 6= det(A) = λ1 λ2 · · · λn , where λ1 , λ2 , . . . , λn are the n eigenval- ues of A. We write λ1 λ2 · · · λn = (λ1 λ2 · · · λk )(λk+1 λk+2 · · · λn ) where λ1 , λ2 , . . . , λk are all the real negative eigenvalues. Observe also that for complex eigenvalues λλ̄ > 0. Hence λk+1 λk+2 · · · λn > 0, sign(det(A)) = (−1)k and the condition that s(A) is even is charac- terized by det(A) > 0. To prove that P (det(A) > 0) = 1/2 note that if B is the matrix obtained by changing the sign of one column of A then det(A) · det(B) < 0 and hence P (det(A) < 0) = P (det(B) > 0). Furthermore, since the entries of A and B have the same distribution we have P (det(B) > 0) = P (det(A) > 0), thus P (det(A) < 0) = P (det(A) > 0) = 1/2. 13 From the above proposition it easily follows: Corollary 4. Consider ẋ = A x, x ∈ Rn with A ∈ Mn×n (R) a random matrix with i.i.d. N(0, 1) entries, let X be the random variable defined above and pk = P (X = k). Then the probabilities pk satisfy all the consequences of Lemma 2. In particular E(X) = n/2. Now we reproduce some experiments to estimate the probabilities pk for low dimensional cases. We apply the Monte Carlo method, that is, for each considered dimension n, we have generated 108 matrices A ∈ Mn×n (R) whose entries are pseudo-random numbers simulat- ing the realizations on n2 independent random variables with N(0, 1) distribution. For each matrix A we have computed the characteristic polynomial, and counted the number of eigen- values with negative real part by using the Routh-Hurwitz zeros counter [10, p. 1076]. We are aware that the stability of the calculation of the coefficients of the characteristic poly- nomial from the entries of a matrix is critical (see [SB, p.378-379] and references therein); however we have only used this calculation for low dimensions, namely n ≤ 4. For n ≥ 5, and in order to decrease the computation time, we have directly computed numerically the eigenvalues of A and counted the number of them with negative real part. 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 with these relations when using the observed frequencies obtained by the Monte Carlo simulation. We give details of one example. Set n = 7, for instance. By Corollary 4 we have 1 p3 = p4 = 2 − p0 − p1 − p2 ; p5 = p2 ; p6 = p1 and p7 = p0 . So, using the notation introduced in Section 3.2, we can write p = Mq + b, where pt = (p0 , . . . , p7 );     1 0 0 0      0 10   0       0 0 1     0        p0  1   −1 −1 −1      M=  ; q= p1    ; and b =  2  .  −1 −1 −1     1  2    p2    0 0 1   0           0 1 0   0  1 0 0 0 The observed relative frequencies in our Monte Carlo simulation are   31643 261137 7124967 1344047 33597117 14248187 1043913 63379 et = p , , , , , , , . 50000000 12500000 50000000 4000000 100000000 100000000 50000000 100000000 By finding the least squares solution of the system (8) ([9, p. 198] or [22, p. 200]), given by 14 (10), we obtain   25333 2088461 28498121 16799573 16799573 28498121 2088461 25333 bt = p , , , , , , , . 40000000 100000000 200000000 50000000 50000000 200000000 100000000 40000000 The other cases follow similarly. We summarize the results of our experiments in the Table 2, 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 predicted absolute error will be of order 10−4 . Observe that in the cases n = 1, 2 the true probabilities are known. We include the results of the Monte Carlo simulations for completeness, but it makes no sense to apply the least squares method. Dimension Observed frequency Least squares Relations (Corol. 4) n=1 pe0 = 0.49996 p0 = 0.5 pe1 = 0.50004 p1 = 0.5 n=2 pe0 = 0.24999 p0 = 0.25 pe1 = 0.50006 p1 = 0.5 pe2 = 0.24995 p2 = 0.25 n=3 pe0 = 0.10447 pb0 = 0.10450 p0 1 pe1 = 0.39542 pb1 = 0.39550 p1 = 2 − p0 1 pe2 = 0.39557 pb2 = 0.39550 p2 = 2 − p0 pe3 = 0.10454 pb3 = 0.10450 p3 = p0 n=4 pe0 = 0.03722 pb0 = 0.03721 p0 1 pe1 = 0.25009 pb1 = 0.25000 p1 = 4 1 pe2 = 0.42556 pb2 = 0.42558 p2 = 2 − 2p0 1 pe3 = 0.24998 pb3 = 0.25000 p3 = 4 pe4 = 0.03715 pb4 = 0.03721 p4 = p0 n=5 pe0 = 0.01126 pb0 = 0.01126 p0 pe1 = 0.13028 pb1 = 0.13024 p1 1 pe2 = 0.35848 pb2 = 0.35850 p2 = 2 − p0 − p1 1 pe3 = 0.35852 pb3 = 0.35850 p3 = 2 − p0 − p1 pe4 = 0.13020 pb4 = 0.13024 p4 = p1 pe5 = 0.01126 pb5 = 0.01126 p5 = p0 15 Dimension Observed frequency Least squares Relations (Corol. 4) n=6 pe0 = 0.00289 pb0 = 0.00288 p0 pe1 = 0.05675 pb1 = 0.05678 p1 1 pe2 = 0.24710 pb2 = 0.24712 p2 = 4 − p0 1 pe3 = 0.38642 pb3 = 0.38644 p3 = 2 − 2p1 1 pe4 = 0.24714 pb4 = 0.24712 p4 = 4 − p0 pe5 = 0.56810 pb5 = 0.05678 p5 = p1 pe6 = 0.00289 pb6 = 0.00288 p6 = p0 n=7 pe0 = 0.00063 pb0 = 0.00063 p0 pe1 = 0.02089 pb1 = 0.02088 p1 pe2 = 0.14250 pb2 = 0.14249 p2 1 pe3 = 0.33601 pb3 = 0.33600 p3 = 2 − p0 − p1 − p2 1 pe4 = 0.33597 pb4 = 0.33600 p4 = 2 − p0 − p1 − p2 pe5 = 0.14248 pb5 = 0.14249 p5 = p2 pe6 = 0.02088 pb6 = 0.02088 p6 = p1 pe7 = 0.00063 pb7 = 0.00063 p7 = p0 n=8 pe0 = 0.00012 pb0 = 0.00012 p0 pe1 = 0.00651 pb1 = 0.00650 p1 pe2 = 0.06948 pb2 = 0.06948 p2 1 pe3 = 0.24356 pb3 = 0.24350 p3 = 4 − p1 1 pe4 = 0.36080 pb4 = 0.36080 p4 = 2 − 2p0 − 2p2 1 pe5 = 0.24346 pb5 = 0.24350 p5 = 4 − p1 pe6 = 0.06946 pb6 = 0.06948 p6 = p2 pe7 = 0.00650 pb7 = 0.00650 p7 = p1 pe8 = 0.00012 pb8 = 0.00012 p8 = p0 n=9 pe0 = 0.00002 pb0 = 0.00002 p0 pe1 = 0.00171 pb1 = 0.00171 p1 pe2 = 0.02880 pb2 = 0.02879 p2 pe3 = 0.14952 pb3 = 0.14955 p3 1 pe4 = 0.31987 pb4 = 0.31993 p4 = 2 − p0 − p1 − p2 − p3 1 pe5 = 0.31999 pb5 = 0.31993 p5 = 2 − p0 − p1 − p2 − p3 pe6 = 0.14958 pb6 = 0.14955 p6 = p3 pe7 = 0.02878 pb7 = 0.02879 p7 = p2 pe8 = 0.00171 pb8 = 0.00171 p8 = p1 pe9 = 0.00002 pb9 = 0.00002 p9 = p0 16 Dimension Observed frequency Least squares Relations (Corol. 4) n = 10 pe0 = 0 pb0 = 0 p0 pe1 = 0.00038 pb1 = 0.00038 p1 pe2 = 0.01015 pb2 = 0.01015 p2 pe3 = 0.07850 pb3 = 0.07850 p3 1 pe4 = 0.23987 pb4 = 0.23985 p4 = 4 − p0 − p2 1 pe5 = 0.34224 pb5 = 0.34224 p5 = 2 − 2p1 − 2p3 1 pe6 = 0.23984 pb6 = 0.23985 p6 = 4 − p0 − p2 pe7 = 0.07849 pb7 = 0.07850 p7 = p3 pe8 = 0.01015 pb8 = 0.01015 p8 = p2 pe9 = 0.00038 pb9 = 0.00038 p9 = p1 pe10 = 0 pb10 = 0 p10 = p0 Table 2. Linear stability indexes for linear random differential systems. 5 Linear random differential equations of order n In this section we consider linear random homogeneous differential equations of order n An x(n) + An−1 x(n−1) + · · · + A2 x00 + A1 x0 + A0 x = 0, (11) where x = x(t), the derivatives are taken in respect to t, and Aj are again i.i.d. random variables with N(0, 1) distribution. To get the stability index for these differential equations, we need we only need to know the probability distributions of the number of roots with negative real part of its associated random characteristic polynomial: Q(λ) = An λn + An−1 λn−1 + · · · + A1 λ + A0 . Let X be the random variable that counts the number of roots of Q(λ) with negative real parts and define pk = P (X = k) for k = 0, 1, . . . , n. Proposition 5. Set pk = P (X = k), where X is the random variable defined above. Then Pn (a) k=0 pk = 1. (b) For all k ∈ {0, 1, . . . , n}, pk = pn−k . P P (c) i odd pi = i even pi = 21 . 17 Proof. The proof of (a) is trivial. To prove (b) consider equation (11) with its characteristic polynomial Q(λ) and also the new differential equation (−1)n An x(n) + (−1)n−1 An−1 x(n−1) + · · · + A2 x00 − A1 x0 + A0 x = 0 (12) with its characteristic polynomial Q∗ (λ) = Q(−λ) = (−1)n An λn + (−1)n−1 An−1 λn−1 + · · · − A1 λ + A0 . Since Q(λ) = 0 if and only if Q∗ (−λ) = 0 we get that pk = p∗n−k where p∗i the probability that Q∗ (λ) has i roots with negative real part. But also pk = p∗k because the coefficients of the equations (11) and (12) have the same distributions. Hence, the result follows. 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 A0 · An < 0, because we can neglect the case of having some roots with zero real part. Since the coefficients of (11) are symmetric independent random variables, the probability that Q(λ) has an odd number of roots with negative real part is 1 1 1 1 1 P ({A0 > 0} ∩ {An < 0}) + P ({A0 < 0} ∩ {An > 0}) = × + × = . 2 2 2 2 2 Corollary 6. Consider the linear random homogeneous differential equation of order n (11), with all Ai being i.i.d. N(0, 1) random variables, let X be defined above, and set pk = P (X = k). Then the probabilities pk satisfy all the conclusions of Lemma 2. In particular E(X) = n/2. For each n, let rn be the probability of the origin to be a global stable attractor (asymp- totically stable equilibrium) for (11), that is rn = pn . By Proposition 5(b) this probability coincides with the probability of being a repeller because pn = p0 . Our results in Proposi- tion 8 seem to indicate that rn decreases with n. Before proving this proposition we need a preliminary result. Lemma 7. Let U, V, S and T be i.i.d. random variables with standard normal distribution. Then p+ := P (U > 0; V > 0; S > 0; T > 0; U T − SV > 0) = 1/32. Proof. Set A± = {U > 0; V > 0; S > 0; T > 0; ±(U T − SV ) > 0}, and A0 = {U > 0; V > 0; S > 0; T > 0; U T −SV = 0}. Denote by p± = P (A± ) and p0 = P (A0 ). Then, since p0 = 0 and A− ∪ A0 ∪ A+ = {U > 0; V > 0; S > 0; T > 0} it holds that p+ + p− = (1/2)4 = 1/16. To end the proof it suffices to show that p+ = p− . Notice first that A+ = {U > 0; V > 0; S > 0; T > 0; U T − SV > 0} = {V > 0; S > 0; T > 0; U T − SV > 0}, A− = {U > 0; V > 0; S > 0; T > 0; U T − SV < 0} = {U > 0; S > 0; T > 0; SV − U T > 0}. 18 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. Proposition 8. With the above notations, rn ≤ 1/2n , for all n ≥ 1. Moreover, r1 = 1/2, r2 = 1/4, r3 = 1/16 and r4 < r3 /2 = 1/32. Proof of Proposition 8. Notice that rn is the probability that the characteristic polynomial Q(λ), associated with the random differential equation (11), is a Hurwitz stable polynomial; that is rn = 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,  A0 , . . . , An such that all roots of P (λ) are in R− ⊂ (n ) (n ) \ [ \ {Ai < 0} {Ai > 0} . (13) i=0 i=0 Since the variables Ai are independent and symmetric n ! n ! \ \ 1 P {Ai < 0} = P {Ai > 0} = . 2n+1 i=0 i=0 As a consequence, n ! n ! \ \ 1 rn ≤ P {Ai < 0} +P {Ai > 0} = , 2n i=0 i=0 and the first statement follows. The equalities r1 = 1/2 and r2 = 1/4 are a simple consequence that for n = 1, 2 the inclusion (13) is an equality. Let us prove that r3 = p3 = 1/16. By using the Routh-Hurwitz criterion [10, p. 1076], it can be seen that a3 λ3 + a2 λ2 + a1 λ + a0 has every root in R− if and only if all its coefficients have the same sign and moreover a1 a2 − a0 a3 > 0. Hence, p3 = p− + 3 + p3 , where p− + 3 := P (A0 < 0; A1 < 0; A2 < 0; A3 < 0; A1 A2 − A0 A3 > 0); and p3 := P (A0 > 0; A1 > 0; A2 > 0; A3 > 0; A1 A2 − A0 A3 > 0), with all the Ai being N(0, 1) distributed and independent. Due to their symmetry, the random variables Ai and −Ai , for i = 0, . . . , 3 19 − have the same distribution and hence p+ + 3 = p3 . Therefore p3 = 2p3 . The result follows now by Lemma 7, which gives p+ 3 = 1/32. Let us prove that r4 < r3 /2. To compare both probabilities, here it will be more conve- nient to write the coefficients of the polynomials with subscripts with increasing ordering, that is qn (x) = a0 xn + a1 xn−1 + · · · + an−1 x + an . With this notation, which also respects the traditional notation when writing the Hurwitz matrices, and when a0 > 0, the Routh- Hurwitz conditions to have stability index n for n = 3, 4 are precisely that the principal minors of the following matrices     a1 a3 0 0 a a 0    1 3   a0 a2 a4 0   a0 a2 0  and  ,    0 a a 0   1 3  0 a1 a3 0 a0 a2 a4 are positive, where the left-hand one corresponds to the case n = 3 and the other to the case when n = 4. Hence, these conditions when a0 > 0 and for n = 3 are: a1 > 0, a1 a2 −a0 a3 > 0 and a3 > 0. Similarly, for n = 4 the conditions are a1 > 0, a1 a2 −a0 a3 > 0, a3 (a1 a2 − a0 a3 ) − a4 a21 > 0 and a4 > 0. Consider now, for n = 3, 4, the random polynomials Qn (x) = Ã0 xn + Ã1 xn−1 + · · · + Ãn−1 x + Ãn , where Ãi ∼ N(0, 1) and are independent (notice that with this notation each coefficient Ãi is the coefficient An−i of the characteristic polynomial). For simplicity we denote with the same name the coefficients of Q3 and Q4 although they are different random variables. As above, r3 = 2p+ + + + 3 , and r4 = 2p4 , where pk = P (Ak ), with A+ 3 = {Ã0 > 0; Ã1 > 0; Ã3 > 0; Ã1 Ã2 − Ã0 Ã3 > 0}, A+ 2 4 = {Ã0 > 0; Ã1 > 0; Ã3 > Ã4 Ã1 /(Ã1 Ã2 − Ã0 Ã3 ); Ã1 Ã2 − Ã0 Ã3 > 0, Ã4 > 0}. Notice that if we define B = {Ã0 > 0; Ã1 > 0; Ã3 > 0; Ã1 Ã2 − Ã0 Ã3 > 0; Ã4 > 0} it is clear that P (B) = p+ + 3 /2 and, moreover A4 ⊂ B3 , with the inclusion being strict. Since the joint density is positive and B ∩ (A+ c 4 ) has positive Lebesgue measure, we have P (B ∩ (A+ c + + + + 4 ) ) > 0. Thus P (A4 ) < P (B), and hence p4 = P (A4 ) < P (B) = p3 /2, and r4 < r3 /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 p0 = p3 = 1/16 and p1 = p2 = 7/16. Proof. By the above proposition, for n = 3, p0 = p3 = r3 = 1/16. Hence, by Proposition 5, p1 = p2 = 7/16. 20 The computations in this case are similar to the ones of the previous section and the obtained results are summarized in Table 3. We only give some comments for the cases b have negative and very small n = 8 and 10, where we have encountered that the vectors p entries. This has occurred because the observed frequencies obtained by the Monte Carlo approach corresponding to these probabilities are not accurate enough. For this reason, we have made a new optimization step. As before, we use the least squares method to obtain a b . However, if negative entries appear in this vector (which is clearly objectionable), vector p we force them to be zero and find a new least squares estimate, which still respects the original linear constraints. We explain this process for the n = 8 order case: The observed relative frequencies vector obtained by the Monte Carlo method is 1 6599 1159359 4996163 45377377 et = p , , , , , 50000000 50000000 50000000 20000000 100000000 ! 4995607 2318357 13497 1 , , , . 20000000 100000000 100000000 100000000 The relations stated in Corollary 6 are p3 = 1/4 − p1 , p4 = 1/2 − 2p0 − 2p2 , p5 = p3 , p6 = p2 , p7 = p1 , p8 = p0 . By solving the system (8) with     1 0 0 0      00 1  0           0 0 1     0       0 −1 0  p0  1       2     1  M= −2 0 −2  ; q=  p1   ; and b =  4       0 −1 0  p2  1     2       0 0 1   0       0 1 0   0      1 0 0 0 we obtain   t 5779 13569 4631293 b = − q , , . 200000000 80000000 200000000 Hence, by (10) we get −5779 13569 4631293 19986431 22687243 b= p , , , , , 200000000 80000000 200000000 80000000 50000000 ! 19986431 4631293 13569 −5779 , , , . 80000000 200000000 80000000 200000000 21 So we impose that p0 = p8 = 0. Thus we have p3 = p5 = 1/4 − p1 , p4 = 1/2 − 2p0 − 2p2 = 1/2 − 2p2 , p6 = p2 and p7 = p1 . We find the least squares solution of the system       pb1 1 0 0        pb2   0 1   0            !  1   pb3   −1 0        pb∗1  4   pb4  =  0 −2  · + 2   1    pb∗2 .      1   pb5   −1 0   4         pb   0 1    0   6     pb7 1 0 0 Using (9) and (10) we obtain   ∗ 13569 13882321 19986431 136117679 19986431 13882321 13569 b = p 0, , , , , , , ,0 80000000 600000000 80000000 300000000 80000000 600000000 80000000 ' (0, 0.00017, 0.02314, 0.24983, 0.45373, 0.24983, 0.02314, 0.00017, 0) . The n = 10 case follows analogously. Dimension Observed frequency Least squares Relations (Corol. 6 and 9) n=1 pe0 = 0.49997 p0 = 0.5 pe1 = 0.50003 p1 = 0.5 n=2 pe0 = 0.24994 p0 = 0.25 pe1 = 0.49999 p1 = 0.5 pe2 = 0.25007 p2 = 0.25 1 n=3 pe0 = 0.06252 p0 = 16 = 0.0625 7 pe1 = 0.43743 p1 = 16 = 0.4375 7 pe2 = 0.43756 p2 = 16 = 0.4375 1 pe3 = 0.06249 p3 = 16 = 0.0625 n=4 pe0 = 0.00928 pb0 = 0.00925 p0 1 pe1 = 0.24998 pb1 = 0.25 p1 = 4 1 pe2 = 0.48152 pb2 = 0.48150 p2 = 2 − 2p0 1 pe3 = 0.24994 pb3 = 0.25 p3 = 4 pe4 = 0.00929 pb4 = 0.00925 p4 = p0 n=5 pe0 = 0.00071 pb0 = 0.00071 p0 pe1 = 0.08405 pb1 = 0.08404 p1 1 pe2 = 0.41526 pb2 = 0.41525 p2 = 2 − p0 − p1 1 pe3 = 0.41523 pb3 = 0.41525 p3 = 2 − p0 − p1 pe4 = 0.08404 pb4 = 0.08404 p4 = p1 pe5 = 0.00071 pb5 = 0.00071 p5 = p0 22 Dimension Observed frequency Least squares Relations (Corol. 6 and 9) n=6 pe0 = 0.00003 pb0 = 0.00005 p0 pe1 = 0.01723 pb1 = 0.01720 p1 1 pe2 = 0.24994 pb2 = 0.24995 p2 = 4 − p0 1 pe3 = 0.46562 pb3 = 0.46560 p3 = 2 − 2p1 1 pe4 = 0.24993 pb4 = 0.24995 p4 = 4 − p0 pe5 = 0.01723 pb5 = 0.01720 p5 = p1 pe6 = 0.00003 pb6 = 0.00005 p6 = p0 n=7 pe0 = 0 pb0 = 0 p0 pe1 = 0.00200 pb1 = 0.00200 p1 pe2 = 0.09571 pb2 = 0.09572 p2 1 pe3 = 0.40224 pb3 = 0.40228 p3 = 2 − p0 − p1 − p2 1 pe4 = 0.40233 pb4 = 0.40228 p4 = 2 − p0 − p1 − p2 pe5 = 0.09573 pb5 = 0.09572 p5 = p2 pe6 = 0.00199 pb6 = 0.00200 p6 = p1 pe7 = 0 pb7 = 0 p7 = p0 n=8 pe0 = 0 pb∗0 =0 p0 pe1 = 0.00013 pb∗1 = 0.00017 p1 pe2 = 0.02319 pb∗2 = 0.02314 p2 1 pe3 = 0.24981 pb∗3 = 0.24983 p3 = 4 − p1 1 pe4 = 0.45377 pb∗4 = 0.45372 p4 = 2 − 2p0 − 2p2 1 pe5 = 0.24978 pb∗5 = 0.24983 p5 = 4 − p1 pe6 = 0.02318 pb∗6 = 0.02314 p6 = p2 pe7 = 0.00013 pb∗7 = 0.00017 p7 = p1 pe8 = 0 pb∗8 =0 p8 = p0 n=9 pe0 = 0 pb0 = 0 p0 pe1 = 0.00001 pb1 = 0.00005 p1 pe2 = 0.00336 pb2 = 0.00337 p2 pe3 = 0.10337 pb3 = 0.10335 p3 1 pe4 = 0.39328 pb4 = 0.39328 p4 = 2 − p0 − p1 − p2 − p3 1 pe5 = 0.39328 pb5 = 0.39328 p5 = 2 − p0 − p1 − p2 − p3 pe6 = 0.10332 pb6 = 0.10335 p6 = p3 pe7 = 0.00338 pb7 = 0.00337 p7 = p2 pe8 = 0 pb8 = 0.00005 p8 = p1 pe9 = 0 pb9 = 0 p9 = p0 23 Dimension Observed frequency Least squares Relations (Corol. 6 and 9) n = 10 pe0 = 0 pb∗0 = 0 p0 pe1 = 0 pb∗1 = 0.00002 p1 pe2 = 0.00030 pb∗2 = 0.00028 p2 pe3 = 0.02784 pb∗3 = 0.02787 p3 1 pe4 = 0.24976 pb∗4 = 0.24972 p4 = 4 − p0 − p2 1 pe5 = 0.44421 pb∗5 = 0.44422 p5 = 2 − 2p1 − 2p3 1 pe6 = 0.24973 pb∗6 = 0.24972 p6 = 4 − p0 − p2 pe7 = 0.02787 pb∗7 = 0.02787 p7 = p3 pe8 = 0.00029 pb∗8 = 0.00028 p8 = p2 pe9 = 0 pb∗9 = 0.00002 p9 = p1 pe10 = 0 pb∗10 = 0 p10 = p0 Table 3. Stability indexes for order n linear random homogeneous differential equations. 6 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 B xk+1 = A xk where x ∈ Rn , (14) where B and each of the n2 entries of the random matrix A are i.i.d. N(0, 1) random variables. Observe that to ensure that the results are invariant under time-scaling, is necessary to add the term B in the left-hand side of Equation (14). Then, given a linear discrete random 1 system (14), its characteristic random polynomial associated with the matrix BA is Q(λ) = Qn λn + Qn−1 λn−1 + · · · + Q1 λ + Q0 where each random variable Qj is a polynomial in the variables 1/B, A1,1 , . . . , An,n which has a complicated distribution function. 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 1 the matrix B A. Also pk denotes the probabilities that X takes the value k. As we will see in the examples, in this case the condition pk = pn−k is no longer satisfied. Among other reasons it happens that the entries of A−1 have complicated distributions. Since we do not know other relations on the probabilities pk apart from the trivial one Pn k=0 pk = 1, and this is directly fulfilled by the observed relative frequencies, in this case we do not perform the least squares refinement. 24 The case n = 1 is the only one that we have been able to solve analytically. Notice that in this situation the only solution of Q(λ) = 0 is λ = A/B, with A and B independent and N(0, 1). Hence p0 = P (|A/B| > 1) and p1 = P (|A/B| < 1) = P (|B/A| > 1). Since A/B and B/A have the same distribution it holds that p0 = p1 = 1/2. As in the other models, for each dimension n ≤ 10, we generate 108 discrete systems 1 of the form (14). For each matrix BA we have computed the characteristic polynomial Q and its associated polynomial Q? (see Equation (7)) 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 4. Dimension Observed frequency Dimension Observed frequency n=1 pe0 = 0.49994 n=6 pe0 = 0.44944 pe1 = 0.50006 pe1 = 0.24419 n=2 pe0 = 0.46348 pe2 = 0.13838 pe1 = 0.27705 pe3 = 0.07536 pe2 = 0.25947 pe4 = 0.04606 n=3 pe0 = 0.45261 pe5 = 0.02449 pe1 = 0.25828 pe6 = 0.02209 pe2 = 0.15351 n=7 pe0 = 0.44937 pe3 = 0.13560 pe1 = 0.24394 n=4 pe0 = 0.45040 pe2 = 0.13723 pe1 = 0.24732 pe3 = 0.07480 pe2 = 0.14799 pe4 = 0.04226 pe3 = 0.08127 pe5 = 0.02636 pe4 = 0.07302 pe6 = 0.01367 n=5 pe0 = 0.44957 pe7 = 0.01236 pe1 = 0.24536 pe2 = 0.13956 pe3 = 0.08116 pe4 = 0.04443 pe5 = 0.03992 25 Dimension Observed frequency Dimension Observed frequency n=8 pe0 = 0.44937 n = 10 pe0 = 0.44934 pe1 = 0.24381 pe1 = 0.24371 pe2 = 0.13702 pe2 = 0.13687 pe3 = 0.07388 pe3 = 0.07358 pe4 = 0.04207 pe4 = 0.04129 pe5 = 0.02394 pe5 = 0.02348 pe6 = 0.01526 pe6 = 0.01388 pe7 = 0.00768 pe7 = 0.00792 pe8 = 0.00698 pe8 = 0.00520 n=9 pe0 = 0.44941 pe9 = 0.00247 pe1 = 0.24374 pe10 = 0.00226 pe2 = 0.13680 pe3 = 0.07371 pe4 = 0.04139 pe5 = 0.02400 pe6 = 0.01374 pe7 = 0.00889 pe8 = 0.00434 pe9 = 0.00397 Table 4. Stability indexes for linear random maps. 7 Linear random difference equations of order n Finally we consider difference equations of order n of type An xk+n + An−1 xk+n−1 + · · · + A1 xk+1 + A0 xk = 0, where all the coefficients are i.i.d. random variables with N(0, 1) distribution. In this situa- tion, the stability index is given by the number of zeros with modulus smaller than 1 of the random characteristic polynomial Q(λ) = An λn + An−1 λn−1 + · · · + A1 λ + A0 . 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 pk = P (X = k) for k = 0, 1, . . . , n. Before proving some relations among the probabilities pk , we give two preliminary lem- Rx 2 mas. Let erf(x) = √2π 0 e−u du be the error function. The following result is stated in [4]. We prove it for the sake of completeness. 26 Lemma 10. For α > 0 and β ∈ R, Z ∞ 2 2 arctan(β/α) F (α, β) := e−α x erf(βx) dx = √ . 0 α π Proof. Fixed α > 0, the function that defines F is absolutely integrable because | erf(x)| ≤ 1. Moreover its partial derivative with respect to β is also absolutely integrable. Hence limβ→0 F (α, β) = F (α, 0) = 0 and Z ∞ Z ∞ ∂F (α, β) ∂  −α2 x2  2 2 2 2 2 = e erf(βx) dx = √ x e−α x e−β x dx ∂β 0 ∂β π 0 Z ∞ 2 2 2 2 1 =√ x e−(α +β ) x dx = 2 √ . π 0 (α + β 2 ) π Therefore Z β Z β ∂F (α, t) 1 arctan(β/α) F (α, β) = F (α, 0) + dt = √ dt = √ , 0 ∂t 0 (α2 2 +t ) π α π as we wanted to prove. The next result is a consequence of the previous lemma. Lemma 11. Let U ∼ N(0, σ 2 ) and V ∼ N(0, ρ2 ) be independent normal random variables. 2 Then P (U 2 − V 2 > 0) = π arctan(σ/ρ). Proof. The joint density function of the random vector (U, V ) is fσ (u)fρ (v), where fs (u) = 2 2 √ e−u /(2s ) /( 2πs). Observe that the points (u, v) ∈ R2 such that u2 − v 2 > 0 is the region where −|u| < v < |u|, hence by symmetry, Z ∞ Z u Z ∞ Z u 4 −u2 /(2σ 2 ) 2 2 2 2 P (U − V > 0) = 4 fσ (u) fρ (v) dv du = e e−v /(2ρ ) dv du 0 0 2πσρ 0 0 2 Z ∞ 2 2  u r π 2 σ  = e−u /(2σ ) erf √ ρ du = arctan , πσρ 0 2ρ 2 π ρ where in the last equality we have used Lemma 10. Notice that with the notation of the above lemma, P (U 2 −V 2 > 0)+P (U 2 −V 2 < 0) = 1. Hence σ  2 ρ 2 P (U 2 − V 2 < 0) = 1 − arctan = arctan , (15) π ρ π σ where we have used the fact that arctan(x)+arctan(1/x) = π/2 or, simply, the same lemma interchanging U and V. Observe also that when σ = ρ, P (U 2 −V 2 > 0) = P (U 2 −V 2 < 0) = 1/2, a result that, in fact, is a straightforward consequence that in this situation U 2 − V 2 and V 2 − U 2 have the same distribution. Proposition 12. With the above notation: 27 Pn (a) k=0 pk = 1. (b) For all k ∈ {0, 1, . . . , n}, pk = pn−k . P P (c) When n is odd, i even pi = i odd pi = 12 . (d) When n = 2k is even,  r r X 2 k+1 X 2  k  pi = arctan and pi = arctan . (16) π k π k+1 i even i odd Proof. The first assertion is obvious. To see the second one we compare the difference equation an xk+n + an−1 xk+n−1 + · · · + a1 xk+1 + a0 xk = 0, with ai ∈ R, i = 0, 1, . . . , n, with characteristic polynomial Q(λ) = an λn +an−1 λn−1 +· · ·+a2 λ2 +a1 λ+a0 with the difference equation an xk + an−1 xk+1 + · · · + a0 xk+n = 0 with characteristic polynomial Q∗ (λ) = an + an−1 λ + · · · + a1 λn−1 + a0 λn . Notice that if Q(λ) has m non-zero roots with modulus smaller than 1 and n − m with modulus bigger than 1, then the converse follows for Q∗ (λ) because Q(λ) = 0 if and only if Q∗ ( λ1 ) = 0. From this result applied to the corresponding random polynomials we get that pk = pn−k , because both have identically distributed coefficients. So we have proved statement (b). To prove items (c) and (d) recall first that it was proved in item (c) of Proposition 5 that a polynomial Q(λ) = an λn + an−1 λn−1 + · · · + a2 λ2 + a1 λ + a0 , without roots with zero real part, has an even number of roots with negative real part if and only if an a0 > 0. By using the polynomial Q? (z) = an (z + 1)n + an−1 (z + 1)n−1 (z − 1) + . . . + a0 (z − 1)n = (an + an−1 + · · · + a1 + a0 )z n + · · · + (an − an−1 + an−2 − · · · + (−1)n a0 ), introduced in Section 3 (Equation (7)) we get that Q(λ), without roots of modulus 1, has an even number (2m) of roots with modulus smaller than 1 if and only if Q? (z) has exactly 2m roots with negative real part and this happens if and only if (an + an−1 + · · · + a1 + a0 ) · (an − an−1 + an−2 − · · · + (−1)n a0 ) > 0. Hence, considering the corresponding random polynomials, we have that X  pi = P (An + An−1 + · · · + A0 ) · (An − An−1 + · · · + (−1)n A0 ) > 0 i even = P (U 2 − V 2 > 0), where U = An + An−2 + An−4 + · · · and V = An−1 + An−3 + An−5 + · · · and the sums end either at A0 or A1 according the parity of n. Since Aj ∼ N(0, 1) and all Aj are independent we get that when n = 2k (resp. n = 2k − 1) then U ∼ N(0, k + 1) (resp. U ∼ N(0.k)) and 28 V ∼ N(0, k) and U and V are independent. Hence, by using Lemma 11, we obtain that when n = 2k − 1, P (U 2 − V 2 > 0) = 1/2 and that when n = 2k, X 2 r k + 1  2 2 pi = P (U − V > 0) = arctan . π k i even The sum of all pi when i is odd can be obtained from the above formula, see also (15). Corollary 13. (i) Consider the linear random homogeneous difference equation of order n, let X be the random variable defined above and pk = P (X = k). Then the probabilities pk satisfy all the conclusions of Lemma 2. In particular E(X) = n/2. (ii) Moreover the new affine relations given in Equations (16) hold. In particular, for √ n = 2, p0 = p2 = π1 arctan( 2) and p1 = π2 arctan( √12 ); and for n = 4, p1 = p3 = q 1 2 π arctan( 3 ). In this case, and for the situations where we have not been able to obtain the exact probabilities we have done similar computations than in the previous section, first with the Monte Carlo method, generating for each order n = 0, . . . , 10, 108 random vectors (A0 , . . . , An ) ∈ Rn+1 whose components are pseudo-random numbers with N(0, 1) distribu- tion. Then, by using the relations in Proposition 12 and Corollary 13 we have performed a least squares refinement. p For instance for n = 4, by Corollary 13 we have p1 = p3 = arctan( 2/3)/π ' 0.217953; p p2 = 2 arctan( 3/2)/π − 2p0 and p4 = p0 . Hence, we fix the values pb1 = p1 and pb3 = p3 and system (8) can be written in the form       1 0 pe0    q     Mbq+b=  p0 ) +   −2  · (b  2 π arctan 3 2  =  pe2  .    1 0 pe4 Hence we can easily find the least squares solution of the above incompatible linear system:       1 0 pe    q    0  (1, −2, 1) ·  −2  · (b   p0 ) +  2  π arctan 3 2    = (1, −2, 1) ·  pe2  ,   1 0 pe4 p  and thus we get the result in Equation (9): pb0 = 16 (e p0 − 2e 2 p2 + pe4 ) + 3π arctan 3/2 , p and therefore pb4 = pb0 and pb2 = 2 arctan( 3/2)/π − 2b p0 . Since our Monte Carlo simulations give   2056203 7169499 10285619 p0 , pe2 , pe4 ) = (e , , ' (0.10281, 0.35847, 0.10286) , 20000000 20000000 100000000 29 p0 , pb2 , pb4 ) ' (0.10282, 0.35846, 0.10282) . the above relations show that (b All our results are collected in Table 5. Order Observed frequency Least squares Relations (Prop. 12 and Cor. 13)) n=1 pe0 = 0.49991 p0 = 0.5 pe1 = 0.50009 p1 = 0.5 1 √ n=2 pe0 = 0.30410 p0 = π arctan( 2) ' 0.304087 2 pe1 = 0.39184 p1 = π arctan( √12 ) ' 0.391826 1 √ pe2 = 0.30406 p2 = π arctan( 2) ' 0.304087 n=3 pe0 = 0.17251 pb0 = 0.17248 p0 1 pe1 = 0.32752 pb1 = 0.32752 p1 = 2 − p0 1 pe2 = 0.32753 pb2 = 0.32752 p2 = 2 − p0 pe3 = 0.17244 pb3 = 0.17248 p3 = p0 n=4 pe0 = 0.10281 pb0 = 0.10282 p0 q 1 pe1 = 0.21792 pb1 = 0.21795 p1 = π arctan( 23 ) ' 0.217953 q 2 pe2 = 0.35847 pb2 = 0.35846 p2 = π arctan( 32 ) − 2p0 q 1 pe3 = 0.21794 pb3 = 0.21795 p3 = π arctan( 23 ) ' 0.217953 pe4 = 0.10286 pb4 = 0.10282 p4 = p0 n=5 pe0 = 0.05909 pb0 = 0.05909 p0 pe1 = 0.15331 pb1 = 0.15333 p1 1 pe2 = 0.28760 pb2 = 0.28758 p2 = 2 − p0 − p1 1 pe3 = 0.28756 pb3 = 0.28758 p3 = 2 − p0 − p1 pe4 = 0.15335 pb4 = 0.15333 p4 = p1 pe5 = 0.05908 pb5 = 0.05909 p5 = p0 n=6 pe0 = 0.03501 pb0 = 0.03502 p0 pe1 = 0.09726 pb1 = 0.09726 p1 q 1 pe2 = 0.23777 pb2 = 0.23779 p2 = π arctan( 43 ) − p0 q 2 pe3 = 0.25985 pb3 = 0.25986 p3 = π arctan( 34 ) − 2p1 q 1 pe4 = 0.23781 pb4 = 0.23779 p4 = π arctan( 43 ) − p0 pe5 = 0.09724 pb5 = 0.09726 p5 = p1 pe6 = 0.03505 pb6 = 0.03502 p6 = p0 30 Order Observed frequency Least squares Relations (Prop. 12 and Cor. 13) n=7 pe0 = 0.02025 pb0 = 0.02025 p0 pe1 = 0.06432 pb1 = 0.06430 p1 pe2 = 0.17174 pb2 = 0.17176 p2 1 pe3 = 0.24376 pb3 = 0.24369 p3 = 2 − p0 − p1 − p2 1 pe4 = 0.24361 pb4 = 0.24369 p4 = 2 − p0 − p1 − p2 pe5 = 0.17177 pb5 = 0.17176 p5 = p2 pe6 = 0.06428 pb6 = 0.06430 p6 = p1 pe7 = 0.02025 pb7 = 0.02025 p7 = p0 n=8 pe0 = 0.01194 pb0 = 0.01196 p0 pe1 = 0.03994 pb1 = 0.03994 p1 pe2 = 0.12272 pb2 = 0.12726 p2 q 1 pe3 = 0.19230 pb3 = 0.19234 p3 = π arctan( 45 ) − p1 q 2 pe4 = 0.25701 pb4 = 0.25700 p4 = π arctan( 54 ) − 2p0 − 2p2 q 1 pe5 = 0.19238 pb5 = 0.19234 p5 = π arctan( 45 ) − p1 pe6 = 0.12724 pb6 = 0.12726 p6 = p2 pe7 = 0.03994 pb7 = 0.03994 p7 = p1 pe8 = 0.01197 pb8 = 0.01196 p8 = p0 n=9 pe0 = 0.00693 pb0 = 0.00693 p0 pe1 = 0.02556 pb1 = 0.02556 p1 pe2 = 0.08711 pb2 = 0.08711 p2 pe3 = 0.15653 pb3 = 0.15653 p3 1 pe4 = 0.22389 pb4 = 0.22386 p4 = 2 − p0 − p1 − p2 − p3 1 pe5 = 0.22382 pb5 = 0.22386 p5 = 2 − p0 − p1 − p2 − p3 pe6 = 0.15654 pb6 = 0.15653 p6 = p3 pe7 = 0.08712 pb7 = 0.08711 p7 = p2 pe8 = 0.02557 pb8 = 0.02556 p8 = p1 pe9 = 0.00693 pb9 = 0.00693 p9 = p0 31 Order Observed frequency Least squares Relations (Prop. 12 and Cor. 13) n = 10 pe0 = 0.00409 pb0 = 0.00411 p0 pe1 = 0.01567 pb1 = 0.01566 p1 pe2 = 0.06089 pb2 = 0.06091 p2 pe3 = 0.11500 pb3 = 0.11497 p3 q 1 pe4 = 0.19950 pb4 = 0.19947 p4 = π arctan( 65 ) − p0 − p2 q 2 pe5 = 0.20978 pb5 = 0.20976 p5 = π arctan( 56 ) − 2p1 − 2p3 q 1 pe6 = 0.19941 pb6 = 0.19947 p6 = π arctan( 65 ) − p0 − p2 pe7 = 0.11499 pb7 = 0.11497 p7 = p3 pe8 = 0.06088 pb8 = 0.06091 p8 = p2 pe9 = 0.01570 pb9 = 0.01566 p9 = p1 pe10 = 0.00408 pb10 = 0.00411 p10 = p0 Table 5. Stability indexes for order n linear random homogeneous difference equations. Acknowledgments We want to thank Professors Maria Jolis and Joan Torregrosa for their helpful indications. We also acknowledge the anonymous reviewer for the careful reading and for addressing to us very interesting suggestions that have allowed us to improve our work. References [1] J.C. Artés, J. Llibre. Statistical measure of quadratic vector fields. Resenhas 6 (2003), 85–97. [2] S. Asmussen, P.W. Glynn. Stochastic simulation: algorithms and analysis. Springer, New York, 2007. [3] P. Bratley, B.L. Fox, L. E. Schrage. A guide to simulation, 2nd ed. Springer, New York, 1987. [4] K. Briggs. Integrals involving erf. http://keithbriggs.info/documents/ erf-integrals.pdf (Accessed February 19, 2021). [5] B.M. Chen-Charpentier, D. Stanescu. Epidemic models with random coefficients. Math. Comput. Modelling 52 (2010), 1004–1010. 32 [6] B.M. Chen-Charpentier, D. Stanescu. Virus propagation with randomness. Math. Com- put. Modelling 57 (2013), 1816–1821. [7] A. Cima, A. Gasull, V. Mañosa. Phase portraits of random planar homogeneous vector fields. Qual. Theory Dyn. Syst. 20, 3 (2021). [8] K. Conrad. Probability distributions and maximum entropy. http://www.math.uconn.edu/~kconrad/blurbs/analysis/entropypost.pdf (Ac- cessed February 19, 2021). [9] G. Dahlquist, Å. Björck. Numerical methods. Prentice-Hall, Englewood Cliffs NJ, 1974. Section 5.7 pp. 196–201. [10] I.S. Gradshteyn, I.M. Ryzhik, Routh-Hurwitz Theorem. §15.715 in Tables of Integrals, Series, and Products, 6th ed. Academic Press, San Diego CA, 2000, p. 1076. [11] M. Grigoriu, T. Soong. Random vibration in mechanical and structural systems. Pren- tice Hall, New Jersey, 1993. [12] F.B. Jones. Connected and disconnected plane sets and the functional equation f (x) + f (y) = f (x + y). Bull. Amer. Math. Soc. 48 (1942), 115–120. [13] C. Lemieux. Monte Carlo and quasi-Monte Carlo sampling. Springer, New York, 2009. [14] H. Lopes, B. Pagnoncelli, C. Palmeira. Coeficientes aleatórios de equações diferenciais ordinárias lineares. Matemática Universitária 14 (2008), 44–50. Translated to English at http://bernardokp.uai.cl/preprint.pdf with title Random linear systems and simulation (Accessed February 19, 2021). [15] G. Marsaglia. Choosing a point from the surface of a sphere. The Annals of Mathe- matical Statistics 43 (1972), 645–646. [16] G. Marsaglia, W. Tsang. The Ziggurat method for generating random variables. Jour- nal of Statistical Software 5(8) (2000), 1–7. [17] M. Matsumoto, T. Nishimura. Mersenne twister: a 623-dimensionally equidistributed uniform pseudo-random number generator. ACM Trans. Model. Comput. Simul. 8, 1 (1998), 3–30. [18] B.J.T. Morgan. Applied stochastic modelling. Arnold Publishers, London, 2000. [19] M.E. Muller. A note on a method for generating points uniformly on N -dimensional spheres. Communications of the ACM 2 (1959), 19–20. 33 [20] J.R. Schott. Matrix analysis for statistics. John Wiley & Sons, New York, 1997. [21] D. Stanescu, B.M. Chen-Charpentier. Random coefficient differential equation models for bacterial growth. Math. Comput. Modelling 50 (2009), 885–895. [22] J. Stoer, R. Bulirsch. Introduction to numerical analysis, 3rd ed. Springer, New York, 2002. [23] S.H. Strogatz. Nonlinear dynamics and chaos. Westview Press, Cambridge MA, 1994. 34