Measuring the criticality of a Hopf bifurcation

Abstract

This work is based on the observation that the first Poincaré–Lyapunov constant is a quadratic function of the coefficients of the two-dimensional vector field at a Hopf bifurcation. From a given parameter point, we find the distance to the “Hopf quadric.” This distance provides a measure of the criticality of the Hopf bifurcation. The viability of the approach is demonstrated through numerical examples.

1 Introduction

An important task in nonlinear science is to determine the nature of the bifurcations the system under study exhibits. A similarly important but less frequently studied problem is the robustness of the system. Quite naturally, the Euclidean distance between the actual parameter values and those of the closest bifurcation can be used to characterize the robustness of the system.

Dobson [1] considers the general dynamical system

$$\begin{aligned} {\dot{\mathbf {x}}}=\mathbf {f}\left( \mathbf {x},{\mathbf {u}}_{0}\right) \end{aligned}$$

(1)

depending on parameters contained in the vector ${\mathbf {u}}_{0}$. The distance-to-bifurcation problem seeks the minimal parameter variation $\gamma =\left\| {\mathbf {u}}_{0}-{\mathbf {u}}\right\|$ for which Eq. (1) has a phase portrait qualitatively different from that of the original system. Dobson [1, 2] proposed an iterative and a direct method of finding the closest bifurcation in parameter space. These initial approaches used numerical methods to find the minimum of $\gamma$, but were computationally demanding and only guaranteed locally optimal solutions. A global search was proposed by Kremer [3]. These methods use the normal vectors at a bifurcation curve in a parameter plane and are suitable for bifurcations of equilibria. Tamba and Lemmon [4] recast the minimum distance-to-bifurcation problem as a sum-of-squares relaxation for nonnegative dynamical systems that have kinetic realizations. Kitajima and Yoshinaga [5] extend Dobson’s method for periodic orbits. These methods and their extensions have been applied to hydraulic systems [3], gene systems [6], power systems [7,8,9] and bifurcations of arbitrary codimension [10].

Measuring the distance from a possible bifurcation in parameter space is also useful in making parametric changes to modify the nature of the bifurcation. Abed and Fu [11] proposed the use of nonlinear feedback control to invert the direction of the bifurcation. Yuen and Bau [12] demonstrated theoretically and experimentally that through the use of a nonlinear feedback controller, one can render a subcritical Hopf bifurcation supercritical and dramatically modify the nature of the flow in a thermal convection loop.

A good measure of parametric robustness is necessary for the design-for-dynamics paradigm. “Designing” the nonlinearities of a system is a growing application in the engineering fields. In particular, we mention nonlinear energy sinks (NESs) that are passive mass–spring–damper systems with negligible linear stiffness term [13, 14]. Undesired vibrational energy of the primary structure flows toward the NES and is locally dissipated there, reducing oscillations over broad frequency and energy ranges due to the absence of a preferential resonance frequency. The basic objective for NES-based vibration reduction/suppression is to provide design rules of the nonlinearities. In O’Neil and Strganac [15], the authors describe a flutter test apparatus designed to permit experimental investigations of prescribed nonlinear response. This is done through a pair of cams shaped to provide a predetermined spring stiffness. Liu et al. [16] utilized a new way to introduce geometrically nonlinear stiffness and damping with the goal of vibration suppression.

Methods of evaluating the distance from a given point ${\mathbf {u}}_{0} \in \mathbb {R}^{n}$ to a bifurcation manifold utilized in the above-cited publications are counterparts of the Optimization Theory algorithms for solving similar problems for arbitrary manifolds in $\mathbb {R}^{n}$ given in parametric or implicit form. Therefore, the conditions for convergence of iterative Newton-type procedures impose the usual restrictions on the properties of the considered manifolds, such as smoothness or convexity. Several specific modifications are developed in Computational Geometry for manifolds like ellipsoids [17] or planar algebraic curves.

It should be noted that the bifurcation conditions arising in applications are frequently represented implicitly via a nonlinear polynomial equation $G({\mathbf {u}})=0$. For this case, an alternative approach for the point-to-manifold distance evaluation problem is based on symbolic algorithms of elimination of variables in the system of algebraic equations. These algorithms implement either the computation of the resultants or, generally, the construction of Gröbner bases [18, 19] and result in an univariate algebraic equation. This variable might be a particular coordinate of the nearest point in the manifold or even the distance to this point. For the case of planar conics, this method can be found in [20, 21], while for the case of general quadric in $\mathbb {R}^{n}$, this ideology was developed in [22, 23] to the concept of distance equation, i.e., the univariate equation whose zero set coincides with the set of critical values of the squared distance function.

The present paper is devoted to application of the distance equation for solving the problem of finding the distance to the Hopf bifurcation manifold in the parameter space. The paper is organized as follows. In Sect. 2, we detail the structure of the first Poincaré–Lyapunov constant. We identify it to be a quadric manifold in the ten-dimensional space of the coefficients of the normal form of a planar dynamical system. In Sect. 3, we briefly outline the general result on the representation of the distance equation for the general quadric in $\mathbb {R}^{n}$ and for finding the coordinates of the nearest point in it (i.e., the footpoint). The cornerstone in this representation is the discriminant of the univariate polynomial, and we exemplify this concept for the particular case of a quartic polynomial. In Sect. 4, the general result is applied to construction of the distance equation for the Hopf quadric. Some illustrative examples are presented in Sect. 5.

Accuracy Although all the approximate computations have been performed within the accuracy $10^{-40}$, the final results are rounded to $10^{-4}$.

2 The Poincaré–Lyapunov constant as a quadric

An oft-studied case of Eq. (1) is the two-dimensional system

