Reducing nonlinear dynamical systems to canonical forms

1. Preamble

Ilya Prigogine was among the first natural scientists in the twentieth century to draw the attention of the scientific community to the hitherto unsuspected complexity of the dynamics generated by systems of nonlinear ordinary differential equations. This complexity was already partially known to mathematicians (H. Poincaré, A. N. Kolmogorov or G. D. Birkhoff), but was largely unfamiliar to the physicists, chemists and biologists.

Emerging phenomena like spatial and/or temporal structures, or, more mathematically, of limit cycles or chaotic regimes are all related to the nonlinearity of the laws that govern many physical systems. For spatially uniform systems, these laws are expressed as systems of nonlinear ordinary differential equations (ODEs). For inhomogeneous systems, diffusion and convection processes represented by spatial derivative operators must be added, thereby transforming the systems of ODEs into systems of nonlinear partial derivative equations.

The attitude of Ilya Prigogine faced with these systems was that of a thermodynamicist and a statistical physicist. He saw them as equations describing the evolution of macroscopic quantities fluctuating out of thermodynamical equilibrium. For small excursions out of that state, the equations are linear in the fluctuations of the macroscopic variables and nothing new appears apart from the equilibrium properties. However, for large fluctuations, these equations become nonlinear. Such large deviations can be sustained in a system by imposing on it external fluxes of matter, momentum and energy. Ilya Prigogine realized that when the values of these fluxes crossed some thresholds, qualitative changes emerged in the spatial and temporal organization of the systems. These transitions were in fact, from a mathematical standpoint, related to bifurcations. Very early, Ilya Prigogine foresaw that the emergence of primitive living organisms from inanimate matter was due to such transitions appearing when some thresholds in the intensities of external fluxes were crossed.

Inspired by this vision he attracted many researchers from a large spectrum of disciplines and engaged them in an extensive and multidisciplinary long-term research programme about the properties of nonlinear dynamical systems. This programme concerned domains as far apart as hydrodynamics, plasma physics or nonlinear optics, and molecular biology, insect societies, plant communities or the growth of cities. This led his group to many important breakthroughs in these fields [1].

Through many discussions, Ilya Prigogine transmitted the virus of nonlinearity to me. I mainly became interested in the mathematical aspects of nonlinear systems. While the main research orientation on these systems was essentially concerned with the topological properties of the phase space portrait, I engaged in a more constructivist approach.

In the following, I present a review of results, some of them quite recent, obtained in this field with several long-term co-workers. Most of these results have already been published in various articles except some of those appearing in §§4–6. In §2, the quasi-polynomial formalism is presented. Section 3 is a brief summary of results on the total or partial integrability of quasi-polynomial systems, their reducibility to lower dimensional systems and their stability properties. Section 4 presents an account of the Lie algebraic structure underlying the quasi-polynomial dynamical systems. Section 5 is dedicated to the solution of quasi-polynomial systems in terms of Taylor series with the analytic form of the general coefficient. Section 6 concerns a recently discovered equivalence between quasi-polynomial dynamical systems and urn stochastic processes.

2. Quasi-polynomial systems and canonical forms

In contrast to linear autonomous differential systems

in which the real independent variable t is time and L is any real and constant matrix, nonlinear differential systems exhibit a great diversity of functional dependance of the functions f_i(x₁,…,x_n). The vector constituted by these functions is called the vector field of the dynamical system (2.2). These functions, depending on the laws governing the phenomena to be modelled, may present any kind of nonlinearity. Thus, the vector field may have components that are combinations of linear functions, polynomials, elementary functions or even special functions. The only mathematical requirement comes from the conditions of the Cauchy–Lipschitz theorem for the existence and uniqueness of the solution of system (2.2). The drawback of this diversity is that it has prevented the construction of a theoretical apparatus as powerful as the one developed from linear algebra for linear dynamical systems of form (2.1).

However, despite their great variety, many nonlinear differential systems can be brought [2] to a unique format, the so-called quasi-polynomial (QP) format, also named Generalized Lotka–Volterra (LV) format [3]. A quasi-polynomial function is defined as a linear combination of monomials in the dependent real variables x₁,…,x_n in which the exponents are real numbers. Any QP system can be written in the following form:

where the numbers N and n are not necessarily equal.

