Reduction of Kinetic Equations to Li\'enard-Levinson-Smith Form: Counting Limit Cycles

TL;DR

This paper introduces a unified method to convert various two-variable systems into Liénard-Levinson-Smith form, deriving conditions for limit cycles and using the Krylov-Bogoliubov method to determine their maximum number, demonstrated on multiple models.

Contribution

It provides a novel unified scheme for transforming systems into LLS form and analyzing limit cycles, including conditions and maximum cycle count, applicable to polynomial damping and restoring forces.

Findings

Derived conditions for limit cycles in LLS systems.

Implemented Krylov-Bogoliubov method for maximum limit cycles.

Validated scheme on multiple model systems.

Abstract

We have presented an unified scheme to express a class of system of equations in two variables into a Li\'enard-Levinson-Smith (LLS) oscillator form. We have derived the condition for limit cycle with special reference to Rayleigh and Li\'enard systems for arbitrary polynomial functions of damping and restoring force. Krylov-Boguliubov (K-B) method is implemented to determine the maximum number of limit cycles admissible for a LLS oscillator atleast in the weak damping limit. Scheme is illustrated by a number of model systems with single cycle as well as the multiple cycle cases.

Equations

$$\frac{dx}{dt}$$

$$\frac{dy}{dt}$$

$$x$$

$$y$$

$$\ddot{\xi }$$

$$={\alpha }_1\{a_0+a_{1}L(\xi ,\dot{\xi })+a_{2}K(\xi ,\dot{\xi })+\phi (\xi ,\dot{\xi })\}+{\alpha }_2\{b_0+b_{1}L(\xi ,\dot{\xi })+b_{2}K(\xi ,\dot{\xi })+\varphi (\xi ,\dot{\xi })\},$$

$$\phi (\xi ,\dot{\xi })=\sum _{n,m=0}^{\infty }{\phi }_{nm}{\xi }^n{\dot{\xi }}^{m}and\varphi (\xi ,\dot{\xi })=\sum _{n,m=0}^{\infty }{\varphi }_{nm}{\xi }^n{\dot{\xi }}^m,$$

$$\ddot{\xi }$$

$$+{\alpha }_2{b}_0+{\alpha }_2{b}_1(c_1\xi +c_2\dot{\xi }+c_L)+{\alpha }_2{b}_2(c_3\xi +c_4\dot{\xi }+c_K),i.e.,$$

$$\ddot{\xi }$$

$$\ddot{\xi }+F(\xi ,\dot{\xi })\dot{\xi }+G(\xi )=0,$$

$$F(\xi ,\dot{\xi })$$

$$G(\xi )$$

$$\ddot{\xi }+F_R(\dot{\xi })\dot{\xi }+G_R(\dot{\xi })\xi =0,$$

$$F_R(\dot{\xi })$$

$$\ddot{\xi }+F_L(\xi )\dot{\xi }+G_L(\xi )=0,$$

$$F_L(\xi )$$

$$F(\xi ,\dot{\xi })$$

$$+A_{03}{\dot{\xi }}^2+A_{13}\xi {\dot{\xi }}^2+A_{23}{\xi }^2{\dot{\xi }}^2+A_{33}{\xi }^3{\dot{\xi }}^2],$$

$$G(\xi )$$

$$\ddot{\xi }+\sigma F_{\sigma }(\xi ,\dot{\xi })\dot{\xi }+G(\xi )=0.$$

$$\ddot{Z}(\tau )+ϵh(Z(\tau ),\dot{Z}(\tau ))+Z(\tau )=0,$$

$$h(Z,\dot{Z})=-\left[\{B_{01}+B_{11}Z+B_{21}{Z}^2+B_{31}{Z}^3+B_{02}\omega \dot{Z}+B_{12}Z\omega \dot{Z}+B_{22}{Z}^2\omega \dot{Z}+B_{32}{Z}^3\omega \dot{Z}+B_{03}{\omega }^2{\dot{Z}}^2+B_{13}Z{\omega }^2{\dot{Z}}^2+B_{23}{Z}^2{\omega }^2{\dot{Z}}^2+B_{33}{Z}^3{\omega }^2{\dot{Z}}^2\}\omega \dot{Z}+B_{20}{Z}^2+B_{30}{Z}^3\right],$$

$$\dot{\overline{r}}$$

$$\dot{\overline{\varphi }}$$