$$\begin{aligned}{}\begin{array}[c]{ll} {\dot{x}}_{1} &{} {=\ \ \omega y}_{2}{+a} _{20}{y}_{1}^{2}{+a}_{11}{y} _{1}{y}_{2}{+a}_{02}{y}_{2}^{2} \\ &{}\quad {+a}_{30}{y}_{1}^{3}{+a}_{21} {y}_{1}^{2}{y}_{2}{+a}_{12}{y} _{1}{y}_{2}^{2}{+a}_{03}{y}_{2} ^{3}+h.o.t\\ {\dot{x}}_{2} &{} {=-\omega y}_{1}{+b} _{20}{y}_{1}^{2}{+b}_{11}{y} _{1}{y}_{2}{+b}_{02}{y}_{2}^{2} {+b}_{30}{y}_{1}^{3}\\ &{}\quad {+b}_{21} {y}_{1}^{2}{y}_{2}{+b}_{12}{y} _{1}{y}_{2}^{2}{+b}_{03}{y}_{2}^{3}+h.o.t \end{array} \end{aligned}$$

(2)

where h.o.t means terms in $x_{1},$ $x_{2}$ of order higher than 3. We assume that the origin $x_{1}=x_{2}=0$ is an isolated equilibrium point. This is a case of a nonhyperbolic equilibrium point whose stability cannot be determined from linear stability analysis. The so-called (first) Poincaré–Lyapunov constant

$$\begin{aligned} L_{1}= & {} \frac{1}{8\omega }\left[ (a_{20}+a_{02})(b_{20}-b_{02}-a_{11} )\, \right. \nonumber \\&+\left. (b_{20}+b_{02})(a_{20}-a_{02}+b_{11})\right] \nonumber \\&+\frac{1}{8}\left( 3a_{30}+a_{12}+b_{21}+3b_{03}\right) \end{aligned}$$

(3)

determines whether the origin is weakly attracting ($L_{1}<0$) or weakly repelling ($L_{1}>0$). In the context of the Hopf bifurcation analysis, Eq. (2) is a parameter-dependent system at a Hopf bifurcation point. The sign of $L_{1}$ is related to the criticality of the Hopf bifurcation. For $L_{1}<0$, a stable limit cycle for (2) appears around the unstable origin when the coefficients in Eq. (2) are varied (supercritical Hopf bifurcation), while for $L_{1}>0$ an unstable limit cycle appears around the stable origin (subcritical Hopf bifurcation).

Even though there are 15 coefficients in Eq. (2) (the angular frequency $\omega$, 7 a’s and 7 b’s), the Poincaré–Lyapunov constant is a function of only 11 of these coefficients, namely $\omega$ and those comprising the ten-dimensional vector

$$\begin{aligned} \mathbf {u}=\left( a_{11},a_{12},a_{20},a_{02},a_{30},b_{11},b_{21} ,b_{20},b_{02},b_{03}\right) ^{\text {T}}. \end{aligned}$$

(4)

Why do we treat $\omega$ and the coefficients in $\mathbf {u}$ separately? In Eq. (2), the angular frequency $\omega$ only appears as a divisor of quadratic terms. With this choice, $L_{1}(\mathbf {u})$ is a quadratic function of the parameter vector $\mathbf {u}$ and thus can be represented as

$$\begin{aligned} L_{1}(\mathbf {u})=\mathbf {u}^{\text {T}}\mathbf {Au}+2\mathbf {b}^{\text {T} }\mathbf {u}, \end{aligned}$$

(5)

where

$$\begin{aligned} \mathbf {A}=\frac{1}{8\omega }\left( \begin{array} [c]{cccccccccc} 0 &{} 0 &{} -\frac{1}{2} &{} -\frac{1}{2} &{} 0 &{} 0 &{} 0 &{} 0 &{} 0 &{} 0\\ 0 &{} 0 &{} 0 &{} 0 &{} 0 &{} 0 &{} 0 &{} 0 &{} 0 &{} 0\\ -\frac{1}{2} &{} 0 &{} 0 &{} 0 &{} 0 &{} 0 &{} 0 &{} 1 &{} 0 &{} 0\\ -\frac{1}{2} &{} 0 &{} 0 &{} 0 &{} 0 &{} 0 &{} 0 &{} 0 &{} -1 &{} 0\\ 0 &{} 0 &{} 0 &{} 0 &{} 0 &{} 0 &{} 0 &{} 0 &{} 0 &{} 0\\ 0 &{} 0 &{} 0 &{} 0 &{} 0 &{} 0 &{} 0 &{} \frac{1}{2} &{} \frac{1}{2} &{} 0\\ 0 &{} 0 &{} 0 &{} 0 &{} 0 &{} 0 &{} 0 &{} 0 &{} 0 &{} 0\\ 0 &{} 0 &{} 1 &{} 0 &{} 0 &{} \frac{1}{2} &{} 0 &{} 0 &{} 0 &{} 0\\ 0 &{} 0 &{} 0 &{} -1 &{} 0 &{} \frac{1}{2} &{} 0 &{} 0 &{} 0 &{} 0\\ 0 &{} 0 &{} 0 &{} 0 &{} 0 &{} 0 &{} 0 &{} 0 &{} 0 &{} 0 \end{array} \right) ,\end{aligned}$$

(6)

$$\begin{aligned} \mathbf {b}=\frac{1}{16}\left( 0,1,0,0,3,0,1,0,0,3\right) ^{\text {T}}. \end{aligned}$$

(7)

      The manifold

$$\begin{aligned} L_{1}(\mathbf {u})=0 \end{aligned}$$

(8)

is a quadric in $\mathbb {R}^{10}$; it will be referred to as the Hopf quadric. The Hopf quadric is degenerate since $\det \mathbf {A}=0$. Its geometry can be established via representation in a canonical form. The eigenvalues of $\mathbf {A}$ are

$$\begin{aligned} \mu _{1,2}=-\frac{1}{8\omega }\sqrt{\frac{3}{2}},\ \mu _{3,4}=\frac{1}{8\omega }\sqrt{\frac{3}{2}},\ \mu _{5,\dots ,10}=0. \end{aligned}$$