As can be seen, in contrast with the linear systems that are characterized by only one square matrix, L, QP systems need two rectangular matrices, A and B to be characterized. The entries of the n×N matrix A are real and constant, and denote the coefficients of the monomials , while the entries of the N×n matrix B, also real and constant, are the exponents of these monomials. To avoid divergences of certain monomials with negative exponents, the system must be such that its solutions stay in an open n-dimensional orthant for initial conditions that are already in it. Let us denote QP(A,B) as such a system (2.3).

Before going further, let us describe the contexts in which QP systems are found. They represent a large subset of the non-equilibrium dynamical systems that appear in various domains and can undergo self-organization transitions [4]. In nonlinear optics, the Lamb equations for a three-mode operating laser cavity are of that form. Electronic circuits obey systems of ODEs with a polynomial vector field. Chemical reactions are governed by rate equations that are systems of ODEs with polynomial nonlinearities. The dynamics of interacting populations also obeys the system of polynomial ODEs. In mechanics, systems with polynomial interaction potentials are represented by systems of ODEs with a polynomial vector field and, consequently, can be cast in QP form (2.3). Even in cosmology such systems appear under the name of generalized LV systems [5]. All these systems can be directly written in the QP form (2.3) without increasing the number of dependent variables.

Moreover, dynamical systems involving non-polynomial functions can sometimes be transformed to the QP format. For instance, biochemical kinetic equations for reactions involving enzymatic catalysis contain the so-called Michaelis–Menten terms. These terms are rational functions in dependent variables. In such cases, the QP format is recovered at the price of increasing the dimensions of the ODE system by introducing new variables that precisely are these terms. More generally, all ODE systems involving non-polynomial functions of the dependent variables in their vector field can be cast in the QP form if these functions obey ODEs with polynomial nonlinearities in terms of these functions [2]. The QP format (2.3) offers a unifying tool that turns out to be quite effective and used, among others, in the fields of biological networks [6] and of feedback control [7] of biological and mechanical systems.

The recognition of the quasi-polynomial notation (2.3) of polynomial vector fields with two matrices A and B paves the way to a novel approach of all the dynamical systems that can be cast in the form QP(A,B). Indeed, the specific notation (2.3) lends itself to identify a group of mappings of the phase-space that preserve the QP format. These transformations read

where the matrix C can be any real and invertible n×n matrix. This mapping transforms equation (2.3) into with and

A consequence of equations (2.6) and (2.7) is

As can be seen, these transformations (2.4) form a group of diffeomorphisms that conserve the QP format (2.3). Moreover, the matrix product BA is an invariant of these mappings. This means that all systems QP(A,B) with the same matrix product BA belong to the same equivalence class under these transformations. Moreover, the phase portrait is preserved by these transformations.

These properties have been found independently several times [8–10] and led their authors to the discovery of two fundamental canonical forms for each equivalence class. The proof goes as follows. Let us replace the dependent variables x₁,…,x_n of the QP system (2.3) by the N monomials appearing in it

with j=1,…,N.

A simple calculation, using equation (2.3), of the time-derivative of these new variables gives

where j=1,…,N and

System (2.10) is the first canonical form. These equations are well known to specialists in population dynamics: they are the LV equations about which there is an extensive number of results. The N×N matrix M is called the population interaction matrix. This terminology is perhaps a bit far fetched for arbitrary matrices A and B; the matrix M does not generally fulfil the conditions that are required by population dynamics. However, properties such as the existence of an explicit Lyapunov function for systems with an interior fixed point and with an admissible matrix M still exist [11]. Another fundamental property is the existence, in an appropriate Poisson structure, of Hamiltonian functions for the N-dimensional LV systems [12].

Equation (2.11) means that all the QP(A,B) systems such that the product BA of their matrices is equal to a given matrix M are equivalent to the system of LV equations (2.10) with that same interaction matrix. They, thus, constitute an equivalence class of label M.

For N>n, obviously the N variables u_j are not all independent. As a consequence, the trajectories of the original system (2.3) are confined to invariant subspaces corresponding to the relations between the variables u_j in the N-dimensional phase-space. In the case n>N, the passage from equation (2.3) to (2.10) corresponds to a reduction of dimensions of the system. More details on these questions can be found in references [3,7,11]. Note that transformation (2.9) is, in general, not bijective except for N=n, in which case it belongs to the group of monomial transformations (2.4). Let us also stress that the LV canonical form is itself a particular N-dimensional QP system with A=M and B=I, where I is the identity matrix. In other words, it is a QP(M,I) system.

