Limit cycle identification in nonlinear polynomial systems

We present a novel formulation, based on the latest advancement in polynomial system solving via linear algebra, for identifying limit cycles in general n -dimensional autonomous nonlinear polynomial systems. The condition for the existence of an algebraic limit cycle is first set up and cast into a Macaulay matrix format whereby polynomials are regarded as coefficient vectors of monomials. This results in a system of polynomial equations whose roots are solved through the null space of another Macaulay matrix. This two-level Macaulay matrix approach relies solely on linear algebra and eigenvalue computation with robust numerical implementation. Furthermore, a state immersion technique further enlarges the scope to cover also non-polynomial (including exponential and logarithmic) limit cycles. Applica-tion examples are given to demonstrate the efficacy of the proposed framework.


Introduction
A limit cycle, which is represented by an isolated closed trajectory in the phase plane, describes a phenomenon of oscillation that is widely observed and studied in various research fields such as electrical circuits and control theory [1], chemistry and medicine [2,3], ecological systems [4], human population [5], etc. Some research focuses on the analysis of properties of limit cycles such as stability and period [6,7]. Besides, a variety of theorems on detecting the existence and the number of limit cycles are also developed by many researchers [8][9][10]. While, beyond these topics, identification of limit cycles remains to be a difficult problem. In previous research, different methods and theories are established for constructing approximated limit cycles (e.g., via harmonic balance, power series) [11,12]. However, little development is achieved in respect of identifying exact analytical formulas of limit cycles for a general system or a family of equations.
In this paper, we propose a new framework for identifying the limit cycles of polynomial systems, inspired by the latest innovation in multivariate polynomial representation and roots finding by formulating the polynomial equations via Macaulay matrices [13,14], thereby turn-ing the problem into an eigenvalue computation problem.
In the proposed framework, an equation of semi-invariant that captures limit cycles is formulated via linear algebraic representation, by assuming that the limit cycle is represented as a multivariate polynomial. Afterwards, a two-level Macaulay matrix approach is applied to solve the set of equations. First, the semi-invariant equation is reformulated by representing multivariate polynomials and their multiplications with coefficient vectors and monomial vectors, based on the concept of Macaulay matrix. A set of polynomial equations are then constructed with respect to the coefficients of the limit cycle, whose roots are found through eigenvalue computation. This framework has a general efficacy for polynomial systems of arbitrary orders and dimensions.
Moreover, our method is applicable not only to systems with polynomial limit cycles, but also to systems with non-polynomial limit cycles containing common terms such as exponential and logarithmic kernels. We achieved this by exploring and flexibly using immersion, a method for eliminating non-polynomial terms, to enlarge the scope of application. Therefore, the proposed framework provides an innovative method for limit cycle identification for polynomial systems.
The remainder of this paper is organised as follows.
The concepts related to limit cycles and polynomials are introduced in Section 2. Then we present in Section 3 a two-level approach based on Macaulay matrix to solve for the limit cycle. In Section 4, we review the concept of immersion and its novel use to extend our framework to non-polynomial limit cycle identification. Section 5 demonstrates the entire procedure with examples. Finally, Section 6 concludes this paper.

Stable Limit Cycles and Semi-Invariants
Given an autonomous system, n n x f x x If one of its trajectories traces out an isolated closed curve, the curve is called the limit cycle of the system. If all trajectories in its neighbourhood spiral towards the limit cycle (outer trajectories shrink and inner trajectories spread) as shown in Figure 1(a), then it is a stable limit cycle.
The LaSalle's theorem in particular [15] can be used to capture stable limit cycles: LaSalle's Theorem: Let be a compact set, positively invariant with respect to system (1). Let be a continuously differentiable function such that in . Let with h, the quotient of divided by , being a polynomial as well. Such that satisfies (2) is studied in a variety of topics, with names such as semi-invariants, second integrals, eigenpolynomials etc. [16]. In this paper, the semi-invariant associated with system (1) plays the role of limit cycle. For (2), notice that results in a special case that . Such is known as a first integral, which expresses a specific kind of stable limit cycles, the neutrally-stable limit cycles. In this case, the system has infinite non-isolated closed trajectories as shown in Figure 1(b).

Representation of Multivariate Polynomials
Any multivariate polynomial can be represented with a coefficient vector multiplied by a vector containing all monomials up to a certain degree. Suppose we have a general bivariate polynomial of degree two, It can be expressed as, Now, consider the multiplication of p and another general polynomial q. Since q can be regarded as a linear combination of monomials, according to the representation above, it can be seen that the coefficient vector of is also a linear combination of coefficient vectors of p multiplied by each monomials in q. By stacking these coefficient vectors, a truncated Macaulay matrix can be constructed, with coefficient vectors of p and its shifted versions as rows. If polynomial of a certain degree, all polynomials and their operations in (2) are represented in the way described in Section 2.2. By equating the coefficient vectors on both sides of (2), a set of polynomial equations are set up whose variables are coefficients of all unknown polynomials in the equation, including those of the limit cycle. Therefore it is sufficient to obtain the accurate expression of the limit cycle by solving those polynomial equations. Generally, a polynomial limit cycle can be of any order. Therefore, with initial trial of a certain order, we iteratively increase the order of the limit cycle. At each iteration the solved coefficients indicate the termination: the recursion stops when a set of valid coefficients of a limit cycle are obtained. We illustrate the idea through a walkthrough example. Given a third-order bivariate polynomial system with a polynomial limit cycle of unknown order, by assuming the limit cycle is of a certain degree, say degree 2, the general form of such limit cycle is

Two-Level Macaulay Matrix Approach
In this section, we present how the two-level Macaulay matrix approach works on finding limit cycles of polynomial form. The entire procedure is divided into two stages. First, with the assumption that the limit cycle is a The corresponding coefficients of each monomial in on both sides must be equal, which generates a set of 15 polynomial equations in 12 variables.
Notice that further investigation based on observation of the numerical simulated limit cycle makes the equations more compact: firstly, the origin does not lie on the limit cycle, which forces 1 nonzero. We might as well set 1 since the limit cycle still holds; moreover, the limit cycle has two intersection points on each of has two roots of 2 . Therefore 4 and 6 are both nonzero. Accordingly, it reduces (4) to a set of 9 equations in 7 variables, as below, Next, these polynomial equations are manipulated by putting coefficients in a Macaulay matrix and unknowns in a monomial vector, such that the problem is converted to a linear algebraic one. Given a polynomial system , the approach at this stage is accomplished in three steps [13,17] x f where for each i f , monomials from degree 0 up to   n M d contains the coefficients of (6). In other word, expressions in (6) can be represented by Whereas K is not directly available, the basis of the nullspace, Z is first computed by SVD. Then K is proposed to be obtained by K ZV  where V is a nonsingular matrix. Let 1 and 2 be two row selection matrices, and S S D  a freely chosen diagonal matrix, such that 1 S KD  . Replacing K with ZV on both sides, where is the pseudo-inverse of  since 1 is not necessarily square. Therefore, V is obtained by computing the eigenvectors of . Once V is available, K is solved with With first element normalized to 1, columns of K represent monomial vectors containing monomials valued at affine roots of the system. Hence, all affine roots are retrieved.
As described above, the two-level Macaulay matrix approach obtains a set of possible coefficients of the limit cycle. This process iterates by increasing the order of the limit cycle and inspecting the validity of the coefficients solved at each iteration, until a valid limit cycle is obtained, if there exists one. For equations in (5), initially of size is constructed. By inspecting , the Macaulay matrix sufficient for solving affine roots is . Through eigenvalue computation, the affine root is obtained as