Integrability analysis of chaotic and hyperchaotic finance systems

In this paper, we consider two chaotic finance models recently studied in the literature. The first one, introduced by Huang and Li, has a form of three first-order nonlinear differential equations $$\begin{aligned} \dot{x}=z+(y-a)x,\,\, \dot{y}=1-b y-x^2,\,\, \dot{z}=-\,x-c z. \end{aligned}$$x˙=z+(y-a)x,y˙=1-by-x2,z˙=-x-cz.The second system, called a hyperchaotic finance model, is defined by $$\begin{aligned} \begin{aligned}&\dot{x}=z+(y-a)x+u, \\&\dot{z}=-\,x-c z, \end{aligned}&\begin{aligned}&\dot{y}=1-b y-x^2,\\&\dot{u}=-\,d xy-k u. \end{aligned} \end{aligned}$$x˙=z+(y-a)x+u,z˙=-x-cz,y˙=1-by-x2,u˙=-dxy-ku.In both models, (a, b, c, d, k) are real positive parameters. In order to present the complexity of these systems Poincaré cross sections, bifurcation diagrams, Lyapunov exponents spectrum and the Kaplan–Yorke dimension have been calculated. Moreover, we show that the Huang–Li system is not integrable in a class of functions meromorphic in variables (x, y, z),  for all real values of parameters (a, b, c), while the hyperchaotic system is not integrable in the case when $$k=c$$k=c and $$\Delta :=1+d(a+d-c)>0$$Δ:=1+d(a+d-c)>0. We give analytic proofs of these facts analyzing properties of the differential Galois groups of variational equations along certain particular solutions. On the other hand, we show that for certain sets of the parameters (a, b, c, d, k), when $$\Delta \le 0$$Δ≤0, the hyperchaotic system possesses a polynomial first integral, which can be used to reduce the dimension of the system by one.


Introduction
Chaos is a very interesting and common phenomena in nature. Having been discovered originally by Lorenz [1] within the context of atmospheric convection, it has found practical applications in various fields of science, in both theoretical and practical point of view. Therefore, for the last five decades researchers have made a great effort for constructing new chaotic models for chaos-needed applications. Among others, let us mention classical ones: the logistic and the Hénon maps [2,3], Chua's circuits [4], the Lorenz-like systems introduced by Chen and Lü in [5,6], and more recent systems, studied by Qi and Li in [7][8][9], the 4D chaotic Duffing system [10], the Chameleon model [11], and many more [12][13][14][15][16][17][18][19].
Among complex dynamical systems exceptional role play hyperchaotic ones. These systems, reported first time by Rössler [20], possess very peculiar prop-erties. For instance, they have at least two positive Lyapunov exponents and thus, the dynamics of hyperchaotic attractor expands in more than one direction. Since hyperchaotic systems possess usually richer dynamics than chaotic ones, researchers have seen its great potential for application in engineering, cryptography and technology to improve the security of communication. Over the last years, hyperchaotic generators have been a topic of active search. Except Rössler's system, many hyperchaotic models have been found, both in theoretical and in practical ways; see, for instance, recent papers [21][22][23][24][25][26][27][28], and references therein.
Early period linear models of economic cycles were criticized for their incompatibility with complex reality. Their dynamics was too regular to assume that they are able to describe properly the complex nature of real economic systems. It was only 1975 when May and Beddington [29,30], as the first, suggested the application of nonlinear and chaotic dynamical systems in economics and financial markets. In contrary to linear and stochastic models, nonlinear deterministic models explain aperiodic business cycles and can give rise to chaotic behaviors. Nowadays, in the era of computers, researchers have found chaotic behaviors in various known models and new ones with very complex dynamic are being proposed. For instance, let us mention the forced van der Pol model [31,32], the Kaldorian model [33,34], various types of the IS-LM model developed by Hick and Hansen [35][36][37], Goodwin's accelerator model [38,39], and recently, the Keynesian model [40], the Cournot-Puu system [41], the Huang-Li model [42,43] and its generalizations to 4D hyperchaotic systems [22,24,28], to cite just a few. We also highly recommend this year special issue, entitled Dynamics Models in Economics and Finance, published in [44].
The presence of parameters is typical for many models of economic processes. For example, in economic growth models, they may represent tools for influencing the economy, while the aim of the analysis is to find such quantities that would lead to the optimal path of growth. However, if the analyzed model has chaotic dynamics, the matter is essentially complicated. The high sensitivity of chaotic system to a change the initial conditions makes it impossible to predict the effects of economic decisions in a long time scale. Therefore, it is important to ask whether there exists any set of values of parameters for which dynamics is regular and a considered system is integrable. It is the significant question from the economical point of view to avoid the chaotic behaviors and make the precise economic prediction possible. In order to get a quick insight into the dynamics of a considered system, we can make a numerical analysis. Nonetheless, numerical methods and techniques, such as Poincaré cross sections, bifurcations diagrams, Lyapunov exponents, power spectra, can be made only for fixed values of parameters describing a system. For other values, results can be completely different and, for particular their sets, the system can be even integrable. However, it is technically and theoretically impossible to make a numerical analysis for all values of parameters. To find integrable cases, one needs a powerful method that enables to distinguish values of parameters for a considered system is suspected to be integrable. The aim of this paper is to present how, in an effective way, we can perform an integrability analysis of complex dynamical systems. For this purpose, we will use differential Galois theory. This approach allows, in an analytical way, to prove the non-integrability of a considered system in a wide class of functions and to distinguish values of parameters for which it is suspected to be integrable. We focus our attention mainly on two examples of chaotic and hyperchaotic finance systems, recently studied in the literature.
The rest of this paper is organized as follows. In Sect. 2, descriptions of systems and their dynamics analysis for specified values of parameters are given. For this purpose, we calculate Poincaré cross sections, bifurcation diagrams, Lyapunov exponents and fractal dimension. The integrability analysis of investigated models by means of the differential Galois approach is made in Sect. 3. In Sect. 4, final conclusions are drawn. In order to make this article self-contained, we have added "Appendix" section with the basic definitions and facts concerning integrability, differential Galois theory and its application to integrability investigations.