$$\ddot{\xi }+\left[(1+a+3b^2)-2b\xi -2bk-3b\dot{\xi }+\xi \dot{\xi }+k\dot{\xi }+{\dot{\xi }}^2\right]\dot{\xi }+(a+b^2)\xi =0;a,b>0,k=b+\frac{b}{a+b^2}$$

$$\ddot{\xi }+\left[-\frac{2a_1\xi }{\alpha }-b+\frac{a_1^2}{{\alpha }^2}+\alpha -\frac{2a_1\dot{\xi }}{{\alpha }^2}+\frac{b\dot{\xi }}{a_1}+\frac{{\dot{\xi }}^2}{{\alpha }^2}+\frac{\xi \dot{\xi }}{\alpha }\right]\dot{\xi }+\frac{a_1^2\xi }{\alpha }=0;a_1,b,\alpha >0$$

$$\ddot{x}-\mu (1-x^2+\alpha x^4-\beta x^6)\dot{x}+x=E\cos \Omega t.$$

$$\ddot{x}-\left({\mu }_1+{\mu }_2\dot{x}+{\mu }_3{\dot{x}}^2+\cdots +{\mu }_{2k}{\dot{x}}^{2k-1}+{\mu }_{2k+1}{\dot{x}}^{2k}\right)\dot{x}+x$$

$$\ddot{z}$$

Full text

1 Introduction

Various open kinetic systemsmurraybio ; murraynld ; epstein ; strogatz ; goldbook ; remickens in physics, chemistry and biology, are generically described by a minimal model of autonomous coupled differential equationsslross ; arnold ; birkhoff ; jjstoker ; smale1967 ; len0 of two variables. They exhibit self sustained oscillation in the form of stable limit cycle in a phase plane in many examples, such as, chemical reactionsepstein ; brusselator2009 ; len4 , biological rhythmsmurraynld ; murraybio ; kaiser83 ; kaiser91 ; k-dsr ; strogatz , vibrations in mechanicalrand2012 , optical system and musical instrumentstrogatz ; dsrrayleigh , to name a few. A Rayleighdsrrayleigh equation in violin string and van der Pol oscillation in electric circuit are the classic examples, in this context. More generally, Liénardperko ; len4 ; strogatz ; lls1 ; lls2 ; lls3 equation underlines the concrete criteria for the existence of at least one limit cycle for a general class of such systems of which van der Pol is a special case of the form $\ddot{x}+f(x)\dot{x}+x=0$ where $f(x)=\epsilon(x^{2}-1)$ and Liénard transformation is $\dot{x}=y-F(x)$ and $\dot{y}=-x$ with $F(x)=\int_{0}^{x}f(\tau)d\tau$. A further generalisation of Liénard equation is the LLS equationlls1 ; lls2 ; lls3 ; remickens , $\ddot{x}+F(x,\dot{x})\dot{x}+G(x)=0$, sometimes called the generalised Liénard equation. Casting a general system of kinetic equations in two variables which describe a variety of scenarios in physical, chemical, bio-chemical and ecological sciences into LLS formremickens is often not straight-forwardremickens ; len4 ; limiso ; strogatz . To this end we have provided a scheme for a wide class of open nonlinear equations, cast in the LLS form so that the later becomes amenable to several techniques in nonlinear dynamics.

Our next objective is to find the nature and the number of limit cycles for a given LLS equation thereby addressing the second part of the Hilbert’s $16^{th}$ problem. The problem of counting limit cycle has a long legacy since Hilbert, Smale and many others and still continues it without complete understandingbirkhoff ; jjstoker ; smale1998 ; perko ; gaiko2008 ; countinglcjkb ; infdampinglcjkb ; lcbounestjkb . Our scheme is based on the KB method of averagingkbbook ; strogatz ; slross ; len0 , a variant of multi-scale perturbation techniquestrogatz ; chen1 ; chen2 ; sarkar2012 to derive amplitude equation with considering the polynomial forms of the nonlinear damping and restoring force functions. We have illustrated our results on a variety of known model systemsperko ; gaiko2008 ; giacomini with single and multiple limit cyclesperko ; gaiko2008 ; giacomini .

2 Reduction of Kinetic Equations to Liénard - Levinson - Smith(LLS) form: Conditions for limit cycle

We consider here a set of two dimensional autonomous kinetic equations for an open system. Our aim is to cast the equations into a form of a variant of LLS oscillatorstrogatz ; len4 ; limiso or LLS oscillatorstrogatz ; len4 ; limiso ; len0 ; remickens which can further be reduced to Rayleigh and Liénard form. Let us begin with the system of autonomous kinetic equations