(9)

Introducing the new column vector $\mathbf {v}$ via the transformation

$$\begin{aligned} \mathbf {u}=\mathbf {S}\mathbf {v} \end{aligned}$$

(10)

with the corresponding orthogonal matrix

$$\begin{aligned} \mathbf {S}=\left( \begin{array} [c]{rrrrrrrrrr} 1/\sqrt{6} &{} 0 &{} 1/\sqrt{6} &{} 0 &{} -\sqrt{2/3} &{} 0 &{} 0 &{} 0 &{} 0 &{} 0\\ 0 &{} 0 &{} 0 &{} 0 &{} 0 &{} 0 &{} 0 &{} 0 &{} 0 &{} 1\\ 1/2 &{} 1/\sqrt{6} &{} -1/2 &{} 1/\sqrt{6} &{} 0 &{} -1/\sqrt{6} &{} 0 &{} 0 &{} 0 &{} 0\\ 1/2 &{} -1/\sqrt{6} &{} -1/2 &{} -1/\sqrt{6} &{} 0 &{} 1/\sqrt{6} &{} 0 &{} 0 &{} 0 &{} 0\\ 0 &{} 0 &{} 0 &{} 0 &{} 0 &{} 0 &{} 0 &{} 0 &{} 1 &{} 0\\ 0 &{} 1/\sqrt{6} &{} 0 &{} 1/\sqrt{6} &{} 0 &{} \sqrt{2/3} &{} 0 &{} 0 &{} 0 &{} 0\\ 0 &{} 0 &{} 0 &{} 0 &{} 0 &{} 0 &{} 0 &{} 1 &{} 0 &{} 0\\ -1/\sqrt{6} &{} -1/2 &{} -1/\sqrt{6} &{} 1/2 &{} -1/\sqrt{6} &{} 0 &{} 0 &{} 0 &{} 0 &{} 0\\ 1/\sqrt{6} &{} -1/2 &{} 1/\sqrt{6} &{} 1/2 &{} 1/\sqrt{6} &{} 0 &{} 0 &{} 0 &{} 0 &{} 0\\ 0 &{} 0 &{} 0 &{} 0 &{} 0 &{} 0 &{} 1 &{} 0 &{} 0 &{} 0 \end{array} \right) \end{aligned}$$

(11)

results in

$$\begin{aligned} \mathbf {v}^{\text {T}}\widetilde{\mathbf {A}}\mathbf {v}+2\widetilde{\mathbf {b} }^{\text {T}}\mathbf {v}=0 \end{aligned}$$

(12)

where

$$\begin{aligned} \widetilde{\mathbf {A}}= & {} \mathbf {S}^{\text {T}}{\mathbf {A}}\mathbf {S}=\frac{1}{8\omega }\sqrt{\frac{3}{2}}\text {diag}\left( -1,-1,1,1,0,0,0,0,0,0\right) ,\nonumber \\ \end{aligned}$$

(13)

$$\begin{aligned} \widetilde{\mathbf {b}}^{\text {T}}= & {} \mathbf {b}^{\text {T}}\mathbf {S}=\frac{1}{16}\left( 0,0,0,0,0,0,3,1,3,1\right) .\nonumber \\ \end{aligned}$$

(14)

Since $\widetilde{\mathbf {b}}$ is in the null space of $\widetilde{\mathbf {A}}$ (i.e. $\widetilde{\mathbf {A}}\widetilde{\mathbf {b}}=\mathbf {0}$), the linear term cannot be eliminated. Therefore, manifold (8) should be treated as a hyperbolic paraboloid in $\mathbb {R}^{10}$.

We are now ready to state the problem: for a given parameter point $\mathbf {u}_{0}\in \mathbb {R}^{10}$, find the distance from $\mathbf {u}_{0}$ to the Hopf quadric (8). This distance can be used as a measure of the “criticality” of the bifurcation. We are also interested in the coordinates of the nearest to $\mathbf {u}_{0}$ bifurcation point $\mathbf {u}^{*}$ in the parameter space.

In the next section, we formulate the main result on computing these metrics for the case of a general quadric in $\mathbb {R}^{n}$.

3 Distance from a point to a quadric

Consider a general quadric in $\mathbb {R}^{n},$ $n\ge 2,$ specified by

$$\begin{aligned} \mathbf {u}^{\text {T}}\mathbf {Au}+2\mathbf {b}^{\text {T}}\mathbf {u}=0, \end{aligned}$$

(15)

where $\mathbf {A}\in \mathbb {R}^{n\times n}$ is a symmetric matrix, $\mathbf {b}\in \mathbb {R}^{n}$ is a column vector. It turns out that the critical values of the squared distance function from a point $\mathbf {u} _{0}\in \mathbb {R}^{n}$ to quadric (15) are the real zeros of a univariate algebraic equation [22, 23]. This algebraic equation can be effectively constructed with the aid of a special algebraic function of the coefficients of a polynomial known as the discriminant.

Formally, the discriminant of a polynomial

$$\begin{aligned} f(\mu )&=A_{0}\mu ^{n}+A_{1}\mu ^{n-1}+\dots +A_{n}\in \mathbb {C}[\mu ],\nonumber \\&\quad A_{0}\ne 0,n\ge 2 \end{aligned}$$

(16)

is defined as the function of its zeros $\mu _{1},\dots ,\mu _{n}$, namely

$$\begin{aligned} \mathcal {D}_{\mu }(f):=A_{0}^{2n-2}\prod _{1\le j<k\le n}(\mu _{k}-\mu _{j} )^{2}\equiv A_0 ^{n-1} \prod _{j=1}^n f^{\prime } (\mu _j), \end{aligned}$$

(17)

Being a symmetric polynomial of the zeros, the discriminant admits polynomial representation in terms of the coefficients of $f(\mu )$:

$$\begin{aligned} \mathcal {D}_{\mu }(f)\equiv \mathfrak {D}(A_{0},A_{1},\dots ,A_{n}). \end{aligned}$$

(18)

For the detailed treatment of this theory, we refer to [18, 19], and here we exemplify this statement via consideration of a quartic polynomial

$$\begin{aligned} f(\mu )=A_{0}\mu ^{4}+A_{1}\mu ^{3}+A_{2}\mu ^{2}+A_{3}\mu +A_{4}\in \mathbb {C} [\mu ],\text { }A_{0}\ne 0 \end{aligned}$$

(19)

(this case will cover our further needs). Its discriminant equals

$$\begin{aligned} \mathcal {D}_{\mu }(f(\mu ))=4\,\mathfrak {I}_{2}^{3}-27\,\mathfrak {I}_{3}^{2}, \end{aligned}$$

(20)

where

$$\begin{aligned}&\mathfrak {I}_{2}:=4\,A_{0}A_{4}-A_{1}A_{3}+\frac{1}{3}A_{2}^{2},\\&\mathfrak {I}_{3}:=-A_{0}A_{3}^{2}-A_{1}^{2}A_{4}+\frac{8}{3}A_{0}A_{2} A_{4}+\frac{1}{3}A_{1}A_{2}A_{3}-\frac{2}{27}A_{2}^{3}. \end{aligned}$$

Due to definition (17), the condition $\mathcal {D}_{\mu }(f(\mu ))=0$ is necessary and sufficient for the coincidence of two zeros for $f(\mu )$, or, in other words, for the existence of a multiple zero for $f(\mu )$.

Theorem 1

(Jacobi) If the polynomial (16) possesses a unique multiple zero $\mu ^{*}$, and its multiplicity equals 2, then the following ratio is valid:

$$ 1 : \mu ^{*} : (\mu ^{*})^{2} : \dots = \frac{\partial \mathfrak {D} }{\partial A_{n}} : \frac{\partial \mathfrak {D}}{\partial A_{n-1}} : \frac{\partial \mathfrak {D}}{\partial A_{n-2}} : \dots $$

From this theorem, it follows that a double zero of a polynomial (in the case of its uniqueness) can be represented as a rational function over $\mathbb {Z}$ in the coefficients of this polynomial. For the quartic polynomial (19), such a formula can be obtained from representation (20 ):

$$\begin{aligned} \mu ^{*}=\frac{2\,A_{1}\mathfrak {I}_{2}^{2}+(3\,A_{1}A_{2}-18\,A_{0} A_{3})\mathfrak {I}_{3}}{(24\,A_{0}A_{2}-9\,A_{1}^{2})\mathfrak {I}_{3} -8\,A_{0}\mathfrak {I}_{2}^{2}}. \end{aligned}$$

(21)

We now turn to the application of this algebraic machinery to our metric problem.

Theorem 2

[22] The square of the distance between quadric (15) and a point $\mathbf {u}_{0}$ outside quadric (15) is the minimal positive zero $z^{*}$ of the distance equation

$$\begin{aligned} \mathcal {F}(z):=\mathcal {D}_{\mu }(\Phi (\mu ,z))=0, \end{aligned}$$

(22)

provided that this zero is not a multiple one. Here,

$$\begin{aligned} \Phi (\mu ,z):=\det \left( \left( \begin{array} [c]{cc} \mathbf {A} &{} \mathbf {b}\\ \mathbf {b}^{\mathrm {T}} &{} 0 \end{array} \right) +\mu \left( \begin{array} [c]{cc} -\mathbf {I} &{} \mathbf {u}_{0}\\ \mathbf {u}_{0}^{\mathrm {T}} &{} z-\mathbf {u}_{0}^{\mathrm {T}}\mathbf {u}_{0} \end{array} \right) \right) \end{aligned}$$

(23)

and $\mathbf {I}$ stands for the nth-order identity matrix. The coordinates of the point in the quadric nearest to $\mathbf {u}_{0}$ are:

$$\begin{aligned} \mathbf {u}^{*}=(\mu ^{*}\mathbf {I}-\mathbf {A})^{-1}(\mathbf {b}+\mu ^{*}\mathbf {u}_{0}). \end{aligned}$$

(24)

Here, $\mu ^{*}$ stands for the location of the multiple zero of the polynomial $\Phi (\mu ,z^{*})$.

Theorems 1 and 2 claim that the coordinates of the nearest to $\mathbf {u}_{0}$ point in the quadric can generically be expressed as rational functions of the squared distance. Nevertheless, one should pay attention to an exceptional case corresponding to the violation of the condition of Theorem 2. Indeed, if the minimal positive zero of the distance equation (22) happens to be a multiple one, then the corresponding point on the quadric might not be unique¹. More troublesome is the case where the minimal positive zero of (22) of an even multiplicity is generated by the pairs of imaginary (complex-conjugate) points in the quadric. This case is exceptional, but it should not be ignored for the case of the Hopf quadric.

4 Distance to the Hopf quadric

Generically, the degree of the distance equation (22) is equal to 2n [23]. However, for the particular case of the Hopf quadric (5), the degree of this equation decreases due to degeneracy of matrix (6). We first present the expression for the polynomial² (23)

$$\begin{aligned} \Phi (\mu ,z)\equiv \frac{1}{2^{20}\omega ^{4}}\mu ^{5}\left( 128\mu ^{2}\omega ^{2}-3\right) \Phi _{1}(\mu ,z). \end{aligned}$$

(25)

Here,

$$\begin{aligned} \Phi _{1}(\mu ,z):=A_{0}\left( z\right) \mu ^{4}+A_{1}\mu ^{3}+A_{2}\left( z\right) \mu ^{2}+A_{3}\mu +A_{4}, \end{aligned}$$

(26)

where

