First and Second Integrals of Hopf–Langford Type Systems

1. Introduction

The qualitative study of turbulence (spatial complex and temporally “random” fluid motion) using dynamical systems theory was first applied by Hopf [1]. In more detail, he reduced the dynamics of the Navier–Stokes equations for turbulent flow to those on an attractor, and after that studied their dynamics (bifurcations) in the context of nonlinear finite-dimensional flows and maps. Later on, in 1971, Ruelle and Takens [2] suggested that the mechanism of a finite sequence of bifurcations is that underlies turbulence. It is intuitively obvious that, although a low-dimensional model simplifies and significantly restricts the full behaviour of a flow with more modes, the low-dimensional model can provide only qualitative information on key mechanisms and the key features of a complex turbulent flow. Thus, the following questions naturally arise: 1) How well does this low-dimensional model approximate the attractor of the original system? and 2) Can we obtain systems that are as small as possible, but will still yield relevant results? The relative simplicity of this kind of model explains why its use continues to this day. To proceed with these crucially important questions, an interaction between physics, mathematics, and computer science as well as a comprehensive approach to the problem is essential. The earliest turbulence model proposed by Hopf in [1] was modified by Langford first in private communication (see [3] p. 106]) and then in a journal article [4]. This modified model predicts regular oscillations due to a Poincaré–Andronov–Hopf bifurcation [3] and chaotic behaviour for a specific choice of the parameters [5,6,7]. However, addressing this issue requires a further extension of the model to the so-called generalized Hopf–Langford (GHL) systems. Such extended models have been proposed recently in [8,9,10,11] and are being analytically and numerically studied, see [8,9,10,11,12,13,14]. The results obtained within the GHL systems provide cues for the transition from simple to complex oscillatory phenomena. Thus, from a more global perspective, the GHL systems underline the links between similar dynamic phenomena occurring in different settings — physical, chemical and biological.

In the present paper, slightly extending the Hopf–Langford type systems studied in [9,10,11], we consider the cubic polynomial differential system