[TABLE]

where $x(t)$ and $y(t)$ are, for example, field variables or populations of species of chemical, biological or ecological process strogatz ; goldbook ; murraynld ; murraybio with $a_{i},b_{i}$ for $i=0,1,2$ are all real parameters expressed in terms of the appropriate kinetic constants. Let, ($x_{s},y_{s}$) be the fixed point of the system and $f(x,y)$ and $g(x,y)$ are the non-linear functions of $x$ and $y$. A first step is shifting the steady state ($x_{s},y_{s}$) to the origin ($0,0$) with the help of a linear transformation as LLS system is a second order homogeneous ordinary differential equation.

The linear transformation can be chosen by introducing a new pair of variables $(\xi,u)$, both of which are functions of $x$ and $y$ where $\xi=\beta_{0}+\beta_{1}x+\beta_{2}y$ with $\beta_{0}=-(\beta_{1}x_{s}+\beta_{2}y_{s})$ i.e. $\xi=\beta_{1}(x-x_{s})+\beta_{2}(y-y_{s})$ such that $\dot{\xi}=u$. $\beta_{1},\beta_{2}$ are weighted constants such that it makes the new steady state at the origin, $\xi_{s}=0$, $u_{s}=0$. $u$ is expressed as $u=\alpha_{0}+\alpha_{1}x+\alpha_{2}y$, with $\beta_{i}$, $\alpha_{i}$ for $i=0,1,2$ are all real constants which can be expressed in terms of system parameters. From the inverse transformation we can easily obtain the expressions for $x$ and $y$ as given by

[TABLE]

provided that $\alpha_{1}\beta_{2}-\alpha_{2}\beta_{1}\neq 0$. Differentiating again, $\dot{\xi}=u$ with respect to the independent variable $t$ we get,

[TABLE]

where, $L(\xi,\dot{\xi})=c_{1}\xi+c_{2}\dot{\xi}+c_{L}$ and $K(\xi,\dot{\xi})=c_{3}\xi+c_{4}\dot{\xi}+c_{K}$ with $\begin{bmatrix}c_{1}&c_{2}&c_{L}\\ c_{3}&c_{4}&c_{K}\end{bmatrix}$ =$\frac{1}{\alpha_{1}\beta_{2}-\alpha_{2}\beta_{1}}$ $\begin{bmatrix}-\alpha_{2}&\beta_{2}&\alpha_{2}\beta_{0}-\alpha_{0}\beta_{2}\\ \alpha_{1}&-\beta_{1}&\alpha_{0}\beta_{1}-\alpha_{1}\beta_{0}\end{bmatrix}$. The functions $\varphi$ and $\phi$ can be expressed as a power series expansion as,

[TABLE]

with, $\phi(\xi,\dot{\xi})=\mu\varphi(\xi,\dot{\xi})$, as the functions $f$ and $g$ are related through $\mu$ by $g=\mu f$, $\mu\in\mathbb{R}$. So, after putting the above form in equation (3) one can find,

[TABLE]

where, $\alpha_{1}a_{0}+\alpha_{2}b_{0}+(\alpha_{1}+\mu\alpha_{2})\varphi_{00}+(\alpha_{1}a_{1}+\alpha_{2}b_{1})c_{L}+(\alpha_{1}a_{2}+\alpha_{2}b_{2})c_{K}=A_{00}=0$ (by definition of a zero fixed point of $\xi$), $A_{10}=\alpha_{1}(a_{1}c_{1}+a_{2}c_{3})+\alpha_{2}(b_{1}c_{1}+b_{2}c_{3})+(\alpha_{1}+\mu\alpha_{2})\varphi_{10}$, $A_{01}=\alpha_{1}(a_{1}c_{2}+a_{2}c_{4})+\alpha_{2}(b_{1}c_{2}+b_{2}c_{4})+(\alpha_{1}+\mu\alpha_{2})\varphi_{01}$, $A_{n0}=(\alpha_{1}+\mu\alpha_{2})\varphi_{n0}$, $A_{n1}=(\alpha_{1}+\mu\alpha_{2})\varphi_{n1}$ and $A_{nm}=(\alpha_{1}+\mu\alpha_{2})\varphi_{nm}$, where indices follow the values as given in the summation over $m,n\in\mathbb{Z}^{+}$. Finally, the above equation looks like,

[TABLE]

where, the functions $F(\xi,\dot{\xi})$ and $G(\xi)$ are given by

[TABLE]