With regard to the second canonical form, it appears as follows. As the LV form (2.10) is an N-dimensional QP system, and provided the matrix M is invertible, one can apply to this system the following monomial transformation:

for k=1,…,N.

This is a transformation of type (2.4) with C=M. Now, using relations (2.6) and (2.7) one readily finds

with j=1,…,N. This is the second canonical form. We shall call it the monomial canonical form. Unlike the first canonical form (2.10) which, in addition to population dynamics, is well known in several domains, such as evolution theory [13], plasma physics [14], cosmology [5], this second canonical form is almost unknown. Only recently and by accident, did we find an occurrence of this equation in relation to the so-called urn stochastic processes. This led us to find a deep equivalence relationship between the latter and the whole class of dynamical systems that are reducible to the QP format, as we shall see later.

Four more comments must be made on the second canonical form. First, in the case where the matrix M is not invertible but A is of maximal rank, it is possible to derive the second canonical form directly [3] from the QP(A,B) system (2.3). A sketch of the proof goes as follows. For N>n, one adds N−n auxiliary dependent variables obeying differential equations that contain only the monomials already present in the original system. This amounts to the original QP system being embedded into a new QP system in N dimensions for which the new matrix corresponds to the old rectangular matrix completed by N−n lines of arbitrary but non-vanishing entries. As for the new matrix , it is the matrix B completed by N−n columns of null entries. The two new matrices and are now square and if the original rectangular A was of maximal rank, that is n, then the square matrix is invertible. One, thus, can apply to the embedded system a monomial transformation of type (2.4) with an N×N matrix . Applying rules (2.6) and (2.7), this leads to a new QP system with a and whose form is purely monomial. In the case N<n, one shows [15] that the original QP system can be reduced to a square QP system of N dimensions on which the same transformations as above can be performed.

Second, the same comments as those already discussed for the LV canonical form can be made about the cases when N>n and N<n. Third, the monomial canonical form is a particular N-dimensional QP system with A=I and B=M, that is it is a QP(I,M) system.

A fourth comment must be made for both the canonical forms. It is current in models to find a diagonal linear term in the vector field of the system. Hence, one must consider systems of the form

where the coefficients a_i are real. Now, this system can be written in the following QP form: with an n×(N+1) matrix such that, for all i=1,…,n, one has for j=1,…,N and , and a (N+1)×n matrix such that, for all k=1,…,n, one has for j=1,…,N and .

Let us summarize the above results in the following

Statement 1.

The set of all the differential dynamical systems that can be cast into the QP format , for any couple of real constant matrices A and B, and whose solutions x_i(t) remain forever in the same open orthant, is divided into equivalence classes labelled by the matrix M=BA. In each equivalence class, there are two canonical forms and , j=1,…,N. Moreover, the solutions to the QP systems that belong to a same equivalence class are exchanged under the group of monomial diffeomorphisms , i=1,…,N, where C is any matrix belonging to , and can be retrieved from the solutions of each canonical form.

Before leaving this chapter, let me give a simple example of the above transformations. Consider the following three-dimensional system where the nonlinearity is not quadratic

The form (2.14) of this system has N=2 and n=3. With the above reasoning, one thus expects a reduction of the dimension from 3 to 2 of the canonical forms. The matrices A and B of form (2.14) are

and while the coefficients of the linear terms are a₁=−1, a₂=3 and a₃=2.

The matrices and of the corresponding QP form (2.15) are

and

The matrix that appears in the two canonical forms of the QP system is

The LV canonical form, hence, is

The last equation has a constant solution u₃=c that can be found in terms of the initial conditions on the x_i(0), i=1,…,3. The above system of equations, thus, becomes

and

This is a two-dimensional quadratic system of LV form.

The monomial canonical form of the QP system is

The last equation can be solved yielding v₃=c′e^t, where c′ is a constant. This solution may be inserted in the two previous equations leading to a reduced system of two monomial ODEs.

The solution of the original system (2.16)–(2.18) can easily be found from the solution of either the LV form or the monomial form that are given above, if these equations can be solved.

The above results lead to many consequences that are exposed in the next chapter. Some of them cannot be detailed here by lack of place and we refer the reader to the publications. Others are presented in more details.

3. Applications

4. Lie algebraic structure of quasi-polynomial systems

5. General solution as Taylor series

