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

Abstract

The behavior and bifurcations of solutions to three-dimensional (three-phase) quadratic polynomial dynamical systems (DSs) are considered. The integrability in elementary functions is proved for a class of autonomous polynomial DSs. The occurrence of bifurcations of the type-twisted fold is discovered on the basis and within the frames of the elements of the developed DS qualitative theory. The discriminant criterion applied originally to two-phase quadratic polynomial DSs is extended to three-phase DSs investigated in terms of their coefficient matrices. Specific classes of D- and S-vectors are introduced and a complete description of the symmetry relations inherent to the DS coefficient matrices is performed using the discriminant criterion.

1. Introduction

Analysis of bifurcations in polynomial dynamical systems (DSs) constitutes 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 non-autonomous 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.

We identify and consider an integrable family of three-dimensional (three-phase) polynomial quadratic DSs. We show, as an extension of a similar study performed in [6] for two-dimensional polynomial quadratic DSs, that one can carry out a 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 characteristic 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 criterion. 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 non-autonomous 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, an unstable front area 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 in 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 oil producing wells and, consequently, to the reduction in oil recovery and oil production stimulation.

Creating a rigorous mathematical model of this process remains an urgent task of oil science and beyond.

In the proposed setting (a version of the growth model) simulating the development of an oil or gas field, the considered three-dimensional autonomous polynomial DSs (system (42) in the text below) describe [8,10,13] temporal flow variations corresponding to the three phases: oil ($D_o$), gas ($D_g$), and water ($D_w$). The discriminants of the growth model for each phase 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 huge amounts of the registered, experimental, and measurement data [7,11,14].

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 the catastrophe theory [14,15,16,17]. 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 for predicting the breakthrough of water and gas into oil producing wells.

When performing qualitative analysis of the three-dimensional DSs addressed in the paper we describe specific modeling scenarios. The treatment is based on processing the data collected as matrices using a technique based upon the notions of D- and S-vectors. This approach enables one to identify all possible types of stable and unstable solutions and their singularities governed in the end by the signs of the polynomials’ discriminants.

A D-vector associated with the coefficient matrix of a three-dimensional DS is a vector-valued function with the domain formed by the matrix rows; for a quadratic polynomial DS, the D-vector components are the row discriminants (and as such, they are considered as functions of the three coefficients of the polynomial forming the right-hand side of every row equation) and the S-vector components are their signs. An S-vector may be represented as an ordered set of three numbers: $-1,0,1$. Thus, S-vectors provide a kind of informative description of a set of D-vectors and may be used to define specific equivalence classes of these vectors.

A more detailed analysis and necessary definitions of D- and S-vectors are given in the subsequent sections of the paper.

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, identify the place of the approach developed in this study, we consider the second- and third-order (two- and three-dimensional) polynomial DSs which 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_{8}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 a dot denotes differentiation with respect to time.

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

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

(3)

Different sets of real coefficients $a,b$ or $a,b,c$ that enter matrices (3) govern the qualitative behavior 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.\mathrm{or}\dot{\mathbf{x}}=\mathrm{A}\mathbf{x},\mathrm{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 corresponding to the following settings: ${\lambda }_1>0$ and ${\lambda }_2>0$, ${\lambda }_1<0$ and ${\lambda }_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 non-zero 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, and saddle).