The Huang-Li financial model
Recently, the Huang-Li system, known also as a chaotic finance model, has been investigated inten-

Chaotic motion
Quasi-periodic motion Periodic motion sively. It is given by the following equations where (x, y, z) are time-dependent variables and (a, b, c) are real nonnegative parameters. Here, x represents the interest rate, y is the investment demand, and z is the price index. Parameters (a, b, c) denote the saving amount, the cost per investment and the elasticity demand of commercial markets, respectively. The complex behavior of system (1) was first noted by Ma and Chen in 2001 [45,46]. Then many papers were published where the dynamics of this model has been investigated by means of various methods and techniques, such as Lyapunov exponents and bifurcation diagrams [22,47]; synchronizations with linear and nonlinear feedbacks [48,49], adaptive [50], sliding mode [51] and passive [51] control methods [51]; control via linear, speed and time-delay feedbacks [52][53][54].

Poincaré sections and phase portraits
To get an idea about complexity of the system dynamics, we made two Poincaré cross sections, which are shown in Figs. 1 and 3. Points creating presented patterns were obtained as traces of intersections of trajectories calculated numerically with suitably chosen plane of section y = 0. The coordinates on this plane are (x, z). The first Poincaré section, visible in Fig. 1, was made for zeroth values of the parameters. As we can see, even in this simple case, the system reveals very complex dynamics. In fact, for a = b = c = 0 system (1) converts into the famous Nosé-Hoover oscillator [55,56]. Poincaré section visible in Fig. 1 is similar to that for conservative systems, i.e., we can observe the coexistence of periodic, quasi-periodic and chaotic motions. For easier understanding particular orbits corresponding to these different types of motion have been marked at the Poincaré section plane, while their corresponding phase curves, with depicted points of intersections with Poincaré cross-sectional plane, are visible in Fig. 2. In particular, we do not observe appearance of an attractor, which is common in the case of systems with dissipation. This is due to the fact that for zeroth values of the parameters, the divergence of (1) is just y and in the considered example its mean value is almost zero to a high precision of numerical calculations. Thus, it is natural to ask whether system (1) has an invariant measure. However, direct calculations have not yield any results, so let us leave this question as the open problem. Poincaré cross section completely changes when we chose nonzero values of the parameters. Figure 3 corresponds to dynamics of system (1) with a = 0.1, b = 0.01, c = 0.67. Since the divergence is y − a − b − c = 0, volume of the phase space is contracted, and we got the expected shapely elegant strange attractor of the  Fig. 3 showing fractals structure of system (1) fractal structure for dissipative system. Indeed, Fig. 4 shows successive magnifications of the small selected region in Fig. 3, which displays the hidden layers preserving self-similarity.
There is more detailed explanation about the differences between the Poincaré cross sections of conservative and dissipative systems the interested reader can find, for example, in the book [57].