Equation (6) is a well known equation of generalised Liénard form called LLS equation. The condition for existence of having at least a locally stable limit cycle of the dynamical system is $F(0,0)<0\implies A_{01}>0$. It can be shown from the linear stability analysis that there is a relation between $F(0,0)$ and eigenvalues ($\lambda_{\pm}$) with, $F(0,0)=-2\ Re(\lambda_{\pm})$. For a LLS system, there are six conditions to have a limit cycle are given in lls2 ; lls3 ; remickens ; len4 . Out of these six conditions, the condition $F(0,0)$ plays an important role to have a locally stable or unstable limit cycle for such kind of systemlen4 ; limiso ; len0 ; remickens depending upon the sign of $F(0,0)$ is $<0$ or $>0$, respectively. In particular, two situations may arise:

I: For $A_{nm}=0$, with $n\geq 2,\forall m$ i.e. there be an unique steady state ($\xi_{s}=0$) with restoring force linear in $\xi$, then the above form of (6) looks like

[TABLE]

where,

[TABLE]

which is in the form of generalised Rayleigh oscillatordsrrayleigh , the limit cycle condition modifies to, $F_{R}(0)<0$.

II: For $A_{nm}=0$, with $m\geq 2,\forall n$, which corresponds to Liénard equation with an unique steady state ($\xi_{s}=0$). This is of the form

[TABLE]

where,

[TABLE]

where the limit cycle condition is $F_{L}(0)<0$. We know that, for a Liénard system, the damping force function, $F_{L}(\xi)$ and the restoring force function, $G_{L}(\xi)$ are even and odd functions of $\xi$, respectively.

However, for generalised Liénard or LLS system the odd-even properties of $G(\xi)$ and $F(\xi,\dot{\xi})$ have complex ramificationsremickens for practical systems. Here, we have examined the properties with the help of Krylov-Boguliubov averaging method.

3 Maximum Number of Limit Cycles

We now restrict ourselves to the case of LLS systems where the $F(\xi,\dot{\xi})$ and $G(\xi)$ are the polynomial functions of $\xi$ and $\dot{\xi}$. It is well known that linear functional forms of $F$ and $G$ preclude the existence of limit cycle. This can be readily seen by considering the typical examples, e.g., a Harmonic oscillator or a weakly nonlinear oscillator with a potential $\frac{1}{2}\omega_{0}^{2}x^{2}+\frac{1}{3}\lambda x^{4}$, $0<\lambda<1$ or a Lotka-Volterra model, where one encounters a center. We therefore consider the polynomial form of nonlinear damping function $F(\xi,\dot{\xi})$ and restoring force function $G(\xi)$ for our analysis of limit cycle. In what follows we employ K-B method of averaging to show that the characteristic even/odd powers of polynomials play crucial role in determining the behaviour of the associated amplitude and phase equations.

To begin with we consider some fixed values of $m,n$ of equation (7) to truncate the series at $M,N$, for the highest power of $\dot{\xi}$ and $\xi$, respectively. For explicit structure of a prototypical example of an amplitude equation we choose upto $M=N=3$ for illustration. This includes all possible cases for the even and odd nature of $F(\xi,\dot{\xi})$ and $G(\xi)$, respectively. Then the above form of $F(\xi,\dot{\xi})$ and $G(\xi)$ will be in the following reduced forms,

[TABLE]

Let us take $|F(0,0)|=\sigma\in\mathbb{R}^{+}$, an arbitrary constant with $F(\xi,\dot{\xi})=\sigma F_{\sigma}(\xi,\dot{\xi})$. Then the LLS equation can be rewritten as

[TABLE]

Therefore the final equation takes the form of a non-linear oscillator after rescaling $t$ by $\tau$ taking, $\omega t\rightarrow\tau$ as

[TABLE]

where, $0<\epsilon=\frac{\sigma}{\omega^{2}}\ll 1$, $\omega^{2}=-A_{10}>0$ and $Z(\tau)\equiv\xi(t)$ and $\omega\dot{Z}(\tau)\equiv\dot{\xi}(t)$. Equation (14) is now ready for the treatment using K-B method with

[TABLE]