In this section, we deduce the Taylor series expansion of the solution to the LV canonical form (4.1) and obtain the analytic expression of its general coefficient.

The starting point is the Carleman embedding equation (4.3). Let us rewrite the action (4.4) of the displacement operators E_i, i=1,…,N, as follows:

with r=(r₁,…,r_N) and where eⁱ is the unit vector along the coordinate axis i. The solution to equation (4.3) or (4.8) is with H given by formula (4.7).

We now express the exponential of the operator H by its Taylor series

The kth power of the Hamiltonian H on the initial function is readily obtained by successive iterations

Introducing in equation (5.4) the relations and (eⁱ.M)_j=M_ij, and inserting it into the Taylor series (5.3), we get

and using relation (4.6), x_i(t)=ψ(0,…,r_i=1,…,0,t)=ψ(eⁱ,t), the Taylor series for the solution of the LV system is obtained

This series in which the general coefficient is analytically known defines a new class of special functions. This class presents a certain analogy with the class of generalized hypergeometric functions. As for these functions, the differential equations they solve as well as their Taylor series, along with the analytic expression of their general coefficient, are known. The knowledge of the general term of series (5.6) allows for

1. evaluating the radius of convergence of the Taylor series

2. performing analytic continuation

3. evaluating the error at each truncation order of the series by calculating the remainder of the series

4. calculating Padé approximants

5. calculating the analytic expression of the solution at various approximation orders via computer algebra languages like MAPLE or Mathematica and construct a computer algebra software for solving QP systems [25].

One must now build the Taylor series that solves the original QP(A,B) systems that are equivalent to the LV canonical form (4.1). In order to do this, one has to resort to an arbitrary monomial transformation on the dependent variables x_i of the LV system

Indeed, we know from relations (2.6) and (2.7) that this transforms the LV equation into

with

For systems with rectangular matrix , the above transformation can be performed after an embedding of that matrix into a square matrix and can be performed only at the condition that is of maximal rank.

All systems (5.8) belong to the same equivalence class with corresponding to the LV canonical form.

This transformation induces the following transformation of the Carleman ψ function

which, using equations (5.5) and (5.9), yields

The solution of the QP systems of form (5.8) is, then, obtained by putting r=eⁱ in the above function

for i=1,…,N and with M=BA.

Here, again, the knowledge of the analytic form of the coefficient of the Taylor series brings the same benefits as enumerated above for the series solving the LV form. It defines also a new special function for any forms of the two rectangular matrices A and B.

One should stress here that the cases of the QP systems (2.3) with n≠N are included in the above result. Indeed, for n<N as well as for n>N, the QP system can be, respectively, extended or reduced to an N-dimensional QP system that transforms into the LV form under a monomial transformation. More details on these question can be found in [3].

One can conclude this chapter by the following statement:

Statement 4.

The Taylor series of the solution to the QP system , for the initial conditions x_i(0)=x_0i, is given by

This series defines an infinite class of special functions labelled by the two rectangular matrices A and B. These functions obey the above QP(A,B) system.

As said above, the knowledge of the general term of the Taylor series gives information on the convergence radius. But it also provides an evaluation of the rest of the series when the latter is truncated at a given order. Moreover, it offers the possibility of calculating the analytic continuation in order to extend the time domain on which the dynamics of the system is known.

6. Equivalence between quasi-polynomial systems and urn stochastic processes

We report here a recent and quite unforeseeable advance. It establishes an equivalence between stochastic processes called urn processes and a sub-class of the QP systems [26].

Urn processes are stochastic processes that offer powerful modelling tools in various scientific fields such as statistical physics, population genetics or epidemiology [27,28]. An urn process consists of a box, the urn, containing at the initial time a distribution of balls of N different colours, and an elementary process of replacement of these balls that we now describe. The initial composition of the urn is given by the vector U⁰=(u₁₀,…,u_N0), where u_i0 is the initial number of balls of colour i in the urn.

The process starts as soon as a ball is picked from the urn with uniform probability, its colour is noted and the ball is re-introduced in the urn. If the colour of the drawn ball was i, i=1,…,N, then specified numbers M_ij of balls of the different colours j=1,…,N are taken from an infinite reservoir of balls and introduced in the urn. If the integer number M_ij is negative, then that number of balls of colour j must be extracted from the urn and put in the reservoir.