A hyperchaotic finance model
One of generalizations of the Huang-Li system is called a hyperchaotic finance model. This system introduced in [22] is described by a set of four nonlinear differential equations where (a, b, c, d, k) are real nonnegative parameters.
Here, x, y, z have the same meaning as in the Huang-Li model, while u is an additional state variable that represents the average profit margin. Since the system is four dimensional, it is rather difficult to interpret its Poincaré cross sections. Therefore, in order to present complexity of this model one can study its bifurcation diagrams, Lyapunov exponents' spectrum and the Kaplan-Yorke dimension.

Bifurcation diagram
Bifurcation diagrams provide an important view to a system dynamics by plotting the state variable as a function of a bifurcation parameter [58]. More precisely, a bifurcation diagram of a dynamical systeṁ x = v(x, k), where k is a control parameter, illustrates how the limit set (fixed points, periodic trajectories, chaotic attractors) changes with a change of k.  2] affects the dynamics of the system. As we see, the periodic motion diverges and splits as k decreases and finally becomes chaotic. Indeed, Fig. 5b shows magnification of the lower right area of the bifurcation diagram, in which we can observe first successive bifurcations that lead the system to chaotic behavior through the period-doubling route.
In the case of continuous dynamical systems, the period-doubling route starts with limit cycle behavior of the system. As a control parameter is changed, a periodic solution may become unstable and gives birth to a period two cycle. In Fig. 5b, we pointed the values of parameter k (to 4 significant figures) at which the subsequent bifurcations appeared.
The period-doubling mechanism, that is one of the scenarios of transition to chaos, was intensively studied by physicist and mathematician M. J. Feigenbaum.
He found that for many systems, which exhibit bifurcations, the intervals between successive thresholds can be determined from the formula where k n is the value for which n bifurcations appears, and k n+1 is the value for which the next bifurcation takes place. The number is the universal constant called Feigenbaum's number.
Relation (3) implies that difference between successive bifurcations (k n+1 − k n ) tends asymptotically to zero, and hence, value of the control parameter k approaches a finite limit. For our system, this critical value is found to be This means that for k > k c the motion never repeats itself and is chaotic.

Lyapunov exponents' spectrum and the Kaplan-Yorke dimension
Although the period-doubling route to chaos is seen at bifurcation diagram in Fig. 5a, the hyperchaotic nature of this system can be recognized only by their Lyapunov exponents. According to the chaos theory, Lyapunov exponents measure the exponential rates of propagation of nearby trajectories (orbits) in the phase space. For four-dimensional autonomous systems, chaos appears when at least one Lyapunov exponent is positive, whereas hyperchaos is manifested by at least two positive Lyapunov exponents. Figure 6 presents a spectrum of the Lyapunov exponents' of system (2) versus parameter k calculated by Wolf's algorithm [59], over interval k ∈ [0, 2]. At first sight, we observe that the Lyapunov exponents' spectrum corresponds well to the bifurcation diagram. Let λ i denote the Lyapunov exponents sorted in the decreasing order. Then, the dynamics of (2) can be summarized in the following way.
1. For periodic motion: λ 1,2,3 = 0, λ 4 < 0. 2. For quasi-periodic motion: λ 1,2 = 0, λ 3,4 < 0.  While Lyapunov exponents measure sensitivity of a dynamical system to initial conditions, the Kaplan-Yorke dimension (or the Lyapunov dimension) of an attractor measures its complexity [60]. The dimension of an attractor is defined in terms of Lyapunov exponents λ i , by formula Figure 7 shows a global spectrum of the Kaplan-Yorke dimension of system (2), versus parameter k ∈ [0, 2], and its magnification over the interval, where the hyperchaotic behavior takes place. Using quantities (5), or directly from Fig. 7b, we can conclude that a maximum value of the Kaplan-Yorke dimension of the system is D KY ∼ = 3.067. For comparison, the famous hyperchaotic Rössler system has an attractor with D KY ∼ = 3.005, whereas a hyperchaotic attractor found in [61] possesses a fractal dimension D KY ∼ = 3.198.
As we can notice, the hyperchaotic behavior of system (2) is not so noticeable as the one studied in [61], where the maximal Lyapunov exponent has been found to be λ max ∼ = 29.79 and D KY ∼ = 3.198. However, a model studied in the mentioned article was especially created to exhibit very strong hyperchaotic behavior and it has no economic interpretation.