where $B_{ij}=\frac{A_{ij}}{\sigma}$, $i,j=0,1,2,3$ with $B_{00}=0$ and $B_{01}$ will take the fixed value, -1, 0, or 1 depending upon the nature of the fixed point is stable focus, center/center-type or limit cycle, respectively. Now choosing, $Z(\tau)\approx r(\tau)\hskip 2.84544pt\cos(\tau+\phi(\tau))$ as a solution of eq. (14) we have $\dot{Z}(\tau)\approx-r(\tau)\hskip 2.84544pt\sin(\tau+\phi(\tau))$ with slowly varying radius $r(\tau)=\sqrt{Z^{2}+\dot{Z}^{2}}$ and phase $\phi(\tau)=-\tau+tan^{-1}(-\frac{\dot{Z}}{Z})$. The function $h(Z,\dot{Z})$ contains all the non-linear terms and $\epsilon$ is the non-linearity controlling parameter i.e. one has to satisfy $0<\sigma\ll\omega^{2}$. Then one can obtain $\dot{r}(\tau)=\epsilon h\sin(\tau+\phi(\tau))$ and $\dot{\phi}(\tau)=\frac{\epsilon h}{r(\tau)}\cos(\tau+\phi(\tau))$ i.e. the time derivative of amplitude and phase are of $O(\epsilon)$. So, after taking a running averagestrogatz ; remickens ; slross of a time dependent function $U$ defined as, $\overline{U}(\tau)=\frac{1}{2\pi}\int_{0}^{2\pi}U(s)ds$, one finds, $\dot{\overline{r}}=\langle\epsilon h\sin(\tau+\phi(\tau))\rangle_{\tau}$ and $\dot{\overline{\phi}}=\langle\frac{\epsilon h}{r(\tau)}\cos(\tau+\phi(\tau))\rangle_{\tau}$, which gives,

[TABLE]

Now from a close look at the equation for $\dot{\overline{r}}$, it is apparent that only even elements of $F(\xi,\dot{\xi})$ appears but none of any elements of $G(\xi)$ is present due to the zero averages of $\sin^{\mu}\cos^{\nu}$ terms with $\mu=1$ and $\nu\in\mathbb{Z}$. The non-zero averages arise only when $\mu,\nu$ both are even i.e. $\mu=2\eta_{1},\nu=2\eta_{2};\eta_{1},\eta_{2}\in\mathbb{Z}$. Thus, the effect in $\dot{\overline{r}}$ appears only through the even coefficients of $F(\xi,\dot{\xi})$ i.e. by examining the respective variables in the $\dot{\overline{r}}$ equation, we find that only some even coefficients appear for the first order correction. On the other hand $\dot{\overline{\phi}}$ contains only even coefficients of $F(\xi,\dot{\xi})$ which are not in amplitude equation along with odd coefficients of $G(\xi)$ which shows that only odd $G(\xi)$ plays a role here.

So, from the equation of $\dot{\overline{r}}$, one finds that there exist at most $4$ non-zero values of $\overline{r}$. If out of the four roots every pair appears as conjugate then there are three possibilities. The cases are, (i) two different sets of complex conjugate roots giving an asymptotically stable solution, (ii) one pair of complex conjugate roots and two real roots of equal magnitude with opposite sign implying a limit cycle solution having only one cycle and (iii) either four real roots of equal magnitude with opposite sign having double multiplicity gives a limit cycle solution with only one cycle or two different sets of real roots of equal magnitude with opposite sign, may give limit cycle solution with two different cycles of different radius. The unique zero values of the roots of $\overline{r}$ gives a center or center-typelen3.5 situation. So, in short, the existence of a nonzero real root will provide the radius of the cycle which will be stable or unstable depending on the $-ve$ or $+ve$ sign of $\frac{d\dot{\overline{r}}}{d\overline{r}}$, at $\overline{r}=\overline{r}_{ss}$ and at $\overline{r}_{ss}=0$ $\frac{d\dot{\overline{r}}}{d\overline{r}}>0$ or $<0$ gives the nature of the fixed point.

As an example, for Kaiser modelkaiser83 ; kaiser91 ; k-y2007a ; k-y2007b ; k-y2007c ; k-y2010 ; k-dsr ; k-y2012 , there exist three limit cycles for a certain range of $\alpha,\beta$. So, if we choose the parameters, $\alpha$ and $\beta$ from the three limit cycle zone then there exist six real roots with three different pairs i.e., three different radii exist according to three cycles. But, slightly away from the three limit cycle zone, there will exist only a pair of real roots with the same magnitude and other four will appear as a complex conjugate pairs and together produces only a stable limit cycle.

Note that, to have a stable limit cycle solution, one condition must be satisfied i.e. $F(0,0)<0$. But, it fails to give how many limit cycles the system can admit. According to the root finding algorithm one can guess the maximum number of cycles of a LLS system. The condition $F(0,0)<0$ plays an important role as a check for the existence of atleast one stable limit cycle. But for 2-cycle situations one can have at first the locally unstable limit cycle before locating the outer stable limit cycle and in this situation $F(0,0)>0$.

