This is a preprint of: “Probability of occurrence of some planar random quasi-homogeneous vector fields”, Bartomeu Coll, Armengol Gasull, Rafel Prohens, Mediterranean J. Math., vol. 19, Article number 278, 2022. DOI: [10.1007/s00009-022-02198-w] PROBABILITY OF OCCURRENCE OF SOME PLANAR RANDOM QUASI-HOMOGENEOUS VECTOR FIELDS B. COLL† , A. GASULL‡ AND R. PROHENS∗ Abstract. The objective of this work is the study of the proba- bility of occurrence of phase portraits in a family of planar quasi- homogeneous vector fields of quasi degree q, that is a natural exten- sion of planar linear vector fields, which correspond to q = 1. We obtain the exact values of the corresponding probabilities in terms of a simple one-variable definite integral that only depends on q. This integral is explicitly computable in the linear case, recovering known results, and it can be expressed in terms of either complete elliptic integrals or of generalized hypergeometric functions in the non-linear one. Moreover, it appears a remarkable phenomenon when q is even: the probability to have a center is positive, in con- trast with what happens in the linear case, or also when q is odd, where this probability is zero. 1. Introduction and main results In this work, a random planar vector field will be a polynomial vector field whose coefficients are normally distributed independent random variables with zero mean and standard deviation one. This is the natu- ral distribution when one is worried about the probability of appearance of some phase portrait for given family of planar vector fields, see for instance [13, 15] or [4, Thm 2.1] to have more details. Some related papers that use a similar approach and also study planar random sys- tems are [2, 3, 16]. In [2] the quadratic systems are considered, in [3] the authors study the homogeneous systems of degree 1,2 and 3 and in [16] the linear systems. Since we will deal with a class of quasi-homogeneous planar polyno- mial vector fields, in next section we will recall some usual notations and definitions and some of their properties. One of the most impor- tant to our interests is that their local phase portraits at the origin determine their global phase portraits. It is also worth to comment that here are many mathematical mod- els involving differential equations where coefficients and parameters come from sources with uncertainty. For instance, this is the case of some epidemiological or viral expansion models. This variability may Key words and phrases. Ordinary differential equations with random coefficients; planar quasi-homogeneous vector fields; index; phase portraits. 2020 Mathematics Subject Classification: 37H10, 34F05 . 1 2 B. COLL, A. GASULL AND R. PROHENS be due to errors in measurements, virus replication and mutation, and many others. A way to tackle this uncertainty, i.e. a way to establish with accuracy the parameters, is by assuming them to be random vari- ables. As we have already said, under this assumption, it is commonly supposed that parameters follow a Gaussian distribution. In this way, diverse phenomena can be approached by using different mathematical approaches. See [5, 6, 7], for instance. In this paper we consider the family of random quasi-homogeneous vector fields ẋ = Axq + By, ẏ = Cx2q−1 + Dxq−1 y, (1) where q is an integer positive number, and A, B, C and D independent and identically distributed (iid) random variables with normal distri- bution N (0, 1), and we study the probability of appearance of each of its phase portraits. When we consider a realization of system (1), we obtain a deterministic equation given by ẋ = axq + by, ẏ = cx2q−1 + dxq−1 y, (2) where a, b, c, d ∈ R. Notice that phase portraits of deterministic sys- tems, characterized by equalities among the parameters a, b, c and d, will have probability zero when they are regarded in the probabilistic setting due to the fact that A, B, C, D are iid random variables with a continuous distribution. For this reason these situations will be disre- garded in our study. As we will prove, the only phase portraits of system (1) with positive probability, depend on the parity of q and correspond to the following cases: saddle, elliptic+hyperbolic, node, center and focus. Moreover, all these phase portraits will be not only local, but global, due the the quasi homogeneity property. Observe that the saddle has index −1 while all the above other type of critical points have index +1. Recall that a critical point is said to be monodromic if given a transversal section through it, the flow defines a return map on this section. In particular, the only monodromic points in the above list are the center and the focus. By way of notation, the subscripts m and nm in our main result stand for “monodromic” and “non-monodromic” critical points of index +1. Theorem 1.1. Consider the quasi-homogeneous random vector field (1), which has weight exponents s1 = 1 and s2 = q and weight degree q. Their global phase portraits with positive probability are: • Saddle, node or focus, when q is odd. • Saddle, elliptic+hyperbolic or center, when q is even. PLANAR RANDOM QUASI-HOMOGENEOUS VECTOR FIELDS 3 Moreover, if Ps = P (saddle), ( P (node) when q is odd, Pnm = Pnm (q) = P (elliptic+hyperbolic) when q is even, ( P (focus) when q is odd, Pm = Pm (q) = P (center) when q is even, it holds that Ps = 1/2, Pnm (q) = 1/2 − Pm (q), and Z s 1 π/2 sin(x) 1 + q2 Pm (q) = dx, where Q = . (3) π 0 sin(x) + Q 2q It is a remarkable and somehow surprising fact that when q is even the probability to have a center is positive, in contrast with what hap- pens in the linear case, or also when q is odd, where this probability is zero. As a corollary of our result we recover a well-known result for the linear case q = 1 : √ √ 1 2−1 2 P (saddle) = , P (node) = , P (focus) = 1 − . 2 2 2 See [3, Thm. 1], for instance. We remark that only for q = 1 the value Pm (q) √ can be obtained explicitly, giving that Pm (1) = P (focus) = 1 − 2/2, see the beginning of Section 4 for a proof. For q > 1, in that section we express Pm (q) in terms of complete elliptic integrals or of generalized hypergeometric functions. We also develop some tools to approach this value. With this objective in mind we study properties of the function q → Pm (q). For instance we prove that it is a decreasing √ convex function or that Pm (q) ∼ M/ q, at q = ∞, for a given explicit M. For completeness we also approximate Pm (q), for some values of q, by using Monte Carlo method. The study of the deterministic system (2) is done in Section 2. Sec- tion 3 is devoted to study the random system (1) and to prove Theo- rem 1.1. In this Theorem one of the relevant results is the computation of the probability Pm (q). This value is obtained in two different ways, firstly simply by reducing a 3-dimensional integral to a 1-variable def- inite integral, and secondly by relating it with the number of expected real roots of a random polynomial that controls the number of quasi- homogenenous invariant curves of (1). This expected value is then computed by using the nice Edelman-Kostlan formula ([10, 11]). 2. Results for the deterministic case In this section we recall some general properties of planar quasi- homogeneous vector fields and we study the global phase portraits of the deterministic differential system (2). 4 B. COLL, A. GASULL AND R. PROHENS A polynomial vector field X = (P, Q) is (s1 , s2 ) quasi-homogeneous, where (s1 , s2 ) ∈ N2 , if P (λs1 x, λs2 y) = λs1 +r−1 P (x, y), Q(λs1 x, λs2 y) = λs2 +r−1 Q(x, y), (4) for all λ ∈ R+ = {a ∈ R : a > 0} and some non negative integer r. We call s1 and s2 the weight exponents of the vector field X and r the weight degree (or quasi-degree) with respect to the weight exponents s1 and s2 . A well-known nice dynamical feature of such systems is that the local behaviour near the origin gives also its global behaviour. For completeness we state and this fact in next proposition. Proposition 2.1. Consider the differential system ẋ = P (x, y), ẏ = Q(x, y), with (P, Q) is a quasi-homogeneous polynomial vector field satisfying (4). If (x(t), y(t)) is one of its solutions, then for all λ ∈ R+ ,  (x̄(t), ȳ(t)) = λs1 x(λr−1 t), λs2 y(λr−1 t) is another solution. As a consequence, the behaviour of all its solutions is controlled by the behaviour of its solutions in any neighbourhood of the origin. Moreover it has no limit cycles. Proof. The proof of the stated property can be done by simple compu- tations. For instance,  ˙ x̄(t) = λs1 +r−1 ẋ(λr−1 t) = λs1 +r−1 P x(λr−1 t), y(λr−1 t)   s1 +r−1 x̄(t) ȳ(t) =λ P , = P (x̄(t), ȳ(t)). λs1 λs2 Then the first result follows because, for a suitable λ, the weighted homothety (x, y) → (λs1 x, λs2 y) transforms any small neighbourhood of the origin to any arbitrarily large neighbourhood of the origin. The non existence of limit cycles is also a consequence of the same property because if the system has a periodic solution then all solutions must be periodic as well and it has a global center.  Next result gives, in the generic cases, the global phase portrait of a realization of system (1). Its proof is based on the classical work of A.F. Andreev about the study of nilpotent critical points, see [1], and on Proposition 2.1, that allows to transform local results into global ones. Proposition 2.2. Consider the quasi-homogeneous vector field ẋ = axq + by, ẏ = cx2q−1 + dxq−1 y, where a, b, c, d ∈ R and define φ = ad − bc and δ = (qa − d)2 + 4qbc = (qa + d)2 − 4qφ. (5) PLANAR RANDOM QUASI-HOMOGENEOUS VECTOR FIELDS 5 Assume that bφ 6= 0. Then the origin is its unique singularity and it is a global saddle when φ < 0 (index -1) and a point of index +1 when φ > 0. Moreover, in this latter case, (i) If δ ≥ 0, the the origin is a global non-monodromic point and: (i.1) if q is even then the phase portrait is globally formed by the union of a hyperbolic and an elliptic sector, (i.2) if q is odd then the origin is a global node, stable if qa+d < 0 and unstable if qa + d > 0. (ii) If δ < 0 then the origin is a global monodromic point and: (ii.1) if q is odd then it is a global focus, stable if qa + d < 0 and unstable if qa + d > 0; and a global center if qa + d = 0, (ii.2) if q is even then it is a global center. Proof. The case q = 1 corresponds to the linear case and the stated results are well-known. From now one, we will concentrate on the case q > 1. Recall that (2) is a quasi-homogeneous vector field with weight ex- ponents s1 = 1 and s2 = q and weight degree q. By Proposition 2.1 to get its global phase portrait it suffices to study its local behaviour at the origin. In [1] there is an explicit result that allows to know the local be- haviour of any planar isolated nilpotent singularity for any analytic vector field, modulus the so called center-focus problem (that is the distinction in the monodromic case between center and focus). An- dreev’s result is usually stated for systems in the normal form ẋ = y + X(x, y), ẏ = Y (x, y), where X and Y have expansions at the origin that begin at least with second order terms in x and y. Then the type of critical points (mod- ulus the center-focus problem) depends on some properties of the two functions Y (x, h(x)) = uxα (1 + O(x)),   ∂X ∂Y Ψ(x) = + (x, h(x)) = vxβ (1 + O(x)), ∂x ∂y where y = h(x) is the analytic function satisfying h(0) = 0 and such that h(x) + X(x, h(x)) = 0 and u 6= 0, 2 ≤ α ∈ N and 1 ≤ β ∈ N and v 6= 0 are well defined, unless Ψ = 0. In fact, this behaviour only depends on α, β, u, v when Ψ 6= 0 and on α and u when Ψ = 0. See [1, 8] for more details. For simplicity, before applying the above result, we start writing (2) in a more suitable form. Straightforward computations give that it is equivalent to the Liénard equation ẍ + f (x)ẋ + g(x) = 0, where f (x) = −(qa + d)xq−1 and g(x) = φx2q−1 . When b 6= 0, this equation 6 B. COLL, A. GASULL AND R. PROHENS leads us to the generalized Liénard system ẋ = z, ż = −g(x) − f (x)z = −φx2q−1 + (qa + d)xq−1 z, (6) that is already in Andreev’s form. Then, under condition bδ 6= 0, the origin is an isolated singularity for (6). Moreover h = 0 and Y (x, h(x)) = Y (x, 0) = −g(x) = −φx2q−1 , Ψ(x) = −f (x) = (qa + d)xq−1 . Hence u = −φ, α = 2q − 1 and either Ψ = 0, or v = qa + d and β = q − 1. In short, either Ψ = 0 or α = 2β + 1, and we can apply Andreev’s approach (we skip the details). We obtain the list of cases given in the statement, where in item (iii) his result only ensures that the origin is a monodromic critical point and we yet have the center- focus disjunctive. To resolve it we study separately the cases q odd and q even. In the first case, q odd, the divergence of the vector field X associated to (6) is div(X) = −f (x) = (qa + d)xq−1 . Hence, when qa + d 6= 0 it does not change sign (only vanishes on the straight line x = 0). Hence, by the divergence theorem, the system has not periodic orbits and as a consequence the origin is a global focus. Clearly, its stability is given the the sign of qa + d. When qa + d = 0, then f = 0 and system (6) is easily integrable and its solutions are contained in the closed curves qz 2 + δx2q = k, for 0 < k ∈ R. Hence a global center arises. When q is even, system (6) is invariant by the change of variables and time (x, z, t) −→ (−x, z, −t). Moreover, by Andreev’s approach we know that it is monodromic. Therefore, the Poincaré reversibility criterion implies that the origin is a global center.  Remark 2.3. (i) In the case q = 1, we get back to the linear case and this classification coincides with the non-degenerate linear one, where φ is the determinant of the matrix and δ is the discriminant of the characteristic polynomial. (ii) When for system (2), bφ = 0, the origin is no more an isolated critical point. The phase portraits for these cases are much easier to be obtained. We will not describe them because in the probabilistic case they will have probability zero. (iii) In this work we are not interested on the behaviour of the or- bits of system (2) on the Poincaré disk. In any case, it is worth to know that near infinity on the Poincaré-Lyapunov disk, the behaviour of the equivalent Liénard system (6) can be easily obtained by using the classification given in the work of F. Dumortier and C. Herssens [9]. Trying to understand why in the non-linear case the determinant φ = ad − bc plays a role similar to its role in the linear situation we have obtained an alternative way for proving Proposition 2.1. The PLANAR RANDOM QUASI-HOMOGENEOUS VECTOR FIELDS 7 point is that system (2) in the new variables ds X = xq , Y =y and = xq−1 (7) dt writes as X 0 = qaX + qbY, Y 0 = cX + dY. (8) Therefore, when q is odd, (7) is an actual change of variables, the new jacobian is simply q(ad − bc), with the same sign that ad − bc, and the new trace is qa + d. Hence the same type of critical points that when q = 1 appear. On the other hand, when q is even, (7) is no more a change of variables, but a folding. Then the phase portrait of (8) on the plane X > 0 appears diffeomorphically in the phase portrait of (2) on x > 0 and as its specular image through a mirror located in x = 0 when x < 0. In short, focus and nodes for system (8) go to centers, saddles go to saddles and nodes are transformed into points with one elliptic and one hyperbolic sector. 3. The probabilistic case: proof of Theorem 1.1 In this section, we prove Theorem 1.1 that gives the probability of the different phase portraits described in Proposition 2.2. Inspired by the results of that proposition we define the new random variables Φ and ∆, Φ = AD − BC and ∆ = (qA − D)2 + 4qBC = (qA + D)2 − 4qΦ. Hence the desired probabilities can be computed as follows Ps = P {ω | Φ(ω) < 0}, Pnm = P {ω | Φ(ω) > 0, ∆(ω) > 0}, Pm = P {ω | Φ(ω) > 0, ∆(ω) < 0} = P {ω | ∆(ω) < 0}, where notice that last equality holds because if ∆ < 0 then Φ > 0 due to the equality ∆ = (qA + D)2 − 4qΦ. Let us prove first that Ps = P (Φ < 0) = 1/2. From the fact that the variables A, B, C and D are continuous, independent and identically distributed, the variables Φ = AD − BC and Ψ = −AD + BC are continuous and identically distributed. As a consequence, P (Φ > 0) = P (Ψ > 0) = P (Φ < 0). Since P (Φ = 0) = 0 it follows that P (Φ > 0) = P (Φ < 0) = Ps = 1/2. Note also that 1 = P (Φ > 0) = P (Φ > 0, ∆ > 0) + P (Φ > 0, ∆ < 0) = Pnm + Pm . 2 The final step is to compute Pm = P (Φ > 0, ∆ < 0). If we de- note by Z = qA − D, then Z p is a normal random variable with zero mean and standard deviation 1 + q 2 . We consider the random vec- tor (X, Y, Z) = (B, C, qA − D), where X, Y, Z are independent normal 8 B. COLL, A. GASULL AND R. PROHENS p variables with zero mean, σX = σY = 1 and σZ = 1 + q 2 . Thus the joint density function is   1 (x2 + y 2 )(1 + q 2 ) + z 2 Ψ(x, y, z) = p exp − . (2π)3/2 1 + q 2 2(1 + q 2 ) Hence Z 2 Pm = P (∆ < 0) = P (Z + 4qXY < 0) = Ψ(x, y, z) dx dy dz, K where K = {(x, y, z) ∈ R : z + 4qxy < 0}. To compute this integral, 3 2 following the same ideas that in [3], we perform the change of variables p x = r sin(t) cos(s), y = r sin(t) sin(s), z = 1 + q 2 r cos(t), where t ∈ (− π2 , π2 ), s ∈ (0, 2π) and r > 0. Since the determinant of the p Jacobian of the change is 1 + q 2 sin(t)r2 , we have Z Z +∞ 1 1 2 P (∆ < 0) = 3/2 | sin(t)|r2 e− 2 r dr ds dt (2π) K̃ 0 where K̃ = {(s, t) : (1 + q 2 ) cos2 (t) + 2q sin2 (t) sin(2s) < 0}. Recall that Z +∞ r 2 − 21 r2 π r e dr = . 0 2 Hence Z 1 P (∆ < 0) = | sin(t)| ds dt. 4π K̃ We consider the 2 qcurve (1 + q )  cos2 (t) + 2q sin2 (t) sin(2s) = 0, that is, 1+q 2 √ 1 t = ± arctan 2q =: ±g(s). Taking into account the − sin(2s) symmetries we obtain that Z Z Z 1 1 π π/2 | sin(t)| ds dt = sin(t) dt ds. 4π K̃ π π/2 g(s) √ By using that cos(arctan(α)) = 1/ α2 + 1, and doing u = 2s, we get Z π Z π/2 Z s 1 2π sin(u) sin(t) dt ds = du, π/2 g(s) 2 π sin(u) − Q 2 where Q = 1+q 2q . Note that Q ≥ 1. Finally, by doing the change of variable x = u − π and taking into account the symmetry of sin(x), Z 2π s Z πs Z π/2 s sin(u) sin(x) sin(x) du = dx = 2 dx π sin(u) − Q 0 sin(x) + Q 0 sin(x) + Q At this point, from all previous calculations we get Z π/2 s 1 sin(x) Pm = Pm (q) = P (∆ < 0) = dx, π 0 sin(x) + Q PLANAR RANDOM QUASI-HOMOGENEOUS VECTOR FIELDS 9 as we wanted to prove. There is a alternative way for obtaining Pm (q) which uses a beautiful result od Edelman and Kostlan for computing the expected number of real zeros of a family of polynomials whose coefficients are Gaussian random variables, see [10, 11] Alternative way for obtaining Pm (q). This different approach has two steps: First, we relate Pm with the expected number of invariant curves of the form αy + βxq = 0 by the flow of the random system (1). Sec- ondly, we compute this expected number. For linear vector fields the existence of invariant lines αx+βy = 0 and the flow over them allows to determine the type of phase portrait. For instance, when a linear planar systems has not invariant straight lines then the origin is monodromic and it is either a focus or a center. In the same way, the existence of invariant curves of the form αxq +βy = 0 and the flow over them plays a similar role to know the global phase portrait of a quasi-homogeneous system (2). Since the line x = 0 is invariant only when b = 0 we can skip this case (this is so, because it corresponds to the event B = 0 for the random vector field (1) and therefore it has probability of appearance 0). Hence, we restrict our attention to find conditions for (2) to have an invariant curve of the form y = λxq , where λ = −α/β. To get these conditions on λ notice that d  0= (y − λxq )|y=λx = ẏ − λqxq−1 ẋ y=λx dt  = (cx2q−1 + dxq−1 y) − λqxq−1 (axq + by) y=λx  = c + (d − qa)λ − qbλ2 x2q−1 . Hence the conditions on λ for the random system (1) to have an in- variant curve y = λxq is that T (λ) = C2 λ2 + C1 λ + C0 = 0, where T is a random polynomial and C2 = −qB, C1 = D − qA and C0 = C. Moreover, since A, B, C, D have N (0, 1) distribution and are independent, that implies C0 is N (0, 1), C1 is N (0, 1 + q 2 ) and C2 is N (0, q 2 ) and they are also independent. In a while we will compute the number of expected roots E(q) of this random polynomial. Let as see, firstly how this value is related with Pm (q). The point is that the following results hold: • If a realization of the polynomial T has two real roots then the corresponding of system (1) has either a saddle or a non- monodromic critical point of index +1. • If a realization of the polynomial T has not real roots then the corresponding of system (1) has a monodromic critical point. 10 B. COLL, A. GASULL AND R. PROHENS • That a realization of the polynomial T has exactly one real root happens with probability zero. As a consequence, the expectation number of real roots of T is E(q) =2 · P (T has 2 real roots) + 1 · P (T has 1 real root) + 0 · P (T has 0 real roots) = 2 · (Ps + Pnm ) + 1 · 0 + 0 · Pm . By using that Ps = 1/2 and that Pnm = 1/2 − Pm we obtain that E(q) Pm (q) = 1 − , (9) 2 which is the desired relation. To end the proof, let us compute E(q). The coefficients of T (λ), Ci , i = 1, 2, 3, are independent normal random variables with zero mean and covariance matrix M where   1 p 0 0 M 1/2 = 0 1 + q 2 0 . 0 0 q Following [10, Thm 3.1], if we define w(t) = M 1/2 ·(1, t, t2 )T and w(t) = w(t)/||ω(t)||, then the expected number of real zeros of T is given by the Edelman-Kostlan formula: Z 1 ∞ E(q) = ||w0 (t)|| dt. π −∞ By straightforward computations, we obtain 1 p  w(t) = p 1, q 2 + 1 t, qt2 q 2 t4 + (q 2 + 1)t2 + 1 and as a consequence Zp ∞ 1 (q 4 + q 2 )t4 + 4q 2 t2 + q 2 + 1 E(q) = dt (10) π −∞ (t2 + 1)(q 2 t2 + 1) Z p 1 ∞ (q 4 + q 2 )s2 + 4q 2 s + q 2 + 1 = √ ds. (11) π 0 (s + 1)(q 2 s + 1) s where this last equality follows by taking s = t2 . Finally, by using (9) we obtain a new expression of Pm (q).  4. Equivalent expressions and properties of Pm We start giving a different expression of the integral given in (3) that is suitable for obtaining Pm (1). Similarly, (11) could be used. In this integral we do the change of variables y = 1/ sin(x). We obtain that Z 1 ∞ 1 Pm (q) = p dy. (12) π 1 y (Qy + 1)(y 2 − 1) PLANAR RANDOM QUASI-HOMOGENEOUS VECTOR FIELDS 11 If √ we substitute when q = 1, then Q = 1, and by introducing w = y − 1 we obtain that Z Z 1 ∞ 1 1 ∞ 2 Pm (1) = √ dy = dw π 1 y(y + 1) y − 1 π 0 (1 + w2 )(2 + w2 ) Z   1 ∞ 2 2 = − dw π 0 1 + w2 2 + w2 ∞ √  ∞ √ 2 2 w 2 = arctan(w) − arctan √ =1− . π 0 π 2 0 2 Another similar expression for Pm (q) follows, again starting from ex- pression (3), by using the related change sin(x) = y. We arrive to Z r 1 1 y Pm (q) = dy. (13) π 0 (y + Q)(1 − y 2 ) Next we present two expressions of Pm (q) in terms of some classical transcendental function obtained by using Maple and Mathematica. First expression, starting from (3):      1 2 3 1 3 3 1 5 1 Pm (q) = √ 24Q Γ 3 F2 , , ; , ; 12π 2πQ3/2 4 4 4 4 2 4 Q2     2 1 3 5 5 3 7 1 −Γ 3 F2 , , ; , ; 4 4 4 4 2 4 Q2 where Γ is the gamma function and 3 F2 is the generalized hypergeo- metric function. Second expression, starting from (12): √ s ! s !! 2 1 Q−1 Q−1 Pm (q) = √ (1 + Q) Π − , − QK π Q Q 2Q 2Q where K( · ), the EllipticK function, is the complete elliptic integral of the first kind while Π( · , · ), the EllipticPi function, is the complete elliptic integral of the third kind. 4.1. Some properties of Pm (q). Recall that from (3) it holds that  2  Z π/2 s 1+q 1 sin(x) Pm (q) = H , where H(Q) = dx. (14) 2q π 0 sin(x) + Q Next result collect several properties of H that can be used to approach Pm (q). Proposition 4.1. Let H : (0, ∞) → R be the function defined in (14). The following holds: (i) It is completely monotone, that is, (−1)n H (n) (Q) > 0, for all n ≥ 0. 12 B. COLL, A. GASULL AND R. PROHENS (ii) For all n ≥ 0, it holds that Z ∞ (n) (−1)n (2n − 1)!! (1 + w2 )n−1 H (1) = dw 2n−1 π 0 (2 + w2 )n+1 and all these values can be computed explicitly in terms of ele- mentary functions. (iii) Let Tm (Q) denote the Taylor polynomial of degree m at Q = 1, P H (n) (1) that is, Tm (Q) = m n=0 n! (Q − 1)n . Then, for all j, k ≥ 1, and all Q > 1, T2j−1 (Q) < H(Q) < T2k (Q). (15) (iv) For all Q > 0, L L √ ≤ H(Q) ≤ √ , (16) Q+1 Q where Γ (3/4) L= √ ≈ 0.38138. 2 π Γ (5/4) √ and hence, H(Q) ∼ L/ Q at Q = ∞. (v) For Q > 1 it holds that X∞   −1/2 Γ (3/4 + n/2) 1 H(Q) = √ n √ . n=0 n 2 π Γ (5/4 + n/2) Q Q q sin(x) Proof. (i) Notice first that for Q > 0 the function sin(x)+Q and all its derivatives with respect to Q are integrable in [0, π/2]. Therefore, for all n ≥ 0, s Z π/2 n s n Z π/2 1 d sin(x) 1 ∂ sin(x) H (n) (Q) = dx = dx π dQn 0 sin(x) + Q π 0 ∂Qn sin(x) + Q Z π/2 s n (−1) (2n − 1)!! 1 sin(x) = n n dx. 2 π 0 (sin(x) + Q) sin(x) + Q Then, clearly (−1)n H (n) (Q) > 0, as we wanted to show. We remark that complete monotone function are a subject of interest by them- selves, see for instance [17]. (ii) By replacing Q = 1 above and introducing the new variable y = 1/ sin(x), and then w2 = y − 1, as in the beginning of this section to obtain Pm (1), we get that Z (n) (−1)n (2n − 1)!! ∞ y n−1 dy H (1) = n n+1 √ 2 π 1 (1 + y) y−1 n Z ∞ 2 n−1 (−1) (2n − 1)!! (1 + w ) = dw, 2n−1 π 0 (2 + w2 )n+1 PLANAR RANDOM QUASI-HOMOGENEOUS VECTOR FIELDS 13 as desired. For each given n the above functions have a primitive in terms of elementary functions and the values H (n) (1) can be explicitly obtained. For instance, √ √ √ 0 2 00 15 2 000 195 2 H (1) = − , H (1) = , H (1) = − . 16 256 2048 (iii) Both inequalities are consequence of Taylor’s formula and the results of item (i). (iv) Clearly, for all 0 ≤ x ≤ π/2, Q ≤ sin(x) + Q ≤ 1 + Q. Hence p s p sin(x) 1 sin(x) sin(x) √ ≤ ≤ √ . π 1+Q π sin(x) + Q π Q By√integrating the above√inequalities between 0 and π/2 we get that L/ 1 + Q ≤ Pm (q) ≤ L/ Q, where Z 1 π/2 p Γ (3/4) L= sin(x) dx = √ . π 0 2 π Γ (5/4) We have used the equality Z π/2 √ α Γ (1/2 + α/2) π sin (x) dx = , (17) 0 Γ (1 + α/2) 2 valid for all α > −1, when α = 1/2. √ (v) If we introduce the new variable U = 1/ Q > 0 we have that H(Q) = G(U ), where G is the smooth function at U = 0, Z s Z s 1 π/2 U 2 sin(x) U π/2 sin(x) G(U ) = 2 dx = 2 dx π 0 U sin(x) + 1 π 0 U sin(x) + 1 Z ∞   U π/2 X p −1/2 = sin(x) U 2n sinn (x) dx π 0 n=0 n ∞   Z π/2 ! U X −1/2 = sinn+1/2 (x) dx U 2n , (18) π n=0 n 0 where we have used the uniform convergence of the series to interchange it with the integral. By using again (17), with√ α = n + 1/2, the final expression of G follows. By replacing U by 1/ Q we get the desired result.  From results (ii) and (iii) of the above proposition we can obtain good approximations of Pm (q) for q not big. This is so just by recalling that, for big values of q, the remainder term of the Taylor’s formula can give a big error. On the contrary, since the boundedness given by expression (16) is better the higher the q value is taken, the results of item (iv) are suitable for big q. Finally, although the results of item (v) apply for all q > 1, its convergence is faster when 1/Q is not big, hence 14 B. COLL, A. GASULL AND R. PROHENS we can apply it for medium and large values of q. Let us see some examples where we approach Pm (2), Pm (10) and Pm (100). An approximation of Pm (2). If q = 2, then Q = 5/4. By item (iii) we know for instance that       5 5 5 T3 1, H(Q) = limk→∞ Hk (Q), where Xk   −1/2 Γ (3/4 + n/2) 1 Hk (Q) = √ n √ . n=0 n 2 π Γ (5/4 + n/2) Q Q Hence Pm (10) = limk→∞ Hk (101/20). If we compute Hk (101/20) for k = 2, 3, 4, 5, 6 we get 0.15895 . . . , 0.158736 . . . , 0.1587702 . . . , 0.1587648 . . . , 0.1587657 . . . providing good approximations of Pm (10) ≈ 0.158766, see again Ta- ble 1. For bigger values of Q the convergence is faster. For instance, taking q = 100, Q = 10001/200 and Pm (100) ≈ H3 (10001/200) = 0.053543958 . . . PLANAR RANDOM QUASI-HOMOGENEOUS VECTOR FIELDS 15 4.2. Some numerical simulations. Although we have a closed form in terms of a defined integral for Pm (q) that allows to compute it simply by computing numerically this integral, in this section we show how to approximate Pm (q) for some values of q, by using the celebrated Monte Carlo method. In our setting we simply will take N = 106 samples of the random vector (A, B, C, D) where the four variables are iid, with distribution N (0, 1), and check how many of them, say J, satisfy ∆ = (qA − D)2 + 4qBC < 0. Then simply Pm (q) ≈ J/N. Due to the law of large numbers and the law of iterated logarithm it is known that this approach gives an absolute error of order O(((log log N )/N )1/2 ), which in practice behaves as O(N −1/2 ), see [12, 15]. Hence, since N = 106 this absolute error is expected to be of order O(10−3 ) and in Table 1 we only show 4 digits of the Monte Carlo results. Notice that comparing these values with the ones obtained by approximating the 1-variable definite integral, given in the second row of that table, the actual differences are the expected ones. To see more details about the probability that |Pm (q) − J/N | is big see the discussion in [4, Sec. 3.2]. q 1 2 5 10 100 Pm (q) 0.292893 0.273080 0.209525 0.158766 0.053544 Monte Carlo 0.2920 0.2732 0.2091 0.1589 0.0535 Table 1. Some approximated values of Pm (q) obtained computing numerically the expression (3) and with Monte Carlo method taking 106 random systems. Acknowledgements The authors are supported by Ministry of Economy, Industry and Competitiveness-State Research Agency of the Spanish Government through grants (MINECO/AEI/FEDER, UE) PID2020-118726GB-I00, first and third authors; and by PID2019-104658GB-I00, second au- thor. The second author is partially supported by the grants 2017- SGR-1617 from AGAUR, Generalitat de Catalunya and Severo Ochoa and Marı́a de Maeztu Program for Centers and Units of Excellence in R&D (CEX2020-001084-M) References [1] A.F. Andreev. Investigation of the behaviour of the integral curve of a system of two differential equations in the neighbourhood of a singular point, Translation of AMS 8 (1958), 183–207. 16 B. COLL, A. GASULL AND R. PROHENS [2] J.C. Artés, J. Llibre, Statistical measure of quadratic vector fields, Resenhas 6 (2003), 85–97. [3] A. Cima, A. Gasull, V. Mañosa. Phase portraits of random planar homogeneous vector fiels. Qual. Theory Dyn. Syst. 20 (2021), paper No 3, 27 pp. [4] A. Cima, A. Gasull, V. Mañosa. Stability index of linear random dynamical sys- tems. Elec. J. of Qual. Theory of Differential Equations, Paper No. 15 (2021), 27 pp. [5] T.P. Chang. Dynamic response of a simplified nonlinear fluid model for vis- coelastic materials under random parameters. Math. Comput. Modelling 54 (2011), 2587–2596. [6] B.M. Chen-Charpentier, D. Stanescu. Epidemic models with random coeffi- cients. Math. Comput. Modelling 52 (2010), 1004–1010. [7] B.M. Chen-Charpentier, D. Stanescu. Virus propagation with randomness. Math. Comput. Modelling 57 (2013), 1816–1821. [8] F. Dumortier, J.C. Artés, J. Llibre. “Qualitative Theory of planar differential systems,”Springer, Berlin, 2006. [9] F. Dumortier, C. Herssens. Polynomial Liénard equations near infinity. J. Dif- ferential Equations 153 (1999), 1–29. [10] A. Edelman, E. Kostlan, How many zeros of a random polynomial are real? Bull. Am. Math. Soc., New Ser. 32, (1995), 1–37. [11] A. Edelman, E. Kostlan, Erratum to “How many zeros of a random polynomial are real?” Bull. Am. Math. Soc., New Ser. 33, (1996), 325. [12] C. Lemieux, “Monte Carlo and quasi-Monte Carlo sampling,”(English) Springer Series in Statistics, New York, NY 2009. [13] G. Marsaglia. Choosing a point from the surface of a sphere. Ann. Math. Statist. 43 (1972) 645–646. [14] B.J.T. Morgan, “Applied stochastic modelling,”2nd ed. Chapman & Hall/CRC Texts in Statistical Science Series. Boca Raton, FL, 2009. [15] M.E. Muller, A note on a method for generating points uniformly on N - dimensional spheres. Commun. ACM 2 (1959), 19–20. [16] B.K. Pagnoncelli, H. Lopes and C.F.B. Palmeira. Sampling linear ODE. (in Portuguese). Matemática Universitária 45 (2009), 44–50. [17] R. Schilling, R. Song, Z. Vondraček, “Bernstein functions. theory and applica- tions,”De Gruyter Berlin, 2002. †∗ Dept. de Matemàtiques i Informàtica, IAC3 Institute of Applied Computing & Community Code, Universitat de les Illes Balears, 07122 Palma de Mallorca, Illes Balears, Spain Email address: † tomeu.coll@uib.cat, ∗ rafel.prohens@uib.cat ‡ Dept. de Matemàtiques, Universitat Autònoma de Barcelona, Ed- ifici C, 08193 Cerdanyola del Vallès (Barcelona), Spain. Centre de Recerca Matemàtica, Edifici Cc, Campus de Bellaterra, 08193 Cer- danyola del Vallès (Barcelona), Spain Email address: gasull@mat.uab.cat