Integrability analysis of the models
The complex behaviors of models (1)-(2) apparent via the Poincaré cross sections, the bifurcation diagram and Lyapunov's exponents spectrum, suggest their nonintegrability. But these numerical evidences of nonintegrability were obtained for chosen values of parameters. For other values, plots can be completely different and, for particular their sets, these systems can be even integrable. However, it is technically and theoretically impossible to make a numerical analysis for all values of parameters. To find integrable cases, one needs a powerful tool that enables us to select these values of parameters for which the considered systems are suspected to be integrable. For Hamiltonian systems for which we have a precise notion of integrability in the sense of Liouville, there are many approaches to the integrability studies: Hamilton-Jacobi theory and separation of variables, perturbation theory, splitting of separatrices and Melnikov's method, Birkhoff's normalization, and recently, Morales-Ramis theory based on the differential Galois approach. In contrast, for non-Hamiltonian systems even there is no commonly accepted definition of integrability. It seems that quite useful is the so-called B-integrability. We say that a given ndimensional systemẋ is , which pairwise commute [u i , u j ] = 0, and k functions F 1 , . . . , F k , which are common first integrals of all vector fields u 1 , . . . , u n−k .
It can be shown that if system (6) is B-integrable, then it is integrable by quadrature. That is, all its solutions can be obtained by means of finite sequence of algebraic operations, superposition and inversion of functions and by calculation of integrals. Basic facts concerning the integrability and in particular a more formal definition of the B-integrability can be found in "Appendix" section.
The aim of this letter is to check whether there exist any sets of values of the parameters for which systems (1) and (2) are B-integrable. The main results are formulated in the following two theorems.

Theorem 1 The Huang-Li nonlinear financial model defined by (1) is not B-integrable in a class of meromorphic functions of variables (x, y, z) for all real values of parameters (a, b, c).
Theorem 2 The hyperchaotic finance system defined by (2) is not B-integrable in a class of meromorphic functions of variables (x, y, z, u) for values of parameters (a, b, c, d), satisfying k = c and To prove these theorems, we investigate variational equations of the considered systems along certain nonstationary solutions and we study their differential Galois groups. This approach of finding necessary conditions for integrability in a framework of differential Galois theory was mostly used for Hamiltonian systems, and it is frequently called Morales-Ramis theory. Thanks to this approach, new integrable cases have been found; see for instance [63][64][65].
If we restrict ourself to the B-integrability, then we have an elegant generalization of Morales-Ramis theory for non-Hamiltonian systems. Namely, with system (6) we can associate its cotangent lift, i.e., a Hamiltonian system defined by the following Hamiltonian function where x = (x 1 , . . . , x n ) and p = ( p 1 , . . . , p n ) are canonical coordinates on a symplectic manifold M = C 2n . In paper [66] Ayoul and Zung showed that if original system (6) is B-integrable, then a lifted system generated by Hamiltonian function (7) is integrable in the Liouville sense. Hence, for both Hamiltonian and non-Hamiltonian systems, we have the same necessary integrability condition, i.e., the identity component of the differential Galois group of variational equations must be Abelian. We summarize the above facts by the following theorem that gives necessary integrability conditions for non-Hamiltonian systems.
Theorem 4 (Ayoul-Zung [66]) Assume that system (6) is meromorphically B-integrable, then the identity component of the differential Galois group of variational equations along a particular non-equilibrium solution ϕ(t) is Abelian.
A more detailed description of differential Galois theory and its application to the integrability studies can be found in "Appendix" section and in cited references therein. In particular, for first applications of differential Galois theory to non-Hamiltonian systems, please consult papers [67,68].