Based on these considerations we have prepared a table(Table-I) illustrating the possible cases for the maximum number of non-zero real roots or the limit cycles.

Table-I: Maxumum number of limit cycles for LLS system

$N$ $M$ $N+M$ Max. No. of Non-zero Max. No. of

Real Roots (Even) Cycle(s)

Even Even Even $N+M-2=(N-1)+(M-1)$ $\frac{N+M}{2}-1$

Even Odd Odd $N+M-1=(N)+(M-1)$ $\frac{N+M-1}{2}$

Odd Even Odd $N+M-3=(N-2)+(M-1)$ $\frac{N+M-3}{2}$

Odd Odd Even $N+M-2=(N-1)+(M-1)$ $\frac{N+M}{2}-1$

Now if we denote the non-zero real values of $\overline{r}$ as an existence of limit cycles as $\overline{r}$ gives the radius of the cycle where at the same time a pair of conjugate (one $+ve$ and one $-ve$) roots of equal magnitude exists for such kind of LLS systems then out of these two roots, radius will be measured by the magnitude and each distinct magnitude counts the number of cycles. For example, if there exists six roots, say, $(p,-p)$ occurring twice and $(q,-q)$ occurring once then the number of cycles will be $2$ of radius $p$ and $q$, respectively. So, if there are all real roots occurring once, then the number of cycles will be atmost $\frac{N+M-2}{2}$ or $\frac{N+M-1}{2}$ or $\frac{N+M-3}{2}$. For LLS equation with $N,M$ are the maximum power of $\xi$ and $\dot{\xi}$ respectively, we have performed the K-B analysis numerically for $N=10,M=10$. The result is given in table-II. For Rayleigh system with $N=1$, for all $M\geq 1$, maximum number of limit cycle will be $\frac{M-1}{2}$ or $\frac{M-2}{2}$ for odd or even $M$, respectively. For Liénard system with $M=1$ for all $N\geq 1$, the maximum number of limit cycle becomes $\frac{N-1}{2}$ or $\frac{N}{2}$ for odd or even $N$, respectively. The above table is valid for an arbitrary finite polynomials of $F$ and $G$. For the case of arbitrary infinite polynomialcountinglcjkb ; infdampinglcjkb ; lcbounestjkb cases maximum number of limit cycles can be stated for finite truncation.

4 Applications to Some Model Systems

Here we have examined three classes of physical models where the analysis of the maximum number of limit cycles holds. This connection with the general model system is discussed with polynomial damping and restoring force function.

4.1 One-Cycle Cases: van der Pol Oscillator, Simple Glycolytic Oscillator, Modified Brusselator Model

Considering the van der Pol oscillatorstrogatz ; limiso ; slross ; kaiser83 ; epstein ; murraynld ; powerlaw with equation, $\ddot{x}+\epsilon(x^{2}-1)\dot{x}+x=0$ having the weak nonlinearity for $0<\epsilon\ll 1$ produces a locally stable limit cycle with $F(0,0)<0$. So, if we compare with the general table we have $N=2$ and $M=1$. This gives a condition for a unique stable limit cycle.

Next considering the Liénard formlimiso ; len3.5 of simple Glycolytic oscillatorstrogatz ; limiso ; len3.5 ; epstein ; murraynld ; goldbook as,

[TABLE]

, has a unique stable limit cycle with $F(0,0)<0$limiso having $N=1$ and $M=3$. This gives one limit cycle.

Furthermore, considering the Modified Brusselator model having the Liénard formlen4 ; limiso ; goldbook ,

[TABLE]

gives a unique stable limit cycle with $F(0,0)<0$len4 ; limiso , where $N=1$ and $M=3$ again giving rise to the same situation.

4.2 Two-Cycle Cases