One can create a 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 the 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, examples of their so-called canonical forms [5] with ’mixed’ variables accepted are illustrated 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 {\mathrm{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 {\mathrm{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 {\mathrm{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],{\mathrm{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 entering the systems, denoted by letters $a,b,c$ with different lower indices characterized each by specific properties of bufircations.

The volume of research and quantity of publications dealing with qualitative analysis of quadratic two-dimensional polynomial DSs may be characterized as enormous; however, many problems remain unsolved. Exemplify the classification of 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. The references to the definitions of antisaddle can be found in [2].

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 so-called T-systems having 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 {\mathrm{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 [18] where the notion of the T-system is properly defined (we also refer to the references therein), where, among all, the pitchfork and Hopf bifurcations occurring in the T-system are reported. Seemingly simple, these systems require the elaboration of advanced mathematical tools within the frames of the singular perturbation and geometric singular perturbation theory [19].

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 {\mathrm{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 performed 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 {\mathrm{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 use a similar matrix notation 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 the qualitative behavior of solutions; however, as we show, many important properties are actually governed by 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, which are the discriminants of the polynomial-entering DSs using the discriminant criterion [6,7].

The proposed matrix representation of polynomial DSs is, among all, a definitely 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 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 the 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 row vectors

$$\begin{array}{c}\hfill {\mathrm{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 the 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({\mathrm{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({\mathrm{A}}^{\left(3\right)}\right)=\left[\mathrm{sing}{D}_1,\mathrm{sing}{D}_2,\mathrm{sing}{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$; they are each considered 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 yielding 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${}^3$ 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${}^3$, 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${}^3$. Namely, the first, second, etc., octants of R${}^3$ correspond to S-vectors $\left[+++\right]$, $\left[-++\right]$, etc.; the coordinate planes in R${}^3$ 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${}^3$ correspond to S-vectors $\left[\pm 00\right]$, $\left[0\pm 0\right]$, $\left[00\pm \right]$, etc.

We can now define 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) \mathrm{is} \mathrm{the} \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-vectors} \left(20\right)<{\mathrm{R}}^3\to {\mathrm{S}}^3>\mathrm{Set} \mathrm{of} \mathrm{S-vectors} \left(20\right).\end{array}$$

(22)

The 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 {\mathrm{A}}^{\left(3\right)}={\mathrm{I}}^{\left(3\right)}=\left[\begin{array}{ccc}1& 0& 0\\ 0& 1& 0\\ 0& 0& 1\end{array}\right],{{\mathrm{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({\mathrm{I}}^{\left(3\right)}\right)={\mathbf{D}}^{\left(3\right)}\left({{\mathrm{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 {\mathrm{A}}^{\left(3\right)}=\left[\begin{array}{ccc}a& 0& 0\\ 0& b& 0\\ 0& 0& c\end{array}\right],{{\mathrm{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({\mathrm{A}}^{\left(3\right)}\right)={\mathbf{D}}^{\left(3\right)}\left({{\mathrm{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 and lower triangular matrix:

$$\begin{array}{c}\hfill {\mathrm{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 {\mathrm{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 {\mathrm{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) and (24), the D- and S-vectors are invariant w.r.t. the interchange (permutation) of the first and third rows, 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${}^3$ 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 the 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 the coefficient matrix (17) as a triple (row) of column vectors,

$$\begin{array}{c}\hfill {\mathrm{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

$$\begin{array}{ccc}\hfill \mathrm{sign}{\mathbf{a}}_1{\mathbf{a}}_3<0\mathrm{or}\mathrm{sign}{\mathbf{a}}_1{\mathbf{a}}_3>0\mathrm{and}{\mathbf{a}}_2^2& >& 4{\mathbf{a}}_1{\mathbf{a}}_3\to {\mathbf{S}}^{\left(3\right)}=\left[+++\right],\hfill \end{array}$$

(30)

$$\begin{array}{ccc}\hfill \mathrm{sign}{\mathbf{a}}_1{\mathbf{a}}_3>0\mathrm{and}{\mathbf{a}}_2^2& <& 4{\mathbf{a}}_1{\mathbf{a}}_3\to {\mathbf{S}}^{\left(3\right)}=\left[---\right],\hfill \end{array}$$

(31)

$$\begin{array}{ccc}\hfill {\mathbf{a}}_2^2& =& 4{\mathbf{a}}_1{\mathbf{a}}_3\to {\mathbf{S}}^{\left(3\right)}=\left[000\right],\hfill \end{array}$$

(32)

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 in 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 rows, 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 {\mathrm{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 non-zero 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 specify the subsets of the equivalence classes of matrices (34):

$$\begin{array}{ccc}\hfill \mathrm{sign}{\mathbf{a}}_1<0\mathrm{or}\mathrm{sign}{\mathbf{a}}_1>0\mathrm{and}{\mathbf{a}}_2^2>{\mathbf{a}}_1& \to & {\mathbf{S}}^{\left(3\right)}=\left[+++\right],\hfill \end{array}$$

(36)

$$\begin{array}{ccc}\hfill \mathrm{sign}{\mathbf{a}}_1>0\mathrm{and}{\mathbf{a}}_2^2<{\mathbf{a}}_1& \to & {\mathbf{S}}^{\left(3\right)}=\left[---\right],\hfill \end{array}$$

(37)

$$\begin{array}{ccc}\hfill {\mathbf{a}}_2^2={\mathbf{a}}_1& \to & {\mathbf{S}}^{\left(3\right)}=\left[000\right].\hfill \end{array}$$

(38)

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]\mathrm{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 the integrability of (15) for the arbitrary set of its nine parameters (coefficients), denoted by letters $a,b,c$ with lower indices, 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], we 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 with specific vector functions $\mathbf{F}$ 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{\mathrm{A}}_k^{ij}{x}_i\left(t\right)x_j\left(t\right),k=1,2,\dots ,n,$$

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

There are several families of algebraically integrable autonomous polynomial DSs mentioned in the Introduction section. In general, however, autonomous polynomial DSs are not integrable. Note that such DSs containing solely equal-index quadratic terms with the non-zero $A_k^{kk}=a_{1k}$ ($k=1,2,\dots ,n)$ and the rest of $A_k^{ij}=0$ constitute 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 non-symmetric and where (each) ith component of $\mathbf{F}$ is an $N_i$-degree polynomial in $x_i,i=1,2,\dots ,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,\dots ,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$. Autonomous n-of dimensional polynomial DSs (41) are algebraically integrable in elementary functions; the corresponding statements are proven in Appendix A.

DSs (41) of the dimension $n=3$ with $\overline{N}=\left[2,2,2\right]$ and $N=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$), which are generally non-symmetric.

As follows from Theorem A1 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.

On the next steps, 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 the quadratic three-dimensional DSs: solutions $\mathbf{x}=\mathbf{x}\left(t\right)$ to these DSs, and particularly to (42) are the curves in ${\mathrm{R}}^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

List all necessary solutions $x\left(t\right)$ to the single autonomous polynomial equation $\dot{x}=f\left(x\right),f\left(x\right)=a{x}^2+bx+c$ which will be used in the subsequent sections. By D, we will denote discriminants of the involved quadratic polynomials.

A detailed description of the properties of the addressed solutions can be found in [6].

• (a) $a=0$; $x\left(t\right)={\frac{1}{b}(-c+Ce}^t),b\ne 0.$
  In what follows, we will consider the equations with $a\ne 0.$
• (b) D > 0; there are three solution families

  $${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,$$

  where the − and + signs correspond, respectively, to $x_1$ and $x_2$:
  Family U, with $C>0$: solutions $x\left(t;C\right)=x_{\mathrm{U}}\left(t;C\right)$ are not stable, since there is a “movable” singular point with $t=t^{*}\left(C\right)=\frac{lnC}{\sqrt{D}}$; next, they are (i) monotonically increasing because $\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; (ii) 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$.
  Family S, with $C>0$: solutions $x\left(t;C\right)=x_{\mathrm{S}}\left(t;C\right)$ are stable and $x_1<x\left(t;C\right)<x_2$.
  Family T, : $x=x_1$ and $x=x_2$ are time-independent solutions such that $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}}}$. These stationary solutions $x=x_1$ and $x=x_2$ are ’nonisolated’: in every neighborhood, there is an infinite number of ‘regular’ solutions $x_{\mathrm{U}}\left(t;C\right)$ or $x_{\mathrm{S}}\left(t;C\right)$.
• (c), $D=0$: all the corresponding solutions

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

  are not stable.
• (d), $D<0$: the corresponding solutions

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

  are not stable.

Note again all properties of the general solutions to autonomous second-order polynomial equation are fully described in [6].

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 non-zero discriminant signs. We identify as well 27 solutions to DS (42) of dimension three as sets D${}^{pmz}$, including all the above-defined eight solution sets, 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 non-zero 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 component-wise 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

$$\begin{array}{l}{\mathbf{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})],\\ \dots \dots ,\\ {\mathbf{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})],\end{array}$$

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)

$$\dots \dots$$

$${\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$.

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 eight possible solution sets of system (42) corresponding to all possible combinations of the solution components having non-zero determinants and characterized by their S-vectors. Empty entries may be replaced by zeros.

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 non-zero and zero determinants and characterized by their S-vectors. Empty entries in (48) may be replaced by zeros.

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,$ where polynomial $P_1\left(x\right)$ is given by (42), are

$$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,\dots ,$$

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.$

4. Analysis of Bifurcations

Families of vector functions ${\mathbf{x}}^{...}(t;\mathbf{C};\mathbf{D})$ exemplified in (44) with $\mathbf{C}=\left[C_1,C_2,C_3\right]\in \mathrm{C}\subset$R${}^3$ and $\mathbf{D}=\left[D_1,D_2,D_3\right]\in \mathrm{D}\subset$R${}^3$, where C, D, and ${}^{...}$ denote, respectively, certain 3D sets (or domains) and any combination of three signs +, − and 0, each specify 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\}$; an example of such surface is presented in Figure 1.

We particularly investigate 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.

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)

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

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 the variation in 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 vectors $\mathbf{C}=\left[C_1,C_2,C_3\right]$ and $\mathbf{D}=\left[D_1,D_2,D_3\right]$ and 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]$.

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 in Figure 3, Figure 5, Figure 9, and Figure 10.

$$\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 11 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 signs as $D_1<0$ increases, which causes a 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 signs and the solution curves in both sets D${}^{-0-}$ and D${}^{+0-}$ preserve the same shape.

5. Conclusions

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 on 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 behavior 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.

Based on the results of this research and a big amount of computations, we can formulate a conjecture (to be proved) that the only type of bifurcations that may occur in 3D quadratic polynomial DSs investigated in this paper is the discovered ‘twisted fold’. This is a remarkable feature of this class of DSs. However, a complete proof of this statement goes far beyond the scope of this paper and may be technically rather difficult (as suggested by similar proofs and findings concerning the classification of bifurcations and critical points published by many authors, revised briefly in Section 2.1). In addition to this, we cannot provide a substantiated report concerning the type of singular points and bifurcations that are observed for each of the identified equivalence classes. We will continue our studies in this direction using the accumulated technical data and will report the results in the following publications.

Indicate in this respect main directions of the planned future research:

-

To develop the method of S- and D-vectors and the discriminant criterion to the polynomial DSs of higher dimensions and the order of the involved polynomials.

-

To clarify that the only type of bifurcations that may occur in quadratic polynomial DSs investigated in this paper is the discovered ‘twisted fold’.

-

To investigate the relations between the described symmetries of the D- and S-vectors and the possible symmetries of solutions to the polynomial DSs.

-

To find the symmetry-breaking bifurcations characteristic to the polynomial DSs under study.

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 unstable behavior of the displacement front.

Using the results of this work, engineers and practitioners will be able to predict the blunt variations in water saturation and well exploitation data by applying the developed growth model that employs three-phase quadratic polynomial DSs.

The results and findings reported in this study make it possible to carry out systematic improvement in handling real-time waterflooding and to control a decrease in the water volume that is not properly put into action and pumped out.