Proof of Theorem 1
System (1) possesses the invariant manifold Indeed, Eq. (1) restricted to N reads as follows: Solving these equations, we obtain our particular solution ϕ(t) = (0, y(t), 0). Let X = [X, Y, Z ] T be the variations of x = [x, y, z] T ; then, the first-order variational equations along ϕ(t) take the form where the explicit form of matrix M(t) is given by and y = y(t) satisfies (9). Since the particular solution corresponds to a motion along y-axis, the equations for X and Z form a subsystem of normal variational equations that can be expressed as one second-order differential equation Next, by means of the change of the independent variable Explicit forms of the coefficients P(s) and Q(s) are the following The classical change of the dependent variable transforms (14) into its reduced form where d := ab − bc − 1.
Here, we should underline one significant fact. Namely, respective transformations (13) and (16) change in general a whole differential Galois group of variational equations (10). However, the clue is that they do not affect the identity component of the group; see, e.g., [62]. Thus, following Theorem 4, in order to prove the non-integrability of our original nonlinear system (1), it is sufficient to show that the identity component of the differential Galois group G of variational equations (10), and thus, its rationalized-reduced form (17) is not Abelian. In order to check whether it is Abelian or not, we usually apply the so-called Kovacic algorithm [69]. This algorithm classifies possible solution forms of second-order differential equations with rational coefficients and their corresponding differential Galois groups.
However, in the considered case, the situation is much simpler. In Eq. (17), we immediately recognize the Whittaker equation with parameters This equation has one regular singularity at z = 0 and one irregular singularity at z = ∞. The differential Galois group of the Whittaker equation is described by the following.

Theorem 5 (Morales-Ruiz [62]) The identity component of the Galois group of equation (18) is Abelian if and only if numbers ( p, q) defined by
Let us assume that the considered system is integrable. Then, according to Theorem 5 the product p · q should always be negative. In our case, however, for (κ, μ) given in (19), we obtain the following equality Thus, p · q > 0 for all real b > 0, and this ends the proof.

Proof of Theorem 2
In order to simplify further calculations, let us write z = v − u, where v is a new variable. System (2) after this linear change of the variable is given by We restrict further analysis to case k = c, in order to find a particular solution in an explicit form. Indeed, in this case system (22) has the two-dimensional invariant manifold and its restriction to N is given bẏ Solution of (24) determines our particular solution ϕ(t)) = (0, y(t), 0, u(t)). The variational equations along ϕ(t) take the following form ⎡ where [X, Y, V, U ] T are the variations of [x, y, v, u] T and y = y(t) satisfies (24). The normal part of this system containing variables (X, V ) can be written as follows Next, we perform change of independent variable (13) that transforms (26) to rational equation (14), with coefficients Then, by means of change of dependent variable (16), we write it in the form of the Whittaker equation (18) with parameters κ and μ defined by Application of Theorem 5 gives the following equality where := 1 + d(a + d − c), has been introduced. From this, we conclude that as long as > 0 the identity component of the differential Galois group of variational equations (26) is not Abelian, and thus, the necessary integrability condition is not satisfied. This ends the proof.

First integrals of the hyperchaotic finance model
In the previous subsection, we have shown that the hyperchaotic finance model is not B-integrable for k = c and = 1 + d(a + d − c) > 0. On the other hand, if such inequality is not satisfied, then there is no integrability obstacles and system (2) can possesses first integrals at least for certain values of parameters  (a, b, c, d, k). So let us look for such first integrals.
Since the right-hand sides of (2) are polynomial in variables (x, y, z, u), we search for a polynomial first integral F(x, y, z, u) ∈ R[(x, y, z, u)] with undetermined coefficients. Postulating the degree of F, computing the Lie derivative L v (F) that is a polynomial of (x, y, z, u) and equaling to zero its all coefficients, we obtain the system of equations for unknown coefficients of the requested first integral F. Its solutions give restrictions for the values of parameters and determine the forms of first integrals which are given below.
In the case when system (2) possesses a linear first integral of the form For the values k = c and it has the cubic first integral Whereas, for k = c and it possesses the quartic first integral It is easy to verify that for all sets of the parameters (30), (32) and (34) with k = c the value of is equal to zero. Thus, the necessary integrability condition is satisfied, as it should be.

