Mediterr. J. Math. (2022) 19:278 https://doi.org/10.1007/s00009-022-02198-w 1660-5446/22/060001-16 published online November 6, 2022 c The Author(s) 2022 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 probability of occurrence of phase portraits in a family of planar quasi-homogeneous vector fields of quasi degree q, that is a natural extension 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 hypergeomet- ric 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 contrast with what happens in the linear case, or also when q is odd, where this probability is zero. Mathematics Subject Classification. 37H10, 34F05. Keywords. Ordinary differential equations with random coefficients, Pla- nar quasi-homogeneous vector fields, Critical point index, Phase por- traits. 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 natural distribution when one is worried about the probability of appearance of some phase por- trait for given family of planar vector fields, see for instance [15,16] or [4, Thm 2.1] to have more details. Some related papers that use a similar ap- proach and also study planar random systems are [2,3,17]. In [2] the quadratic systems are considered, in [3] the authors study the homogeneous systems of degree 1,2 and 3 and in [17] the linear systems. Since we will deal with a class of quasi-homogeneous planar polynomial vector fields, in next section we will recall some usual notations and definitions and some of their properties. One of the most important to our interests is 278 Page 2 of 16 B. Coll et al. MJOM that their local phase portraits at the origin determine their global phase portraits. It is also worth to comment that here are many mathematical models involving differential equations where coefficients and parameters come from sources with uncertainty. For instance, this is the case of some epidemio- logical or viral expansion models. This variability may 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 variables. As we have already said, under this assumption, it is commonly supposed that parameters follow a Gauss- ian distribution. In this way, diverse phenomena can be approached by using different mathematical approaches. See [5–7], for instance. A similar analysis can be performed in the case of linear systems of differential and difference equations of arbitrary dimension and not necessary Gaussian random variables. See [4,8,9], for instance. In this paper we consider the family of random quasi-homogeneous vec- tor 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 distribution N (0, 1), and we study the probability of appearance of each of its phase por- traits. 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 systems, characterized by equalities among the parameters a, b, c and d, will have prob- ability 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 disregarded 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. MJOM Planar Random Quasi-homogeneous Vector Fields Page 3 of 16 278 • Saddle, elliptic+hyperbolic or center, when q is even. 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   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 happens 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 Sect. 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 Sect. 2. Section 3 is devoted to study the random system (1) and to prove Theorem 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 definite integral, and secondly by relating it with the number of expected real roots of a random polynomial that controls the number of quasi-homogeneous invariant curves of (1). This expected value is then computed by using the nice Edelman-Kostlan formula ( [12,13]). 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). 278 Page 4 of 16 B. Coll et al. MJOM 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 computations. For instance,   ˙ x̄(t) = λs1 +r−1 ẋ(λr−1 t) = λs1 +r−1 P x(λr−1 t), y(λr−1 t)   x̄(t) ȳ(t) = λs1 +r−1 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 real- ization 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) Assume that bφ = 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: MJOM Planar Random Quasi-homogeneous Vector Fields Page 5 of 16 278 (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 expo- nents 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 behaviour 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). Andreev’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 (modulus 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 = 0, 2 ≤ α ∈ N and 1 ≤ β ∈ N and v = 0 are well defined, unless Ψ = 0.. This behaviour only depends on α, β, u, v when Ψ = 0 and on α and u when Ψ = 0. See [1,11] 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 = 0, this equation 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δ = 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 278 Page 6 of 16 B. Coll et al. MJOM 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 = 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 orbits 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 using the classification given in the work of F. Dumortier and C. Herssens [10]. 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 point is that system (2) in the new variables ds X = xq , Y = y and = xq−1 (7) dt writes as X  = qaX + qbY, Y  = 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. MJOM Planar Random Quasi-homogeneous Vector Fields Page 7 of 16 278 3. The Probabilistic Case: Proof of Theorem 1.1 In this section, we prove Theorem 1.1 that gives the probability of the dif- ferent 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 vari- ables A, B, C and D are continuous, independent and identically distributed, the variables Φ = AD − BC and Ψ = −AD + BC are continuous and iden- tically 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) 2 = Pnm + Pm . The final step is to compute Pm = P (Φ > 0, Δ < 0). If we denote by Z = qA − D, then Z is a normal random variable with zero mean and standard deviation 1 + q 2 . We consider the random vector (X, Y, Z) = (B, C, qA − D), where X, Y, Z are independent normal 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) = exp − . (2π)3/2 1 + q 2 2(1 + q 2 ) Hence  2 Pm = P (Δ < 0) = P (Z + 4qXY < 0) = Ψ(x, y, z) dx dy dz, K where K = {(x, y, z) ∈ R3 : z 2 + 4qxy < 0}. To compute this integral, following the same ideas that in [3], we perform the change of variables 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 Jacobian of the change is 1 + q 2 sin(t)r2 , we have   +∞ 1 1 2 P (Δ < 0) = | sin(t)|r2 e− 2 r dr ds dt (2π)3/2 K̃ 0 where K̃ = {(s, t) : (1 + q 2 ) cos2 (t) + 2q sin2 (t) sin(2s) < 0}. Recall that  +∞ 1 2 π r2 e− 2 r dr = . 0 2 278 Page 8 of 16 B. Coll et al. MJOM Hence  1 P (Δ < 0) = | sin(t)| ds dt. 4π K̃ 2 2 2  the curve (1 +q ) cos (t) + 2q sin (t) sin(2s) = 0, that is, t = We consider 1+q 2 √ 1 ± arctan 2q =: ±g(s). Taking into account the symmetries − sin(2s) we obtain that   π  π/2 1 1 | sin(t)| ds dt = sin(t) dt ds. 4π K̃ π π/2 g(s) √ Using that cos(arctan(α)) = 1/ α2 + 1, and doing u = 2s, we get  π  π/2   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),  2π   π  π/2  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   1 π/2 sin(x) Pm = Pm (q) = P (Δ < 0) = dx, π 0 sin(x) + Q 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 re- al zeros of a family of polynomials whose coefficients are Gaussian random variables, see [12,13] Proof. (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). Secondly, 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 MJOM Planar Random Quasi-homogeneous Vector Fields Page 9 of 16 278 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 invariant 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 corre- sponding 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 corre- sponding of system (1) has a monodromic critical point. • 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 . 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 0 0 M 1/2 = ⎝0 1 + q 2 0⎠ . 0 0 q Following [12, 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:  1 ∞ E(q) = ||w (t)|| dt. π −∞ 278 Page 10 of 16 B. Coll et al. MJOM By straightforward computations, we obtain 1   w(t) = 1, q 2 + 1 t, qt2 2 4 2 2 q t + (q + 1)t + 1 and as a consequence  ∞ 1 (q 4 + q 2 )t4 + 4q 2 t2 + q 2 + 1 (10) E(q) = dt π −∞ (t2 + 1)(q 2 t2 + 1) ∞ 1 (q 4 + q 2 )s2 + 4q 2 s + q 2 + 1 (11) = √ ds. π 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  1 ∞ 1 Pm (q) = dy. (12) π 1 y (Qy + 1)(y 2 − 1) √ If we substitute when q = 1, then Q = 1, and by introducing w = y − 1 we obtain that   1 ∞ 1 1 ∞ 2 Pm (1) = √ dy = dw π 1 y(y + 1) y − 1 π 0 (1 + w2 )(2 + w2 )    1 ∞ 2 2 = − dw π 0 1 + w2 2 + w2 ∞ √  ∞ √ 2 2 w 2 = arctan(w) − arctan √ =1− . π 0 π 2 2 0 Another similar expression for Pm (q) follows, again starting from expres- sion (3), using the related change sin(x) = y. We arrive to  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 using Maple and Mathematica. First expression, starting from (3):      1 3 1 3 3 1 5 1 Pm (q) = √ 24Q Γ 2 3 F2 , , ; , ; 12π 2πQ3/2 4 4 4 4 2 4 Q2     1 3 5 5 3 7 1 − Γ2 3 F2 , , ; , ; 4 4 4 4 2 4 Q2 where Γ is the gamma function and 3 F2 is the generalized hypergeometric function. MJOM Planar Random Quasi-homogeneous Vector Fields Page 11 of 16 278 Second expression, starting from (12): √       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     1 + q2 1 π/2 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. (ii) For all n ≥ 0, it holds that  (−1)n (2n − 1)!! ∞ (1 + w2 )n−1 H (n) (1) = dw 2n−1 π 0 (2 + w2 )n+1 and all these values can be computed explicitly in terms of elementary functions. (iii) Let Tm (Q) denote the Taylor polynomial of degree m at Q = 1, that is, m (n) Tm (Q) = n=0 H n!(1) (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 ∞   −1/2 Γ (3/4 + n/2) 1 H(Q) = √ √ . n Q n=0 n 2 π Γ (5/4 + n/2) 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 278 Page 12 of 16 B. Coll et al. MJOM n ≥ 0,     π/2 1 dn sin(x) 1 π/2 ∂ n sin(x) H (n) (Q) = dx = dx π dQn 0 sin(x) + Q π 0 ∂Q n sin(x) + Q   (−1)n (2n − 1)!! π/2 1 sin(x) = dx. 2n π 0 (sin(x) + Q)n 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 themselves, see for instance [18]. (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  (−1)n (2n − 1)!! ∞ y n−1 dy H (n) (1) = n n+1 √ 2 π 1 (1 + y) y −1 n  ∞ 2 n−1 (−1) (2n − 1)!! (1 + w ) = dw, 2n−1 π 0 (2 + w2 )n+1 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, √ √ √  2  15 2  195 2 H (1) = − , H (1) = , H (1) = − . 16 256 2048 (iii) Both inequalities are consequence of Taylor’s formula and the re- sults of item (i). (iv) Clearly, for all 0 ≤ x ≤ π/2, Q ≤ sin(x) + Q ≤ 1 + Q. Hence  sin(x) 1 sin(x) sin(x) √ ≤ ≤ √ . π 1+Q π sin(x) + Q π Q √ √ the above inequalities between 0 and π/2 we get that L/ 1 + Q ≤ By integrating Pm (q) ≤ L/ Q, where  π/2 1 Γ (3/4) L= sin(x) dx = √ . π 0 2 π Γ (5/4) We have used the equality  π/2 √ Γ (1/2 + α/2) π sinα (x) dx = , (17) 0 Γ (1 + α/2) 2 valid for all α > −1, when α = 1/2. MJOM Planar Random Quasi-homogeneous Vector Fields Page 13 of 16 278 √ (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,     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  ∞   U π/2  −1/2 = sin(x) U 2n sinn (x) dx π 0 n=0 n ∞    π/2  U −1/2 (18) = sinn+1/2 (x) dx U 2n , π n=0 n 0 where we have used the uniform convergence of the series to interchange it with the integral. 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 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 k   −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 Table 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 . . . 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, 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 [14,16]. 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]. Acknowledgements Our manuscript has no associate data. Data sharing not applicable to this ar- ticle as no datasets were generated or analysed during the current study. The authors are supported by Ministry of Economy, Industry and Competitiveness- State Research Agency of the Spanish Government through grants MJOM Planar Random Quasi-homogeneous Vector Fields Page 15 of 16 278 (MINECO/AEI/FEDER, UE) PID2020-118726GB-I00, first and third au- thors; and by PID2019-104658GB-I00, second author. 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) Funding Information Open Access funding provided thanks to the CRUE- CSIC agreement with Springer Nature. Open Access. This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and re- production in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party ma- terial in this article are included in the article’s Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article’s Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permis- sion directly from the copyright holder. To view a copy of this licence, visit http:// creativecommons.org/licenses/by/4.0/. Publisher’s Note Springer Nature remains neutral with regard to jurisdic- tional claims in published maps and institutional affiliations. References [1] Andreev, A.F.: Investigation of the behaviour of the integral curve of a system of two differential equations in the neighbourhood of a singular point. Transl. AMS 8, 183–207 (1958) [2] Artés, J.C., Llibre, J.: Statistical measure of quadratic vector fields. Resenhas 6, 85–97 (2003) [3] Cima, A., Gasull, A., Mañosa, V.: Phase portraits of random planar homoge- neous vector fiels. Qual. Theory Dyn. Syst. 20, paper No 3, 27 pp (2021) [4] Cima, A., Gasull, A., Mañosa, V.: Stability index of linear random dynamical systems. Elec. J. of Qual. Theory Differ. Equ. 15, 27 (2021) [5] Chang, T.P.: Dynamic response of a simplified nonlinear fluid model for vis- coelastic materials under random parameters. Math. Comput. Model. 54, 2587– 2596 (2011) [6] Chen-Charpentier, B.M., Stanescu, D.: Epidemic models with random coeffi- cients. Math. Comput. Model. 52, 1004–1010 (2010) [7] Chen-Charpentier, B.M., Stanescu, D.: Virus propagation with randomness. Math. Comput. Model. 57, 1816–1821 (2013) [8] Cortés, J.-C., Navarro-Quiles, A., Romero, J.-V., Roselló, M.-D.: Probabilistic solution of random autonomous first-order linear systems of ordinary differen- tial equations. Romanian Rep. Phys. 68, 1397–1406 (2016) [9] Cortés, J.-C., Navarro-Quiles, A., Romero, J.-V., Roselló, M.-D.: Full solution of random autonomous first-order linear systems of difference equations. Ap- plication to construct random phase portrait for planar systems. Appl. Math. Lett. 68, 150–156 (2017) 278 Page 16 of 16 B. Coll et al. MJOM [10] Dumortier, F., Herssens, C.: Polynomial Liénard equations near infinity. J. Differ. Equ. 153, 1–29 (1999) [11] Dumortier, F., Artés, J.C., Llibre, J.: Qualitative Theory of planar differential systems. Springer, Berlin (2006) [12] Edelman, A., Kostlan, E.: How many zeros of a random polynomial are real? Bull. Am. Math. Soc. New Ser. 32, 1–37 (1995) [13] Edelman, A., Kostlan, E.: Erratum to How many zeros of a random polynomial are real? Bull. Am. Math. Soc. New Ser. 33, 325 (1996) [14] Lemieux, C.: Monte Carlo and quasi-Monte Carlo sampling, (English) Springer Series in Statistics, New York, NY (2009) [15] Marsaglia, G.: Choosing a point from the surface of a sphere. Ann. Math. Statist. 43, 645–646 (1972) [16] Muller, M.E.: A note on a method for generating points uniformly on N- dimensional spheres. Commun. ACM 2, 19–20 (1959) [17] Pagnoncelli, B.K., Lopes, H., Palmeira, C.F.B.: Sampling linear ODE. Matemática Universitária (in Portuguese) 45, 44–50 (2009) [18] Schilling, R., Song, R., Vondraček, Z.: Bernstein functions. Theory and appli- cations. De Gruyter, Berlin (2002) B. Coll and R. Prohens Dept. de Matemàtiques i Informàtica, IAC3 Institute of Applied Computing and Community Code Universitat de les Illes Balears 07122 Palma Illes Balears Spain e-mail: rafel.prohens@uib.cat B. Coll e-mail: tomeu.coll@uib.cat A. Gasull Dept. de Matemàtiques Universitat Autònoma de Barcelona Edifici C, 08193 Cerdanyola del Vallès Barcelona Spain e-mail: gasull@mat.uab.cat and Centre de Recerca Matemàtica, Edifici Cc, Campus de Bellaterra, 08193 Cerdanyola del Vallès Barcelona Spain Received: January 15, 2022. Revised: April 29, 2022. Accepted: October 18, 2022.