$$\begin{aligned} A_{0}\left( z\right)= & {} 8192\,\omega ^{2}z,\text { }\\ A_{1}= & {} 8192\omega ^{2}L_{1}(\mathbf {u}_{0}),\\ A_{2}\left( z\right)= & {} 32\big [20\omega ^{2}-6z+3(b_{20}+b_{02})^{2} +3(a_{20}+a_{02})^{2}\\&+2(a_{20}-a_{02}+b_{11})^{2} +2(b_{20}-b_{02}+a_{11} )^{2}\big ],\\ A_{3}= & {} -24(3a_{30}+a_{12}+b_{21}+3b_{03}),\text { }\\ A_{4}= & {} -15 , \end{aligned}$$

and $L_{1}(\mathbf {u}_{0})$ is the first Poincaré–Lyapunov constant defined in (3). Utilization of Theorem 2 is subject to an additional condition.

Theorem 3

If the point $\mathbf {u}_{0}\in \mathbb {R}^{10}$ does not lie in either of the eigenspaces of matrix (6) or a sum of any pairs of these eigenspaces, then the square of the distance from $\mathbf {u}_{0}$ to the Hopf quadric (8) equals the minimal positive zero of the degree 5 equation

$$\begin{aligned} \mathcal {F}(z)= & {} \mathcal {D}_{\mu }(\Phi _{1}(\mu ,z))=-3\cdot 2^{36}\omega ^{2}(C_{0}z^{5}\nonumber \\ \quad&+C_{1}z^{4}+C_{2}z^{3}+C_{3}z^{4}+C_{4}z+C_{5})=0.\nonumber \\ \end{aligned}$$

(27)

Coefficients $C_{j}$ are polynomials of degree 2j in $\omega$ and components of the vector $\mathbf {u}_{0}$:

$$\begin{aligned} C_{0}= & {} 12960,\nonumber \\ C_{1}= & {} 216\left[ 800\omega ^{2}-3(3a_{30}+a_{12}+b_{21}+3b_{03})^{2}\right. \nonumber \\&\left. -200((a_{20}+a_{02})^{2}+(b_{20}+b_{02})^{2})\right. \nonumber \\&\left. +80((b_{11} +2a_{20})(2a_{02}-b_{11}) \right. \nonumber \\&\left. +(2b_{20}-a_{11})(a_{11}+2b_{02}))\right] , \end{aligned}$$

(28)

while the expressions for the other coefficients are rather cumbersome, and we note only that

$$\begin{aligned} C_{5}=\omega ^{2}L_{1}^{2}\left( \mathbf {u}_{0}\right) P_{6}(\omega ,\mathbf {u}_{0}) \end{aligned}$$

(29)

where $P_{6}(\omega ,\mathbf {u}_{0})$ is a degree 6 polynomial.

Proof

On elimination of the extraneous factor $\mu ^5$ and numerical factor $1/(2^{20} \omega ^{4})$ from (25), one has to find the discriminant of the product of two polynomials in $\mu$, namely (26) and $128\mu ^{2}\omega ^{2}-3 \equiv 128(\mu -\mu _1)(\mu -\mu _3)$; here, $\mu _{1}$ and $\mu _3$ are the eigenvalues (9) of matrix (6). Utilizing definition (17) of the discriminant, one gets

$$\begin{aligned}&\mathcal {D}_{\mu } \left( (128\mu ^{2}\omega ^{2}-3) \Phi _{1}(\mu ,z) \right) \nonumber \\&\quad \equiv 128^2 (\mu _1-\mu _3)^2 \mathcal {D}_{\mu }(\Phi _{1}(\mu ,z))\Phi _{1}^2(\mu _1,z) \Phi _{1}^2(\mu _3,z) \, .\nonumber \\ \end{aligned}$$

(30)

The first factor in the right-hand side does not vanish. Two last factors happen to be independent of z but still be dependent on the components of the vector $\mathbf {u}_{0}$. For instance, the following identity is valid

$$\begin{aligned} \Phi _{1}(\mu _3,z)\equiv \frac{9}{\omega ^2} \left\{ \left( \mathbf {u}_{0}^{\mathrm {T}} \mathbf {S}_{[3]} \right) ^2+ \left( \mathbf {u}_{0}^{\mathrm {T}} \mathbf {S}_{[4]} \right) ^2 \right\} \end{aligned}$$

where $\mathbf {S}_{[j]}$ stands for the jth column of matrix (11). This value becomes zero iff the vector $\mathbf {u}_{0}$ is normal to the eigenvectors $\mathbf {S}_{[3]}$ and $\mathbf {S}_{[4]}$ of the symmetric matrix (6), or, equivalently, iff $\mathbf {u}_{0}$ lies in the subspace of $\mathbb R^{10}$ spanned by the remained columns of (11). Therefore, for such points of $\mathbb R^{10}$, expression (30), as a polynomial in z, vanishes identically and the distance equation cannot be constructed via the formula of the theorem. Similar statement is valid for the expression $\Phi _{1}^2(\mu _1,z)$, and the vectors $\mathbf {u}_{0}$ normal to the eigenvectors $\mathbf {S}_{[1]}$ and $\mathbf {S}_{[2]}$.

If both conditions $\Phi _{1}(\mu _1,z) \ne 0$ and $\Phi _{1}(\mu _3,z) \ne 0$ are satisfied, then the distance equation for the Hopf quadric coincides, up to a numerical factor, with $\mathcal {D}_{\mu }(\Phi _{1}(\mu ,z))$. Since $\deg _{\mu } \Phi _{1}(\mu ,z) = 4$, application of formula (20) together with some simplifications with the aid of Maple results in (27). $\square$

5 Examples

To demonstrate the use of the theorems, we now present some examples.

Example 1

Find the distance from $\mathbf {u}_{0} =(5,-2,4,-6,-3,-7,2,4,-2,7)^{\mathrm {T}}$ to the Hopf quadric.

Solution One has

