Bifurcations Associated with Three-Phase Polynomial Dynamical Systems and Complete Description of Symmetry Relations Using Discriminant Criterion

1. Introduction

Analysis of bifurcations in polynomial dynamical systems (DSs) consitutes a classical problem of the qualitative theory of DSs and differential equations and catastrophe theory and has been a subject of intense studies [1,2,3,4,5]. The main attention has been paid to nonautonomous polynomial DSs where many characteristic bifurcation types have been identified including equilibriums, sinks, saddles, limit cycles. Autonomous polynomial DSs investigated in these studies, especially the case of two dimensions, are non-integrable and characterized by ’mixed’ right-hand sides involving second-order polynomials in two variables.

In this regard, a complete investigation of particular families of DSs, and first of all quadratic polynomial DSs, which reveals all possible singular points and solutions is of definite interest, both theoretically and from the viewpoint of various applications. A specific interest here is towards identifying particular parameter DS families that admit complete description of all their singular points and bifurcations.

In this work, an integrable family of three-dimensional (three-phase) polynomial quadratic DSs is identified and considered. We show, as an extension of a similar study performed in [6] for two-dimensional polynomial quadratic DSs, that one can carry out complete investigation of singularities and in this manner create a complete qualitative theory in the three-phase case as well. Two-dimensional integrable polynomial quadratic DSs considered in [6,7] preserve certain clear types of symmetry; all characterisic cases are summarized in [6] as collections of curves on the phase planes. The symmetry is governed by the occurrence of a finite number of integrable combinations of the DS solutions; all of them are described and vizualized in [6]. We extend in the paper this finding to the 3D-case and demonstrate that there is finitely many integrable combinations of the three-phase DS solutions which can be fully classified by three parameters applying the so-called discriminant criterium. However, this number (which we have determined) is huge, much greater than in the 2D-case. In spite of a broad variety of such combinations, we state the absence of many bifurcations and critical modes known for autonomous and nonautonomous polynomial DSs, like limit circles. On the other hand, we report a specific type of bifurcation characteristic for the DSs under consideration, stating a conjecture that the three-dimensional polynomial quadratic DSs have only this very bifurcation type.

As it is known from the theory of two-phase filtration, when oil is displaced by water in reservoir conditions, an unstable front zone is formed, at which water saturation has a triple value and changes discontinuously, which in turn leads to water breakthrough to oil producing wells. The model and solution of this problem is proposed in [8,9,10,11,12,13] as well as various applications of the proposed technique. In the case of three-phase filtration, the gas phase is involved in the process, which significantly complicates the solution of the problem, since gas behaves the same way as water and breaks through to the producing wells even more actively. As a result, all this leads to premature breakthrough of gas and water to the producing oil wells and, consequently, to the reduction of oil recovery and oil production stimulation. At present there are no theoretical mathematical solutions to this complex and important problem of oil science.

In a model of regulating and monitoring the process of waterflooding and the development of an oil or gas field, the considered three-dimensional autonomous polynomial DSs (system (42) in the tex below) describe [8,10,13] the change in time, at a given interval $t\in \left(t_0,t_1\right),t_0\ge 0,$ and under predetermined constant values of the problem parameters $a_{ij}$, $i,j=1,2,3,$the accumulated oil flow rate, x(t), water, y(t), and gas, z(t) corresponding to the three phases: oil, water, and gas. In practice, the accumulated and current well rates are known, which makes it possible to determine the coefficients of the polynomials entering DSs as constant values over a given interval. The discriminants of the growth model for each phase: for oil ($D_o$), gas ($D_g$), and water ($D_w$), are used as appropriate control parameters.

The particular form of three-dimensional autonomous polynomial DSs with ’separated’ variables (as (42)) is dictated by the results of long-standing observations and monitoring of oil and gas fields and huge amount of the registered, experimental and measurement data [7,11,18].

The present paper proposes a new approach to create qualitative theory and describe singularities and bifurcations of the DSs under study on the basis of the growth model from catastrophe theory [18,19,20,21]. The model has been equipped with new tools that enable one to classify all possible solution types into equivalence classes defined in terms of the so-called discriminant (D-) vectors.

With a certain combination of values and signs of discriminants for oil, water and gas and using the technique developed in this work, it is possible to formulate criteria that allow predicting the breakthrough of water and gas into oil producing wells.

Performing qualitative analysis of three-dimensional DSs addressed in the paper we describe specific scenarios for modeling of the concrete processes under study. The results obtained are summarized in a compact form of a matrix and a criterion that employs D- and S-vectors and enables one to identify all possible different types of stable and unstable solutions and their singularities governed in the end by the signs of the polynomials’ discriminants.

2. A Family of Polynomial Dynamical Systems

2.1. Matrix Representation of Polynomial Dynamical Systems

In order to propose and apply a general matrix framework for polynomial DSs and in this manner to identify the place of the anzats treated in this study consider the second- and third-order (two- and three-dimensional) polynomial DSs that inolve quasi-polynomials (QPs); the DSs of this class can be written in the general form containing all quadratic terms