We rewrite the Liénard form according to Ref.perko ; giacomini as, $\dot{x(t)}=y(t)-F(x(t)),\dot{y(t)}=-x(t)$, where $F(x(t))$ is an odd polynomial. After taking derivative it takes the following form $\ddot{x}+F^{\prime}(x)\dot{x}+x=0$, where $F^{\prime}(x)=\frac{\partial F(x)}{\partial x}$ now becomes the form of an even polynomial. For $F(x)=a_{1}x+a_{2}x^{2}+a_{3}x^{3}$, it has been shownperko ; giacomini ; lins that the system allows a unique limit cycle if $a_{1}a_{3}<0$, which will be stable if $a_{1}<0$ and unstable if $a_{1}>0$. This corresponds to the table $N=2,M=1$. Further extension by Rychkovrychkov shows that for $F(x)=(a_{1}x+a_{3}x^{3}+a_{5}x^{5})$ the number of limit cycle is atmost two. Numerical simulation corroborates this observation when $F(x)$ is chosen as in Ref.perko ; giacomini ; lins , $F(x)=0.32x^{5}-\frac{4}{3}x^{3}+0.8x$. For this case the inner one is unstable limit cycle as $F(0,0)=0.8>0$ i.e. a stable fixed point but the outer one is a stable limit cycle. Here, as per Table-I we have $N=4$ and $M=1$. The above table thus gives the strategies to find out the number of limit cycles(both stable and unstable) a system can have. On the other hand our analysis by K-B method provides a hint towards a choice of the parameter space for search of real roots of the radial equation.

4.3 Three-Cycles Case: Kaiser Bi-rythmicity Model

Extending van der Pol oscillator model with a nonlinear function of higher order polynomial, Kaiserkaiser83 ; kaiser91 ; k-y2007a ; k-y2007b ; k-y2007c ; k-y2010 ; k-dsr ; k-y2012 has described bi-rythmicity with the nonlinear equation,

[TABLE]

Here, $\alpha,\beta,\mu>0$ tune the non linearity. This is a prototype self-sustained oscillatory system in absence of $E$ and $\Omega$ which are the amplitude and the frequency of the external excitation, respectively. The model exhibits an extremely rich bifurcation behaviour and the system actually produces bi-rhythmicity. It has been emphasized that in the undriven case, the model is a multi-limit cycle oscillator and has three limit cycles, two of them are stable and between the two stable limit cycles there is an unstable one which divides the basins of attraction of the two stable cycles. In presence of $E$ and $\Omega$, the above system exhibits some interesting featureskaiser83 ; kaiser91 ; k-y2007a ; k-y2007b ; k-y2007c ; k-y2010 ; k-y2012 . From table-I one has $N=6$ and $M=1$ with even-odd sub cases, while $E=0$. Thus, there may have 6 roots for the radial equation if $\mu>0$ and $\alpha,\beta$(controlling parameters of the radii) are chosen from three limit cycle zone($\alpha=0.144,\beta=0.005$) and finally, the number of distinct values will be 3 which implies that the system can have atmost three limit cycles (but here it is exactly 3). Further, if we choose $\beta=0$ for the above undriven Kaiser model with $\alpha=0.1$, one can have two limit cycles with radii $\approx 2.35$ and $\approx 3.80$, respectively, of which the smaller one will be stable and the larger one will be unstable.

Note that, if there are odd number of limit cycles, say $l$, then out of the $l-$cycles, $\frac{l+1}{2}$ will be stable limit cycles and the remaining $\frac{l-1}{2}$ will be unstable limit cycles. For example, for the van der Pol oscillator, Glycolytic oscillator, Modified Brusselator model etc. only one limit cycle exists which is stable. For Kaiser model, $l=3$ and one can observe the situations accordingly. So, for odd number of cycles innermost one will be locally stable.

4.4 k-Cycle Cases:

4.4.1 A Model With N=1 and M=2k+1

For counting the number of limit cycles Gaikogaiko2008 has shown, for a Liénard-type system i.e., LLS equation having the form,

[TABLE]

can have atmost $k$ limit cycles if and only if, $\mu_{1}>0$. The resultgaiko2008 correlates with our result. For any value of $k\in\mathbb{Z}$ it fits the odd-odd case of the general table and accordingly, $M$ and $N$ are $2k+1$ and $1$, respectively, and finally the number of cycles will be atmost $\frac{N+M}{2}-1=k$.

4.4.2 A Model With N=2k and M=1

Blows and Lloydblows ; perko have stated that “For the Liénard or LLS system $\dot{x}=y-F(x),\dot{y}=-g(x)$ with $g(x)=x$ and $F(x)=a_{1}x+a_{2}x^{2}+\dots+a_{2k+1}x^{2k+1}$ has at most $k$ local limit cycles and there are coefficients with $a_{1},a_{3},\dots,a_{2k+1}$ altering in sign”. This can be found from the table-I with $N=2k$ and $M=1$ to give the condition of atmost $k$ limit cycles. For example, taking $k=3$ with $F(x)=-\epsilon(72x-\frac{392}{3}x^{3}+\frac{224}{5}x^{5}-\frac{128}{35}x^{7})$ has exactly three limit cycles for sufficiently small $\epsilon\neq 0$ which are circles with radii $1,2$ and $3$. The above statement nicely corresponds to the Theorem-6, pp-260perko .