$$\begin{aligned} \Phi _{1}(\mu ,z)= & {} 8192\mu ^{4}\omega ^{2}z+(12288\omega ^{2}+4096\omega )\mu ^{3}\nonumber \\&+(640\omega ^{2}-192z+1408)\mu ^{2}-288\mu -15,\nonumber \\ \end{aligned}$$

(31)

and expression (27) for the polynomial $\mathcal {F}(z)$ is

$$\begin{aligned}&\frac{\mathcal {F}(z)}{-3\cdot 2^{36}\omega ^{2}}=12960\,z^{5}+3456(50\omega ^{2}-137)z^{4}\nonumber \\&\quad +1152(750\omega ^{4}-4340\omega ^{2}+1980\omega +5442)z^{3} \nonumber \\&\quad +3072(625\omega ^{6}-6050\omega ^{4}-2475\omega ^{3}\nonumber \\&\quad -5608\omega ^{2} -16245\omega -11884)z^{2}\nonumber \\&\quad +512(3125\omega ^{8}-53500\omega ^{6}-103500\omega ^{5}\nonumber \\&\quad +725610\omega ^{4}+212040\omega ^{3}\nonumber \\&\quad +1655540\omega ^{2}+528084\omega +163097)z\nonumber \\&\quad -2048(625\omega ^{6}-6675\omega ^{4}-16200\omega ^{3}\nonumber \\&\quad +111081\omega ^{2} +33372\omega +13189)(3\omega +1)^{2}. \end{aligned}$$

(32)

For different values of $\omega$, the distance equation possesses 1, 3 or 5 real zeros; in any case, all of them are positive. For $\omega =2$, the distance equation possesses the single real zero $z^{*}\approx 6.1994$. Therefore, the distance from $\mathbf {u}_{0}$ to the Hopf quadric equals $\sqrt{z^{*} }\approx 2.4898$. To find the nearest point on the quadric, we first evaluate the multiple zero for the polynomial $\Phi _{1}(\mu ,z_{*})$ (defined by (26)) via formula (21) to get $\mu ^{*}\approx -0.1508$. Then, we apply formula (24):

$$\begin{aligned} \mathbf {u}^{*}\approx (4.3258,-2.4145,3.1624,-6.4152,-4.2435,\nonumber \\ -7.2112,1.5855,4.1837,-3.1646,5.7565)^{\mathrm {T}}.\nonumber \\ \end{aligned}$$

(33)

Check $\Vert \mathbf {u}_{0}-\mathbf {u}^{*}\Vert =\sqrt{z^{*}}$, $L_{1}(\mathbf {u}^{*})=0$ and $\mathbf {u}_{0}-\mathbf {u}^{*}$ is normal to the Hopf quadric at $\mathbf {u}=\mathbf {u}^{*}$, i.e., it is parallel to the gradient to this quadric at $\mathbf {u}^{*}$:

$$\begin{aligned}&\det \left( [\mathbf {u}_{0}-\mathbf {u}^{*},\mathbf {A}\mathbf {u}^{*}+\mathbf {b}]^{\mathrm {T}}[\mathbf {u}_{0}-\mathbf {u}^{*},\mathbf {A} \mathbf {u}^{*}+\mathbf {b}]\right) \nonumber \\&\quad \approx \left| \begin{array} [c]{cc} z^{*} &{} 0.9347\\ 0.9348 &{} 0.1409 \end{array} \right| \approx 0. \end{aligned}$$

(34)

Therefore, $\mathbf {u}^{_{*}}$ is indeed a stationary point for the distance function from point $\mathbf {u}_{0}$ to the points in the quadric.

For $\omega =1/2$, the distance equation has three real zeros, namely

$$\begin{aligned} z^{*}\approx 2.0166,14.0618,14.2527. \end{aligned}$$

(35)

In this case, the distance to the Hopf quadric equals $\sqrt{z^{*}} \approx 1.4201$. The multiple zero for the polynomial $\Phi _{1}(\mu ,z^{*})$ is $\mu ^{*}\approx -0.6831$. The nearest point on the quadric is

$$\begin{aligned} \mathbf {u}^{*}\approx (4.4582,-2.0914,3.3131,-6.2738,-3.2745,\nonumber \\ -7.2066, 1.9085,4.1061,-2.9772,6.7255)^{\mathrm {T}}.\nonumber \\ \end{aligned}$$

(36)

As mentioned in Sect. 3, the case of existence of a multiple zero for the distance equation (22) may cause a problem in evaluating the distance for a general quadric (15). This indeed happens for the Hopf quadric.

Theorem 4

If $\mathbf {u}_{0}$ belongs to any of the eigenspaces of matrix (6), then polynomial $\mathcal {F}(z)$ given by Eq. (27) possesses a multiple zero.

Proof

Let us take $\mathbf {u}_{0}$ in the eigenspace associated with the eigenvalue $\mu _{1}=-\frac{1}{8\omega }\sqrt{\frac{3}{2}}$. This point can be represented as

$$\begin{aligned} \mathbf {u}_{0}=\tau _{1}\mathbf {S}_{[1]}+\tau _{2}\mathbf {S}_{[2]} \end{aligned}$$

(37)

for some real scalars $\tau _{1},\tau _{2},$ where $\mathbf {S}_{[j]}$ stands for the jth column of matrix (11). Substitution of such $\mathbf {u}_{0}$ into the polynomial $\mathcal {F}(z)$ results in

$$\begin{aligned} \mathcal {F}(z)=-(6\,z-3\tau +20\,\omega ^{2})^{2}\mathcal {F}_{1}(z) \end{aligned}$$

(38)

where $\tau :=\tau _{1}^{2}+\tau _{2}^{2}$ and