Reduction of the hyperchaotic finance model
In the previous subsection, we have shown that for certain sets of parameters (a, b, c, d, k), the hyperchaotic system possesses polynomial first integral. Thus, for prescribed values of the parameters we can use it to reduce the dimension of the system by one. So, let us In order to get a quick insight into the dynamics of this three-dimensional system, we make the bifurcation diagram and the Poincaré cross section, which are shown in Figs. 8 and 9. As we can notice, for parameters b = 0.01, c = 0.71, k = 0.02, f = 0, system (36) still reaches very complex dynamics. Like in the Huang-Li system, we obtain shapely elegant strange attractor of the fractal structure; see especially Fig. 9b showing hidden layers. Indeed, if we put a = − 1/c in (1), and f = k = 0 in (36), then the reduced hyperchaotic model (36) coincides with chaotic one (1). Similar reductions can be made for other sets of the parameters for which there exist cubic (33) and quartic (35) first integrals, respectively. However, its final form will be more complicated because of nonlinear dependence of all variables in first integrals.
It is worth mentioning that system (36) on level f = 0 carries particular solution (9). Making the analogues calculations as in Sects. 3.1-3.2, we conclude that as long as 1 − k/c > 0, system (36) is not B-

Conclusions
It seems that analysis of a differential Galois group of variational equations gives one of the strongest known obstructions to integrability. This method only requires knowledge of a particular solution that is not an equilibrium position. Indeed, in order to find such a solution for the hyperchaotic finance model we put k = c. Even though it has some limitation, it is weaker than many of assumptions required in other methods and the result of application of this method is the non-integrability proof of a considered system in a wide class of functions meromorphic in variables. Furthermore, if a system depends on certain parameters (e.g., as the models investigated in this paper), then we can usually prove its non-integrability for almost all values of parameters except some finite set of values. In this way, new integrable cases can be found.
We showed the application of this approach as exemplified by intensively studied nonlinear models. We proved that the Huang-Li finance model is not integrable in a class of functions meromorphic in variables (x, y, z) for all real values of parameters (a, b, c). In the case of its generalization, called hyperchaotic finance model (2), we proved its non-integrability for k = c and = 1 + d(a + d − c) > 0. Moreover, for ≤ 0, we found linear, cubic and quartic polynomial first integrals for certain sets of parameters (a, b, c, d, k).
Although some values of the parameters have no direct economic interpretation (a < 0), it is very surprising that the four-dimensional hyperchaotic financial model possesses certain first integrals, which allow to reduce it to a three-dimensional chaotic system. Since in the direct search method we restricted ourselves to the polynomial first integrals, it is natural to ask about the existence of Darboux's polynomials, which can be very useful to construct rational first integrals, and thus to further reduction of the hyperchaotic finance model. Let us leave this question as the open problem.

Appendix: Integrability and basics facts from general theory
In this section, we give basic facts concerning the notion of integrability, differential Galois theory and its application to investigations of integrability. A more detailed introduction to this theory can be found, e.g., in [62,70,71].