$$\left\{\begin{array}{c}\dot{x}=a_0+a_{1}x+a_{2}y+a_3{x}^2+a_4{y}^2+a_{5}xy,\\ \dot{y}=b_0+b_{1}x+b_{2}y+b_3{x}^2+b_4{y}^2+b_{5}xy,\end{array}\right.$$

(1)

or

$$\left\{\begin{array}{c}\dot{x}=a_0+a_{1}x+a_{2}y+a_{3}z+a_4{x}^2+a_5{y}^2+a_6{z}^2+a_{7}xy+a_{7}xz+a_{9}yz,\\ \dot{y}=b_0+b_{1}x+b_{2}y+b_{3}z+b_4{x}^2+b_5{y}^2+b_6{z}^2+b_{7}xy+b_{8}xz+b_{9}yz,\\ \dot{z}=c_0+c_{1}x+c_{2}y+c_{3}z+c_4{x}^2+c_5{y}^2+c_6{z}^2+c_{7}xy+c_{8}xz+c_{9}yz,\end{array}\right.$$

(2)

where dot denotes differentiation with respect to time.

In view of the technique developed in this work, it is convenient to introduce the coefficient 2 x 5 or 3 x 9 matrices associated with DS (1) or (2),

$$\begin{array}{c}\hfill A^{\left(2\right)}=\left[\begin{array}{cccc}{a}_0& a_1& \cdots & a_5\\ b_0& b_1& \cdots & b_5\end{array}\right],A^{\left(3\right)}=\left[\begin{array}{cccc}{a}_0& a_1& \cdots & a_9\\ b_0& b_1& \cdots & b_9\\ c_0& c_1& \cdots & c_9\end{array}\right].\end{array}$$

(3)

Different sets of real coefficients $a,b$ or $a,b,c$ that enter matrices (3) govern qualitative behaviour of solutions and the occurrence and character of singularities and bifurcations. Note that generally all these DSs are non-integrable.

A well-known particular case of (1) constitutes linear homogeneous two-dimensional DSs that admit a matrix representation employing a square matrix

$$\left\{\begin{array}{c}\dot{x}=a_{1}x+a_{2}y,\\ \dot{y}=b_{1}x+b_{2}y,\end{array}\right.or\dot{\mathbf{x}}=A\mathbf{x},A=\left[\begin{array}{cc}{a}_1& a_2\\ b_1& b_2\end{array}\right]$$

(4)

(because the number of unknowns equals the number of equations). All their bifurcation types (critical modes) at the origin are identified and classified in terms of the eigenvalues ${\lambda }_1,{\lambda }_2$ of the DS coefficient (real-valued) matrix A; more precisely, by the combination of their signs (or of the signs of their imaginary parts if they are complex) and zero values ${\lambda }_{1,2}>0$, ${\lambda }_{1,2}<0$, ${\lambda }_1>0,{\lambda }_2<0$, ${\lambda }_1=0,{\lambda }_2\ne 0$ etc (node (source, sink), saddle, focus, center etc). Every such combination naturally gives rise to a class of equivalence on the set $R^4$ of square $2\times 2$ matrices with real entries characterized by the unique combination of signs of eigenvalues or their real and imaginary parts. Limiting ourselves to the set $R^{4,re}$ of matrices having real nonzero eigenvalues one may identify three such classes (without separating the cases of double roots) denoting them by the ’sign’ vectors ${\mathbf{s}}_1=\left[++\right]$, ${\mathbf{s}}_2=\left[--\right]$, and ${\mathbf{s}}_3=\left[-+\right]$. Each class has its particular type of bifurcation at the origin (source, sink, saddle).

One can create general qualitative theory for the DSs described by (4). However, it is not possible to create any general qualitative theory for the whole DS class described by (1) or (2). Therefore, a common practice is that researchers (beginning from Hilbert and Poincare) specify certain sub-classes of quadratic polynomial DSs, mainly in 2D, which are amenable to qualitative analysis. The present study is on this track and picks up a certain sub-class which can be fully investigated as a family of integrable DSs.

In [2,5] a detailed qualitative investigation is performed of two-dimensional polynomial DSs and a long list of the relevant publications can be found. A global bifurcation theory of such systems is presented, including particularly the issues connected with solution to the famous Hilbert’s Sixteenth Problem concerning the determination of the maximum number and relative position of limit cycles. In this respect, certain specific families of two-dimensional polynomial DSs (1) are identified. Particularly, give examples of their so-called canonical forms [5] with ’mixed’ variables accepted in the literature

$$\left\{\begin{array}{c}\dot{x}=-y-\alpha y^2-xy,\\ \dot{y}=x+(\lambda +\beta +\gamma )y+a{x}^2+c\gamma y^2+(\alpha +\beta +\gamma )xy,\end{array}\right.$$

(5)

with

$$\begin{array}{c}{a}_0=a_1=a_3=0,a_2=-1,a_4=-\alpha ,a_5=-1,\\ b_0=0,b_1=1,b_2=\lambda +\beta +\gamma ,b_3=a,b_4=c\gamma ,b_5=\alpha +\beta +\gamma ,\end{array}$$

(6)

and the coefficient matrix

$$\begin{array}{c}\hfill A^{\left(2\right)}=\left[\begin{array}{cccccc}0& 0& -1& 0& -\alpha & -1\\ 0& 1& \lambda +\beta +\gamma & a& c\gamma & \alpha +\beta +\gamma \end{array}\right],\end{array}$$

(7)

or

$$\left\{\begin{array}{c}\dot{x}=-y-\nu y^2,\nu =0,1,\\ \dot{y}=x+(\lambda +\beta +\gamma )y+a{x}^2+c\gamma y^2+(\beta +\gamma )xy,\end{array}\right.$$

(8)

in (1) with

$$\begin{array}{c}{a}_0=a_1=a_3=a_5=0,a_2=-1,a_4=-\nu ,\\ b_0=0,b_1=1,b_2=\lambda +\beta +\gamma ,b_3=a,b_4=c\gamma ,b_5=\beta +\gamma ,\end{array}$$

(9)

in (1) and the coefficient matrix

$$\begin{array}{c}\hfill A^{\left(2\right)}=\left[\begin{array}{cccccc}0& 0& -1& 0& -\nu & 0\\ 0& 1& \lambda +\beta +\gamma & a& c\gamma & \beta +\gamma \end{array}\right],\end{array}$$

(10)

or

$$\left\{\begin{array}{c}\dot{x}=-y+\nu y^2,\nu =0,1,\\ \dot{y}=x+(\lambda +\beta )y+a{x}^2+c{y}^2+\beta xy\end{array}\right.$$

(11)

in (1) and the coefficient matrices

$$\begin{array}{c}\hfill A_{\nu =0}^{\left(2\right)}=\left[\begin{array}{cccccc}0& 0& -1& 0& 0& 0\\ 0& 1& \lambda +\beta & a& c& \beta \end{array}\right],A_{\nu =1}^{\left(2\right)}=\left[\begin{array}{cccccc}0& 0& -1& 0& 1& 0\\ 0& 1& \lambda +\beta & a& c& \beta \end{array}\right].\end{array}$$

(12)

There is a number of other forms employing different sets of real coefficients $a,b,c$ characterized each by specific properties of bufircations.

The volume of research and quantity of publication dealing with qualitative analysis of quadratic two-dimensional polynomial DSs may be characterized as enormous; however, many problems remain unsolved. Exemplify the classification of the number and character of finite singularities reported for these systems; particularly for (11) it is established that there may be one saddle and three antisaddles, three saddles and one antisaddle, two saddles and two antisaddles etc depending on the parameter sets involved.

Qualitative analysis of quadratic three-dimensional polynomial DSs (2) is less developed; the amount of the obtained results could be hardly compared with that achieved for two-dimensional polynomial DSs. Note in this way a class of the T-systems of the form

$$\left\{\begin{array}{c}\dot{x}=y-x,\\ \dot{y}=mx-xz,\\ \dot{z}=-nz+xy\end{array}\right.$$

(13)

and the coefficient matrix composed according to (3)

$$\begin{array}{c}\hfill A^{\left(3\right)}=\left[\begin{array}{ccccccccccc}0& 1& -1& 0& 0& 0& 0& 0& 0& 0& \\ 0& m& 0& 0& 0& 0& 0& 0& -1& 0& \\ 0& 0& 0& -n& 0& 0& 0& 1& 0& 0\end{array}\right].\end{array}$$

(14)

The DSs belonging to the class of T-systems are studied in [22] (we also refer to the references therein), where, among all, the pitchfork and Hopf bifurcations occurring in the T-system are reported as well as the influence of the symmetry of the system on its dynamics. Seemingly simple, these systems require the elaboration of advanced mathematical tools within the frames the singular perturbation and geometric singular perturbation theory [23].

The latter has become a driving force: in our studies, we address, as well as in this one, a particular integrable family of quadratic autonomous three-dimensional polynomial DSs with ’separated’ variables and free terms

$$\left\{\begin{array}{c}\dot{x}=a_0+a_{1}x+a_4{x}^2,\\ \dot{y}=b_0+b_{2}y+b_5{y}^2,\\ \dot{z}=c_0+c_{3}z+c_6{z}^2.\end{array}\right.$$

(15)

with the coefficient matrix (3)

$$\begin{array}{c}\hfill A^{\left(3\right)}=\left[\begin{array}{ccccccccccc}{a}_0& a_1& 0& 0& a_4& 0& 0& 0& 0& 0& \\ b_0& 0& b_2& 0& 0& b_5& 0& 0& 0& 0& \\ c_0& 0& 0& c_3& 0& 0& c_6& 0& 0& 0\end{array}\right].\end{array}$$

(16)

Representation (15) where the number of terms on the right-hand side equals the number of dependent variables enables one to write the corresponding DS in a compact quadratic-matrix form which will be done below in this section. In view of this, it is convenient to rearrange the coefficient matrix of (15) to the symmetric $3\times 3$ form

$$\begin{array}{c}\hfill A^{\left(3\right)}=\left[\begin{array}{ccc}{a}_0& a_1& a_4\\ b_0& b_2& b_5\\ c_0& c_3& c_6\end{array}\right]={∥a_{kj}∥}_{k,j=1,2,3}.\end{array}$$

(17)

Here we see a formal similarity between the DS classes described by (4) and (15), (17) (partly because in both cases the number of the equation terms equals the number of equations yielding square DS coefficient matrices). In our analysis we extend this similarity and introduce and describe the equivalence classes using the ’sign’ S-vectors on the set of $3\times 3$ coefficient matrices.

Investigation of the DS (15) may be, among all, the first step towards analysis of more complicated polynomial DSs with ’mixed’ variables as in (2). Here, the sets of nine real coefficients $a,b,c$ in (15) govern qualitative behaviour of solutions; however, as we show, many important properties are governed by actually three decisive parameters. Namely, there is a finite number of integrable combinations of solutions to (15) and, remarkably, they all can be fully classified by three parameters, the discriminants of the polynomial entering DSs using the discriminant criterium [6,7].

The proposed matrix representation of polynomial DSs is among all definitely a useful tool of saving information about the DS in a compact form. Indeed, as we see from (5)–(12) and (13), many known kinds of DSs can be described using this compact matrix anzats.

A conjecture is that there are deep specific relations between the structure of the coefficient matrix (its rank, symmetry) and qualitative properties of a polynomial DS. In this study we confirm this statement for a specific integrable family of such DSs.

2.2. Discriminant Criterion and Matrix Representation of 3D Polynomial Dynamical Systems

The discriminant criterion employs introduction of the D-vectors which is based on the following general anzats: A DS coefficient matrix can be represented, as well as any $3\times 3$ matrix, in the form of a column of raw vectors

$$\begin{array}{c}\hfill A^{\left(3\right)}=\left[\begin{array}{ccc}{a}_0& a_1& a_4\\ b_0& b_2& b_5\\ c_0& c_3& c_6\end{array}\right]={∥a_{kj}∥}_{k,j=1,2,3}=\left[\begin{array}{c}{{\mathbf{a}}^{\prime }}_1\\ {{\mathbf{a}}^{\prime }}_2\\ {{\mathbf{a}}^{\prime }}_3\end{array}\right].\end{array}$$

(18)

Let $f(x,y,z):$$R^3\to$R denote a smooth real-valued function of three variables. An f-vector ${\mathbf{F}}^{\left(3\right)}$ associated with function f and matrix (18) can be defined as a ’3D vector-valued functional’ according to

$$\begin{array}{c}\hfill {\mathbf{F}}^{\left(3\right)}=\left[\begin{array}{c}f\left({{\mathbf{a}}^{\prime }}_1\right)\\ f\left({{\mathbf{a}}^{\prime }}_2\right)\\ f\left({{\mathbf{a}}^{\prime }}_3\right)\end{array}\right].\end{array}$$

(19)

A concrete form of f is dictated by specific needs of analysis.

Particularly, we define f-vectors of the DS coefficient matrices as the discriminant (D-) column vector; simultaneously we define the (S-) row vector of the discriminant signs associated with the DS coefficient matrix (17):

$$\begin{array}{c}\hfill {\mathbf{D}}^{\left(3\right)}={\mathbf{D}}^{\left(3\right)}\left(A^{\left(3\right)}\right)=\left[\begin{array}{c}{D}_1\\ D_2\\ D_3\end{array}\right],{\mathbf{S}}^{\left(3\right)}={\mathbf{S}}^{\left(3\right)}\left(A^{\left(3\right)}\right)=\left[sign{D}_1,sign{D}_2,sign{D}_3\right].\end{array}$$

(20)

S-vectors may be equally represented as sets of the three ordered numbers $-1,0,1$, so that, for example,

$$\begin{array}{c}\hfill {\mathbf{S}}^{\left(3\right)}=\left[+0-\right]=\left[10-1\right]if{D}_1>0,D_2=0,D_3<0.\end{array}$$

The quantities ${D=D}_i={a_{i2}}^2-4a_{i1}{a}_{i3}$ in (20) specifying the concrete form of f are the discriminants of the quadratic trinomials $P_i$entering the ith row of system (42); they are considered each as a real-valued function of three variables $a_{i1},a_{i2},a_{i3}$, $i=1,2,3$, with the range being the set of all real numbers.

Generally an f-vector (19) and particularly a D-vector (20) is an aggregate quantity (a ’vector-valued functional’) describing in a compact form certain important properties of a $3\times 3$ matrix or particularly the DS coefficient matrix and in this manner the DS itself, including its symmetries, bifurcations and singular points. The corresponding S-vector in its turn is an informative three-symbol description of a (set of) D-vector specifying characteristic classes of these vectors and yelding much less information than the D-vector: the range of D-vectors is the same as that of three-dimensional vectors with real components while the range of S-vectors is the set of three symbols (or three numbers $-1,0,1$). S-vectors form a set denoted by S³ which contains a finite number of elements (this number is determined in the next section).

D-vectors (20) may be naturally considered as elements of a three-dimensional space R³ and a subset of D-vectors corresponding to an S-vector in (20) composed by a particular triple of signs will be a particular set of R³. Namely: the first, second, etc octants of R³ correspond to S-vectors $\left[+++\right]$, $\left[-++\right]$, etc; the coordinate planes in R³ correspond to S-vectors $\left[\pm \pm 0\right]$, $\left[\pm 0\pm \right]$, $\left[0\pm \pm \right]$; and the coordinate axes in R³ correspond to S-vectors $\left[\pm 00\right]$, $\left[0\pm 0\right]$, $\left[00\pm \right]$ etc.

We can define now the following sets and relations (mappings, denoted by $<>$) that couple these sets:

$$\begin{array}{c}\hfill \mathrm{Set} \mathrm{of} \mathrm{polynomial} \mathrm{DSs}\left(15\right)\sim \mathrm{Set}{\mathrm{R}}^9\mathrm{of} \mathrm{their} \mathrm{coefficient} \mathrm{matrices} \left(17\right)<{\mathrm{R}}^9\to {\mathrm{R}}^3>\end{array}$$

(21)

$$\begin{array}{c}\hfill \mathrm{Set} {\mathrm{R}}^3\mathrm{of} \mathrm{D}-\mathrm{vectors}\left(20\right)<{\mathrm{R}}^3\to {\mathrm{S}}^3>\mathrm{Set} \mathrm{of} \mathrm{S}-\mathrm{vectors}\left(20\right).\end{array}$$

(22)

Application of the discriminant criterion means that, on the first step, we assign to a coefficient matrix (17) the discriminant D-vector and to the latter the S-vector of the discriminant signs (20). On the second step, one establishes classes of equivalence (invariance, symmetry) on the sets of the DS coefficient matrices and DS general solutions and the D-vectors in terms, respectively, of the D- or S-vectors (this issue is addressed in the next section). On the third step, analysis is performed of the DS qualitative behavior within a chosen equivalence class including the occurrence and character of bifurcations, as well as specific analysis of the transitions between classes (when one of the discriminants changes the sign).

Exemplify the relations between particular coefficient matrix families and its D- and S-vectors which, as will be shown below, govern qualitative properties of the corresponding DS solutions.

Identity and ’anti-identity’ matrices:

$$\begin{array}{c}\hfill A^{\left(3\right)}=I^{\left(3\right)}=\left[\begin{array}{ccc}1& 0& 0\\ 0& 1& 0\\ 0& 0& 1\end{array}\right],{I^{\prime }}^{\left(3\right)}=\left[\begin{array}{ccc}0& 0& 1\\ 0& 1& 0\\ 1& 0& 0\end{array}\right],\end{array}$$

(23)

$$\begin{array}{c}\hfill {\mathbf{D}}^{\left(3\right)}\left(I^{\left(3\right)}\right)={\mathbf{D}}^{\left(3\right)}\left({I^{\prime }}^{\left(3\right)}\right)=\left[\begin{array}{c}0\\ 1\\ 0\end{array}\right],{\mathbf{S}}^{\left(3\right)}=\left[0+0\right],\end{array}$$

Diagonal and ’anti-diagonal’ matrices:

$$\begin{array}{c}\hfill A^{\left(3\right)}=\left[\begin{array}{ccc}a& 0& 0\\ 0& b& 0\\ 0& 0& c\end{array}\right],{A^{\prime }}^{\left(3\right)}=\left[\begin{array}{ccc}0& 0& a\\ 0& b& 0\\ c& 0& 0\end{array}\right],\end{array}$$

(24)

$$\begin{array}{c}\hfill {\mathbf{D}}^{\left(3\right)}\left(A^{\left(3\right)}\right)={\mathbf{D}}^{\left(3\right)}\left({A^{\prime }}^{\left(3\right)}\right)=\left[\begin{array}{c}0\\ b^2\\ 0\end{array}\right],{\mathbf{S}}^{\left(3\right)}=\left[0+0\right].\end{array}$$

Upper- (lower-) triangular matrix:

$$\begin{array}{c}\hfill A^{\left(3\right)}=\left[\begin{array}{ccc}a& b& c\\ 0& d& e\\ 0& 0& f\end{array}\right],{\mathbf{D}}^{\left(3\right)}=\left[\begin{array}{c}{b}^2-4ac\\ d^2\\ 0\end{array}\right],{\mathbf{S}}^{\left(3\right)}=\left[\pm +0\right];\end{array}$$

(25)

$$\begin{array}{c}\hfill A^{\left(3\right)}=\left[\begin{array}{ccc}f& 0& 0\\ e& d& 0\\ a& b& c\end{array}\right],{\mathbf{D}}^{\left(3\right)}=\left[\begin{array}{c}0\\ d^2\\ b^2-4ac\end{array}\right],{\mathbf{S}}^{\left(3\right)}=\left[0+\pm \right].\end{array}$$

(26)

Symmetric matrix:

$$\begin{array}{c}\hfill A^{\left(3\right)}=\left[\begin{array}{ccc}a& b& c\\ b& d& e\\ c& e& f\end{array}\right],{\mathbf{D}}^{\left(3\right)}=\left[\begin{array}{c}{b}^2-4ac\\ d^2-4be\\ e^2-4cf\end{array}\right],{\mathbf{S}}^{\left(3\right)}=\left[\pm \pm \pm \right].\end{array}$$

(27)

One may continue this list for other particular matrix families.

We see that for matrices (23)–(24), the D- and S-vectors are invariant w.r.t. interchange (permutation) of the first and third raws, creating in this manner a definite symmetry (and the invariance or equivalence classes). This issue is discussed in more detail in the next section.

Next, according to relations (21) and (22) and between subsets of the space R³ of D- and S-vectors, the D- and S-vectors may be treated as quantities establishing classes of equivalence (invariance, symmetry) on the sets of (i) the DS coefficient matrices and DS general solutions in terms of the D- or S-vectors (when each D-vector corresponds to a particular subset of the DS coefficient matrices and each S-vector corresponds to a particular subset of D-vectors) and (ii) the D-vectors in terms of the S-vectors using relations (22). There are finitely many equivalence classes specified by condition (ii) and they are described in the next section.

Whether the DS coefficient matrix (17) belongs to a certain equivalence class determines the presence and nature (type) of singular points and bifurcations of the corresponding polynomial DS (15). This is a crucial reason to introduce D- and S-vectors and investigate quadratic polynomial DSs in terms of the equivalence classes (symmetry relations) defined using these vector quantities.

2.3. Symmetry Relations on the Sets of Coefficient Matrices and D-Vectors

D-vectors possess definite symmetry relations; namely, representing coefficient matrix (17) as a triple (raw) of column vectors,

$$\begin{array}{c}\hfill A^{\left(3\right)}={∥a_{kj}∥}_{k,j=1,2,3}=\left[{\mathbf{a}}_1{\mathbf{a}}_2{\mathbf{a}}_3\right]\end{array}$$

(28)

we can write a D-vector with components (20) as

$$\begin{array}{c}\hfill {\mathbf{D}}^{\left(3\right)}={\mathbf{D}}^{\left(3\right)}({\mathbf{a}}_1,{\mathbf{a}}_2,{\mathbf{a}}_3)={∥{a_{i2}}^2-4a_{i1}{a}_{i3}∥}_{i=1,2,3}={\mathbf{a}}_2^2-4{\mathbf{a}}_1{\mathbf{a}}_3\end{array}$$

(29)

and deduce that

and the following symmetry relations hold

$$\begin{array}{c}\hfill {\mathbf{D}}^{\left(3\right)}({\mathbf{a}}_1,{\mathbf{a}}_2,{\mathbf{a}}_3)={\mathbf{D}}^{\left(3\right)}({\mathbf{a}}_1,-{\mathbf{a}}_2,{\mathbf{a}}_3),{\mathbf{D}}^{\left(3\right)}({\mathbf{a}}_1,{\mathbf{a}}_2,{\mathbf{a}}_3)={\mathbf{D}}^{\left(3\right)}({\mathbf{a}}_3,{\mathbf{a}}_2,{\mathbf{a}}_1).\end{array}$$

(33)

The same symmetry relations (33) hold for S-vectors (D- and S-vectors are invariant. w.r.t. the sign of the central column and interchange of the side columns).

We see that for the matrices satisfying (30)–(32), the D- and S-vectors are the same; i.e. they are invariant w.r.t. any interchange (permutation) of the matrix raws creating in this manner definite symmetries (and the invariance or equivalence classes).

It is reasonable to introduce and consider a subset of coefficient matrices (17) of the form

$$\begin{array}{c}\hfill A^{\left(3\right)}=\left[\begin{array}{ccc}a& 2b& 1\\ c& 2d& 1\\ e& 2f& 1\end{array}\right]=\left[{\mathbf{a}}_{1}2{\mathbf{a}}_2\mathbf{1}\right],\end{array}$$

(34)

with the D-vectors

$$\begin{array}{c}\hfill {\mathbf{D}}^{\left(3\right)}={\mathbf{D}}^{\left(3\right)}({\mathbf{a}}_1,{\mathbf{a}}_2,\mathbf{1})={\mathbf{a}}_2^2-{\mathbf{a}}_1=\left[\begin{array}{c}{b}^2-a\\ d^2-c\\ f^2-e\end{array}\right]\end{array}$$

(35)

corresponding to DSs (15) with reduced polynomials having nonzero quadratic terms. Unlike (17) which is a nine-parameter matrix set, coefficient matrices (34) constitute a six-parameter subset of $3\times 3$ matrices. The following conditions hold that specify the subsets of the equivalence classes of matrices (34):

A case which may important is the set of degenerate coefficient matrices (34) having the forms

$$\begin{array}{c}\hfill \left[\alpha ·\mathbf{1}2{\mathbf{a}}_2\mathbf{1}\right]=\left[\begin{array}{ccc}\alpha & 2b& 1\\ \alpha & 2d& 1\\ \alpha & 2f& 1\end{array}\right],\left[{\mathbf{a}}_1\beta ·\mathbf{1}\mathbf{1}\right]=\left[\begin{array}{ccc}a& \beta & 1\\ c& \beta & 1\\ e& \beta & 1\end{array}\right]or\left[{\mathbf{a}}_1\gamma ·{\mathbf{a}}_1\mathbf{1}\right]=\left[\begin{array}{ccc}a& \gamma a& 1\\ c& \gamma c& 1\\ e& \gamma e& 1\end{array}\right].\end{array}$$

(39)

For these coefficient matrices the conditions (36) -(38) that specify belonging to the equivalence classes governed by the corresponding S-vectors may be greatly simplified; e.g.

$$\begin{array}{c}\hfill {\mathbf{D}}^{\left(3\right)}(\alpha ·\mathbf{1},2{\mathbf{a}}_2,\mathbf{1}):\alpha <0\to {\mathbf{S}}^{\left(3\right)}=\left[+++\right].\end{array}$$

(40)

Many more symmetry (equivalence, invariance) relations may be discovered and described for the set of D-vectors. Their complete description goes far beyond the scope of the present study and may be a subject of a great number of future works performed by other researchers.

3. Classification of Solutions to Autonomous Polynomial Equations

3.1. Representations of Autonomous and Integrable Polynomial Dynamical Systems

In view of integrability of (15) for arbitray set of its nine parameters (coefficients) $a,b,c$ we propose and develop a method of analysis of singularities and bifurcations based on direct investigation of its explicitly obtained general solutions and solutions to initial-value problems (rather than on the methods employing general bifurcation theory, characteristic DS polynomials and the like). This technique of the ’direct’ DS qualitative analysis has proved to be a good complement to general methods of analysis.

Following [6,7] present a brief description of a family of autonomous DSs written in the compact vector form as

$$\dot{\mathbf{x}}=\mathbf{F}\left(\mathbf{x}\right),\mathbf{x},\mathbf{F}\in {\mathrm{R}}^n,$$

where ${\mathrm{R}}^n$denotes n-dimensional vector space. A sub-family of autonomous DSs for which vector-functions $\mathbf{F}$ are QPs, namely, finite sums of monomials containing powers of the dependent variables that can be real numbers are exemplified in (1)–(15).

Autonomous symmetric n-dimensional polynomial DSs form a class of autonomous DSs,

$$\frac{d{x}_k}{dt}=\sum _{i,j=1}^n{A}_k^{ij}{x}_i\left(t\right)x_j\left(t\right),k=1,2,\cdots ,n,$$

involving symmetric polynomials in n variables of degree two with equal-index quadratic terms ($i=j$), where $A_k^{ij}=A_k^{ji}$are given real numbers and $\mathbf{x}=\left[x_1\left(t\right),x_2\left(t\right),\cdots ,x_n\left(t\right)\right];$ below, when convenient, we will a use the notation $\mathbf{x}=\left[x,y,z\right]$ in the 3D case under an assumption that ${\mathrm{R}}^{\mathbf{3}}$is equipped with a Cartesian coordinate system.

There are several families of algebraically integrable autonomous polynomial DSs mentioned in Introduction. In general, however, autonomous polynomial DSs are not integrable. Note that such DSs containing solely equal-index quadratic terms with the nonzero $A_k^{kk}=a_{1k}$ ($k=1,2,\cdots ,n)$and the rest of $A_k^{ij}=0$ consitute a sub-family of degenerate autonomous polynomial DSs as in (15) and are algebraically integrable in elementary functions.

Autonomous n-dimensional polynomial DSs which are generally nonsymmetric and where (each) ith component of $\mathbf{F}$ is a polynomial of degree $N_i$ in one variable$x_i,i=1,2,\cdots ,n,$ with real coefficients can be represented as

$$\dot{\mathbf{x}}=\mathbf{P}\left(\mathbf{x}\right),\mathrm{or}\frac{d{x}_i}{dt}=P_i\left(x_i\right),P_i\left(s\right)=\sum _{j=1}^{N_i+1}{a}_{ji}{s}^{N_i-j+1},i=1,2,\cdots ,n,$$

(41)

where $a_{ji}$are given real numbers. To every DS of the type (41) the vector of the polynomial dimensions $\overline{N}={\left\{N_i\right\}}_{i=1}^n$ can be assigned as well as the number $N={\sum }_{i=1}^n{N}_i$, the DS total order. Autonomous n-of dimensional polynomial DSs (41) are algebraically integrable in elementary functions; the corresponding statements are proved in Appendix A.

DSs (41) of the dimension $n=3$ with $\overline{N}=\left[2,2,2\right]$ and total order 6,

$$\dot{\mathbf{x}}=\mathbf{P}\left(\mathbf{x}\right),\mathrm{or}\frac{d{x}_k}{dt}=P_k\left(x_k\right),P_k\left(s\right)=\sum _{j=1}^3{a}_{kj}{s}^{j-1},k=1,2,3$$

(42)

(this form is equivalent to (15) with matrix (17) having the entries $a_{11}=a_0$, $a_{12}=a_1$, $a_{13}=a_4$, $a_{21}=b_0$, $a_{22}=b_2$, $a_{23}=b_5$, and $a_{31}=c_0$, $a_{32}=c_3$, $a_{33}=c_6$) are generally nonsymmetric, each component of $\mathbf{F}$ is a polynomial of degree two (with real coefficients) in one variable$x_j,j=1,2,3$.

As follows from Theorem 1 in Appendix A, DSs (15) are algebraically integrable in elementary functions. In other words, DSs (2) with coefficient matrices (16) are algebraically integrable.

The sets of solutions to (42) (or (15)) are described in the next section.

Obtaining general solution to DSs (42) is based on

(a) factorization of the polynomials $P_i\left(s\right)$on the right-hand side of DS (15), (41), or (42) to a product of irreducible polynomials performed in Appendix A,

(b) subsequent representation of functions ${P_i\left(s\right)}^{-1}$ as partial fraction decomposition, and

(c) integration of the obtained decomposition.

All steps (a)—(c) can be accomplished explicitly for polynomials of degree two ending up with general solution to the corresponding three-dimensional polynomial DSs.

In this paper we investigate qualitatively the behavior of the solutions to DS (42). With this purpose we make use of the DS coefficient matrices and divide the DS family under study into the equivalence classes using the specially defined discriminant vectors. Based on this framework, we introduce and identify the particular bifurcations using the proposed methods of analysis of polynomial quadratic DSs. We employ the characteristic property inherent just to quadratic three-dimensional DSs: solutions $\mathbf{x}=\mathbf{x}\left(t\right)$ to these DSs and particularly to (42) are curves in ${\mathrm{R}}^{\mathbf{3}}$parametrized by time t.

We consider the properties of the general solution of each of the equations in the system using the example of the equation $\dot{x}=f\left(x\right),$ and solutions of the Cauchy problem $\dot{x}=f\left(x\right)$, $x\left(0\right)=q_0$(of the quadratic polynomial DS) on the phase plane x, t depending on the problem parameters.

3.2. General Solutions to Autonomous Second-Order Polynomial Equations

The general solutions of the autonomous polynomial equation $\dot{x}=f\left(x\right),f\left(x\right)=a{x}^2+bx+c,$ is obtained for all possible factorization cases (with respect to the signs of the discriminants) and combinations of the polynomial coefficients in the system. The discriminants of polynomials (namely, their signs) are the defining parameters that control the essential properties of solutions specify their complete classification.

• (a)$a=0$; $\dot{x}=f\left(x\right)$ is a first-order linear equation with constant coefficients and its solutions $x\left(t\right)={\frac{1}{b}(-c+Ce}^t),b\ne 0,$ are differentiable monotone functions on the whole line.
  In what follows we will consider the equations with $a\ne 0.$
• (b)D> 0; there are three families of the general solution

  $${x\left(t;C\right)=x}_1+\frac{\sqrt{D}}{a}\frac{1}{-1+C{e}^{-\sqrt{D}t}},x_{1,2}=-x_0\pm \frac{\sqrt{D}}{2a},x_0=\frac{b}{2a},D=b^2-4ac,$$

  of the equation $\dot{x}=f\left(x\right)$:
  Family U, $C>0$: solutions $x\left(t;C\right)=x_{\mathrm{U}}\left(t;C\right)$ are unstable, since they
  (i) have a “movable” singularity with $t=t^{*}\left(C\right)=\frac{lnC}{\sqrt{D}}$,
  (ii) are monotonically increasing functions on the whole line for $t\ne t^{*}$ since $\frac{dx\left(t;C\right)}{dt}=\sqrt{D}\frac{C{e}^{-\sqrt{D}t}}{{(C{e}^{-\sqrt{D}t}-1)}^2}>0$for C> 0,
  (iii) satisfy the condition $x\left(t;C\right)>x_2,x<t^{*},x\left(t;C\right)<x_1,x>t^{*},$ and
  (iv) have two horizontal asymptotes $x=x_1\mathrm{or}x=x_2$; these solutions are not bounded on the whole line (the latter means instability).
  Family S, $C>0$: solutions $x\left(t;C\right)=x_{\mathrm{S}}\left(t;C\right)$are stable since they have no singularities, are monotonically decreasing $\left(\frac{dx\left(t;C\right)}{dt}<0\right)$and bounded functions on the whole line, satisfying the condition $x_1<x\left(t;C\right)<x_2$, and have two asymptotes $x=x_1\mathrm{and}x=x_2$ (the last condition written in the form of limits (A1) means stability).
  Family T, : $x=x_1$ and $x=x_2$are stationary solutions of the equation $\dot{x}=f\left(x\right)$satisfying the initial conditions $x\left(0\right)=x_{1}orx\left(0\right)=x_2,$and the first corresponds to $C=0$ in ${\displaystyle {x\left(t;C\right)=x}_1+\frac{\sqrt{D}}{a}\frac{1}{-1+C{e}^{-\sqrt{D}t}}}$.
  Stationary solutions $x=x_1$ and $x=x_2$ are ’nonisolated’, so that in every neighborhood of each of these solutions there are infinitely many ’non-singular’ solutions $x_{\mathrm{U}}\left(t;C\right)$ or $x_{\mathrm{S}}\left(t;C\right)$from families U or S:

  $$\forall \epsilon >0\exists B>0,C\ne 0:\left|x_{\mathrm{U}}\left(t;C\right)-x_{1,2}\right|<\epsilon ,\left|t\right|>B\mathrm{or}\left|x_{\mathrm{S}}\left(t;C\right)-x_{1,2}\right|<\epsilon ,\left|t\right|>B,$$

  and an arbitrarily small change of the parameter $q_0$ in the initial condition $x\left(0\right)=q_0=x_1$or $x\left(0\right)=q_0=x_2$ translates the stationary solution $x=x_1$ or $x=x_2$to one of the stable or unstable solutions from the families S or U.
• (c), $D=0$: all solutions

  $$x\left(t;C\right)=-\frac{1}{a}\left(\frac{b}{2}+\frac{1}{t-C}\right)$$

  are unstable since they have a “moving” singularity at $t=C,$ are monotonically increasing or decreasing (depending on the sign of a) functions for $t\ne C$ are not bounded on the whole line (the latter means instability) and have the asymptote $x={-x}_0$. Stationary solutions $x={-x}_0$ are special in the sense that in any neighborhood there are infinitely many ’non-singular’ solutions $x\left(t;C\right)$.
• (d), $D<0$: all solutions

  $$x\left(t;C\right)=-x_0+\frac{\sqrt{-D}}{2a}tan\frac{\sqrt{-D}}{2}\left(t+C\right),D<0,$$

  are unstable since they have “movable” singularities for $t=t_n^{*}\left(C\right)=\frac{2}{\sqrt{-D}}\left(\frac{\pi }{2}\left(2n+1\right)-C\right),$$n=0,\pm 1,\pm 2,\cdots ,$are monotonically increasing or decreasing (depending on the sign of parameter a) functions for $t\ne t_n^{*}\left(C\right)$ because, for example,

  $$\frac{dx\left(t;C\right)}{dt}=\frac{(-D)}{4a}\frac{1}{{\left(cos\frac{\sqrt{-D}}{2}(t+C)\right)}^2}>0$$

  for $C>0$, $D<0$, and $a>0$ and are not bounded on the whole line (the latter means instability).

In case (d) there are no solutions that are special in the same sense as stationary solutions in case (b).

3.3. Equivalence Classes of D-Vectors and General Solutions to Autonomous Polynomial Equation Systems

Summarize all possible combinations of general solutions to DS (42) of dimension three identifying eight solution sets $D^{+++}$, $D^{++-}$, $D^{+-+}$, $D^{+--}$, $D^{-++}$, $D^{-+-}$, $D^{--+}$, and $D^{---}$ governed by all possible combinations of the nonzero discriminant signs. We identify as well 27 solution to DS (42) of dimension three as sets D^pmz, where p, m and z denote plus, minus, or zero and $pmz$ is any of their combinations including those with two or three repetitions, corresponding to all possible combinations of the zero and nonzero discriminant signs.

A set D^pmz with a fixed combination $pmz$ of signs is the equivalence class (defined in terms of S-vectors) in the sense that all its elements have the same S-vector S^pmz.

These solution sets are presented below in the convenient form of matrices (47) and (48).

In what follows we will consider the equations in system (42) with $a_{i1}\ne 0$, $i=1,2,3.$ and make use of the function families

$$g^{+}\left(t;D,C,a\right)=\frac{\sqrt{D}}{a}\frac{1}{-1+C{e}^{-\sqrt{D}t}},D>0,$$

$$g^0\left(t;C,a,b\right)=-\frac{1}{a}\left(\frac{b}{2}+\frac{1}{t-C}\right),D=0,$$

(43)

$$g^{-}\left(t;D,C,a\right)=\frac{\sqrt{-D}}{2a}tan\frac{\sqrt{-D}}{2}\left(t+C\right),D<0,$$

where C is an arbitrary constant, specifying componentwise the general solutions of (42) for the respective signs of the discriminants ${D=D}_i={a_{i2}}^2-4a_{i1}{a}_{i3}$ of polynomials $P_i$ in system (42), $i=1,2,3$.

3.4. Description of All Possible solution Combinations in Terms of Discriminants

In view of (43) and the relation between the coefficient matrix (17) and vectors (20), introduce the eight vector-functions

g⁺⁺⁺=[$g^{+}(t;D_1,C_1,a_{11}),g^{+}(t;D_2,C_2,a_{21}),g^{+}(t;D_3,C_3,a_{31})],$

⋯⋯,

g⁻⁻⁻=[$g^{-}(t;D_1,C_1,a_{11}),g^{-}(t;D_2,C_2,a_{21}),g^{-}(t;D_3,C_3,a_{31})],$corresponding to the respective solution sets D⁺⁺⁺, ..., D⁻⁻⁻, so that the upper indices are the S-vectors that indicate the respective equivalence classes; the solution sets are formed by the vector-functions

$${\mathbf{x}}^{+++}={\mathbf{x}}_1+{\mathbf{g}}^{+++},\mathrm{or}\left\{\begin{array}{c}{x\left(t;C_1\right)=x}_1+\frac{\sqrt{D_1}}{a_{11}}\frac{1}{-1+C_1{e}^{-\sqrt{D_1}t}},\\ {y\left(t;C_2\right)=y}_1+\frac{\sqrt{D_2}}{a_{21}}\frac{1}{-1+C_2{e}^{-\sqrt{D_2}t}},\\ {z\left(t;C_3\right)=z}_1+\frac{\sqrt{D_3}}{a_{31}}\frac{1}{-1+C_3{e}^{-\sqrt{D_3}t}},\end{array}\right.$$

(44)

$$\cdots \cdots$$

$${\mathbf{x}}^{---}={\mathbf{x}}_1+{\mathbf{g}}^{---},\mathrm{or}\left\{\begin{array}{c}{x\left(t;C_1\right)=x}_1+\frac{\sqrt{-D_1}}{2a}tan\frac{\sqrt{-D_1}}{2}\left(t+C_1\right),\\ {y\left(t;C_2\right)=y}_1+\frac{\sqrt{-D_2}}{2a}tan\frac{\sqrt{-D_2}}{2}\left(t+C_2\right),\\ {z\left(t;C_3\right)=z}_1+\frac{\sqrt{-D_3}}{2a}tan\frac{\sqrt{-D_3}}{2}\left(t+C_3\right);\end{array}\right.$$

(45)

here

$$\begin{array}{ccc}\hfill x_{1,2}\left(D_1\right)& =& -x_0\pm \frac{\sqrt{D_1}}{2a_{11}},x_0=\frac{a_{12}}{2a_{11}},\hfill \\ \hfill y_{1,2}\left(D_2\right)& =& -y_0\pm \frac{\sqrt{D_2}}{2a_{21}},y_0=\frac{a_{22}}{2a_{21}},\hfill \\ \hfill z_{1,2}\left(D_3\right)& =& -z_0\pm \frac{\sqrt{D_3}}{2a_{31}},z_0=\frac{a_{32}}{2a_{31}}\hfill \end{array}$$

(46)

$a_{i1}\ne 0$ and $D_i>0$, $i=1,2,3$.

Figure 1. Surface ${\mathbf{X}}_{D_1}^{-0-}\cup {\mathbf{X}}_{D_1}^{+0-}$ formed by the parametrized curves $\mathbf{x}(t;D_1)={\mathbf{x}}_{\mathbf{1}}+{\mathbf{g}}^{-0-}$ (black), $C_1=-1$, $-4<D_1<0$, and $\mathbf{x}(t;D_1)={\mathbf{x}}_{\mathbf{1}}+{\mathbf{g}}^{+0-},0\le t\le 1.3$ (green), $0<D_1<4$, with $a_{11}=a_{13}=a_{21}=1$ and $a_{22}=2$ given by (43)–(46) corresponding to the solutions from the sets $D^{-0-}$ and $D^{+0-}$.The explicit formulas (A22)–(A26) for the parametrized components are in Appendix A.

Figure 1. Surface ${\mathbf{X}}_{D_1}^{-0-}\cup {\mathbf{X}}_{D_1}^{+0-}$ formed by the parametrized curves $\mathbf{x}(t;D_1)={\mathbf{x}}_{\mathbf{1}}+{\mathbf{g}}^{-0-}$ (black), $C_1=-1$, $-4<D_1<0$, and $\mathbf{x}(t;D_1)={\mathbf{x}}_{\mathbf{1}}+{\mathbf{g}}^{+0-},0\le t\le 1.3$ (green), $0<D_1<4$, with $a_{11}=a_{13}=a_{21}=1$ and $a_{22}=2$ given by (43)–(46) corresponding to the solutions from the sets $D^{-0-}$ and $D^{+0-}$.The explicit formulas (A22)–(A26) for the parametrized components are in Appendix A.

The matrix

$$\begin{array}{c}\hfill {\mathbf{G}}^{(+-)}=\left[\begin{array}{ccc}{\mathbf{g}}^{+++}& {\mathbf{g}}^{+-+}& {\mathbf{g}}^{++-}\\ {\mathbf{g}}^{+--}& & {\mathbf{g}}^{---}\\ {\mathbf{g}}^{-++}& {\mathbf{g}}^{--+}& {\mathbf{g}}^{-+-}\end{array}\right]\end{array}$$

(47)

associated with DS (42) comprises all possible eight solution sets of system (42) corresponding to all possible combinations of the solution components having nonzero determinants and characterized by their S-vectors. Empty entries may be replaced by zeros.

Figure 2. Surface ${\mathbf{X}}_{D_1}^{-0-}\cup {\mathbf{X}}_{D_1}^{+0-}$ formed by the parametrized curves $\mathbf{x}(t;D_1)={\mathbf{x}}_{\mathbf{1}}+{\mathbf{g}}^{-0-}$ (right, black), $-4<D_1<0$, and $\mathbf{x}(t;D_1)={\mathbf{x}}_{\mathbf{1}}+{\mathbf{g}}^{+0-},0\le t\le 1.3$ (left, blue), $0<D_1<4$, with $a_{11}=a_{13}=a_{21}=1$ and $a_{22}=2$ given by (43)–(46) corresponding to the solutions from the sets $D^{-0-}$ and $D^{+0-}$. The explicit formulas (A16)–(A21) for the parametrized components are in Appendix A.

Figure 2. Surface ${\mathbf{X}}_{D_1}^{-0-}\cup {\mathbf{X}}_{D_1}^{+0-}$ formed by the parametrized curves $\mathbf{x}(t;D_1)={\mathbf{x}}_{\mathbf{1}}+{\mathbf{g}}^{-0-}$ (right, black), $-4<D_1<0$, and $\mathbf{x}(t;D_1)={\mathbf{x}}_{\mathbf{1}}+{\mathbf{g}}^{+0-},0\le t\le 1.3$ (left, blue), $0<D_1<4$, with $a_{11}=a_{13}=a_{21}=1$ and $a_{22}=2$ given by (43)–(46) corresponding to the solutions from the sets $D^{-0-}$ and $D^{+0-}$. The explicit formulas (A16)–(A21) for the parametrized components are in Appendix A.

The matrix

$$\begin{array}{c}\hfill {\mathbf{G}}^{(+-0)}=\left[\begin{array}{cccccc}{\mathbf{g}}^{+++}& {\mathbf{g}}^{++-}& {\mathbf{g}}^{+--}& {\mathbf{g}}^{+-+}& {\mathbf{g}}^{++0}& {\mathbf{g}}^{+-0}\\ {\mathbf{g}}^{-++}& {\mathbf{g}}^{-+-}& {\mathbf{g}}^{---}& {\mathbf{g}}^{--+}& {\mathbf{g}}^{-+0}& {\mathbf{g}}^{--0}\\ {\mathbf{g}}^{00+}& {\mathbf{g}}^{00-}& {\mathbf{g}}^{000}& & & \\ {\mathbf{g}}^{+00}& {\mathbf{g}}^{-00}& {\mathbf{g}}^{\mathbf{0}+0}& & & \\ {\mathbf{g}}^{0-0}& {\mathbf{g}}^{+0+}& {\mathbf{g}}^{\mathbf{0}++}& & & \\ {\mathbf{g}}^{-0+}& {\mathbf{g}}^{0+-}& {\mathbf{g}}^{+0-}& {\mathbf{g}}^{0-+}& {\mathbf{g}}^{-0-}& {\mathbf{g}}^{0--}\end{array}\right]\end{array}$$

(48)

comprises all possible 27 solution sets of system (42) corresponding to all possible combinations with nonzero and zero determinants and characterized by their S-vectors. Empty entries in (48) may be replaced by zeros.

Figure 3. Surface ${\mathbf{X}}_{D_1}^{-0-}\cup {\mathbf{X}}_{D_1}^{+0-}$ formed by the parametrized curves $\mathbf{x}(t;D_1)={\mathbf{x}}_{\mathbf{1}}+{\mathbf{g}}^{-0-}$ (black), $C_1=-1.3$, $-4<D_1<0$, and $\mathbf{x}(t;D_1)={\mathbf{x}}_{\mathbf{1}}+{\mathbf{g}}^{+0-},0\le t\le 1.3$ (green), $0<D_1<4$, with $a_{11}=a_{13}=a_{21}=1$ and $a_{22}=2$ given by (43)–(46) corresponding to the solutions from the sets $D^{-0-}$ and $D^{+0-}$.

Figure 3. Surface ${\mathbf{X}}_{D_1}^{-0-}\cup {\mathbf{X}}_{D_1}^{+0-}$ formed by the parametrized curves $\mathbf{x}(t;D_1)={\mathbf{x}}_{\mathbf{1}}+{\mathbf{g}}^{-0-}$ (black), $C_1=-1.3$, $-4<D_1<0$, and $\mathbf{x}(t;D_1)={\mathbf{x}}_{\mathbf{1}}+{\mathbf{g}}^{+0-},0\le t\le 1.3$ (green), $0<D_1<4$, with $a_{11}=a_{13}=a_{21}=1$ and $a_{22}=2$ given by (43)–(46) corresponding to the solutions from the sets $D^{-0-}$ and $D^{+0-}$.

3.5. Analysis of Solutions to Cauchy Problems

Solutions to the Cauchy problems, e.g. to $\dot{x}=P_1\left(x\right)$, $x\left(0\right)=q_0,$ is

$$x\left(t;C_1^{0,+}\right)=x_1+g^{+}\left(t;D_1,C_1^{0,+},a_{11}\right)=x_1+\frac{\sqrt{D_1}}{a_{11}}\frac{1}{-1+C_1^{0,+}{e}^{-\sqrt{D_1}t}},$$

$$C_1^{0,+}=1+\frac{a_{11}}{\sqrt{D_1}}\frac{1}{q_0-x_1}$$

for $D_1>0$ and

$$x\left(t;C_1^{0,-}\right)=-x_0+g^{-}\left(t;D_1,C_1^{0,-},a_{11}\right)=-x_0+\frac{\sqrt{-D_1}}{2a_{11}}tan\frac{\sqrt{-D_1}}{2}\left(t+C_1^{0,-}\right),$$

$$C_1^{0,-}=\frac{2}{\sqrt{-D_1}}arctan\left[\frac{2a_{11}}{\sqrt{-D_1}}\left(q_0+x_0\right)\right],C_1^{0,-}\ne \frac{\pi }{2}+\pi n,n=0,\pm 1,\pm 2,\cdots ,$$

for $D_1<0,$so that the solution to $\dot{\mathbf{x}}=\mathbf{P}\left(\mathbf{x}\right),\mathbf{x}\left(0\right)={\mathbf{q}}^0,{\mathbf{q}}^0=\left[q_1^0,q_2^0,q_3^0\right],$ can be written for different solution sets D⁺⁺⁺, ..., D⁻⁻⁻, as

$$\mathbf{x}={\mathbf{x}}_{\mathbf{1}}+{\mathbf{g}}^{+++},\mathrm{or}\left\{\begin{array}{c}x\left(t;C_1^{0,+}\right)=x_1+\frac{\sqrt{D_1}}{a_{11}}\frac{1}{-1+C_1^{0,+}{e}^{-\sqrt{D_1}t}},\\ y\left(t;C_2^{0,+}\right)=y_1+\frac{\sqrt{D_2}}{a_{21}}\frac{1}{-1+C_2^{0,+}{e}^{-\sqrt{D_2}t}},\\ z\left(t;C_3^{0,+}\right)=z_1+\frac{\sqrt{D_3}}{a_{31}}\frac{1}{-1+C_3^{0,+}{e}^{-\sqrt{D_3}t}},\end{array}\right.$$

$C_i^{0,+}=1+\frac{a_{i1}}{\sqrt{D_i}}\frac{1}{q_i^0-x^{\left(i\right)}},$ for D⁺⁺⁺ etc, where $i=1,2,3,$${x^{\left(1\right)}=x}_1,{x^{\left(2\right)}=y}_1,{x^{\left(3\right)}=z}_1.$

Figure 4. Surface ${\mathbf{X}}_{D_1}^{0--}\cup {\mathbf{X}}_{D_1}^{+0-}$ formed by the parametrized curves $\mathbf{x}(t;D_1)={\mathbf{x}}_{\mathbf{1}}+{\mathbf{g}}^{\mathbf{0}--}$ (black), $C_1=-1$, $-4<D_1<0$, and $\mathbf{x}(t;D_1)={\mathbf{x}}_{\mathbf{1}}+{\mathbf{g}}^{+0-},0\le t\le 1.3$ (green), $0<D_1<4$, with $a_{11}=a_{13}=a_{21}=1$ and $a_{22}=2$ given by (43)–(46) corresponding to the solutions from the sets $D^{0--}$ and $D^{+0-}$.

Figure 4. Surface ${\mathbf{X}}_{D_1}^{0--}\cup {\mathbf{X}}_{D_1}^{+0-}$ formed by the parametrized curves $\mathbf{x}(t;D_1)={\mathbf{x}}_{\mathbf{1}}+{\mathbf{g}}^{\mathbf{0}--}$ (black), $C_1=-1$, $-4<D_1<0$, and $\mathbf{x}(t;D_1)={\mathbf{x}}_{\mathbf{1}}+{\mathbf{g}}^{+0-},0\le t\le 1.3$ (green), $0<D_1<4$, with $a_{11}=a_{13}=a_{21}=1$ and $a_{22}=2$ given by (43)–(46) corresponding to the solutions from the sets $D^{0--}$ and $D^{+0-}$.

4. Analysis of Bifurcations

Families of vector-functions (44) ${\mathbf{x}}^{...}(t;\mathbf{C};\mathbf{D})$ of the type (44) with $\mathbf{C}=\left[C_1,C_2,C_3\right]\in \mathrm{C}\subset$R³ and $\mathbf{D}=\left[D_1,D_2,D_3\right]\in \mathrm{D}\subset$R³, where C, D, and ... denote, respectively, certain 3D sets (or domains) and any combination of three signs +, − and 0, specify each a 3D-surface ${\mathbf{X}}_{\Psi }^{...}$ formed by the corresponding curves ${\mathbf{x}}^{...}(t;.)$ parametrized by $t\in \left[T_0,T_1\right],T_0\ge 0$ and constructed for different values of one (the rest being fixed) real parameter $\Psi =D_j$ or $\Psi =C_j$, $j=1,2,3$: ${\mathbf{X}}_{\Psi }^{...}=\left\{{\mathbf{x}}^{...}(t;\Psi ),t\in \left[T_0,T_1\right],\Psi \in \left[{\Psi }_0,{\Psi }_1\right]\right\}$.

Figure 5. Surface ${\mathbf{X}}_{D_1}^{-0-}$ formed by the parametrized curves $\mathbf{x}(t;D_1)={\mathbf{x}}_{\mathbf{1}}+{\mathbf{g}}^{-0-}$, $0\le t\le 1.45,-2\le D_1<0$, $a_{11}=a_{13}=1$ given by (43)–(46) corresponding to the solutions from the set $D^{-0-}$.

Figure 5. Surface ${\mathbf{X}}_{D_1}^{-0-}$ formed by the parametrized curves $\mathbf{x}(t;D_1)={\mathbf{x}}_{\mathbf{1}}+{\mathbf{g}}^{-0-}$, $0\le t\le 1.45,-2\le D_1<0$, $a_{11}=a_{13}=1$ given by (43)–(46) corresponding to the solutions from the set $D^{-0-}$.

We investigate particularly the occurrence of bifurcations for the combined surfaces ${\mathbf{X}}_{\Psi =D_1}^{-0-}\cup {\mathbf{X}}_{\Psi =D_1}^{+-0}$ when the discriminant of one (or two) component varies.

Figure 6. Enlarged fragment of the projection on the $(x,y)$-plane of the surface presented in Figure 4 close to the point of bifuraction associated with the stationary solution $y=-1$ of the second component in system (51) (graphs of the parametrized curves $\mathbf{x}\left(t\right)={\mathbf{x}}_{\mathbf{1}}+{\mathbf{g}}^{\mathbf{0}--}$ (black), $C_1=-1$, $-4<D_1<0$, and $\mathbf{x}\left(t\right)={\mathbf{x}}_{\mathbf{1}}+{\mathbf{g}}^{+0-},0\le t\le 1.3$ (green), $0<D_1<4$, with $a_{11}=a_{13}=a_{21}=1$ and $a_{22}=2$ given by (43)–(46) corresponding to the solutions from the sets $D^{0--}$ and $D^{+0-}$).

Figure 6. Enlarged fragment of the projection on the $(x,y)$-plane of the surface presented in Figure 4 close to the point of bifuraction associated with the stationary solution $y=-1$ of the second component in system (51) (graphs of the parametrized curves $\mathbf{x}\left(t\right)={\mathbf{x}}_{\mathbf{1}}+{\mathbf{g}}^{\mathbf{0}--}$ (black), $C_1=-1$, $-4<D_1<0$, and $\mathbf{x}\left(t\right)={\mathbf{x}}_{\mathbf{1}}+{\mathbf{g}}^{+0-},0\le t\le 1.3$ (green), $0<D_1<4$, with $a_{11}=a_{13}=a_{21}=1$ and $a_{22}=2$ given by (43)–(46) corresponding to the solutions from the sets $D^{0--}$ and $D^{+0-}$).

Figure 7. Graphs of the parametrized curves $\mathbf{x}(t;D_1)={\mathbf{x}}_{\mathbf{1}}+{\mathbf{g}}^{-0-}$ (black), $C_1=-1$, $C_2=-0.04$, $-4<D_1<0$, and $\mathbf{x}(t;D_1)={\mathbf{x}}_{\mathbf{1}}+{\mathbf{g}}^{+0-},0\le t\le 1.4$ (yellow), $0<D_1<4$, with $a_{11}=a_{13}=a_{21}=1$ and $a_{22}=2$ given by (43)–(46) corresponding to the solutions from the sets $D^{-0-}$ and $D^{+0-}$.

Figure 7. Graphs of the parametrized curves $\mathbf{x}(t;D_1)={\mathbf{x}}_{\mathbf{1}}+{\mathbf{g}}^{-0-}$ (black), $C_1=-1$, $C_2=-0.04$, $-4<D_1<0$, and $\mathbf{x}(t;D_1)={\mathbf{x}}_{\mathbf{1}}+{\mathbf{g}}^{+0-},0\le t\le 1.4$ (yellow), $0<D_1<4$, with $a_{11}=a_{13}=a_{21}=1$ and $a_{22}=2$ given by (43)–(46) corresponding to the solutions from the sets $D^{-0-}$ and $D^{+0-}$.

We have, according to (43),

$$\begin{array}{ccc}\hfill {lim}_{D_1\to +0}\left[-x_0\pm \frac{\sqrt{D_1}}{2a_{11}}+g^{+}\left(t;D_1,C,a_{11}\right)\right]& =& -x_0,\hfill \end{array}$$

(49)

$$\begin{array}{ccc}\hfill {lim}_{D_1\to -0}\left[-x_0+g^{-}\left(t;D_1,C,a_{11}\right)\right]& =& -x_0,\hfill \end{array}$$

(50)

uniformly with respect to $t\in \left[t_0,T_0\right],t_0\ge 0$. (49) and (50) mean that there is no bifurcation associated with the transition $D^{-0-}\to$$D^{+0-}$ as illustrated by Figure 2.

Figure 8. Enlarged fragment of the surface close to the point of bifuraction associated with the stationary solution $y=-1$ of the second and third components in system (51) (graphs of the parametrized curves $\mathbf{x}(t;D_1)={\mathbf{x}}_{\mathbf{1}}+{\mathbf{g}}^{-00}$ (black), $C_1=-1$, $-4<D_1<0$, $C_2=1.6$, $C_3=2$, and $\mathbf{x}(t;D_1)={\mathbf{x}}_{\mathbf{1}}+{\mathbf{g}}^{+00},0\le t\le 1.5$ (yellow), $0<D_1<4$, with $a_{11}=a_{13}=a_{21}=1$ and $a_{22}=2$ given by (43)–(46) corresponding to the solutions from the sets $D^{-00}$ and $D^{+00}$).

Figure 8. Enlarged fragment of the surface close to the point of bifuraction associated with the stationary solution $y=-1$ of the second and third components in system (51) (graphs of the parametrized curves $\mathbf{x}(t;D_1)={\mathbf{x}}_{\mathbf{1}}+{\mathbf{g}}^{-00}$ (black), $C_1=-1$, $-4<D_1<0$, $C_2=1.6$, $C_3=2$, and $\mathbf{x}(t;D_1)={\mathbf{x}}_{\mathbf{1}}+{\mathbf{g}}^{+00},0\le t\le 1.5$ (yellow), $0<D_1<4$, with $a_{11}=a_{13}=a_{21}=1$ and $a_{22}=2$ given by (43)–(46) corresponding to the solutions from the sets $D^{-00}$ and $D^{+00}$).

Figure 2, Figure 3, Figure 4, Figure 5, Figure 6, Figure 7 and Figure 8 illustrate the formation of a bifuraction of the type twisted fold (a ’cusp’). This bifuraction occurs with variation of discriminant $D_1$ of a single (first or second) solution component within single family $D^{-0-}$ or $D^{0--}$ at certain values of constant $C_1$ (namely, at $C_1=-1,-1.3$) when the first and third (or the second and third) components (independent of variable $D_1$) are bounded in a closed interval $t\in \left[0,T_0\right]$ where the parametrization is applied.

A twisted fold manifests itself as a ’rotation’ caused by a varied $D_1\in (-4,0)$ of a 3D-surface formed by the curves ${\mathbf{x}}^{-0-}(t;D_1)$ or ${\mathbf{x}}^{0--}(t;D_1)$ parametrized by $t\in \left[T_0,T_1\right],T_0\ge 0$, with fixed vector $\mathbf{C}=\left[C_1,C_2,C_3\right]$ and $\mathbf{D}=\left[D_1,D_2,D_3\right]$ with constant $D_2$ and $D_3$ and varied $D_1$. The ’rotation’ takes place around a point of equilibrium where all the curves ${\mathbf{x}}^{-0-}(t;D_1)$ or ${\mathbf{x}}^{0--}(t;D_1)$ with different $D_1\in (-4,0)$ merge at a certain $t\in \left[T_0,T_1\right]$.

Figure 9. Graph of the initial condition $x\left(0\right)=\frac{\sqrt{D_1+4}}{2}-\frac{\sqrt{-D_1}}{2}tan\frac{\sqrt{-D_1}}{2}$ vs $D_1\in \left(-4,0\right]$ of the first component of the solution (A22), (A23) from the set $D^{-0-}$ presented in Figure 1.

Figure 9. Graph of the initial condition $x\left(0\right)=\frac{\sqrt{D_1+4}}{2}-\frac{\sqrt{-D_1}}{2}tan\frac{\sqrt{-D_1}}{2}$ vs $D_1\in \left(-4,0\right]$ of the first component of the solution (A22), (A23) from the set $D^{-0-}$ presented in Figure 1.

Figure 10. Graphs of the parametrized curves $\mathbf{x}\left(t\right)={\mathbf{x}}_{\mathbf{1}}+{\mathbf{g}}^{---}$, $0\le t\le 2,-4\le D_1,D_3<0$, $a_{11}=a_{13}=1$, $a_{22}=a_{32}=2$, $C_j=-1$, $j=1,2,3$, given by (43)–(46) corresponding to the solutions from the set $D^{---}$ with two variable components.

Figure 10. Graphs of the parametrized curves $\mathbf{x}\left(t\right)={\mathbf{x}}_{\mathbf{1}}+{\mathbf{g}}^{---}$, $0\le t\le 2,-4\le D_1,D_3<0$, $a_{11}=a_{13}=1$, $a_{22}=a_{32}=2$, $C_j=-1$, $j=1,2,3$, given by (43)–(46) corresponding to the solutions from the set $D^{---}$ with two variable components.

To understand the nature of the observed bifurcation consider system (42) corresponding to families $D^{\mp 0-}$ and the cases involving the specific values of parameters (coefficients of the trinomials in (42)) illustrated by Figure 3 and Figure 5

$$\left\{\begin{array}{c}\dot{x}={(x-\frac{1}{2}\sqrt{D_1+4})}^2-\frac{1}{4}{D}_1,\\ \dot{y}={(y+1)}^2,\\ \dot{z}={(z-\frac{1}{2})}^2+\frac{3}{4},\end{array}\right.$$

(51)

where $D_1<0$ and $D_1>0$ are, respectively, for $D^{-0-}$ and $D^{+0-}$. The second component in (51) has the degenerate stationary solution $y=-1$ while the first and third components preserve the (positive) sign when parameter $D_1$ varies. Figure 9 displays the initial condition $x\left(0\right)=\frac{\sqrt{D_1+4}}{2}-\frac{\sqrt{-D_1}}{2}tan\frac{\sqrt{-D_1}}{2}$ vs $D_1\in \left(-4,0\right]$ of the first component of the solution (A22), (A23) from the set$D^{-0-}$ presented in Figure 1. We see that the initial condition changes sign as $D_1<0$ increases which causes twist of the curves, contrary to the surface ${\mathbf{X}}_{D_1}^{-0-}\cup {\mathbf{X}}_{D_1}^{+0-}$ in Figure 2 where the solutions from the sets $D^{-0-}$ and $D^{+0-}$ are given by (A16)–(A21) and the initial condition $x\left(0\right)=\frac{\sqrt{D_1+4}}{2}$ for the first component considered vs $D_1\in \left(-4,0\right]$ does not change sign and the solution curves in both sets $D^{-0-}$ and $D^{+0-}$ preserve the same shape.

Figure 11. Graphs of the parametrized curves $\mathbf{x}\left(t\right)={\mathbf{x}}_{\mathbf{1}}+{\mathbf{g}}^{-+-}$ (right) and $\mathbf{x}\left(t\right)={\mathbf{x}}_{\mathbf{1}}+{\mathbf{g}}^{++-},0\le t\le 1.3,-4\le D_1\le 6$, $a_{11}=a_{13}=1$ given by (43)–(46) corresponding to the solutions from the sets $D^{-+-}$ and $D^{++-}$.

Figure 11. Graphs of the parametrized curves $\mathbf{x}\left(t\right)={\mathbf{x}}_{\mathbf{1}}+{\mathbf{g}}^{-+-}$ (right) and $\mathbf{x}\left(t\right)={\mathbf{x}}_{\mathbf{1}}+{\mathbf{g}}^{++-},0\le t\le 1.3,-4\le D_1\le 6$, $a_{11}=a_{13}=1$ given by (43)–(46) corresponding to the solutions from the sets $D^{-+-}$ and $D^{++-}$.

5. Conclusion

We have extended the discriminant criterion applied originally to two-phase quadratic polynomial DSs to three-phase DSs investigated in terms of their coefficient matrices. Based in this criterion and the established integrability in elementary functions for the considered class of autonomous polynomial DSs, we have proposed a general approach to qualitative analysis of these DSs. Specific classes of D- and S-vectors have been introduced which enable one to classify the DS solutions depending on the discriminant signs. Complete description of the symmetry relations inherent to the DS coefficient matrices has been performed using the discriminant criterion. The DS solution properties have been investigated in terms of two-parameter surfaces. The behaviour and bifurcations of solutions to three-dimensional quadratic polynomial DSs have been considered. The occurrence of bifurcations of the type twisted fold has been discovered on the basis and within the frames of the developed DS qualitative theory.

Based on the results of the study, we have made the following conclusion: whether the DS coefficient matrix belongs to a particular equivalence class determines the presence and nature (type) of singular points and bifurcations of the corresponding polynomial DS.

As far as the application of the developed approach and technique to regulating and monitoring the waterflooding and the development of oil or gas fields is concerned, the proposed approach enables one to simulate the consequences of instable behavior of the displacement front. Using the results of this work, engineers and practitioners will be able to predict natural abrupt variations in water saturation and dependent well exploatation parameters using the analysis of the growth model employing three-phase quadratic polynomial DSs.

The results and findings reported in this study make it possibile to carry out systematic optimization of non-stationary waterflooding and the prospect of enhanced oil recovery, as well as to control reduction in the amount of water that is not efficiently injected and withdrawn.

Appendix A.1. Autonomous Polynomial Dynamical Systems Integrable in Elementary Functions

A polynomial $P_i\left(x\right)$ of degree $N_i\ge 2$ ($i=1,2,\cdots ,n$) with real coefficients can be uniquely factorized and represented as

$$P_i\left(x\right)={\Pi }_{s=1}^{m_i}{(x-x_s^{\left(i\right)})}^{k_s^{\left(i\right)}}{\Pi }_{t=1}^{n_i}{\left(r_t^{\left(i\right)}\right)}^{p_t^{\left(i\right)}},$$

(A1)

where ${\displaystyle \sum _{s=1}^{m_i}{k}_s^{\left(i\right)}+2\sum _{t=1}^{n_i}{p}_t^{\left(i\right)}=N_i}$, $m_i,n_i\ge 0$, and

$$r_s^{\left(i\right)}\left(x\right)=a_{s,1}^{\left(i\right)}{x}^2+a_{s,2}^{\left(i\right)}x+a_{s,3}^{\left(i\right)}=a_{s,1}^{\left(i\right)}\left[{(x+x_{s,0}^{\left(i\right)})}^2-\frac{D_s^{\left(i\right)}}{4{\left(a_{s,1}^{\left(i\right)}\right)}^2}\right](a_{s,1}^{\left(i\right)}\ne 0)$$

(A2)

are irreducible quadratic trinomials with ${\displaystyle x_{s,0}^{\left(i\right)}=\frac{a_{s,2}^{\left(i\right)}}{2a_{s,1}^{\left(i\right)}}}$ and the discriminants

$$D_s^{\left(i\right)}={\left(a_{s,2}^{\left(i\right)}\right)}^2-4a_{s,1}^{\left(i\right)}{a}_{s,3}^{\left(i\right)}<0,s=1,2,\cdots ,n_i,i=1,2,\cdots ,n.$$

(A3)

A rational function $Q_i\left(x\right)=P_i^{-1}\left(x\right)$ can be uniquely decomposed [17] into partial fractions,

$$Q_i\left(x\right)=\sum _{s=1}^{m_i}\sum _{j=1}^{k_s^{\left(i\right)}}\frac{A_{sj}^{\left(i\right)}}{{(x-x_s^{\left(i\right)})}^j}+\sum _{t=1}^{n_i}\sum _{j=1}^{p_t^{\left(i\right)}}\frac{B_{tj}^{\left(i\right)}x+C_{tj}^{\left(i\right)}}{{\left(r_t^{\left(i\right)}\right)}^j}$$

(A4)

where $A_{sj}^{\left(i\right)}$, $B_{tj}^{\left(i\right)}$, and $C_{tj}^{\left(i\right)}$ are uniquely determined real numbers.

Lemma A1.

Each term in (A4) is integrable in elementary functions

Proof. 

Verification of the statement is obvious for the terms entering the first double sum in (A4).

Consider the terms entering the second double sum in (A4). To prove the required integrability it is sufficient to calculate two indefinite integrals

$$I_j^{\left(1\right)}\left(x\right)=\int \frac{xdx}{r^j\left(x\right)}and{I}_j^{\left(2\right)}\left(x\right)=\int \frac{dx}{r^j\left(x\right)},r\left(x\right)=a{x}^2+bx+c,$$

(A5)

with integer $j\ge 1$ and $D=b^2-4ac<0$. For the first integral, we have

$$\begin{array}{ccc}\hfill I_j^{\left(1\right)}\left(x\right)& =& \int \frac{xdx}{{\left[a\left({(x+x_0)}^2-\frac{D}{4a^2}\right)\right]}^j}=\hfill \\ & =& -\frac{1}{2a^j(j-1)}\frac{1}{{\left[{(x+x_0)}^2-\frac{D}{4a^2}\right]}^{j-1}}-\frac{x_0}{a^j}{I}_j^{\left(2\right)}\left(x\right)\hfill \end{array}$$

(A6)

for $j\ge 2$ and

$$\begin{array}{ccc}\hfill I_1^{\left(1\right)}\left(x\right)& =& \frac{1}{a}\int \frac{xdx}{\left[{(x+x_0)}^2-\frac{D}{4a^2}\right]}=\hfill \end{array}$$

(A7)

$$\begin{array}{ccc}& =& \frac{1}{2a}ln\left[{(x+x_0)}^2-\frac{D}{4a^2}\right]-\frac{x_0}{a}{I}_1^{\left(2\right)}\left(x\right)=\hfill \\ & =& \frac{1}{2a}ln\left[{(x+x_0)}^2-\frac{D}{4a^2}\right]-\frac{2x_0}{\sqrt{-D}}{tan}^{-1}\left[\frac{2a(x+x_0)}{\sqrt{-D}}\right]+C\hfill \end{array}$$

(A8)

for $j=1$. For the second integral, we have for $j\ge 2$,

$$\begin{array}{c}\hfill I_j^{\left(2\right)}\left(x\right)=-\frac{2ax+b}{(j-1)D{r}^{j-1}\left(x\right)}-\frac{2a(2j-3)}{(j-1)D}{I}_{j-1}^{\left(2\right)}\left(x\right).\end{array}$$

(A9)

This is actually a first-order inhomogeneous difference equation which can be written in a compact form

$$\begin{array}{c}\hfill I_j^{\left(2\right)}={\alpha }_{j-1}{I}_{j-1}^{\left(2\right)}+{\beta }_{j-1},j=2,3,\cdots ,I_1^{\left(2\right)}=\frac{2a}{\sqrt{-D}}{tan}^{-1}\left[\frac{2a(x+x_0)}{\sqrt{-D}}\right],\end{array}$$

(A10)

The solution to (A10) is

$$\begin{array}{c}\hfill I_j^{\left(2\right)}=\sum _{s=1}^j{\beta }_{j-s}{\Pi }_{t=1}^s{\alpha }_{j-t+1}+I_1^{\left(2\right)}{\Pi }_{t=1}^j{\alpha }_t,\end{array}$$

(A11)

where

$$\begin{array}{c}\hfill {\alpha }_{j-1}=-\frac{2ax+b}{(j-1)D{r}^{j-1}\left(x\right)},{\beta }_{j-1}=-\frac{2a(2j-3)}{(j-1)D}.\end{array}$$

(A12)

Thus, both $I_j^{\left(1\right)}$ and $I_j^{\left(2\right)}$ are expressed in elementary functions using (A11) and (A12). The indefinite integral

$$\int Q_i\left(x\right)dx=\sum _{s=1}^{m_i}\sum _{j=1}^{k_s^{\left(i\right)}}{A}_{sj}^{\left(i\right)}{q}_{j,s}^{\left(i\right)}\left(x\right)+T_i\left(x\right),$$

(A13)

where

$$q_{j,s}^{\left(i\right)}\left(x\right)=\left\{\begin{array}{ccc}\hfill & -\frac{1}{(j-1){(x-x_s^{\left(i\right)})}^{j-1}},\hfill & \hfill j\ge 2,\\ \hfill & ln\left|x-x_s^{\left(i\right)}\right|,\hfill & \hfill j=1,\end{array}\right.$$

(A14)

and $T_i\left(x\right)$ is a linear combination of $I_j^{\left(1\right)}$ and $I_j^{\left(2\right)}$ which proves the lemma. □

Theorem A1.

DSs (41) (and (42)) are integrable in elementary functions.

Proof. 

Each particular row equation in (41) written as

$$\begin{array}{ccc}\hfill \frac{dx}{dt}& =& P_i\left(x\right),P_i\left(x\right)=\sum _{j=1}^{N_i+1}{a}_{ji}{x}^{N_i-j+1}=\hfill \\ \hfill & =& a_{1i}{x}^{N_i}+a_{2i}{x}^{N_i-1}+\cdots +a_{N_i,i}x+a_{N_i+1,i}(i=1,2,\cdots ,n)\hfill \end{array}$$

(A15)

yields a rational function $Q_i\left(x\right)=P_i^{-1}\left(x\right)$ which is uniquely decomposed into partial fractions (A4) and is therefore integrable in elementary functions according to Lemma A1.

□

Appendix A.2. Examples with Bifurcations

Parametrized curves $\mathbf{x}(t;D_1)={\mathbf{x}}_{\mathbf{1}}+{\mathbf{g}}^{-0-}$, $-4<D_1<0$, and $\mathbf{x}(t;D_1)={\mathbf{x}}_{\mathbf{1}}+{\mathbf{g}}^{+0-}$, $0<D_1<4$ for $0\le t\le 1.3$, forming surface ${\mathbf{X}}_{D_1}^{-0-}\cup {\mathbf{X}}_{D_1}^{+0-}$ corresponding to the solutions from the sets $D^{-0-}$ and $D^{+0-}$ displayed in Figure 2

$${\mathbf{x}}^{-0-}={\mathbf{x}}_1+{\mathbf{g}}^{-0-}=$$

(A16)

$$\left\{\begin{array}{c}x\left(t;C_1=0\right)=-x_0+g^{-}\left(t;D_1,C_1=0,1\right)\\ y\left(t;C_2\right)=g^0\left(t;C_2=-1,a_{21},a_{22}\right),\\ z\left(t;C_3\right)=-z_0+g^{-}\left(t;D_3=-3,C_3=0,1\right),\end{array}\right.=$$

(A17)

$$=\left\{\begin{array}{c}\frac{\sqrt{D_1+4a_{11}{a}_{13}}}{2a_{11}}+\frac{\sqrt{-D_1}}{2a_{11}}tan\left(\frac{\sqrt{-D_1}}{2}t\right),\\ -\frac{1}{a_{21}}\left(\frac{a_{22}}{2}+\frac{1}{t-C_2}\right),\\ \frac{\sqrt{D_3+4a_{31}{a}_{33}}}{2a_{31}}+\frac{\sqrt{-D_3}}{2a_{31}}tan\left(\frac{\sqrt{-D_3}}{2}t\right),\end{array}\right.=\left\{\begin{array}{c}\frac{1}{2}\sqrt{D_1+4}+\frac{\sqrt{-D_1}}{2}tan\left(\frac{\sqrt{-D_1}}{2}t\right),\\ -1-\frac{1}{t+1},\\ \frac{1}{2}+\frac{\sqrt{3}}{2}tan\left(\frac{\sqrt{3}}{2}t\right).\end{array}\right.$$

(A18)

$${\mathbf{x}}^{+0-}={\mathbf{x}}_1+{\mathbf{g}}^{+0-}=$$

(A19)

$$\left\{\begin{array}{c}x\left(t;C_1=-1\right)=x_1+\frac{\sqrt{D_1}}{a_{11}}\frac{1}{-1+C_1{e}^{-\sqrt{D_1}t}},\\ y\left(t;C_2\right)=g^0\left(t;C_2=-1,a_{21},a_{22}\right),\\ z\left(t;C_3\right)=-z_0+g^{-}\left(t;D_3=-3,C_3=0,1\right),\end{array}\right.=$$

(A20)

$$=\left\{\begin{array}{c}\frac{\sqrt{D_1+4a_{11}{a}_{13}}}{2a_{11}}+\frac{\sqrt{D_1}}{2a_{11}}+\frac{\sqrt{D_1}}{a_{11}}\frac{1}{-1+C_1{e}^{-\sqrt{D_1}t}},\\ -\frac{1}{a_{21}}\left(\frac{a_{22}}{2}+\frac{1}{t-C_2}\right),\\ \frac{\sqrt{D_3+4a_{31}{a}_{33}}}{2a_{31}}+\frac{\sqrt{-D_3}}{2a_{31}}tan\left(\frac{\sqrt{-D_3}}{2}t\right),\end{array}\right.=\left\{\begin{array}{c}\frac{1}{2}\sqrt{D_1+4}-\frac{\sqrt{D_1}}{1+e^{-\sqrt{D_1}t}},\\ -1-\frac{1}{t+1},\\ \frac{1}{2}+\frac{\sqrt{3}}{2}tan\left(\frac{\sqrt{3}}{2}t\right).\end{array}\right.$$

(A21)

Parametrized curves $\mathbf{x}(t;D_1)={\mathbf{x}}_{\mathbf{1}}+{\mathbf{g}}^{-0-}$, $C_1=-1$, $-4<D_1<0$, and $\mathbf{x}(t;D_1)={\mathbf{x}}_{\mathbf{1}}+{\mathbf{g}}^{+0-}$, $0<D_1<4$, for $0\le t\le 1.3$ with $a_{11}=a_{13}=a_{21}=1$ and $a_{22}=2$ forming surface ${\mathbf{X}}_{D_1}^{-0-}\cup {\mathbf{X}}_{D_1}^{+0-}$ corresponding to the solutions from the sets $D^{-0-}$ and $D^{+0-}$ displayed in Figure 1

$${\mathbf{x}}^{-0-}={\mathbf{x}}_1+{\mathbf{g}}^{-0-}=$$

(A22)

$$\left\{\begin{array}{c}x\left(t;C_1=-1\right)=-x_0+g^{-}\left(t;D_1,C_1=-1,1\right)\\ y\left(t;C_2\right)=g^0\left(t;C_2=-1,a_{21},a_{22}\right),\\ z\left(t;C_3\right)=-z_0+g^{-}\left(t;D_3=-3,C_3=0,1\right),\end{array}\right.=$$

(A23)

$$=\left\{\begin{array}{c}\frac{\sqrt{D_1+4a_{11}{a}_{13}}}{2a_{11}}+\frac{\sqrt{-D_1}}{2a_{11}}tan\left(\frac{\sqrt{-D_1}}{2}(t-1)\right),\\ -\frac{1}{a_{21}}\left(\frac{a_{22}}{2}+\frac{1}{t-C_2}\right),\\ \frac{\sqrt{D_3+4a_{31}{a}_{33}}}{2a_{31}}+\frac{\sqrt{-D_3}}{2a_{31}}tan\left(\frac{\sqrt{-D_3}}{2}t\right),\end{array}\right.=\left\{\begin{array}{c}\frac{1}{2}\sqrt{D_1+4}+\frac{\sqrt{-D_1}}{2}tan\left(\frac{\sqrt{-D_1}}{2}(t-1)\right),\\ -1-\frac{1}{t+1},\\ \frac{1}{2}+\frac{\sqrt{3}}{2}tan\left(\frac{\sqrt{3}}{2}t\right),\end{array}\right.$$

$${\mathbf{x}}^{+0-}={\mathbf{x}}_1+{\mathbf{g}}^{+0-}=$$

(A24)

$$=\left\{\begin{array}{c}x\left(t;C_1=-1\right)=x_1+\frac{\sqrt{D_1}}{a_{11}}\frac{1}{-1+C_1{e}^{-\sqrt{D_1}t}},\\ y\left(t;C_2\right)=g^0\left(t;C_2=-1,a_{21},a_{22}\right),\\ z\left(t;C_3\right)=-z_0+g^{-}\left(t;D_3=-3,C_3=0,1\right),\end{array}\right.=$$

(A25)

$$=\left\{\begin{array}{c}\frac{\sqrt{D_1+4a_{11}{a}_{13}}}{2a_{11}}+\frac{\sqrt{D_1}}{2a_{11}}+\frac{\sqrt{D_1}}{a_{11}}\frac{1}{-1+C_1{e}^{-\sqrt{D_1}t}},\\ -\frac{1}{a_{21}}\left(\frac{a_{22}}{2}+\frac{1}{t-C_2}\right),\\ \frac{\sqrt{D_3+4a_{31}{a}_{33}}}{2a_{31}}+\frac{\sqrt{-D_3}}{2a_{31}}tan\left(\frac{\sqrt{-D_3}}{2}t\right),\end{array}\right.=\left\{\begin{array}{c}\frac{1}{2}\sqrt{D_1+4}-\frac{\sqrt{D_1}}{1+e^{-\sqrt{D_1}t}},\\ -1-\frac{1}{t+1},\\ \frac{1}{2}+\frac{\sqrt{3}}{2}tan\left(\frac{\sqrt{3}}{2}t\right).\end{array}\right.$$

(A26)