We shall assume here that the process is balanced, that is at each step the total variation of the number of balls in the urn is a constant σ. This means that . At each step of the process, this procedure is repeated with the same replacement N×N matrix M. After n successive replacement steps, the initial composition vector U₀ evolves into a final vector U=(u₁,…,u_N).

One would like to calculate the probability of such a transition in n steps between given vectors U₀ and U. We adopt here the combinatorial approach developed in [29]. The trajectory of the n-step process in the N-dimensional composition space is called a history of length n. Let us introduce the number of possible histories of length n linking U₀ to U, H_n(u₁₀,…,u_N0,u₁,…,u_N). We also define the number of n-steps histories starting from U₀ irrespective of the final composition vector, H_n(u₁₀,…,u_N0).

As the process is balanced, all the histories in n steps are equiprobable. Hence, the probability we are looking for is given by

However, the most adapted tools for our purpose here are not the number of length n histories but their counting generating function defined as

In terms of that function the searched probability reads

where the notation [x^m]S(x) for a power series S(x) represents, as usual, the coefficient of x^m in that series.

Now, there is a connection between the above generating function and the monomial canonical form (2.13) of an equivalence class of QP systems. Ph. Flajeolet, Ph. Dumas and V. Puyhaubert established the following deep property [29]: any N-colours balanced urn process of replacement matrix M is isomorphic to the system of ordinary differential equation

with i=1,…,N. This isomorphism amounts to the following equation relating the history counting generating function of the urn process and the solution of the above system where X_i(z), i=1,…,N, is the solution at time t=z of system (6.4) with initial conditions X_i(0)=x_i0.

Apparently, the authors of this theorem did not know the theory of QP systems, otherwise they would immediately have recognized system (6.4) as one of the two canonical forms of that theory.

Using our statement 1, this identification allows us to make the following new statement:

Statement 5.

To any balanced urn process of replacement matrix M, one can associate an infinite equivalence class of QP(A,B) systems, , i=1,…,n, with BA=M. The counting generating function of this urn process is obtained via the relation where the X_i(z) are the solutions of the monomial canonical differential system associated with the above-mentioned equivalence class: , i= 1,…,N for t=z and for the initial conditions X_i(0)=x_i0. Conversely, the solution of this monomial differential system for the initial conditions X_i(0)=x_i0 can be calculated from the counting generating function of the associated urn process by: X_i(z)=H(x₁₀,…,x_N0,0,…,u_i0=1,…0,z), for z=t. As a consequence, the solution of any QP(A,B) system with BA=M can be calculated from the previous solution by applying to it the monomial transformations for any invertible matrix C such that A=C⁻¹M and B=C.

A promising application of this result would be to use the above result to build an entirely new approach to computer simulation of dynamical systems. One could, indeed, compute the solutions of the monomial canonical form with initial conditions X_i(0)=x_i0 by simulating the associated urn process and taking advantage of the relation X_i(z)=H(x₁₀,…,x_N0,0,…,u_i0=1,…0,z). The coding in a computer simulation program of an elementary urn replacement process of coloured balls should be easy to implement and its computation, fast. The method would be limited, of course, to solving QP(A,B) systems such that the product matrix BA has only integer entries.

7. Conclusion

I hope to have shown that theory of quasi-polynomial systems is interesting. What started as merely a new notation of differential dynamical systems led to many applications. These were based on the property of systems with a large spectrum of nonlinearities to be transformed into two canonical forms, the LV and the monomial forms. The first form contains the simplest type of non-trivial quadratic nonlinearity and has been the object of a huge number of publications. The second, though, much less known, possesses many potentialities that should be envisaged. Future results on these two types of equations are waiting to be translated into new properties for QP systems belonging to their equivalence classes.

Prospectively, the study of the relation between urn processes and dynamical processes prefigures interesting and new interactions between nonlinear dynamics and the theory of stochastic processes. Moreover, the study of the special functions corresponding to the Taylor series derived in this article could be of interest for special functions theory. The study of the asymptotic behaviour in time of these functions could shed some new light on the properties of chaotic regimes, for instance.

From the computing viewpoint, the equivalence urn processes-dynamical systems as well as the availability of Taylor series solutions with analytic general coefficient offer new perspectives. They provide radically novel tools for the numerical and symbolical computation of the solutions of these systems.

Data accessibility

This article has no additional data.

Competing interests

I have no competing interest.

Funding

No funding has been received for this article.