$$\begin{aligned} \mathcal {F}_{1}(z):= & {} 18\,z^{3}+(-54\tau +120\,\omega ^{2})z^{2}\nonumber \\&+(300\,\tau \omega ^{2}+54\tau ^{2}+200\,\nonumber \\&\omega ^{4})z-3\tau ^{2}(6\tau +5\omega ^{2}). \end{aligned}$$

(39)

Thus, $\mathcal {F}(z)$ indeed has a multiple zero

$$\begin{aligned} z_{1}=\frac{1}{2}\tau -\frac{10}{3}\omega ^{2}. \end{aligned}$$

(40)

If $z_{1}>0$, then it is the minimal positive zero for this polynomial. Indeed,

$$\begin{aligned} \mathcal {F}_{1}(0)<0,\text { }\mathcal {F}_{1}(z_{1})= & {} -\frac{1}{4}\tau (3\tau -80\omega ^{2})^{2}\le 0,\nonumber \\&\mathcal {F}_{1}(+\infty )>0 \end{aligned}$$

(41)

and

$$\begin{aligned} \mathcal {F}_{1}^{\prime }(z)=54(z-\tau )^{2}+240\omega ^{2}z+300\omega ^{2} \tau +200\,\omega ^{4}>0 \end{aligned}$$

(42)

for $z>0$. Therefore, $\mathcal {F}_{1}(z)$ increases as $z\rightarrow +\infty$ and therefore possesses a unique positive zero greater than $z_{1}$.

For $z=z_{1}$, polynomial (26) considered as a polynomial in $\mu$ possesses a multiple zero:

$$\begin{aligned} \Phi _{1}(\mu ,z_{1})&=\frac{1}{6}(16\omega \mu -\sqrt{6})^{2}(96\mu ^{2} \tau \nonumber \\&\quad -640\omega ^{2}\mu ^{2}-80\sqrt{6}\omega \mu -15). \end{aligned}$$

(43)

This zero coincides with the eigenvalue $\mu _{3}$ of matrix (6). However, substitution of this value into polynomial (25) results in its identical vanishing.

Similar arguments are valid for the eigenspaces associated with other eigenvalues (9).

The scenario with the identical vanishing of polynomial (25) which occurs in the proof of the previous theorem might signal the existence of infinite number of points in the Hopf quadric that are “nearest” to $\mathbf {u}_{0}$.

Example 2

Find the distance from

$$\begin{aligned} \mathbf {u}_{0}&= \left( \frac{7}{\sqrt{6}},0,\frac{7}{2}+\frac{5}{\sqrt{6} },\frac{7}{2}-\frac{5}{\sqrt{6}},0,\frac{5}{\sqrt{6}},0, \right. \nonumber \\&\quad \left. -\frac{7}{6}\sqrt{6}-\frac{5}{2},\frac{7}{\sqrt{6}}-\frac{5}{2},0\right) ^{\mathrm {T}} \end{aligned}$$

(44)

to the Hopf quadric specified for $\omega =3$.

Solution Point $\mathbf {u}_{0}$ lies in manifold (37) ($\tau _{1}=7,\tau _{2}=5$). The minimal positive zero of (38) is evaluated via Eq. (40) and equals 7. It is a double zero, and the problem is to discover whether it is generated by the real points of the Hopf quadric or the imaginary ones. To determine this, we slightly disturb the coordinates of point $\mathbf {u}_{0}$ pushing it out from the eigenspace. Since the multiple zero of a polynomial is unstable under perturbations, one may expect its splitting into simple ones. For the disturbance $10^{-10}$ in $a_{20}$, this is indeed what happens and both zeros become imaginary: $z_{1,2}\approx 7.00000005\pm 0.00000008i$. For these values of z, polynomial (26) possesses a multiple imaginary zero $\widetilde{\mu }_{3}$ which is very close to the value $\mu _{3}$ defined by (9), i.e., in our example, to $\approx 0.0510$. Substitution $\widetilde{\mu }_{3}$ into (24) yields the point

$$\begin{aligned}&\mathbf {u}^{*}\approx (1.4287-2.0371i,1.2247,2.7706\nonumber \\&\quad +4.1583i,0.7294+0.8317i,\nonumber \\&3.6742,1.0206+1.6633i,1.2247,-2.6789\nonumber \\&\quad +4.0743i,0.1789,3.6742)^{\mathrm {T}}. \end{aligned}$$

(45)

This point satisfies Eq. (5) and the condition $(\mathbf {u}_{0}-\mathbf {u}^{*})^{\mathrm {T}}(\mathbf {u}_{0}-\mathbf {u} ^{*})\approx z_{1,2}$. The imaginary parts of several entries of $\mathbf {u}^{*}$ are far away from zero. 50 random perturbations of the same magnitude $10^{-10}$ in all the components of the column $\mathbf {u}_{0}$ result in approximations $\mathbf {u}^{*}$ with the same real parts but distinct imaginary ones. This allows one to conclude that the double zero $z_{1}=7$ of polynomial (38) is generated by the infinite number of imaginary, complex conjugate pairs of points $\{\mathbf {u}^{*} ,\overline{\mathbf {u}^{*}}\}$. The true distance is given by the unique real zero of polynomial (39), i.e., it equals $\approx \sqrt{17.2072}\approx 4.1481$.

$\square$

Example 3

The essential advantage of an analytical solution over numerical one is that they provide an opportunity to trace the dynamics of the result under variation in parameters. We will illuminate this statement with the reference to the delayed Liénard equation

$$\begin{aligned} \ddot{x}(t)+f(x(t))\dot{x}+g(x(t-\tau ))=0. \end{aligned}$$

(46)

Expanding (46) in the neighborhood of the null solution up to third order yields