Integrability
Let us consider a system of differential equations where v(x) = (v 1 (x), . . . , v n (x)) is a smooth vector field in the considered domain U . Evolution of a dynamical system is given by solutions of these equations with certain initial conditions x(t 0 ) = x 0 that can be chosen arbitrarily. However, general solutions are explicitly known only for few systems, such as simple mathematical pendulum, damped harmonic oscillator or standard Kepler problem. As it was observed a long time ago, a significant role in the study of dynamical systems is played by first integrals or other invariant quantities remaining constant on all solutions, which help to find their explicit solutions. Let us recall that a non-constant function F(x) : R n → R is called a first integral of (37) if F(x(t)) = const for all solutions x(t). In other words, F(x) is a conservation law for system (37). The invariance of F(x) with respect to the evolution described by (37) (or, mathematically, with respect to the flow of system (37)) can be expressed for differentiable function in the form where L v is the Lie derivative along vector field v(x). The existence of a first integral implies the relation between coordinates x 1 (0), . . . , x n (0)) = const, which, in principle, allows to eliminate one variable, e.g., x n = ψ(x 1 , . . . , x n−1 ) to reduce a dimension of a considered system by one. If we have F 1 , . . . , F k first integrals, which are functionally independent, then solutions of system (37) lie on n−k-dimensional hyper-surface in R n given by a common level of constant values Assuming, without loss of generality, that we can usẽ x = (x 1 , . . . , x n−k−1 ) as coordinates on the manifold M f 1 ,..., f k , we reduce system (37) to the following onė where i = 1, . . . , n−k−1, and thus x j = x j (x, f ), for j = n − k, . . . , n. In particular, if for n-dimensional vector field we know n − 1 first integrals, which are functionally independent, then we can reduce system (37) to only one first-order equation that can be easily solved by one quadrature and it gives a solution in an implicit form t = γ (x 1 ).
If we assume that a function γ is invertible, then we can determine x 1 = γ −1 (t), and other variables can be determined from the algebraic operations such as addition, subtractions or calculations of roots. In this situation, we say that the considered system is integrable by quadratures.
Suppose that system (37) is a Hamiltonian one. Then, n = 2m and v = v H is a Hamiltonian vector field generated by a smooth scalar function H : R 2m → R, called Hamiltonian. In this context, we consider R 2m as a linear symplectic manifold parametrized by means of canonical coordinates x = (q, p), where q = (q 1 , . . . , q n ) and p = ( p 1 , . . . , p n ). It is equipped with the standard symplectic form and the equations of motion v H generated by H are given by For Hamiltonian systems, we have the well-known notion of integrability in the sense of Liouville. Definition 1 (Liouville) A Hamiltonian system with m degrees of freedom is completely integrable iff it has m first integrals F 1 , . . . , F m , whose Poisson's brackets vanish {F i , F j } = 0, for every i, j = 1, . . . , m, and are functionally independent.
In the case of non-Hamiltonian systems, there is no commonly accepted definition of integrability. Below we give two of them.
Definition 2 (Jacobi) We say that system (37) is integrable in the Jacobi sense iff it admits n−2 functionally independent first integrals and differential m-form such that Equation (45) expresses the fact that ω is an invariant quantity with respect to system (37). In the old literature, the invariant m-form is called the Jacobi last multiplier. This integrability definition is frequently used in non-holonomic mechanics; see, e.g., [68]. Other definition of integrability was proposed by Bogoyavlensky [72,73].

Definition 3 (Bogoyavlensky)
We say that a given ndimensional system (37) is B-integrable if there exist 1 ≤ k ≤ n functionally independent first integrals F 1 , . . . , F k , and n − k symmetries, i.e., vector fields u 1 In this definition, the second condition means that functions F 1 , . . . , F k are common first integrals of all vector fields u 1 , . . . , u n−k . It can be shown that Bintegrability for non-Hamiltonian systems is similar to the integrability in the Liouville sense for Hamiltonian mechanics; see, e.g., [71].

Differential Galois approach to integrability
In order to connect the integrability with properties of differential Galois theory, we have to assume that system (37) is considered over the field of complex functions. More precisely, we assume that the phase space is a complex manifold M = C n and the time is a complex variable, i.e., d dt Most of differential equations describing evolution of physical systems are nonlinear, and thus, it is very difficult to analyze them. However, suppose that we know a particular solution t → ϕ(t), then we can make a lineralization of system (47) around this solution. In other words, we calculate variational equations along ϕ(t), which read as follows: . (48) By means of the variational equations, we investigate the behavior of solutions of nonlinear system (37) around the phase curve corresponding to a particular solution ϕ(t). Since the times of Henri Poincaré, it has been well known that if an original system has an analytic first integral, then variational equations possess a first integral that is a polynomial in variables X. (37), then the first non-vanishing term of its Taylor expansion F(ϕ(t) + X) = f m (X) + · · · , where f m (X) = 0 is a homogeneous polynomial in X of degree m ≥ 0, is a first integral of variational equations (48).