$$\left\{\begin{array}{c}{\dot{x}}_1=\alpha x_1-\beta x_2+\eta x_1{x}_3+\epsilon x_1{x}_3^2,\hfill \\ {\dot{x}}_2=\beta x_1+\alpha x_2+\eta x_2{x}_3+\epsilon x_2{x}_3^2,\hfill \\ {\dot{x}}_3=\mu x_3-\gamma \left(x_1^2+x_2^2\right)-\delta x_3^2,\hfill \end{array}\right.$$

(1)

where, as usual, the overdots denote differentiation with respect to the independent variable t, and $\mu$, $\alpha$, $\beta$, $\gamma$, $\delta$, $\eta$ and $\epsilon$ are real numbers.

The motivation for the current work is to find first and second integrals of the system (1) as their existence and specific forms would provide a better understanding of the dynamics (local and global bifurcation behaviour) and stability of such a nonlinear dynamical system. Moreover, the Darbouxian theory of integrability provides a link between the integrability of polynomial vector fields and the number of invariant algebraic curves that they have.

The paper is organized as follows. In Section 2, by introducing polar coordinates in the phase space $(x_1,x_2,x_3)$ of the system (1), we show that it reduces to a two-dimensional Liénard system, which corresponds to a second-order Liénard equation. Then, in Section 3, we present first and second integrals of this Liénard system identifying the values of the involved parameters for which these integrals exist. In Section 4, we show that a generic Liénard equation is factorizable if and only if the corresponding Liénard system admits a second integral of a special form. In Section 5, using the established first and second integrals, we give some typical examples of phase portraits and particular solution curves of Liénard systems of the form (12). In Section 6, we briefly discuss the obtained results.

2. Dimensional Reduction of the Extended Hopf–Langford System (1)

It is easy to verify that upon the transformation of the variables of the form

$$x_1=rcos\theta ,x_2=rsin\theta ,x_3=x,$$

i.e., upon introducing polar coordinates in the phase space $(x_1,x_2,x_3)$, the system (1) takes the form

$$\left\{\begin{array}{c}\dot{x}=\mu x-\gamma y-\delta x^2,\hfill \\ \dot{y}=2y\left(\alpha +\eta x+\epsilon x^2\right),\hfill \\ \dot{\theta }=\beta ,\hfill \end{array}\right.$$

(2)

where $y:=r^2$. Thus, the regarded system splits into the equation $\dot{\theta }=\beta$ that gives

$$\theta \left(t\right)=\beta t+{\beta }_0,$$

(3)

where ${\beta }_0$ is an arbitrary real constant, and the two-dimensional amplitude system

$$\left\{\begin{array}{c}\dot{x}=\mu x-\gamma y-\delta x^2,\hfill \\ \dot{y}=2y\left(\alpha +\eta x+\epsilon x^2\right)\hfill \end{array}\right.$$

(4)

whose associated vector field reads

$${\mathcal{X}}_1=\left(\mu x-\gamma y-\delta x^2\right){\displaystyle \frac{\partial }{\partial x}}+2y\left(\alpha +\eta x+\epsilon x^2\right){\displaystyle \frac{\partial }{\partial y}}·$$

(5)

In the case $\gamma =0$, the system (4) separates into the master equation

$$\dot{x}=\mu x-\delta x^2$$

(6)

and slave equation

$$\dot{y}=2y\left(\alpha +\eta x+\epsilon x^2\right).$$

(7)

The general solution of the system (6) & (7) is readily obtained to be

$$x\left(t\right)={\displaystyle \frac{\mu e^{(t+t_0)\mu }}{\delta e^{(t+t_0)\mu }-e^{{\sigma }_1}}},y\left(t\right)={\sigma }_2{\left[e^{\mu {\sigma }_1}-\delta e^{\mu \left(t+t_0\right)}\right]}^{{\displaystyle \frac{\delta \eta +\mu \epsilon }{{\delta }^2}}}{e}^{\left[\alpha \left(t+t_0\right)+{\displaystyle \frac{\mu \epsilon e^{\mu {\sigma }_1}}{{\delta }^2\left(e^{\mu {\sigma }_1}-\delta e^{\mu \left(t+t_0\right)}\right)}}\right]},$$

where $t_0$, ${\sigma }_1$ and ${\sigma }_2$ are arbitrary real constants.

So, hereafter we assume $\gamma \ne 0$.

Now, solving the first one of Equation (4) with respect to the variable y, one obtains

$$y={\displaystyle \frac{1}{\gamma }}\left(\mu x-\delta x^2-\dot{x}\right),\dot{y}={\displaystyle \frac{1}{\gamma }}\left(\mu \dot{x}-2\delta x\dot{x}-\ddot{x}\right).$$

(8)

Next, substituting expressions (8) for the variable y and its derivative $\dot{y}$ in the second one of equations (4), one arrives at the polynomial Liénard equation (see [15])

$$\ddot{x}+F\left(x\right)\dot{x}+G\left(x\right)=0,$$

(9)

with

$$F\left(x\right)=-2\alpha -\mu -2(\eta -\delta )x-2\epsilon x^2,$$

(10)

and

$$G\left(x\right)=2x(\mu -\delta x)(\alpha +\eta x+\epsilon x^2),$$

(11)

in the considered case. Equation (9) is equivalent to the following, defined in the phase plane $(x,y)$, $y=\dot{x}$, differential system

$$\left\{\begin{array}{c}\dot{x}=y,\hfill \\ \dot{y}=-F\left(x\right)y-G\left(x\right),\hfill \end{array}\right.$$

(12)

which is frequently called the Liénard system, with

$${\mathcal{X}}_2=y{\displaystyle \frac{\partial }{\partial x}}-\left[F\left(x\right)y+G\left(x\right)\right]{\displaystyle \frac{\partial }{\partial y}}$$

(13)

as the associated vector field. Notice that in the passage from the equation (9) to the system (12), we kept, for the sake of simplicity, the same notation, ”y”, for the newly introduced dependent variable as those used in the system (4). Actually, under the transformation

$$x\to x,y\to -{\displaystyle \frac{1}{\gamma }}\left(y-\mu x+\delta x^2\right)$$

(14)

the system (4) maps into the system (12) in which the functions $F\left(x\right)$ and $G\left(x\right)$ are given by the equations (10) and (11), respectively. Clearly, the inverse transformation is

$$x\to x,y\to -\gamma y+\mu x-\delta x^2.$$

(15)

3. First and Second Integrals of the Liénard System (12)

To specify the terminology, let us first recall the definitions of the first and second integrals of a real planar polynomial differential system

$$\dot{x}=P(x,y),\dot{y}=Q(x,y),P,Q\in \mathbb{R}[x,y],$$

(16)

where $\mathbb{R}[x,y]$ is the ring of polynomials in x and y with real coefficients.

Definition 1. 

A function $I(x,y)$ is said to be first integral of the system (16) if

$$P(x,y){\displaystyle \frac{\partial I(x,y)}{\partial x}}+Q(x,y){\displaystyle \frac{\partial I(x,y)}{\partial y}}=0.$$

(17)

Let $I(x,y)$ be a first integral of a differential system of the form (16). From the condition (17) it is evident that $\dot{I}(x,y)=0$ on the solutions of this system. In other words, each solution of the regarded system lies on a level set $I(x,y)=c$, determined by certain real constant c, of its first integral. In this sense, the phase portrait of such a system may be considered completely determined.

Definition 2. 

A function $J(x,y)$ is said to be second integral of the system (16) if there exist a polynomial $\Lambda (x,y)$, called the cofactor, such that

$$P(x,y){\displaystyle \frac{\partial J(x,y)}{\partial x}}+Q(x,y){\displaystyle \frac{\partial J(x,y)}{\partial y}}=\Lambda (x,y)J(x,y).$$

(18)

In the context of dynamical systems theory, the second integrals are known in the literature under several different names, see the comments on p. 46 in [16]. Here, they will also be referred to as the Darboux polynomials or invariant algebraic curves of system (16).

According to Poincaré, the basic problem of nonlinear dynamics is studying near-integrable differential systems. In general, the dynamics of higher-dimensional systems is difficult to study. But, if a system has functionally independent first (and second) integrals, then its dynamics can be reduced to that of a lower-dimensional system, which is easier to investigate. That is why the problem of the existence and determination in an explicit form of first and second integrals of a given dynamical system is important and attracts the interest of many researchers, see, e.g., [10,11,16,17,18,19] and the references therein.

Consider now a generic Liénard equation

$${\displaystyle \frac{d^{2}x}{d{t}^2}}+f\left(x\right){\displaystyle \frac{dx}{dt}}+g\left(x\right)=0,$$

(19)

together with the corresponding Liénard system consisting of the two coupled equations

$$\left\{\begin{array}{c}\dot{x}=y,\hfill \\ \dot{y}=-f\left(x\right)y-g\left(x\right),\hfill \end{array}\right.$$

(20)

and its associated vector field

$$\mathcal{X}=y{\displaystyle \frac{\partial }{\partial x}}-\left[f\left(x\right)y+g\left(x\right)\right]{\displaystyle \frac{\partial }{\partial y}}.$$

(21)

As usual, this vector field will also be regarded as an operator acting on the smooth functions defined on some open subset of the phase plain $(x,y)$.

It is well-known (cf., e.g., [20]) that each Liénard system (20) is Hamiltonian, with a first integral

$$I_0(x,y)=y^2+2\int g\left(x\right)dx,$$

(22)

under the condition $f\left(x\right)\equiv 0$. Indeed, if $f\left(x\right)\equiv 0$, then it is easy to check that

$$\mathcal{X}{I}_0(x,y)=0.$$

(23)

In the case of the system (12), we end up with the following result.

Proposition 1. 

Each system of the form (12) in which the parameters involved in the functions $F\left(x\right)$ and $G\left(x\right)$ defined by the equations (10) and (11), respectively, obey the conditions

$$2\alpha +\mu =0,\eta =\delta ,\epsilon =0,$$

(24)

admits a first integral $I_0(x,y)$ of the form

$$I_0(x,y)=\left(y+\mu x-\delta x^2\right)\left(y-\mu x+\delta x^2\right).$$

(25)

Proof. 

In this case we have $f\left(x\right)=F\left(x\right)$ and $g\left(x\right)=G\left(x\right)$. Hence, bearing in mind the Equations (10) and (11), it is clear that the condition $F\left(x\right)\equiv 0$ is fulfilled only for those systems of the form (12) for which the conditions (24) hold. Then, from the Equation (22), for each such system one obtains the first integral (25). □

Another first integral arises when the parameters $\delta$ and $\eta$ involved in a system of the form (12) are not equal and $2\alpha +\mu \ne 0$.

Before presenting this result as well as the results concerning the derivation of second integrals, however, let us make a remark. Given a nonlinear differential system, no universal method exists for constructing its first and second integrals. In the present work, we take some ”promising” polynomials with arbitrary coefficients and try to satisfy the conditions (17) or (18) that in the case of a system of the form (12) leads to an over-determined system of algebraic equations whose solutions (if exist) identify those values of the parameters of the regarded system for which first or second integrals exist.

All the results presented below in the Propositions , and are obtained following this approach.

Proposition 2. 

Each system of the form (12) in which the parameters involved in the functions $F\left(x\right)$ and $G\left(x\right)$ defined by the equations (10) and (11), respectively, obey the conditions

$$2\alpha \delta +\mu \eta =0,\epsilon =0,2\alpha +\mu \ne 0,\delta \ne \eta ,$$

(26)

admits a first integral $I_1(x,y)$ of the form

$$I_1(x,y)=J_1(x,y)J_2{(x,y)}^{-{\displaystyle \frac{\mu }{2\alpha }}}$$

(27)

where

$$J_1(x,y)=y-2\alpha x+2{\displaystyle \frac{\alpha \delta }{\mu }}{x}^2,J_2(x,y)=y-\mu x+\delta x^2.$$

(28)

Proof. 

Taking into account the Equations (10), (11) and (13), one can easily verify by direct computation that

$${\mathcal{X}}_2{I}_1(x,y)={\mathcal{X}}_2\left[J_1(x,y)J_2{(x,y)}^{-{\displaystyle \frac{\mu }{2\alpha }}}\right]=0$$

(29)

provided that relations (26) hold and hence, according to Definition 1, $I_1(x,y)$ is a first integral of the corresponding system of the form (12). □

Actually, $J_1(x,y)$ and $J_2(x,y)$ are second integrals (Darboux polynomials) of any system of the form (12) for which the conditions (26) are satisfied. Indeed, if this is the case, then it is easy to check that

$${\mathcal{X}}_2{J}_1=\left(\mu -2\delta x\right)J_1,{\mathcal{X}}_2{J}_2=\left(2\alpha +2\eta x+2\epsilon x^2\right)J_2.$$

(30)

Hence, according to Definition 2, $J_1$ and $J_2$ are second integrals of the respective systems of form (12). Notably, the second one of the above relations is valid even if the conditions (26) are not satisfied meaning that any Liénard system that corresponds to a generalized Hopf–Langford system of the form (1) has a second integral.

Proposition 3. 

Each Liénard system of the form (12) admits the function

$$J_2=y-\mu x+\delta x^2$$

(31)

as a second integral of cofactor

$${\Lambda }_2=2\alpha +2\eta x+2\epsilon x^2.$$

(32)

Proof. 

Again, taking into account the Equations (10), (11) and (13), one can easily verify by direct computation that

$${\mathcal{X}}_2{J}_2=J_2{\Lambda }_2.$$

(33)

without any additional assumptions. Hence, according to Definition 2, $J_2$ is a second integral of the system (12) of the cofactor ${\Lambda }_2$ as in (32). □

Below, a family of Liénard systems of the form (12) with $\epsilon \ne 0$ is specified which admits second integrals in addition to that formulated in Proposition 3.

Proposition 4. 

Any Liénard system of the form (12) for which

$$\delta =-{\displaystyle \frac{2\eta \mu }{4\alpha +3\mu }},\epsilon ={\displaystyle \frac{4\alpha {\eta }^2}{{(4\alpha +3\mu )}^2}}$$

(34)

admits the function

$$J_3=y-{\displaystyle \frac{\alpha (4\alpha +3\mu )}{3\eta }}-(2\alpha +\mu )x+\eta \left({\displaystyle \frac{\mu }{4\alpha +3\mu }}-1\right)x^2-{\displaystyle \frac{8\alpha {\eta }^2}{3{(4\alpha +3\mu )}^2}}{x}^3$$

(35)

as a second integral (Darboux polynomial) of cofactor

$${\Lambda }_3={\displaystyle \frac{6\eta \mu }{4\alpha +3\mu }}x.$$

(36)

Proof. 

Using the Equations (10), (11), (13), (34) and (35) it is easy to verify that

$${\mathcal{X}}_2{J}_3=J_0{\Lambda }_3.$$

(37)

Hence, according to Definition 2, $J_3$ is a second integral of the respective system of the form (12) of cofactor ${\Lambda }_3$. □

In the next Section, we will give an insight into the relation between the existence of second integrals admitted by a Liénard system of the form (20) and the factorizability of the corresponding Liénard Equation (19).

4. Factorization of the Liénard Equation

Definition 3. 

An equation of the form (19) is said to be factorizable if there exist tow functions ${\varphi }_1\left(x\right)$ and ${\varphi }_2\left(x\right)$ of the dependent variable x such that

$${\displaystyle \frac{d^{2}x}{d{t}^2}}+f\left(x\right){\displaystyle \frac{dx}{dt}}+g\left(x\right)=\left[{\displaystyle \frac{d}{dt}}+{\varphi }_2\left(x\right)\right]\left[{\displaystyle \frac{d}{dt}}+{\varphi }_1\left(x\right)\right]x.$$

(38)

Proposition 5. 

A Liénard type equation of the form (19) is factorizable, i.e., it can be represented in the form (38) with the aid of certain two functions ${\varphi }_1\left(x\right)$ and ${\varphi }_2\left(x\right)$, if and only if the functions $f\left(x\right)$ and $g\left(x\right)$ can be represented as follows

$$f\left(x\right)={\displaystyle \frac{d{\varphi }_1\left(x\right)}{dx}}x+{\varphi }_1\left(x\right)+{\varphi }_2\left(x\right),g\left(x\right)={\varphi }_1\left(x\right){\varphi }_2\left(x\right)x.$$

(39)

Proof. 

The proof of this assertion is straightforward. It suffices to perform the differentiations in the right-hand side of the Equation (38) and compare the obtained result with the left-hand side of this relation. □

Proposition 6. 

Let a Liénard type equation of the form (19) be factorizable, i.e., there exist two functions ${\varphi }_1\left(x\right)$ and ${\varphi }_2\left(x\right)$ such that it can be represented in the form (38). Then, the functions

$$J=x{\varphi }_1\left(x\right)+y,\Lambda =-{\varphi }_2\left(x\right)$$

(40)

satisfy the relation

$$\mathcal{X}J=\Lambda J$$

(41)

that is, J is a second integral (Darboux polynomial) of the corresponding Liénard system (20) with cofactor Λ.

Proof. 

Let the functions $f\left(x\right)$ and $g\left(x\right)$ be such that the corresponding Liénard equation (19) is factorizable. Then, according to Proposition 5, there exist two functions ${\varphi }_1\left(x\right)$ and ${\varphi }_2\left(x\right)$ such that the representation (39) holds. Now, taking into account the specific form (21) of the operator X and the assumption (40) one obtains

$$\begin{array}{ccc}\hfill \mathcal{X}J-\Lambda J& =& y{\displaystyle \frac{\partial J}{\partial x}}-\left[f\left(x\right)y+g\left(x\right)x\right]{\displaystyle \frac{\partial J}{\partial y}}-\Lambda J\hfill \\ & =& y{\displaystyle \frac{\partial }{\partial x}}\left[x{\varphi }_1\left(x\right)+y\right]-\left[f\left(x\right)y+g\left(x\right)x\right]{\displaystyle \frac{\partial }{\partial y}}\left[x{\varphi }_1\left(x\right)+y\right]+{\varphi }_2\left(x\right)\left[x{\varphi }_1\left(x\right)+y\right]\hfill \\ & =& \left[{\displaystyle \frac{\partial {\varphi }_1\left(x\right)}{\partial x}}x+{\varphi }_1\left(x\right)+{\varphi }_2\left(x\right)-f\left(x\right)\right]y-\left[g\left(x\right)-{\varphi }_1\left(x\right){\varphi }_2\left(x\right)\right]x.\hfill \end{array}$$

Thus, (41) holds whenever the representation (39) is assumed. □

Proposition 7. 

Let the functions $f\left(x\right)$, $g\left(x\right)$ and $J_0\left(x\right)$ be such that the Liénard system (20) admits a second integral (Darboux polynomial) of the form

$$J(x,y)=J_0\left(x\right)+y$$

(42)

with a given cofactor $\Lambda \left(x\right)$. Then, the corresponding Liénard equation (19) is factorizable.

Proof. 

Let a function $J(x,y)$ given by the expression (42) be a Darboux polynomial with cofactor $\Lambda \left(x\right)$ of a Liénard system (20) determined by the functions $f\left(x\right)$ and $g\left(x\right)$, that is

$$y{\displaystyle \frac{\partial J(x,y)}{\partial x}}-\left(f\left(x\right)y+g\left(x\right)\right){\displaystyle \frac{\partial J(x,y)}{\partial y}}-\Lambda \left(x\right)J(x,y)=0.$$

(43)

Then, substituting (42) into (43) one obtains

$$y\left({\displaystyle \frac{\partial J_0\left(x\right)}{\partial x}}-\Lambda \left(x\right)-f\left(x\right)\right)-\Lambda \left(x\right)J_0\left(x\right)-g\left(x\right)=0,$$

(44)

and hence

$$f\left(x\right)={\displaystyle \frac{\partial J_0\left(x\right)}{\partial x}}-\Lambda \left(x\right),g\left(x\right)=-\Lambda \left(x\right)J_0\left(x\right).$$

(45)

It is easy to see now, that any tow functions of the form (45) can be represented in the form (39) taking

$${\varphi }_1\left(x\right)={\displaystyle \frac{1}{x}}{J}_0\left(x\right),{\varphi }_2\left(x\right)=-\Lambda \left(x\right).$$

(46)

Therefore, according to Proposition 5, the corresponding Liénard equation is factorizable since the identity

$${\displaystyle \frac{d^{2}x}{d{t}^2}}+\left({\displaystyle \frac{\partial J_0\left(x\right)}{\partial x}}-\Lambda \left(x\right)\right){\displaystyle \frac{dx}{dt}}+\lambda \left(x\right)J_0\left(x\right)=\left[{\displaystyle \frac{d}{dt}}-\Lambda \left(x\right)\right]\left[{\displaystyle \frac{d}{dt}}+{\displaystyle \frac{1}{x}}{J}_0\left(x\right)\right]x$$

(47)

holds for any two functions $J_0\left(x\right)$ and $\Lambda \left(x\right)$. □

Combining the above propositions, we can formulate the following theorem.

Theorem 1. 

A Liénard type equation of the form (19) is factorizable if and only if the corresponding Liénard system (20) admits a second integral (Darboux polynomial) of the form (42).

It should be noted that the factorizability of Liénard type equations has been studied in a number of works (see, e.g., [21] and the relevant references therein). However, to the best of our knowledge, for the first time here it is related to the existence of second integrals (Darboux polynomials) of the corresponding Liénard systems.

As a consequence of Theorem 1 and Proposition 3, we obtain the following result.

Corollary 1. 

Each Liénard system of the form (12) is factorizable.

5. Phase Portraits and Particular Solutions of Liénard Systems of the form (12)

In this section, we will give several typical examples of phase portraits and particular solutions of Liénard systems of the form (12) using the established first and second integrals.

First, let us determine the equilibrium points of the systems of form (12). In view of the expressions (10) and (11), the equilibrium points in the phase portrait of the system (12) are:

$${\mathbf{O}}_0=(0,0),{\mathbf{O}}_1=\left({\displaystyle \frac{\mu }{\delta }},0\right),$$

(48)

and

$${\mathbf{O}}_2=\left(-{\displaystyle \frac{\alpha }{\eta }},0\right)$$

(49)

if $\epsilon =0$, or ${\mathbf{O}}_0$, ${\mathbf{O}}_1$,

$${\mathbf{O}}_3=\left(-{\displaystyle \frac{1}{2\epsilon }}\left[\eta +\sqrt{{\eta }^2-4\alpha \epsilon }\right],0\right),{\mathbf{O}}_4=\left(-{\displaystyle \frac{1}{2\epsilon }}\left[\eta -\sqrt{{\eta }^2-4\alpha \epsilon }\right],0\right),$$

(50)

if $\epsilon \ne 0$ and ${\eta }^2-4\alpha \epsilon \ge 0$.

Evidently (see, e.g., Theorem 9.3.2 in [22] Ch. 9, p. 410]), the nonlinear system (12) is locally linear in the neighbourhood of each equilibrium point ${\mathbf{O}}_i$$(i=0,1,2,3,4)$ since the right-hand sides of the equations in (12) are polynomials in x and y and hence they are infinitely smooth, i.e., $C^{\infty }$, functions. Therefore, the behaviour of the solutions to the system (12) near its equilibrium points is characterized by the corresponding Jacobian matrix evaluated at those points.

In the case of the system (12), taking into account the Equations (10) and (11), the Jacobian matrix $\mathbf{J}$ reads

$$\mathbf{J}=\left[\begin{array}{cc}0& 1\\ -2\alpha \mu +4(\alpha \delta -\eta \mu )x+6(\delta \eta -\mu \epsilon )x^2+8\delta \epsilon x^3+2y(2\epsilon x-\delta +\eta )& 2\alpha +\mu +2\left(\eta -\delta \right)x+2\epsilon x^2\end{array}\right].$$

(51)

This matrix evaluated at the equilibrium points ${\mathbf{O}}_i$$(i=0,1,2,3,4)$, takes the forms

$${\mathbf{J}}_{{\mathbf{O}}_0}=\left[\begin{array}{cc}0& 1\\ -2\alpha \mu & 2\alpha +\mu \end{array}\right],$$

(52)

$${\mathbf{J}}_{{\mathbf{O}}_1}=\left[\begin{array}{cc}0& 1\\ 2\alpha \mu +2{\displaystyle \frac{{\mu }^2}{{\delta }^2}}\left(\delta \eta +\mu \epsilon \right)& 2\alpha -\mu +2{\displaystyle \frac{\mu }{{\delta }^2}}\left(\delta \eta +\mu \epsilon \right)\end{array}\right],$$

(53)

$${\mathbf{J}}_{{\mathbf{O}}_2}=\left[\begin{array}{cc}0& 1\\ {\displaystyle \frac{2\alpha \left(\alpha \delta +\eta \mu \right)}{\eta }}& {\displaystyle \frac{2\alpha \delta }{\eta }}+\mu \end{array}\right],$$

(54)

and

$${\mathbf{J}}_{{\mathbf{O}}_{3,4}}=\left[\begin{array}{cc}0& 1\\ -{\displaystyle \frac{\eta }{{\epsilon }^2}}\left(\eta \mp \sqrt{{\eta }^2-4\alpha \epsilon }\right)(\delta \eta +\mu \epsilon )+{\displaystyle \frac{2\alpha \delta }{\epsilon }}\left(2\eta \mp \sqrt{{\eta }^2-4\alpha \epsilon }+2{\displaystyle \frac{\mu \epsilon }{\delta }}\right)& {\displaystyle \frac{\delta }{\epsilon }}\left(\eta \mp \sqrt{{\eta }^2-4\alpha \epsilon }+{\displaystyle \frac{\mu \epsilon }{\delta }}\right)\end{array}\right],$$

(55)

respectively. In the latter expression, the sign minus corresponds to ${\mathbf{O}}_3$, while the sign plus fits the equilibrium point ${\mathbf{O}}_4$.

As is well known (see, e.g., [22] Sec. 9]), the local stability or instability and the type of a point of equilibrium is determined by the eigenvalues of the Jacobian matrix corresponding to the system under consideration evaluated at that particular point.

5.1. Case 1

Consider the systems of the form (12) in which the parameters involved in the functions $F\left(x\right)$ and $G\left(x\right)$ defined by the Equations (10) and (11), respectively, obey the conditions (24), i.e., $2\alpha +\mu =0$, $\eta =\delta$, and $\epsilon =0$. By virtue of Proposition 1, each such system admits a first integral $I_0(x,y)$ of the form (25), i.e., the solution curves are determined by the equation

$$y^2-x^2{\left(\mu -\delta x\right)}^2-c_0=0,$$

(56)

where $c_0$ is an arbitrary real constant.

The phase portrait of such a system of the form (12) in which the parameters are fixed to be $\alpha =1$, $\mu =-2$, $\epsilon =0$, $\eta =\delta =2$ is shown in Figure 1 (left panel). In this case:

• using the Formulas (48) and (49), the equilibrium (fixed) points of the system are calculated to be ${\mathbf{O}}_0=(0,0)$, ${\mathbf{O}}_1=(-1,0)$ and ${\mathbf{O}}_2=(-0.5,0)$;
• according to the Formulas (52), (53) and (54), the Jacobian matrices evaluated at these points become

  $${\mathbf{J}}_{{\mathbf{O}}_0}={\mathbf{J}}_{{\mathbf{O}}_1}=\left[\begin{array}{cc}0& 1\\ 4& 0\end{array}\right],{\mathbf{J}}_{{\mathbf{O}}_2}=\left[\begin{array}{cc}0& 1\\ -2& 0\end{array}\right];$$
• the eigenvalues of these matrices are: ${\lambda }_1=2>0$ and ${\lambda }_2=-2<0$ at the points ${\mathbf{O}}_0$ and ${\mathbf{O}}_1$ as well as ${\lambda }_1=i\sqrt{2}$ and ${\lambda }_2=-i\sqrt{2}$ at the point ${\mathbf{O}}_2$.

Therefore, ${\mathbf{O}}_0$ and ${\mathbf{O}}_1$ are unstable saddle points while ${\mathbf{O}}_2$ is a centre – the phase portrait around it consists of a continuum of concentric closed curves. The special trajectories presented by the thick curves passing through ${\mathbf{O}}_0$ and ${\mathbf{O}}_1$ (excluding the fixed points themselves) are separatrices. The remaining trajectories have the separatrices as asymptotes.

5.2. Case 2

Next, consider the system of the form (12) in which the parameters involved in the functions $F\left(x\right)$ and $G\left(x\right)$ defined by the Equations (10) and (11), respectively, are: $\alpha =-0.5$, $\mu =2$, $\delta =1$, $\epsilon =0$, $\eta =0.5$. Evidently, the conditions (26) are satisfied for these values of the parameters. Therefore, by virtue of Proposition 2, this system has a first integral $I_1(x,y)$ of the form (27) and, consequently, its solution curves are determined by the equation

$${\displaystyle \frac{1}{2}}\left(2y+2x-x^2\right){\left(y-2x+x^2\right)}^2-c_1=0,$$

(57)

where $c_1$ is an arbitrary real constant.

The phase portrait of the regarded system is depicted in Figure 1 (right panel). In this case:

• using the Formulas (48) and (49), the equilibrium (fixed) points of the system are calculated to be ${\mathbf{O}}_0=(0,0)$, ${\mathbf{O}}_1=(2,0)$ and ${\mathbf{O}}_2=(1,0)$;
• according to the Formulas (52), (53) and (54), the Jacobian matrices evaluated at these points become

  $${\mathbf{J}}_{{\mathbf{O}}_0}=\left[\begin{array}{cc}0& 1\\ 2& 1\end{array}\right],{\mathbf{J}}_{{\mathbf{O}}_1}=\left[\begin{array}{cc}0& 1\\ 2& -1\end{array}\right],{\mathbf{J}}_{{\mathbf{O}}_2}=\left[\begin{array}{cc}0& 1\\ -1& 0\end{array}\right];$$
• the eigenvalues of these matrices are: ${\lambda }_1=2>0$ and ${\lambda }_2=-1<0$ at the point ${\mathbf{O}}_0$, ${\lambda }_1=1>0$ and ${\lambda }_2=-2<0$ at ${\mathbf{O}}_1$ and ${\lambda }_1=i$ and ${\lambda }_2=-i$ at the point ${\mathbf{O}}_2$.

Therefore, ${\mathbf{O}}_0$ and ${\mathbf{O}}_1$ are unstable saddle points while ${\mathbf{O}}_2$ is a centre. The special trajectories presented by the thick curves passing through ${\mathbf{O}}_0$ and ${\mathbf{O}}_1$ (excluding the fixed points themselves) are separatrices. The remaining trajectories have the separatrices as asymptotes.

It is seen that the two phase portraits shown in Figure 1 are qualitatively equivalent because the two saddle points have a common separatrix (or saddle connection) and between them a centre point exists.

5.3. Case 3

Finally, we will consider tow systems, (S1) and (S2), of the form (12) in which the parameters involved in the functions $F\left(x\right)$ and $G\left(x\right)$ defined by the Equations (10) and (11), respectively, are:

(S1)

$\alpha =0.3$, $\mu =0.1$, $\delta =0.04$, $\epsilon =0.048$, $\eta =-0.3$;

(S2)

$\alpha =-0.2$, $\mu =0.1$, $\delta =-0.12$, $\epsilon =-0.288$, $\eta =-0.3$.

The parameters of both systems (S1) and (S2) meet the conditions (34). Therefore, by virtue of Proposition 3 and Proposition 4, they both have second integrals $J_2(x,y)$ and $J_3(x,y)$ of the form (31) and (35), respectively. Consequently, the zero level sets of these second integrals provide the solution curves of the regarded systems depicted in Figure 2.

For the system (S1):

• using the Formulas (48) and (50), the equilibrium (fixed) points are calculated to be ${\mathbf{O}}_0=(0,0)$, ${\mathbf{O}}_1=(2.5,0)$, ${\mathbf{O}}_3=(1.25,0)$ and ${\mathbf{O}}_3=(5,0)$;
• according to the Formulas (52), (53) and (55), the Jacobian matrices evaluated at these points become

  $${\mathbf{J}}_{{\mathbf{O}}_0}=\left[\begin{array}{cc}0& 1\\ -0.06& 0.7\end{array}\right],{\mathbf{J}}_{{\mathbf{O}}_1}=\left[\begin{array}{cc}0& 1\\ -0.03& -0.4\end{array}\right],$$

  $${\mathbf{J}}_{{\mathbf{O}}_3}=\left[\begin{array}{cc}0& 1\\ 0.0225& 0\end{array}\right],{\mathbf{J}}_{{\mathbf{O}}_4}=\left[\begin{array}{cc}0& 1\\ 0.18& -0.3\end{array}\right];$$
• the eigenvalues of these matrices are: ${\lambda }_1=0.1>0$ and ${\lambda }_2=0.6>0$ at the point ${\mathbf{O}}_0$, ${\lambda }_1=-0.3<0$ and ${\lambda }_2=-0.1<0$ at the point ${\mathbf{O}}_1$, ${\lambda }_1=-0.15<0$ and ${\lambda }_2=0.15>0$ at the point ${\mathbf{O}}_3$ and ${\lambda }_1=-0.6<0$ and ${\lambda }_2=0.3>0$ at the point ${\mathbf{O}}_4$.

Therefore, ${\mathbf{O}}_0$ is an unstable node, ${\mathbf{O}}_1$ is a stable node, while ${\mathbf{O}}_3$ and ${\mathbf{O}}_4$ are unstable saddle points.

For the system (S2):

• using the Formulas (48) and (50), the equilibrium (fixed) points are calculated to be ${\mathbf{O}}_0=(0,0)$ and ${\mathbf{O}}_1=(-0.833,0)$;
• according to the Formulas (52), (53) and (55), the Jacobian matrices evaluated at these points become

  $${\mathbf{J}}_{{\mathbf{O}}_0}=\left[\begin{array}{cc}0& 1\\ 0.04& -0.3\end{array}\right],{\mathbf{J}}_{{\mathbf{O}}_1}=\left[\begin{array}{cc}0& 1\\ -0.03& -0.4\end{array}\right].$$
• the eigenvalues of these matrices are: ${\lambda }_1=-0.4<0$ and ${\lambda }_2=0.1>0$ at the point ${\mathbf{O}}_0$, and ${\lambda }_1=-0.3<0$ and ${\lambda }_2=-0.1<0$ at the point ${\mathbf{O}}_1$.

These results allow us to conclude that ${\mathbf{O}}_0$ is a saddle point, while ${\mathbf{O}}_1$ is a stable node.

6. Concluding Remarks and Discussion

In this work, we have studied the nonlinear, seven-parameter, three-dimensional, cubic polynomial, autonomous differential system (1), which is proposed here as a generalization the classical [4] and recently studied [9,10,11] quadratic polynomial Hopf–Langford type systems. By introducing polar coordinates in its phase space, it has been shown that the regarded system reduces to the two-dimensional Liénard system (12), see Section 2. We have found and presented in explicit form in Section 3 polynomial first and second integrals (Darboux polynomials) of Liénard systems of the form (12). The values of the parameters for which these integrals exist have been identified.

It has been proved in Section 4 (see Theorem 1) that a generic Liénard equation is factorizable if and only if the corresponding Liénard system admits a second integral of the special form (42). It is worth noting that the factorization of Liénard type equations has been studied in a number of works (see, e.g., [21] and the references therein). However, to the best of our knowledge, here, for the first time, the factorizability of the Liénard equation has been related to the existence of second integrals (Darboux polynomials) of the corresponding Liénard system. It has been established (see Corollary 1) that each Liénard system corresponding to a Hopf–Langford system of the considered type admits a second integral of the required form (42) and hence the respective Liénard equation is factorizable.

Finally, in Section 5, using the established first and second integrals, we give some typical examples of phase portraits of Liénard systems of the form (12).