$$\begin{aligned} \left\{ \begin{aligned} \dot{x}(t)&=y(t)-kx(t)+ax^{2}(t)+bx^{3}(t),\\ \dot{y}(t)&=-x(t-\tau )+cx^{2}(t-\tau )+dx^{3}(t-\tau ). \end{aligned} \right. \end{aligned}$$

(47)

It was shown in [24] that the parametric curve $\left\{ \tau \left( \omega \right) ,k\left( \omega \right) \right\}$, $\omega \in (0,1)$, where

$$\begin{aligned} k\left( \omega \right) =\frac{\sqrt{1-\omega ^{4})}}{\omega },\ \tau \left( \omega \right) =\frac{1}{\omega }\arctan \frac{k}{\omega }, \end{aligned}$$

(48)

separates the linearly stable and unstable regions in the $(\tau ,k)$-parameter plane. Zhao and Kalmár-Nagy performed center manifold reduction of the delay equation (47) on the stability boundary (48), and they derived the Poincaré–Lyapunov constant as [24]

$$\begin{aligned} L_{1}(\tau ,k, \omega )=l_{1}d_{12}+l_{2}d_{22} \, , \end{aligned}$$

(49)

where

$$\begin{aligned} \left\{ \begin{aligned} l_{1}&=\frac{3}{8}d\omega ^{2}+\frac{a^{2}}{4}\omega \zeta (1+4\omega ^{2}-2\omega ^{4})\\&\quad -\frac{ac}{2}k\omega \zeta (1+\omega ^{2}+\omega ^{4})\\&~~~+\frac{c^{2}}{4}\omega \zeta (\frac{11}{2}+k^{2}+2\omega ^{2}+12\omega ^{4}-12\omega ^{6}),\\ l_{2}&=\frac{3}{8}b\omega -\frac{3}{8}dk\omega +\frac{a^{2}}{2}k\omega ^{2}\zeta (1-\omega ^{2})\\&\quad +\frac{ac}{4}\zeta (\frac{7}{2}+\omega ^{2}+10\omega ^{4}-10\omega ^{6})\\&~~~+\frac{c^{2}}{4}k\zeta (-\frac{11}{2}+\omega ^{2}-12\omega ^{4} +12\omega ^{6}),\\ d_{12}&=\frac{2\left( k-\omega ^{2}\tau \right) }{(k-\tau \omega ^{2} )^{2}+(2\omega +k\tau \omega )^{2}},\\ d_{22}&=\frac{2\left( 2\omega +k\tau \omega \right) }{(k-\tau \omega ^{2} )^{2}+(2\omega +k\tau \omega )^{2}},\\ \zeta&=\frac{2\omega }{5+12\omega ^{4}-8\omega ^{6}}. \end{aligned} \right. \end{aligned}$$

(50)

We consider the case where the dependency of distance from $\mathbf {u}_{0}$ depends on three parameters, namely $\tau ,k$ and $\omega \in (0,1)$. Here, the coefficients of (2) contain the parameters in a complicated manner, like, for instance:

$$\begin{aligned} a_{20}= \frac{4a\omega ^2(4+2k\tau )+2c\omega ^4(k-\tau \omega ^2)}{(k-\tau \omega ^2)^2+(2\omega +k\tau \omega )^2} \, . \end{aligned}$$

Fortunately, this dependence involves only rational functions, and this permits one to reduce the distance-to-bifurcation problem to solving an algebraic equation. For instance, if $a=1,b=-2,c=5,d=-2,\tau =3/2,k=8/10$, this distance equals $\approx 0.0015$ for $\omega = 68/100$ and $\approx 4.3515$ for $\omega =99/100$.

This example gives rise to an alternative version of the problem. Instead of its treatment in the $\mathbb R^{10}$ space of vectors (4), it is of practical interest to deal with the space of actual parameters, namely $\tau , k$ and $\omega$.

For any choice of $\omega$, the bifurcation curve $L_1=0$ can be represented as an algebraic one. For example, for $a=1,b=-2,c=5,d=-2,\omega =1/2$, one has:

$$\begin{aligned} \widetilde{L}_1(k,\tau )&:=400\,k^3-2284\,k^2\tau -210 k^2+255k\tau \nonumber \\&\quad -1848\, k-630\tau +405 =0 \, . \end{aligned}$$

(51)

To find the distance to this algebraic surface from the point $(\tau _0,k_0)$, one can extend the ideology of distance equation suggested in Sect. 3 for a quadric. This time the curve is of the order 3, but it is linear in $\tau$:

$$\begin{aligned} \widetilde{L}_1(k,\tau )\equiv A_0(k)\tau + A_1(k) \quad \text{ with } \ \{ A_0(k),A_1(k)\} \subset \mathbb R[k] \, . \end{aligned}$$

For the point $(\tau ,k)= (0,0)$, the distance equation can then be constructed as

$$\begin{aligned} \mathcal F(z):= \mathcal D_{k} ([A_0(k)]^2(z-k^2)- [A_1(k)]^2)=0 \end{aligned}$$

with $\mathcal D$ standing for discriminant (17). Polynomial $\mathcal F(z)$ is of the degree 7 with coefficients of the orders up to $10^{146}$ and with real zeros $z_{*} \approx 0.0415, \approx 1.0609, \approx 1.1426$. Distance from (0, 0) to curve (51) equals $\sqrt{z_{*}} \approx 0.2037$. $\square$

6 Conclusion

The problem of estimating a measure of criticality of the Hopf bifurcation is treated as the distance evaluation between a point and a quadric in a ten-dimensional parameter space. Although the dimension of the parameter space is high, the resulting distance equation is of a sufficiently low degree and its zero set can be easily obtained.

The extension of the suggested approach to other metric bifurcation analysis problems (including those mentioned in Introduction) is possible, at least in principle. As demonstrated in Example 3, the concept of distance equation can be extended for the case of manifolds of arbitrary order in $\mathbb {R}^{n}$ [25]. Unfortunately, the computational difficulties grow tremendously with the increase in the number of variables and order of the manifolds (traditional curse of symbolic algorithms). For instance, in [26], the problem of robust Schur stability of the order 3 matrix family depending on 3 parameters leads to the distance equation of the order 162.