Lemma 1 If F(x) is a first integral of system
Furthermore, it is not difficult to prove that if F(x) is a meromorphic first integral of (47), then the variational equations have a rational and homogeneous in variables X first integral. Moreover, Ziglin's lemma (see [74]) guaranties that if system (47) has k functionally independent and commuting first integrals, then also variational equations (48) have the same number of functionally independent and commuting first integrals. These observations show that integrability of linearized system is related to integrability of an original, nonlinear model. It means that in order to formulate necessary integrability conditions for system (47), we can use its lineralization (48). With linear system (48), we can connect a group. Invariants of this group are related to first integrals of (48). We describe just main ingredients of this relation and its application to integrability study; for more details, see, e.g., books [62,75] and papers [71,76].
We assume that entries of matrix M(t) in equation (48) are elements of a differential field F, e.g., of rational functions F = C(t), or meromorphic functions F = M(t). It is obvious that solutions of (47) are not necessarily elements of F. If we add to the field F, n linearly independent solutions of (48), then the new differential field L is called the Picard-Vessiot extension of F. The differential Galois group G = G(L/F) of linear system (48) is the group of all automorphisms of the Picard-Vessiot extension that leave elements of F fixed and commute with differentiation. It can be shown that G is a linear algebraic subgroup of the group GL(n, C). As a set group G consists of a certain number of separated parts, one of them contains the identity matrix and it is called the identity component G 0 .
If f (g(X)) = f (X) for every element g of the differential Galois group of linear system (48), then we say that f is an invariant of G. The reason why properties of the differential Galois group of variational equations are connected with the integrability is as follows. (47) has k functionally independent meromorphic first integrals in a neighborhood of phase curve , then the differential Galois group G of variational equations (48) along has k functionally independent rational invariants.

Lemma 2 If system
The relation between the presence of first integrals of a nonlinear system and invariants of a certain group related to variational equations was first noticed by Ziglin for monodromy group, see [74]. Later, Morales and Ramis extended this observation to the differential Galois group that contains the monodromy group. For Hamiltonian systems and integrability in the Liouville sense, this relation is particularly elegant. For these systems, G is an algebraic subgroup of the symplectic group Sp(2n, C) and the abelianity of the Lie algebra of first integrals implies abelianity of the Lie algebra of the differential Galois group. This is more precisely stated in the following theorem.
Theorem 6 (Morales-Ramis) Assume that a Hamiltonian system is meromorphically integrable in the Liouville sense in a neighborhood of an analytic phase curve . Then, the identity component of the differential Galois group of the variational equations along is Abelian.
In order to apply this theorem, it is necessary to use effective methods which allow to determine properties of the differential Galois group of linear equations. In most of the investigated systems, variational equations contain a closed smaller s-th dimensional subsystem called normal variational equations and analysis of its differential Galois group is sufficient. Moreover, frequently normal variational equations can be rationalized, i.e., transformed by means of change of independent variable t → f (x(t)), into s-th order differential equation with rational coefficients. In the case when s = 2, the abelianity of the identity component of the differential Galois group can be investigated by means of the Kovacic algorithm [69]. After successful applications of Morales-Ramis theory to various Hamiltonian systems, it was natural to extend its application for non-Hamiltonian systems. As we already mentioned, for such systems, there is no commonly accepted definition of integrability. However, if we restrict ourselves to the B-integrability, then we have an elegant generalization of this theory. Namely, with system (6) we can associate its cotangent lift, i.e., a Hamiltonian system with the following Hamiltonian function where x = (x 1 , . . . , x n ) and p = ( p 1 , . . . , p n ) are canonical coordinates defined in a symplectic manifold M = C 2n . Hence, Hamiltonian vector field is given by In paper [66], Ayoul and Zung showed that if system (47) is B-integrable with k first integrals F 1 , . . . , F k and commuting n − k symmetries u 1 (x) = v(x), u 2 (x), . . . , u n−k (x), then a lifted system generated by function (49) is integrable in the Liouville sense. The above observation was nicely sketched in [71]. Namely, we start with defining the following functions Calculation of the Poisson bracket of these functions with Hamiltonian (49) gives where we used assumption [u j , v] = 0 for 1 ≤ j ≤ n − k. Hence, F j+k are first integrals of system (49). Using the fact that [u j , u i ] = 0, one can also show that F 1 , . . . , F n pairwise commute. Let ϕ(t) = (ϕ 1 (t), . . . , ϕ n (t)) be a particular solution of equations (47). Then t −→ ψ(t) := (ϕ(t), 0) ∈ C 2n (52) is a particular solution of Hamiltonian system (50). The variational equations along this particular solution are given bẏ where M(t) = ∂v ∂ x (ϕ(t)).
Notice that first system of (53) is simply variational equations (48), and the second one is just its adjoint. Hence, if is a fundamental matrix of the systeṁ X = M(t)X, then := ( −1 ) T is a fundamental matrix of the equationsẎ = −M(t) T Y . It implies that the differential Galois group of whole system (53) coincides with the differential Galois group of variational equations (48). Based on this observation, Ayoul and Zung proved the following theorem that gives necessary integrability condition for non-Hamiltonian systems.