Stability index of linear random dynamical systems

TL;DR

This paper investigates the stability index of linear random dynamical systems, providing methods to compute or estimate the probabilities of different stability indices using Monte Carlo and least squares techniques.

Contribution

It introduces a novel approach combining Monte Carlo simulations with least squares to estimate stability index probabilities for random dynamical systems.

Findings

Derived exact or estimated probabilities for stability indices.

Applied methods to linear differential and difference equations.

Analyzed stability properties in high-dimensional systems.

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 $p_k$ with $k=0,1,\ldots,n,$ denote the probabilities that $P(X=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,\ldots,n.$ The particular case of $n$-order homogeneous linear random differential or difference equations is also studied in detail.

Equations

$$\dot{\mathbf{x}}=A\mathbf{x}\text{mbox}or{\mathbf{x}}_{k+1}=A{\mathbf{x}}_k,$$

$$\dot{\mathbf{x}}=A\mathbf{x}+f(\mathbf{x}),\text{mbox}or{\mathbf{x}}_{k+1}=A{\mathbf{x}}_k+f({\mathbf{x}}_k),$$

$$\left(\begin{array}[]{c}\dot{x}\\ \dot{y}\end{array}\right)=\left(\begin{array}[]{cc}A&B\\ C&D\end{array}\right)\left(\begin{array}[]{c}x\\ y\end{array}\right)$$

$$\left(\frac{A}{S},\frac{B}{S},\frac{C}{S},\frac{D}{S}\right),\text{mbox}withS=\sqrt{A^2+B^2+C^2+D^2},$$

$$\left(\dfrac{X_{1}}{S},\dfrac{X_{2}}{S},\ldots,\dfrac{X_{n}}{S}\right),\quad\mbox{with}\quad S=\Big{(}\sum_{i=1}^{n}X_{i}^{2}\Big{)}^{1/2},$$

$$\alpha =P(AD-BC<0)=\frac{1}{4{\pi }^2}{\int }_{\mathcal{U}}{\mathrm{e}}^{-\frac{a^2+b^2+c^2+d^2}{2}}\mathrm{d}a\mathrm{d}b\mathrm{d}c\mathrm{d}d,$$

$$g(x_1,\dots ,x_n)=h(x_1^2+\cdots +x_n^2),$$

$$g_{R,\Theta}(r,\theta)=g(r\cos(\theta_{1}),\ldots,r\sin(\theta_{1})\cdots\sin(\theta_{n-1}))\,r^{n-1}\sin^{n-2}(\theta_{1})\sin^{n-1}(\theta_{2})\cdots\sin(\theta_{n-2})\cdot\mbox{\Large$\chi$}$$

$$\mbox{\Large$\chi$}:=\mbox{\Large$\chi$}_{[0,\infty)}(r)\cdot\,\mbox{\Large$\chi$}_{[0,2\pi)}(\theta_{n-1})\,\cdot\,\prod\limits_{i=1}^{n-2}\mbox{\Large$\chi$}_{[0,\pi)}(\theta_{i}),$$

$$g_{\Theta |R}(r,\theta ):=\frac{g_{R,\Theta }(r,\theta )}{g_R(r)},$$

$$g_R(r):={\int }_0^{\pi }\cdots {\int }_0^{\pi }{\int }_0^{2\pi }g(r\cos ({\theta }_1),\dots ,r\sin ({\theta }_1)\cdots \sin ({\theta }_{n-1}))\mathrm{d}\mathcal{S},$$

$$\mathcal{S}({\Sigma }_r)={\int }_{{\alpha }_1}^{{\beta }_1}\cdots {\int }_{{\alpha }_{n-1}}^{{\beta }_{n-1}}\mathrm{d}\mathcal{S}=r^{n-1}A(\alpha ,\beta )$$

$$A(\alpha ,\beta )={\int }_{{\alpha }_1}^{{\beta }_1}\cdots {\int }_{{\alpha }_{n-1}}^{{\beta }_{n-1}}{\sin }^{n-2}({\theta }_1){\sin }^{n-1}({\theta }_2)\cdots \sin ({\theta }_{n-2})\mathrm{d}{\theta }_{n-1}\cdots \mathrm{d}{\theta }_1$$

$$P({\Sigma }_r)=\Gamma \left(\frac{n}{2}\right)\frac{A(\alpha ,\beta )}{2{\pi }^{\frac{n}{2}}},$$

$$P({\Sigma }_r)={\int }_{{\alpha }_1}^{{\beta }_1}\cdots {\int }_{{\alpha }_{n-1}}^{{\beta }_{n-1}}{g}_{\Theta |R}\mathrm{d}\theta$$

$$\int_{\alpha_{1}}^{\beta_{1}}\!\cdots\!\int_{\alpha_{n-1}}^{\beta_{n-1}}\frac{g(r\cos(\theta_{1}),\ldots,r\sin(\theta_{1})\cdots\sin(\theta_{n-1}))\,r^{n-1}\sin^{n-2}(\theta_{1})\sin^{n-1}(\theta_{2})\cdots\sin(\theta_{n-2})\cdot\mbox{\Large$\chi$}}{g_{R}(r)}\,{\rm{{d}}}\theta$$

$$=\frac{\Gamma \left(\frac{n}{2}\right)}{2{\pi }^{\frac{n}{2}}}{\int }_{{\alpha }_1}^{{\beta }_1}\cdots {\int }_{{\alpha }_{n-1}}^{{\beta }_{n-1}}{\sin }^{n-2}({\theta }_1){\sin }^{n-1}({\theta }_2)\cdots \sin ({\theta }_{n-2})\mathrm{d}\theta ,$$

$$\frac{\Gamma \left(\frac{n}{2}\right)}{2{\pi }^{\frac{n}{2}}}=\frac{g(r\cos ({\theta }_1),\dots ,r\sin ({\theta }_1)\cdots \sin ({\theta }_{n-1}))r^{n-1}}{g_R(r)},$$

$$f(x_1)\cdots f(x_n)=h(x_1^2+\cdots +x_n^2)\text{mbox}forall(x_1,\dots ,x_n)\in {\mathbb{R}}^n$$

$$f(x_1)\dots f(x_n)=\frac{h(x_1^2)}{(f(0){)}^{n-1}}\cdots \frac{h(x_n^2)}{(f(0){)}^{n-1}}=\frac{h(x_1^2)}{(f(0){)}^n}\cdots \frac{h(x_n^2)}{(f(0){)}^n}(f(0){)}^n=h(x_1^2+\cdots +x_n^2).$$

$$\frac{h(x_1^2)}{h(0)}\cdots \frac{h(x_n^2)}{h(0)}=\frac{h(x_1^2+\cdots +x_n^2)}{h(0)}.$$

$$H(u_1)\cdots H(u_n)=H(u_1+\cdots +u_n)\text{mbox}withH(0)=1.$$

$$\phi (u_1)+\cdots +\phi (u_n)=\phi (u_1+\cdots +u_n)\text{mbox}with\phi (0)=0.$$

$$f(x)=\frac{1}{\sqrt{2\pi {\sigma }^2}}{\mathrm{e}}^{-\frac{x^2}{2{\sigma }^2}},$$

$$P(\mathcal{A})=\frac{1}{(\sqrt{2\pi }{)}^m}{\int }_{\mathcal{A}}{\mathrm{e}}^{-||\mathbf{a}|{|}^2/2}\mathrm{d}\mathbf{a},$$

$$Q^{\star }(z)=q_n(z+1{)}^n+q_{n-1}(z+1{)}^{n-1}(z-1)+\dots +q_0(z-1{)}^n.$$

$$E(Y)$$

$$\displaystyle=\frac{n}{2}\big{(}(p_{0}+p_{n})+(p_{1}+p_{n-1})+\cdots+(p_{\frac{n-1}{2}}+p_{\frac{n+1}{2}})\big{)}=\frac{n}{2}.$$

$$\mathrm{Var}(P_M)=\mathrm{Var}\left(\frac{1}{M}\sum _{j=1}^M{B}_j\right)=\frac{1}{M^2}M\mathrm{Var}(B_1)=\frac{p(1-p)}{M}\le \frac{1}{4M},$$

$$P\left(|X-E(X)|<\epsilon \right)\ge 1-\frac{\mathrm{Var}(X)}{{\epsilon }^2}.$$

$$P\left(|P_M-p|<\epsilon \right)\ge 1-\frac{p(1-p)}{M{\epsilon }^2}\ge 1-\frac{1}{4M{\epsilon }^2}.$$

Full text

Mathematics Subject Classification 2010: 37H10, 34F05, 39A25, 37C75.

Keywords: Stability index; random differential equations; random difference equations; random dynamical systems.

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 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 [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,

[TABLE]

where both $\mathbf{x},\mathbf{x}_{k}\in{\mathbb{R}}^{n}$ and $A$ is an $n\times 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 as

[TABLE]

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 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 searched 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 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, 3, 13, 18]. It is the case that using $M$ samples this approach gives the true value with an absolute error of order O\big{(}((\log\log M)/M)^{1/2}\big{)}, 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=10^{8},$ 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†††Concretely, we use the commands RandomVariable(Normal(0,1)) and Sample.. This procedure use the Mersenne Twister method with period $2^{19937}-1$ to generate uniformly-distributed pseudo-random numbers, 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=10^{8}$ is much smaller than the period of the pseudo-random number generator, which is greater that $10^{6001}$.

• •

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 the Monte Carlo’s observed relative frequencies 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 and quasi-Monte Carlo methods [2, 13].

To have a flavour of the type of 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 [math] solution to be a global stable attractor (stability index equals $n$) for them. Then, for differential equations, $r_{n}\leq 1/2^{n}.$ Furthermore, $r_{1}=1/2,$ $r_{2}=1/4,$ $r_{3}=1/16$ and our two step approach gives that $r_{4}\simeq 0.00925,$ $r_{5}\simeq 0.00071,$ and that $r_{k}$ is smaller that $10^{-4}$ for bigger $k.$ In the case of difference equations we prove that $r_{1}=1/2$ and $r_{2}=\frac{1}{\pi}\arctan(\sqrt{2})\simeq 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:

[TABLE]

where $A,B,C,D$ are random variables, so we can set the sample space to be $\Omega={\mathbb{R}}^{4}$. 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 $\mathbb{R}^{4}$. But there is no uniform distribution for unbounded probability spaces. Nevertheless, 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 proportional change in the time scale). This means that in the space of parameters, ${\mathbb{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

[TABLE]

has a uniform distribution on the sphere ${\mathbb{S}}^{3}\subset{\mathbb{R}}^{4}$. This achieves our objective, since ${\mathbb{S}}^{3}$ 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 the next theorem. We remark that the converse result is well–known [15, 19].

Theorem 1**.**

Let $X_{1},X_{2},\ldots,X_{n}$ be i.i.d. one-dimensional random variables with a continuous positive density function $f$. The random vector

[TABLE]

has a uniform distribution in $\mathbb{S}^{n-1}\subset\mathbb{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 suggests to use those distributions that maximize the entropy, i.e. the quantity $h(f)=-\int_{\Omega}f(x)\ln(f(x))\mathrm{d}x$ for any given density $f$. The one-dimensional random variables with continuous probability density function $f$ on $\Omega={\mathbb{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 $(\Omega,\mathcal{F},P)$ where $\Omega={\mathbb{R}}^{4}$, $\mathcal{F}$ is the $\sigma$-algebra generated by the open sets of ${\mathbb{R}}^{4}$ and $P:\mathcal{F}\rightarrow[0,1]$ is the probability function with density $\frac{1}{4\pi^{2}}\mathrm{e}^{-(a^{2}+b^{2}+c^{2}+d^{2})/{2}},$ where for simplicity we take variance 1 in each marginal density function.

For instance, assume that we want to compute the probability $\alpha$ 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 possibility 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 $\alpha$ coincides with the probability of having a saddle (stability index $1$) at the origin, i.e. $AD-BC<0$. Then, the open set $\mathcal{U}:=\{(a,b,c,d)\in{\mathbb{R}}^{4}:ad-bc<0\}$ belongs to $\mathcal{F}$ and

[TABLE]

which is $1/2$ by symmetry, as we will see.

Proof of Theorem 1.

Let $(X_{1},\ldots,X_{n})$ be the random vector associated with the random variables of the statement, with joint continuous density function $g(x_{1},\ldots,x_{n})$. We claim that

[TABLE]

for some continuous function $h$.

Taking spherical coordinates, we consider the new random vector $(R,\Theta)\in\mathbb{R}^{n}$ where $R=(X_{1}^{2}+\cdots+X_{n}^{2})^{1/2}$ and $\Theta=(\Theta_{1},\ldots,\Theta_{n-1}).$ We have $X_{1}=R\cos\Theta_{1}$, $X_{2}=R\sin\Theta_{1}\cos\Theta_{2},\ldots$ $X_{n-1}=R\sin\Theta_{1}\sin\Theta_{2}\cdots\sin\Theta_{n-2}\cos\Theta_{n-1}$ and $X_{n}=R\sin\Theta_{1}\sin\Theta_{2}\cdots\sin\Theta_{n-2}\sin\Theta_{n-1}$. By the change of variables theorem, the joint density function of $(R,\Theta)$ is

[TABLE]

where $\theta=(\theta_{1},\ldots,\theta_{n-1})$, and

[TABLE]

where \mbox{\Large\chi}_{A} stands for the characteristic function of the set $A.$

The density function of $(R,\Theta)$ conditioned to $R$, $g_{\Theta|R},$ is

[TABLE]

where $g_{R}(r)$ is the marginal density of $R$:

[TABLE]

where ${\rm{{d}}}\mathcal{S}=r^{n-1}\sin^{n-2}(\theta_{1})\sin^{n-1}(\theta_{2})\cdots\sin(\theta_{n-2}){\rm{{d}}}\theta_{n-1}\cdots{\rm{{d}}}\theta_{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(x_{1},\ldots,x_{n})$ such that when we fix $R=r$, the probability on the $(n-1)$-dimensional sphere of radius $r$, denoted by $\mathbb{S}^{n-1}(r)$, is uniformly distributed. In such a case the partial spherical segment $\Sigma_{r}=\{R=r,\,\theta_{i}\in[\alpha_{i},\beta_{i}]$ for $i=1,\ldots,n-1\}$ must have probability $P(\Sigma_{r})={\mathcal{S}(\Sigma_{r})}/{\mathcal{S}(\mathbb{S}^{n-1}(r))}$ where $\mathcal{S}$ denotes the surface area. Set $\alpha=(\alpha_{1},\ldots,\alpha_{n-1})$ and $\beta=(\beta_{1},\ldots,\beta_{n-1})$. Notice that

[TABLE]

where

[TABLE]

and $\mathcal{S}(\mathbb{S}^{n-1}(r))=\frac{2\pi^{\frac{n}{2}}}{\Gamma\left(\frac{n}{2}\right)}\,r^{n-1}.$ Hence, on the one hand,

[TABLE]

which does not depend on $r$. On the other hand,

[TABLE]

where ${\rm{{d}}}\theta={\rm{{d}}}\theta_{n-1}\cdots{\rm{{d}}}\theta_{2}{\rm{{d}}}\theta_{1}$. This implies that

[TABLE]

[TABLE]

​​for all $\alpha_{i},\beta_{i}\in[0,\pi)$ for $i=1,\ldots,n-2$ with $\alpha_{i}<\beta_{i}$ and $\alpha_{n-2},\beta_{n-2}\in[0,2\pi)$ with $\alpha_{n-2}<\beta_{n-2}$. This last equality implies that almost everywhere

[TABLE]

and therefore $g(r\cos(\theta_{1}),\ldots,r\sin(\theta_{1})\cdots\sin(\theta_{n-1}))$ is a function that only depends on $r.$ In consequence, writing this fact in Cartesian coordinates, we get that almost everywhere $g(x_{1},\ldots,x_{n})=h(x_{1}^{2}+\cdots+x_{n}^{2}),$ for some continuous function $h$ and the claim (4) follows.

Now we complete the proof. Since $X_{1},\ldots,X_{n}$ are i.i.d. with positive density $f$, we know that $g(x_{1},\cdots,x_{n})=f(x_{1})\ldots f(x_{n}).$ So equation (4) can be expressed as

[TABLE]

where $h$ is a positive function. Taking $x_{2}=\cdots=x_{n}=0$ we have that $f(x_{1})\,f(0)^{n-1}=h(x_{1}^{2})$ and $h(0)=(f(0))^{n}>0.$ Thus,

[TABLE]

Hence, using that $h(0)=(f(0))^{n}>0,$

[TABLE]

Taking $H(\xi):=h(\xi)/h(0),$ and $u_{i}=x_{i}^{2}$, it holds that

[TABLE]

Hence, $\varphi(U)=\log(H(u))$ is a continuous function that satisfies the Cauchy’s functional equation

[TABLE]

It is well–known that all its continuous solutions are $\varphi(x)=ax,$ for some $a\in{\mathbb{R}}.$ Hence all continuous solutions of (5) are $H(x)=\mathrm{e}^{ax}.$

As a consequence, $f(x)=b\,\mathrm{e}^{ax^{2}}$ for some $(a,b)\in{\mathbb{R}}^{2}.$ Since $f$ is a density function, $a<0.$ Moreover, using $\int_{-\infty}^{\infty}b\mathrm{e}^{ax^{2}}{\rm{{d}}}x=b\sqrt{-\pi/a}=1$, and setting $a=-1/(2\sigma^{2})$, we get that

[TABLE]

so each variable $X_{i}$ is a normal random variable $\mathrm{N}(0,\sigma^{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 the 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 $\dot{\mathbf{x}}=A\,\mathbf{x}$ where $\mathbf{x}\in{\mathbb{R}}^{n}$ and $A$ is a real constant $n\times n$ matrix,

• (b)

Homogeneous linear differential equation of order $n$ with constant coefficients: $a_{n}x^{(n)}+a_{n-1}x^{(n-1)}+\cdots+a_{1}x^{\prime}+a_{0}x=0,$

• (c)

Linear discrete systems $b\,\mathbf{x}_{k+1}=A\,\mathbf{x}_{k}$ where $\mathbf{x}_{k}\in{\mathbb{R}}^{n},$ $b\in{\mathbb{R}}$; and $A$ is a real constant $n\times n$ matrix,

• (d)

Linear homogeneous difference equation of order $n$ with constant coefficients $a_{n}x_{k+n}+a_{n-1}x_{k+n-1}+\cdots+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, 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 $(\Omega,\mathcal{F},P)$, where $\Omega={\mathbb{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, $\mathcal{F}$ is the $\sigma$-algebra generated by the open sets and for each $\mathcal{A}\in\mathcal{F}$,

[TABLE]

where $\mathbf{a}=(a_{1},a_{2},\ldots,a_{m}),$ $||\mathbf{a}||^{2}=\sum_{j=1}^{m}a_{j}^{2}$ and ${\rm{{d}}}\mathbf{a}={\rm{{d}}}a_{1}\,{\rm{{d}}}a_{2}\ldots{\rm{{d}}}a_{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(\lambda)=q_{n}\lambda^{n}+q_{n-1}\lambda^{n-1}+\cdots+q_{1}\lambda+q_{0}$ with $q_{j}\in{\mathbb{C}},$ by using the conformal transformation $\lambda=\frac{z+1}{z-1},$ we get the associated polynomial

[TABLE]

It is straightforward to observe that $\lambda_{0}\in\mathbb{C}$ is a root of of $Q(\lambda)$ such that $|\lambda_{0}|<1$ if and only if $z_{0}=(\lambda_{0}+1)/(\lambda_{0}-1)$ is a root of $Q^{\star}(z)$ such that $\mathrm{Re}(z_{0})<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\leq k\leq n,$ is measurable. Hence $\mathcal{A}_{k}:=\{\mathbf{a}\in{\mathbb{R}}^{m}\,:\,X(\mathbf{a})=k\}\in\mathcal{F}$ and its probability $p_{k}:=P(\mathcal{A}_{k})$ 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 $p_{k}$ for $k\leq 10.$

3.1 A preliminary result

In three of the above considered cases we will apply the following auxiliary result:

Lemma 2**.**

Let $(\Omega,\mathcal{F},P)$ be a probability space and let $Y:\Omega\to\mathbb{R}$ be a discrete random variable with image $\mathrm{Im}(Y)=\{0,1,\ldots,n\}$, and probability mass function $p_{k}=P(Y=k)$ such that $p_{k}=p_{n-k}$ for all $k=0,\ldots,n$. Then $E(Y)=\sum_{k=0}^{n}kp_{k}=n/2.$ Moreover

1. (a)

If $n$ is odd then $2\sum_{k=0}^{\frac{n-1}{2}}p_{k}=1.$ In particular, when $n=1,$ $p_{0}=p_{1}=\frac{1}{2}.$ 2. (b)

If $n$ is even and $n\geq 2$ then $2\sum_{k=0}^{\frac{n}{2}-1}p_{k}\,+\,p_{\frac{n}{2}}=1.$

If, additionally, $n$ is even‡‡‡When $n$ is odd the imposed equalities automatically hold.* and $\sum_{i\,\mathrm{odd}}p_{i}=\sum_{i\,\mathrm{even}}p_{i}=\frac{1}{2}$, then*

1. (c)

If $\frac{n}{2}$ is even, then $2\sum_{k=0,\,k\,\mathrm{even}}^{\frac{n}{2}-2}p_{k}+p_{\frac{n}{2}}=\frac{1}{2},$ and $2\sum_{k=1,\,k\,\mathrm{odd}}^{\frac{n}{2}-1}p_{k}=\frac{1}{2}.$ In particular, when $n=4,$ $p_{1}=p_{3}=\frac{1}{4},$ $p_{2}=\frac{1}{2}-2p_{0}$ and $p_{4}=p_{0}.$ 2. (d)

If $\frac{n}{2}$ is odd, then $2\sum_{k=0,\,k\,\mathrm{even}}^{\frac{n}{2}-1}p_{k}=\frac{1}{2},$ and $\sum_{k=1,\,k\,\mathrm{odd}}^{\frac{n}{2}-2}p_{k}+p_{\frac{n}{2}}=\frac{1}{2}.$ In particular, when $n=2,$ $p_{0}=p_{2}=\frac{1}{4}$ and $p_{1}=\frac{1}{2}.$

Proof.

We start proving that $E(Y)=n/2.$ Assume for instance that $n$ is odd. Since $p_{k}=p_{n-k},$ its holds that $kp_{k}+(n-k)p_{n-k}=np_{k},$ for each $k\leq(n-1)/2.$ Hence,

[TABLE]

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 probabilities $p_{k}$ 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=10^{8}$ to obtain an estimation, say $\widetilde{p},$ for a probability $p:=P(\mathcal{A})$ like the one given in equality (6) for different measurable sets $\mathcal{A}.$ Further details for each concrete situation are given in each of the following subsections.

In a few words, recall that $\widetilde{p}$ is given by the proportion of samples that are in $\mathcal{A}$. For studying, for a given $M,$ how close are $p$ and $\widetilde{p},$ let $B_{j},j=1,\ldots,M,$ i.i.d. Bernoulli random variables, where each one of them that takes the value $1$ with probability $p$ and the value [math] with probability $1-p.$

Define $P_{M}=\frac{1}{M}\sum_{j=1}^{M}B_{j}.$ Then, the value obtained for the random variable $P_{M},$ $\widetilde{p}$ is the approximation of $p$ given by Monte Carlo method. Let us see, by using Chebyshev inequality or the Central limit theorem, that with very high probability, $\widetilde{p}$ is a good approximation of $p.$

Notice first that $\operatorname{E}(P_{M})=p$ and due to the independence of the $B_{j},$

[TABLE]

because $p(1-p)\leq 1/4.$ Recall also that for each $\varepsilon>0$ and any random variable $X,$ with $\operatorname{E}(X^{2})<\infty,$ the Chebyshev inequality reads as

[TABLE]

Hence, applying the Chebyshev inequality to $X=P_{M}$ we get that

[TABLE]

Taking $M=10^{8},$ as in our computations, denoting $\widetilde{p}=P_{10^{8}},$ and considering $\varepsilon=10^{-3}$ we get that

[TABLE]

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

[TABLE]

can be practically considered to be a random variable $Z$ with distribution $\mathrm{N}(0,1).$ In fact, in Statistics it is usually imposed that $p(1-p)M>18.$ Hence, unless $p$ is very close to [math] or $1$, the value $M=10^{8}$ is big enough. Hence

[TABLE]

where $\Phi$ is the distribution function of a $\mathrm{N}(0,1).$ Taking again $M=10^{8}$ and $\varepsilon=10^{-3}$ or $\varepsilon=2\times 10^{-4}$ we get

[TABLE]

In fact, for instance looking at the values $p_{k}$ of Table 2 for $n=2$ in Section 4, that can also be obtained analytically, we get that $|\widetilde{p}_{k}-p_{k}|\leq 6\times 10^{-5},$ for $k=0,1,2.$ So, the actual bound is smaller that the bounds obtained above.

Finally, and 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 second order difference equation $A_{2}x_{k+2}+A_{1}x_{k+1}+A_{0}x_{k}=0$ where $A_{i}$ are i.i.d. random variables with $\mathrm{N}(0,1)$ distribution. The stability index is given by the number of zeroes with modulus smaller than $1$ of the characteristic polynomial $Q(\lambda)=A_{2}\lambda^{2}+A_{1}\lambda+A_{0}$. Let $X$ be the random variable that counts the number of roots with modulus smaller than $1$ of $Q(\lambda)$, and $p_{k}=P(X=k)$ for $k=0,1,2$. The true value of the probabilities $p_{k}$ is obtained in Corollary 13. Performing Monte Carlo simulations with $M=10^{m}$ with $m=2,\ldots,10$ we obtain the observed frequencies $\widetilde{p}_{2}(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 $p_{0}=p_{2}=\arctan(\sqrt{2})/\pi$ and, of course, $p_{1}=1-p_{0}-p_{2}={2}\arctan(1/\sqrt{2})/\pi$). For $M=10^{m}$ we denote the absolute error $e_{m}=|\widetilde{p}_{2}(m)-p_{2}|$:

[TABLE]

Table 1. Observed frequency and absolute error of $p_{2}$ for second order difference equations, using that $p_{2}=\arctan(\sqrt{2})/\pi\approx 0.304086723985$.

With the above results, the regression line of $Y=\log(e_{m})$ versus $X=\log(M)=m\log(10)$ is $Y=-0.505\,X-1.260$ with $R^{2}=0.893$. The slope is therefore $-0.505\approx-1/2$ as was expected a priori since, in practice, the absolute error behaves as $O(M^{-1/2})$, 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 $p_{k}$ satisfy some affine relations, like the ones in Lemma 2 or the ones in Proposition 12 below. Then, if we denote $\mathbf{p}=(p_{0},\ldots,p_{n})^{t}\in{\mathbb{R}}^{n+1}$ it is possible to write $\mathbf{p}=\mathrm{M}\mathbf{q}+b$ where $\mathbf{q}\in{\mathbb{R}}^{k}$ with $k\leq n$ is a vector whose components are different elements of ${p_{0},\ldots,p_{n}}$; $\mathrm{M}\in\mathcal{M}_{n\times k}({\mathbb{R}})$; and $b\in{\mathbb{R}}^{k}$. Let $\widetilde{\mathbf{p}}=(\widetilde{p}_{0},\ldots,\widetilde{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 [20, Def. 6.1] of the the system,

[TABLE]

which is

[TABLE]

see [9, p. 198] or [22, p. 200]. So we can find some improved estimations $\widehat{\mathbf{p}}$, via

[TABLE]

Some detailed examples are given in Sections 4, 5 and 7.

4 Linear random differential systems

Consider linear differential systems $\dot{\mathbf{x}}=A\,\mathbf{x}\mbox{ where }\mathbf{x}\in{\mathbb{R}}^{n},\,A\in\mathcal{M}_{n\times n}({\mathbb{R}}),$ where $A$ is a random matrix whose entries are i.i.d. random variables with $\mathrm{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 $p_{k}=P(X=k).$ The following holds:

1. (a)

$\sum_{k=0}^{n}p_{k}=1.$** 2. (b)

For all $k\in\{0,1,\ldots,n\},$ $p_{k}=p_{n-k}.$ 3. (c)

$\sum_{i\,\mathrm{odd}}p_{i}=\sum_{i\,\mathrm{even}}p_{i}=\frac{1}{2}.$**

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 $q_{m}$ the probability that $B$ has $m$ eigenvalues with negative real part, we get that $p_{m}=q_{m}.$ This is so, because if $X\sim\mathrm{N}(0,1)$ then $-X\sim\mathrm{N}(0,1)$ and as a consequence the entries of $A$ and $B$ have the same distribution. Then, $q_{k}=p_{n-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 and, moreover, $P(\det(A)>0)=1/2.$ From this claim we get the result because $\sum_{i\,\text{even}}\,p_{i}$ is the probability of $s(A)$ being even. To prove the first part of the claim notice first that we can assume that $0\neq\det(A)=\lambda_{1}\lambda_{2}\cdots\lambda_{n},$ where $\lambda_{1},\lambda_{2},\ldots,\lambda_{n}$ are the $n$ eigenvalues of $A.$ We write $\lambda_{1}\lambda_{2}\cdots\lambda_{n}=(\lambda_{1}\lambda_{2}\cdots\lambda_{k})(\lambda_{k+1}\lambda_{k+2}\cdots\lambda_{n})$ where $\lambda_{1},\lambda_{2},\ldots,\lambda_{k}$ are all the real negative eigenvalues. Observe also that for complex eigenvalues $\lambda\bar{\lambda}>0.$ Hence $\lambda_{k+1}\lambda_{k+2}\cdots\lambda_{n}>0,$ ${\rm sign}(\det(A))=(-1)^{k}$ and the condition that $s(A)$ is even is characterized 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)\cdot\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.$

From the above proposition it easily follows:

Corollary 4**.**

Consider $\dot{\mathbf{x}}=A\,\mathbf{x},$ $\mathbf{x}\in{\mathbb{R}}^{n}$ with $A\in\mathcal{M}_{n\times n}({\mathbb{R}})$ a random matrix with i.i.d. $\mathrm{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 generated $10^{8}$ matrices $A\in\mathcal{M}_{n\times n}({\mathbb{R}})$ whose entries are pseudo-random numbers simulating the realizations on $n^{2}$ independent random variables with $\mathrm{N}(0,1)$ distribution. For each matrix $A$ we have computed the characteristic polynomial, and counted the number of eigenvalues 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 polynomial 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\leq 4$. For $n\geq 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 $p_{3}=p_{4}=\frac{1}{2}-p_{0}-p_{1}-p_{2}$; $p_{5}=p_{2}$; $p_{6}=p_{1}$ and $p_{7}=p_{0}$. So, using the notation introduced in Section 3.2, we can write $\mathbf{p}=\mathrm{M}\mathbf{q}+b$, where $\mathbf{p}^{t}=(p_{0},\ldots,p_{7})$;

[TABLE]

The observed relative frequencies in our Monte Carlo simulation are

[TABLE]

By finding the least squares solution of the system (8) ([9, p. 198] or [22, p. 200]), given by (10), we obtain

[TABLE]

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.

[TABLE]

[TABLE]

[TABLE]

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$

[TABLE]

where $x=x(t)$, the derivatives are taken in respect to $t$, and $A_{j}$ are again i.i.d. random variables with $\mathrm{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:

[TABLE]

Let $X$ be the random variable that counts the number of roots of $Q(\lambda)$ with negative real parts and define $p_{k}=P(X=k)$ for $k=0,1,\ldots,n.$

Proposition 5**.**

Set $p_{k}=P(X=k),$ where $X$ is the random variable defined above. Then

1. (a)

$\sum_{k=0}^{n}p_{k}=1.$** 2. (b)

For all $k\in\{0,1,\ldots,n\},$ $p_{k}=p_{n-k}.$ 3. (c)

$\sum_{i\,\mathrm{odd}}p_{i}=\sum_{i\,\mathrm{even}}p_{i}=\frac{1}{2}.$**

Proof.

The proof of (a) is trivial. To prove (b) consider equation (11) with its characteristic polynomial $Q(\lambda)$ and also the new differential equation

[TABLE]

with its characteristic polynomial $Q^{*}(\lambda)=Q(-\lambda)=(-1)^{n}A_{n}\lambda^{n}+(-1)^{n-1}A_{n-1}\lambda^{n-1}+\cdots-A_{1}\lambda+A_{0}.$ Since $Q(\lambda)=0$ if and only if $Q^{*}(-\lambda)=0$ we get that $p_{k}=p^{*}_{n-k}$ where $p^{*}_{i}$ the probability that $Q^{*}(\lambda)$ has $i$ roots with negative real part. But also $p_{k}=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(\lambda)$ has an odd number of roots with negative real part if and only if $A_{0}\cdot A_{n}<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(\lambda)$ has an odd number of roots with negative real part is

[TABLE]

Corollary 6**.**

Consider the linear random homogeneous differential equation of order $n$ (11), with all $A_{i}$ being i.i.d. $\mathrm{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 conclusions 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 (11), that is $r_{n}=p_{n}.$ By Proposition 5(b) this probability coincides with the probability 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.

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;UT-SV>0)=1/32$.

Proof.

Set $\mathcal{A}^{\pm}=\{U>0;V>0;S>0;T>0;\pm(UT-SV)>0\},$ and $\mathcal{A}^{0}=\{U>0;V>0;S>0;T>0;UT-SV=0\}.$ Denote by $p^{\pm}=P(\mathcal{A}^{\pm})$ and $p^{0}=P(\mathcal{A}^{0}).$ Then, since $p^{0}=0$ and $\mathcal{A^{-}}\cup\mathcal{A}^{0}\cup\mathcal{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

[TABLE]

This is so, because for instance in the definition of $\mathcal{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 $\mathcal{A}^{+}$ and $\mathcal{A}^{-}$. Since all variables are independent $\mathrm{N}(0,1)$, both sets have the same probability and $p^{+}=p^{-},$ as we wanted to prove.

Proposition 8**.**

With the above notations, $r_{n}\leq{1}/{2^{n}},$ for all $n\geq 1.$ Moreover, $r_{1}=1/2,$ $r_{2}=1/4,$ $r_{3}=1/16$ and $r_{4}<r_{3}/2=1/32.$

Proof of Proposition 8.

Notice that $r_{n}$ is the probability that the characteristic polynomial $Q(\lambda)$, associated with the random differential equation (11), is a Hurwitz stable polynomial; that is $r_{n}=P(\mbox{Every root of }Q(\lambda)\mbox{ belongs to }\mathfrak{R}^{-}),$ where $\mathfrak{R}^{-}=\{z\in\mathbb{C}\mbox{ such that }\mathrm{Re}(z)<0\}$. It is well–known that a necessary condition for a polynomial to have every root in $\mathfrak{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,

[TABLE]

Since the variables $A_{i}$ are independent and symmetric

[TABLE]

As a consequence,

[TABLE]

and the first statement follows.

The equalities $r_{1}=1/2$ and $r_{2}=1/4$ are a simple consequence that for $n=1,2$ the inclusion (13) is an equality.

Let us prove that $r_{3}=p_{3}=1/16.$ By using the Routh-Hurwitz criterion [10, p. 1076], it can be seen that $a_{3}\lambda^{3}+a_{2}\lambda^{2}+a_{1}\lambda+a_{0}$ has every root in $\mathfrak{R}^{-}$ if and only if all its coefficients have the same sign and moreover $a_{1}a_{2}-a_{0}a_{3}>0$. Hence, $p_{3}=p_{3}^{-}+p_{3}^{+},$ where $p_{3}^{-}:=P(A_{0}<0;A_{1}<0;A_{2}<0;A_{3}<0;A_{1}A_{2}-A_{0}A_{3}>0);$ and $p_{3}^{+}:=P(A_{0}>0;A_{1}>0;A_{2}>0;A_{3}>0;A_{1}A_{2}-A_{0}A_{3}>0),$ with all the $A_{i}$ being $\mathrm{N}(0,1)$ distributed and independent. Due to their symmetry, the random variables $A_{i}$ and $-A_{i}$, for $i=0,\ldots,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}+\cdots+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-Hurwitz conditions to have stability index $n$ for $n=3,4$ are precisely that the principal minors of the following matrices

[TABLE]

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 $a_{0}>0$ and for $n=3$ are: $a_{1}>0,a_{1}a_{2}-a_{0}a_{3}>0$ and $a_{3}>0.$ Similarly, for $n=4$ the conditions are $a_{1}>0,a_{1}a_{2}-a_{0}a_{3}>0,$ $a_{3}(a_{1}a_{2}-a_{0}a_{3})-a_{4}a_{1}^{2}>0$ and $a_{4}>0.$

Consider now, for $n=3,4$, the random polynomials $Q_{n}(x)=\tilde{A}_{0}x^{n}+\tilde{A}_{1}x^{n-1}+\cdots+\tilde{A}_{n-1}x+\tilde{A}_{n},$ where $\tilde{A}_{i}\sim\mathrm{N}(0,1)$ and are independent (notice that with this notation each coefficient $\tilde{A}_{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(\mathcal{A}^{+}_{k}),$ with

[TABLE]

Notice that if we define

[TABLE]

it is clear that $P(\mathcal{B})=p_{3}^{+}/2$ and, moreover $\mathcal{A}^{+}_{4}\subset\mathcal{B}_{3},$ with the inclusion being strict. Since the joint density is positive and $\mathcal{B}\cap(\mathcal{A}^{+}_{4})^{c}$ has positive Lebesgue measure, we have $P(\mathcal{B}\cap(\mathcal{A}^{+}_{4})^{c})>0.$ Thus $P(\mathcal{A}^{+}_{4})<P(\mathcal{B})$, and hence $p_{4}^{+}=P(\mathcal{A}^{+}_{4})<P(\mathcal{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, $p_{1}=p_{2}=7/16.$

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 $n=8$ and $10,$ where we have encountered that the vectors $\widehat{\mathbf{p}}$ have negative and very small entries. This has occurred because the observed frequencies obtained by the Monte Carlo approach corresponding to these probabilities are not enough accurate. For this reason, we have made a new optimization step. As before, we use the least squares method to obtain a vector $\widehat{\mathbf{p}}$. However, if negative entries appear in this vector (which is clearly objectionable), 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

[TABLE]

The relations stated in Corollary 6 are $p_{3}=1/4-p_{1}$, $p_{4}=1/2-2p_{0}-2p_{2}$, $p_{5}=p_{3}$, $p_{6}=p_{2}$, $p_{7}=p_{1}$, $p_{8}=p_{0}$. By solving the system (8) with

[TABLE]

we obtain

[TABLE]

Hence, by (10) we get

[TABLE]

So we impose that $p_{0}=p_{8}=0$. Thus we have $p_{3}=p_{5}=1/4-p_{1}$, $p_{4}=1/2-2p_{0}-2p_{2}=1/2-2p_{2}$, $p_{6}=p_{2}$ and $p_{7}=p_{1}$. We find the least squares solution of the system

[TABLE]

Using (9) and (10) we obtain

[TABLE]

[TABLE]

The $n=10$ case follows analogously.

[TABLE]

[TABLE]

[TABLE]

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

[TABLE]

where $\mathcal{B}$ and each of the $n^{2}$ entries of the random matrix $A$ are i.i.d. $\mathrm{N}(0,1)$ random variables. Observe that to ensure that the results are invariant under time-scaling, is necessary to add the term $\mathcal{B}$ in the left-hand side of Equation (14). Then, given a linear discrete random system (14), its characteristic random polynomial associated with the matrix $\frac{1}{\mathcal{B}}A$ is

[TABLE]

where each random variable $Q_{j}$ is a polynomial in the variables $1/\mathcal{B},A_{1,1},\ldots,A_{n,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 the matrix $\frac{1}{\mathcal{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 longer 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}$ apart from the trivial one $\sum_{k=0}^{n}p_{k}=1$. Since this 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 solution of $Q(\lambda)=0$ is $\lambda=A/B,$ with $A$ and $B$ independent and $\mathrm{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\leq 10,$ we generate $10^{8}$ discrete systems of the form (14). For each matrix $\frac{1}{\mathcal{B}}A$ we have computed the characteristic polynomial $Q$ and its associated polynomial $Q^{\star}$ (see Equation (7)) and have counted the number of roots of this last polynomial by using the Routh-Hurwitz zero counter. For $n\geq 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.

[TABLE]

[TABLE]

[TABLE]

[TABLE]

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

[TABLE]

where all the coefficients are i.i.d. random variables with $\mathrm{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(\lambda)=A_{n}\lambda^{n}+A_{n-1}\lambda^{n-1}+\cdots+A_{1}\lambda+A_{0}.$ As in the preceding sections let $X$ be the random variable that counts the number of roots of $Q(\lambda)$ with modulus smaller than 1 and set $p_{k}=P(X=k)$ for $k=0,1,\ldots,n.$

Before proving some relations among the probabilities $p_{k},$ we give two preliminary lemmas. Let $\operatorname{erf}(x)=\frac{2}{\sqrt{\pi}}\int_{0}^{x}{\rm e}^{-u^{2}}{\rm d}u$ be the error function. The following result is stated in [4]. We prove it for the sake of completeness.

Lemma 10**.**

For $\alpha>0$ and $\beta\in{\mathbb{R}},$

[TABLE]

Proof.

Fixed $\alpha>0,$ the function that defines $F$ is absolutely integrable because $|\operatorname{erf}(x)|\leq 1.$ Moreover its partial derivative with respect to $\beta$ is also absolutely integrable. Hence $\lim_{\beta\to 0}F(\alpha,\beta)=F(\alpha,0)=0$ and

[TABLE]

Therefore

[TABLE]

as we wanted to prove.

Next result is a consequence of the previous lemma.

Lemma 11**.**

Let $U\sim\mathrm{N}(0,\sigma^{2})$ and $V\sim\mathrm{N}(0,\rho^{2})$ be independent normal random variables. Then $P(U^{2}-V^{2}>0)=\frac{2}{\pi}\arctan(\sigma/\rho).$

Proof.

The joint density function of the random vector $(U,V)$ is $f_{\sigma}(u)f_{\rho}(v),$ where $f_{s}(u)={\rm e}^{-u^{2}/(2s^{2})}/(\sqrt{2\pi}s)$. Observe that the points $(u,v)\in\mathbb{R}^{2}$ such that $u^{2}-v^{2}>0$ is the region where $-|u|<v<|u|$, hence by symmetry,

[TABLE]

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

[TABLE]

where we have used the fact that $\arctan(x)+\arctan(1/x)=\pi/2$ or, simply, the same lemma interchanging $U$ and $V.$ Observe also that when $\sigma=\rho,$ $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:

1. (a)

$\sum_{k=0}^{n}p_{k}=1.$** 2. (b)

For all $k\in\{0,1,\ldots,n\},$ $p_{k}=p_{n-k}.$ 3. (c)

When $n$ is odd, $\sum_{i\,\,\mathrm{even}}p_{i}=\sum_{i\,\,\mathrm{odd}}p_{i}=\frac{1}{2}.$ 4. (d)

When $n=2k$ is even,

[TABLE]

Proof.

The first assertion is obvious. To see the second one we compare the difference equation $a_{n}x_{k+n}+a_{n-1}x_{k+n-1}+\cdots+a_{1}x_{k+1}+a_{0}x_{k}=0,$ with $a_{i}\in{\mathbb{R}},\,i=0,1,\ldots,n,$ with characteristic polynomial $Q(\lambda)=a_{n}\lambda^{n}+a_{n-1}\lambda^{n-1}+\cdots+a_{2}\lambda^{2}+a_{1}\lambda+a_{0}$ with the difference equation $a_{n}x_{k}+a_{n-1}x_{k+1}+\cdots+a_{0}x_{k+n}=0$ with characteristic polynomial ${Q}^{*}(\lambda)=a_{n}+a_{n-1}\lambda+\cdots+a_{1}\lambda^{n-1}+a_{0}\lambda^{n}.$ Notice that if $Q(\lambda)$ has $m$ non-zero roots with modulus smaller than 1 and $n-m$ with modulus bigger than 1, then the converse follows for $Q^{*}(\lambda)$ because $Q(\lambda)=0$ if and only if $Q^{*}(\frac{1}{\lambda})=0.$ From this result applied to the corresponding random polynomials we get that $p_{k}=p_{n-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(\lambda)=a_{n}\lambda^{n}+a_{n-1}\lambda^{n-1}+\cdots+a_{2}\lambda^{2}+a_{1}\lambda+a_{0}$, without roots with zero real part, has an even number of roots with negative real part if and only if $a_{n}a_{0}>0.$ By using the polynomial

[TABLE]

introduced in Section 3 (Equation (7)) we get that $Q(\lambda),$ without roots of modulus 1, has an even number ($2m$) of roots with modulus smaller than 1 if and only if $Q^{\star}(z)$ has exactly $2m$ roots with negative real part and this happens if and only if $(a_{n}+a_{n-1}+\cdots+a_{1}+a_{0})\cdot(a_{n}-a_{n-1}+a_{n-2}-\cdots+(-1)^{n}a_{0})>0.$ Hence, considering the corresponding random polynomials, we have that

[TABLE]

where $U=A_{n}+A_{n-2}+A_{n-4}+\cdots$ and $V=A_{n-1}+A_{n-3}+A_{n-5}+\cdots$ and the sums end either at $A_{0}$ or $A_{1}$ according the parity of $n.$ Since $A_{j}\sim\mathrm{N}(0,1)$ and all $A_{j}$ are independent we get that when $n=2k$ (resp. $n=2k-1$) then $U\sim\mathrm{N}(0,k+1)$ (resp. $U\sim\mathrm{N}(0.k)$) and $V\sim\mathrm{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,$

[TABLE]

The sum of all $p_{i}$ 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 $p_{k}=P(X=k)$. Then the probabilities $p_{k}$ 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,$ $p_{0}=p_{2}=\frac{1}{\pi}\arctan(\sqrt{2})$ and $p_{1}=\frac{2}{\pi}\arctan(\frac{1}{\sqrt{2}})$; and for $n=4,$ $p_{1}=p_{3}=\frac{1}{\pi}\arctan(\sqrt{\frac{2}{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,\ldots,10,$ $10^{8}$ random vectors $(A_{0},\ldots,A_{n})\in{\mathbb{R}}^{n+1}$ whose components are pseudo-random numbers with $\mathrm{N}(0,1)$ distribution. Then, by using the relations in Proposition 12 and Corollary 13 we have performed a least squares refinement.

For instance for $n=4$, by Corollary 13 we have $p_{1}=p_{3}=\arctan(\sqrt{2/3})/\pi\simeq 0.217953$; $p_{2}=2\arctan(\sqrt{3/2})/\pi-2p_{0}$ and $p_{4}=p_{0}$. Hence, we fix the values $\widehat{p}_{1}=p_{1}$ and $\widehat{p}_{3}=p_{3}$ and system (8) can be written in the form

[TABLE]

Hence we can easily find the least squares solution of the above incompatible linear system:

[TABLE]

and thus we get the result in Equation (9): $\widehat{p}_{0}=\frac{1}{6}\left(\widetilde{p}_{0}-2\widetilde{p}_{2}+\widetilde{p}_{4}\right)+\frac{2}{3\pi}\arctan\left(\sqrt{{3}/{2}}\right),$ and therefore $\widehat{p}_{4}=\widehat{p}_{0}$ and $\widehat{p}_{2}=2\arctan(\sqrt{3/2})/\pi-2\widehat{p}_{0}.$ Since our Monte Carlo simulations give

[TABLE]

the above relations show that $(\widehat{p}_{0},\widehat{p}_{2},\widehat{p}_{4})\simeq\left(0.10282,0.35846,0.10282\right).$

All our results are collected in Table 5.

[TABLE]

[TABLE]

[TABLE]

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 us very interesting suggestions that have allowed us to improve our work.