Table-II: Table for highest degree polynomial $N+M(\oplus)$ for LLS system together with maximum number of distinct conjugate roots$(R)$, with $1\leq N,M\leq 10$

[TABLE]

Counting the number of limit cycles through RG method in first order will give similar result which was done by Das et. al.countinglcjkb ; infdampinglcjkb ; lcbounestjkb for some models. We have verified similar results for (3,3) polynomial cases for ($F$,$G$) functions using RG method which become increasingly very difficult and almost impossible upto (10,10) case than K-B averaging method as tabulated in this work. It is very useful to count the number of limit cycles from the table by just looking at the LLS form. For example, the number of limit cycles of all models in Ref.countinglcjkb ; infdampinglcjkb ; lcbounestjkb along with the models in our work can be estimated from our table. The table-I can also be utilized to prepare a model of a desired number of limit cycles in a systematic way.

5 Conclusions

We have presented a scheme to cast a set of a class of coupled nonlinear equations in two variables into a LLS form. By expressing the nonlinear damping and forcing functions as polynomial we have implemented K-B method of averaging to explore the number of admissible limit cycles of the dynamical systems. Our results can be summarised as follows:

For a LLS system, the number of limit cycles will be atmost $\frac{N+M-2}{2}$ when $N$ and $M$ degree of the polynomials for damping and restoring force both are even or odd. Again, $\frac{N+M-1}{2}$ cycles can be found when $N$ is even and $M$ is odd and finally, $\frac{N+M-3}{2}$ cycles when $N$ is odd and $M$ is even. 2. 2.

For a Liénard system, in particular, the formula of counting the number of limit cycles follows the same with $M=1$ and $N\in\mathbb{Z^{+}}$. Also for the generalised Rayleigh situation there occurs a linear restoring force so that $N$ is $1$ and $M\in\mathbb{Z^{+}}$ . 3. 3.

We have validated our general result with the help of a variety of physical systems with one, two, three upto arbitrary k-cycles. 4. 4.

This method stated in our work can also be utilized to prepare a model of a desired number of limit cycles in a systematic way.

Acknowledgement

Sandip Saha acknowledges RGNF, UGC, India for the partial financial support.

**Compliance with ethical standards

Conflicts of interests**

The authors declare that there is no conflict of interests regarding the publication of this paper.

6 Appendix: Lotka-Volterra System

To obtain the LLS form of Lotka-Volterra Systemepstein ; strogatz ; goldbook , let us set $z=\delta x+\beta y$ then $\dot{z}=\alpha\delta x-\beta\gamma y=u$ $\implies$ $x=\frac{\dot{z}+\gamma z}{(\alpha+\gamma)\delta}$ and $y=\frac{-\dot{z}+\alpha z}{(\alpha+\gamma)\beta}$. After taking $t$ derivative upon $\dot{z}$ one can have,

[TABLE]

The fixed point $(0,0)$ gives a saddle solution which is not of any interest in the present context. Choosing the remaining non-zero fixed point for further investigations, and after taking perturbation $z=\xi+z_{s}$ around the fixed point $z_{s}=\alpha+\gamma=\delta x_{s}+\beta y_{s}\neq 0$, one can get the LLS form with $F(\xi,\dot{\xi})=a_{1}\xi+a_{2}\dot{\xi}$ with $a_{1}=\frac{\alpha-\gamma}{\alpha+\gamma}$ and $a_{2}=-\frac{1}{\alpha+\gamma}$. It is to be noted that $G(\xi)$ contains nonlinearity with $G(\xi)=\omega^{2}\xi+a_{3}\xi^{2}$ where $\omega=\sqrt{\alpha\gamma}=Im(\lambda)$($+ve$ sense) and $a_{3}=\frac{\alpha\gamma}{\alpha+\gamma}$. After introducing a small parameter $\epsilon_{1}$ (say) in the constants, $a_{i},b_{i}$ such that $a_{i}=\epsilon_{1}b_{i},i=1,2,3$ the above equation reduces to $\ddot{\xi}+\epsilon_{1}(b_{1}\xi+b_{2}\dot{\xi})\dot{\xi}+\omega^{2}\xi+\epsilon_{1}b_{3}\xi^{2}=0$.