Qualitative Theory of Two-Dimensional Polynomial Dynamical Systems

: A qualitative theory of two-dimensional quadratic-polynomial integrable dynamical systems (DSs) is constructed on the basis of a discriminant criterion elaborated in the paper. This criterion enables one to pick up a single parameter that makes it possible to identify all feasible solution classes as well as the DS critical and singular points and solutions. The integrability of the considered DS family is established. Nine speciﬁc solution classes are identiﬁed. In each class, clear types of symmetry are determined and visualized and it is discussed how transformations between the solution classes create new types of symmetries. Visualization is performed as series of phase portraits revealing all possible catastrophic scenarios that result from the transition between the solution classes.


Introduction
DSs, in particular two-dimensional ones, are widely used to model various processes described by several dependent quantities varying in time [1,2]. Visualization and analysis of the DS phase portraits, or geometric representations of solutions [3], are in turn one of the most popular tools [4] to create a qualitative theory of DSs, reveal qualitative properties of the process, and apply the results of modeling (www.pks.mpg.de/nonlinear-dynamicsand-time-series-analysis/visualization-of-dynamical-systems/, accessed on: 20 September 2021). The qualitative analysis of DSs employs catastrophe theory [5] and diverse advanced mathematical techniques [6,7].
A mission of particular interest is to distinguish families of DSs, particularly integrable, for which one can construct a complete qualitative theory including all possible solution scenarios and end up with a complete collection (atlas) of phase portraits. This task may be solved by considering autonomous two-dimensional polynomial DSs [8,9] which still produce a variety of unsolved problems [9,10]. One can sort out among them a set of DSs integrable in elementary functions with the right-hand sides comprising quadratic trinomials. This will be a six-parameter DS family with several possibilities of polynomial factorization which generates a rich variety of solutions possessing comprehensive qualitative behavior.
Investigation of this family of polynomial DSs constitutes the main objective of the present work. We prove the integrability of the addressed DS family and show that qualitatively, the solutions of the DSs under study can be divided into nine classes defined using one real parameter and that all the DS evolution scenarios are associated with the transitions between these solution classes. The introduction of this descriptive parameter forms the core of the discriminant criterion proposed in the paper. In each class defined in the paper, specific solution examples are presented and clear types of symmetry are determined and visualized. It is discussed how transformations between the solution classes create new types of symmetries.
Visualization is performed as series of phase portraits revealing all possible catastrophic scenarios that result from the transition between the solution classes.
An atlas of phase portraits is created providing visualization of the general solutions that gives a complete description of the qualitative theory of the DS under study.
The paper contains recommendations that may assist in extending the presented research analysis. In particular, the collection of the transition/transformation pictures illustrating the transition between classes and symmetry issues addressed in the paper can be used by researchers and engineers in future studies when similar DSs are applied. Additionally, in the paper, explicit formulas are presented (together with examples) that enable one to create more visualizations applicable for different parameter sets and in diverse situations and models where the trinomial DSs are applied.

Classification of Autonomous Polynomial Dynamical Systems
Consider autonomous DS (ADS) of two dimensionṡ where x=[x 1 (t), x 2 (t)], R 2 is the two-dimensional vector space, and differentiation with respect to parameter t (time) is denoted by dot.
A subclass of ADS for which vector functions F are componentwise quasipolynomials with respect to all variables, that is, finite sums of monomials including the powers of the dependent real variables, is called polynomial DS (PDS). Autonomous symmetric PDS (ASPDS) of two variables can be represented in general form [8] comprising symmetric second-order polynomials of two variables, where A ij k = A ji k are given real numbers.
Autonomous PDS of two-dimensions, usually nonsymmetric, when (each) jth component of vector function F is a polynomial of degree two of one variable x j , j = 1, 2, with real coefficients have the form [11] where a jk are real numbers. APDS (3) at a 1k = a 2k = 0 are ASPDS containing only quadratic terms with nonzero A kk k = a 0k (k = 1, 2) and A ij k = 0; such ASPDS are integrable in elementary functions.

On Integrability of Polynomial Dynamical Systems
Generally speaking, if equations in (1) are functions of two variables, in particular, second-order polynomials f (x) = Ax 2 (2), a DS is nonintegrable. However (3) contains differential equations with separable variables and, thus, each equation in (3) is integrable in elementary functions [11] since its right-hand side is a quadratic trinomial. The number of possible factorizations of polynomials (four) means that system (3) will have 16 different combinations of them. Together with the presence of six parameters, this means that the picture and variety of qualitative features of the process modeled by DS (3) will be extremely far reaching. On the other hand, the creation of a qualitative theory that includes the description of all singular points and CPs of APDS (3) and the construction of the entire collection of its phase portraits is a difficult task indeed. The latter is the main goal of this work. We will obtain general solutions and solutions to Cauchy problems in all possible classes determined by various types of factorization, which, as will be shown, are determined by a combination of signs of the discriminants of polynomials, and study their most important properties. Solutions x = x(t) of two-dimensional DSs and, in particular, APDS (3) are represented by curves in R 2 parameterized by t (time). This is the most convenient and compact form of solutions that allows one to visualize and investigate their properties taking into account the specific requirements and characteristics of the model in which the APDSs are applied in [11][12][13][14][15][16][17].
Obtaining a general solution to APDS (3) presupposes [11] (i) the expansion of polynomials P i (s) on the right-hand side of DS to a product of irreducible polynomials, (ii) the subsequent representation of functions P i (s) −1 as a sum of elementary (irreducible) fractions, and (iii) computing the integral of the resulting expansions. All steps (i)-(iii) can be performed for polynomials of the second degree, which will lead to the possibility of obtaining a general solution of the corresponding three-dimensional APDS in the form of quadratures using the technique and results [11].
Conducting a qualitative analysis of two-dimensional APDS, we provide a description and visualization of various scenarios of models of specific processes. The results obtained are generalized in a compact matrix and criterial form, which makes it possible to determine various types of stable and unstable solutions and their singular and critical points (CPs) depending on the signs of the discriminants of the polynomials. The latter is of decisive importance for practical applications of the results of the qualitative theory constructed in this work.
In particular, in the regulation and monitoring model, for example, of the process of waterflooding and the development of an oil or gas field [15][16][17], system (3) describes the change in time on a given interval t ∈ (t 0 , t 1 ), t 0 ≥ 0 and at given constants values of the parameters of the problem a ij , i, j = 1, 2, of the accumulated oil, x(t), and water, y(t), production rates.
We carry out a qualitative study of the dynamics of solutions to system (3) using the mathematical apparatus [11] of polynomial quadratic DS and the obtained explicit formulas for general solutions.

Symmetries of Quadratic Polynomial DSs
According to well-accepted terminology, symmetry of a DS is a "transformation that takes solutions to solutions" [12,13]. The DS family considered in the paper (which may describe a physical or biological system) have definite symmetries caused by the modeling assumptions, and, to a certain extent, simplifying normal form transformations [13]. The mentioned transformations are realized explicitly by varying the DS characteristic parameters, namely, the trinomial discriminants. Nine specific solution classes are identified. In each class, clear types of the symmetry are determined and visualized and it is discussed how a transformation from class to class (between solution classes) creates new types of symmetries.
We would like to note that a deep investigation of equivariant quadratic (trinomial) DSs in relation to the discriminant criterion established by the authors goes far beyond the scope of this work and will be a subject of a separate study. An obvious aim of this particular paper, which naturally forms a base for this future study, is to obtain explicitly all general solutions in all nine classes described according to the discriminant criterion, specify examples for each class, and visualize the typical solution behavior, including the cases reflecting main transformations between classes.

Extension to Higher-Degree Polynomials
The methods and analysis of the paper can be extended to DSs of the type (1) in R n , n > 2, involving higher-degree polynomials. In fact, when a polynomial P n of degree n > 2 is factorized and its inverse P −1 n can be expanded in partial (elementary) fractions, then, under certain restrictions, the corresponding polynomial DS can be integrated (in elementary functions) and an approach can be developed similar to that described in the study. The results open a broad field for visualizations and applications. In particular, for even-or odd-degree polynomials (like e.g. ax + bx 3 + cx 5 + . . . or dx 2 + gx 4 + . . . ) one can determine symmetries and associate with them equivariant DSs [13]. However, this analysis is tedious, goes beyond the borders of this paper, and will be a subject of a separate publication.

Classes D of Solutions to Quadratic Polynomial DSs. Discriminant Criterion. Linearization
Consider system (3) and the equation obtained from it, where x = x(t) and y = y(t) are two related time-dependent quantities.
In particular, in a mathematical model of, e.g., waterflooding and development of an oil field, they denote [11] the accumulated volumes of oil and water produced over the current period of time (for example, month), and Equation (4) describes dynamic changes in y depending on x. This dynamics is described as a result of solving (4) under different "initial" conditions y(x 0 ) = y 0 and the corresponding curves y = y(x) on the phase plane (x, y).

Discriminant Criterion
To carry out a complete qualitative study of solutions to an autonomous system, including analysis and visualization of all its possible critical and singular points and solutions, we will consider the behavior of solutions of system (3) and the corresponding Equation (4) on the phase plane (x, y) depending on the problem parameters.
With the aim of constructing a qualitative theory of polynomial quadratic DSs and introducing a discriminant criterion for a qualitative theory, we concretize the definitions of the classes D ++ , D −+ , D +− and D −− as sets { f , h} of pairs of polynomials f (x) = ax 2 + bx + c and h(y) = my 2 + ny + l such that 1.
Factorization of trinomials and successive integration of Equation (4) that leads to explicit solution formulas in classes D ± are performed in [11]. Using the integration methods and formulas [11] (where one can find a comprehensive analysis of the solution properties) we obtain and summarize in Table 1 explicit formulas for the general solutions of Equation (4) (all polynomial quadratic DSs under study) with respect to the phase variables (x, y) for the classes D ± with nonzero discriminants, as well as for the classes D 0± and D ±0 when one of the discriminants of trinomials vanishes.  (4) according to the classical scheme [18] is possible only in the class D ++ when each of trinomials has two different real zeros. As a result of linearization one can describe all four CPs, their character, and location on the phase plane.

Linearization of DS (3) or the corresponding Equation
The following statement holds.

Lemma 1.
In the class D ++ there are four CPs, and in the vicinity of each of them, one can perform linearization by the classical scheme. A CP may be either a saddle (unstable) or node (stable or unstable).
To prove Lemma 1, note first that the following factorization is possible in the class D ++ (with a, m = 0) where The existence and character of the four CPs: saddles (unstable) or nodes (stable or unstable) are determined, after linearization using (6), by the value of the quantity The values of the derivatives of trinomials are actually eigenvalues λ 1 = aX = − √ D 1 and λ 2 = mY = − √ D 2 of the matrix of the linearized system (3). This enables one to introduce the quantity ρ = D 2 D 1 > 0 and consider it as a universal "control parameter", which together with its sign specifies the type (saddles or nodes, stable or unstable) of a CP of Equation (4) in D ++ .
Next, using the identities write the linearized representationF(x, y) of the right-hand side of (4) in a vicinity of a CP A kn = (x k , y n ), k, n = 1, 2, as According to the classical linearization scheme, the combination of signs of eigenvalues λ 1 and λ 2 of the matrix of the linearized system (3) determines the type of each CP. However, one can directly determine the type of each of the four CPs by solving the linearized equation using the variable substitutionỹ = y − y n ,x = x − x k . The general solutions to the linearized equations are Among the four different CPs situated at the vertices of a "critical" rectangle A 11 A 12 A 22 A 21 , two of them, A kk = (x k , y k ), k = 1, 2, are of the type node (on the main diagonal A 11 A 22 ) while two other points A kn = (x k , y n ), k, n = 1, 2, k = n, are of the type saddle (on the other diagonal A 12 A 21 ).
CPs A 22 = (x 2 , y 2 ) = (5, 11) and A 11 = (x 1 , Class D ++ plays a special role in the construction of a qualitative theory of the 6parameter family of Equation (4), since it is in this class that the solutions have the largest number of CPs (four). The general solution of (4) in D ++ has the form [11] y(x; C) = y 1 − Particularly, if D 1 = D 2 then ρ = 1. Note that in the family of functions of the general solution one can single out a special stationary particular solution y(x; 0) = y 2 , which is special in the sense that in each of its neighborhoods there are infinitely many "nonsingular" solutions, and the solutions are continuous and bounded (stable) for all x if C ≤ 0.
For different values of constant C, four different families of curves can be distinguished (analyzed in detail in [11]  Let us illustrate properties of this family considering an example with the parameters D 1 = D 2 = 1, a = m = 1, x 1 = 1, x 2 = 2, y 1 = 1, y 2 = 2. Then h(y) = my 2 + ny + l = y 2 − 3y + 2 = (y − 1)(y − 2), and the general solution of the corresponding equation If C ≤ 0, the general solution has no singularities if we consider its extension y(x; ; C) = 1 + 1 If C > 0, the general solution has movable singularities If C = 1, the general solution has the singular point x * ,2 (C) = 3 2 .

Analysis of Solutions to Cauchy Problems
Solutions to the Cauchy problems y = F(x, y), y(x 0 ) = y 0 are particular solutions to Equation (4); that is, solutions obtained from the general solution family for a specific value of an arbitrary constant that corresponds to the set of initial conditions. Their solution curves pass through a given point (x 0 , y 0 ) of the phase plane (x, y), and for them the local uniqueness holds.
Below we present typical examples of solutions to the Cauchy problems in classes D ++ , D −+ , D +− , and D −− . Their visualization is performed in Section 4.
The graphs in Figures 3 and 4 plot the solutions to the Cauchy problems in class D ++ with a given initial point highlighting their trajectories on the solution curve. For class D +− at a = 1, b = 0, c = −0.25, D 1 = 1, m = 1, n = 2, l = 1.25, D 2 = −1, for every point (x 0 , y 0 ) with x 0 = ±0.5, x = ±0.5 and the function solves the Cauchy problem There are at least two possible scenarios for the behavior of solutions to the Cauchy problem: (a) for −0.5 < x 0 < 0.5 and any y 0 the solutions are monotonically increasing functions that lie in the strip between two special stationary solutions x = −0.5 and x = 0.5; (b) for x 0 > 0.5 or x 0 < −0.5 and any y 0 the solutions are monotonically decreasing unbounded functions that lie outside the strip between two stationary solutions. In the considered case, the entire phase plane is filled with the solution curves.
For class D −+ the function We see that at 0 < y 0 < 1 we have C 0 < 0 and therefore the solutions for any x 0 are bounded monotonic functions situated in a band between two stationary solutions y = 0 and y = 1.

Qualitative Theory of Polynomial Quadratic DSs. Visualization of Transitions
The main property of a DS is its stability, which is used to evaluate how the system responds to outside impacts, assuming that its current state remains unchanged. If disturbances gain momentum over time, this results in a loss of steadiness accompanied by an abrupt (catastrophic) transition to a qualitatively new state.
A systematic study of the properties of solutions to DS (3) or Equation (4) allows us to draw the following conclusions: the formation of catastrophic changes in solutions within the framework of various scenarios for the development of catastrophes occurs when the sign of one of the discriminants of trinomials changes; that is, when passing from one class D ± or D ±0 to another. The latter corresponds to the transition of the system described by the DS from one state to a qualitatively different one. Moreover, each class corresponds to a specific structure of the solution curves. A change in this structure, including a catastrophic one, accompanied by the appearance of new types of particular solutions that are absent in the "previous" class, occurs when one parameter is varied, namely, discriminant D 1 or D 2 of one of the trinomials. This is the defining property of quadratic polynomial DSs (3) which distinguishes them from other types of polynomial systems and allows one to construct a complete qualitative theory for them and conduct an exhaustive visualization of all possible qualitative transitions and catastrophe scenarios.
Class D ++ plays a special role in these evolutions as a "source" of such changes, since, as has been shown, it is in this class that solutions have the largest number of CPs (four).

Analysis of Solutions and Phase Portraits
A series of phase portraits is plotted below using explicit formulas of the general solution to (4) with different constants C visualize a transition D ++ → D 0+ and other characteristic transitions.
The source transitions from D ++ to D 0+ or from D ++ to D +0 take place when the discriminants . The latter means that that the distance between the xor y-coordinates of CPs tends to 0, they approach each other pairwise along the sides of the "critical" rectangle in Figure 1 and merge.
By analyzing the behavior of the solution curves visualized in the section, among all, the following conclusion can be drawn: the order of change in the CP character when they merge (which occurs when passing from class D ++ to another class), depending on the degree of preservation of their properties and eigenvalues of the matrix of the linearized system, is as follows: Stable node → Unstable node → Saddle where λ 1 and λ 2 are eigenvalues of the matrix of linearized system (3). After merging: from a stable node and an unstable node, a stable node is formed; from an unstable node and a saddle, an unstable node is formed; from a stable node and a saddle, a stable node is formed.
Ultimately, for any set of trinomial coefficients in class D ++ , after the successive vanishing of discriminants and the resulting transition to class D 00 , there will always be one singular point, a stable node. The analysis is summarized in Table 2. Table 2. CPs after merging, transition from D ++ to D 00 .

Stable Node
Unstable Node Saddle

Atlas of Phase Portraits
The atlas of phase portraits is a collection of plots of general solutions of (4) in all nine solution classes illustrating all possible transitions between classes. These transitions reveal in turn all possible catastrophic scenarios of the process modeled by DS (3). Below, we analyze the visualized data presented as series of phase portraits.
Approaching the transition from the D ++ class to the D −+ class through the D 0+ class, that is, the convergence and merging of nodal and saddle CPs in Figure 2 which are located on the upper and lower sides of the rectangle (top left) occurs with decreasing and tending to 0 of the discriminant D 1 (D 1 = 1, 0.64, 0.36, 0.16, 0.04, first and second rows left to right in Figure 2, respectively), so that the distance between the CP x-coordinates tends to 0. For D 1 = 0, in the class D 0+ , there is a merging (a) of two nodal and saddle singular points 1 and 2 as a result of which one point 1 is formed, an unstable "one-sided" node, and (b) two nodal and saddle singular points 3 and 4 as a result of which one point 2 is formed, a stable "one-sided" node (second row, extreme right in Figure 2). A rupture is also formed at x = 0 and a discontinuous "fold" (second row, extreme right in Figure 2).
The transition from D 0+ to D −+ in Figure Figure 2, respectively). For D 1 < 0, in class D −+ , the solutions have no singularities, so that the discontinuity, characteristic to class D 0+ disappears and turns into a "fold" which smoothes out as |D 1 | increases. Additionally, it can be noted that the graph was inverted, the red (with positive C values) and blue (with negative C values) lines swapped relative to the asymptotes formed by the stationary solutions y = y 1 and y = y 2 .
Approaching the transition from D −+ to D −− bypassing class D −0 , in Figure 5 occurs with a decrease, vanishing, and then a transition to the region of negative values of the discriminant D 2 (D 2 = 0.01, 0.0064, 0.0036, 0.0016, 0.0004, first and second rows left to right in Figure 5, respectively, and D 2 = 0 highlighted with blue in the second row). When D 2 → 0, the distance in y-coordinate between CPs decreases, the points begin to influence each other more and more and merge into one CP (second row in Figure 5). Having passed the zero value, for D 2 < 0 in D −− close to the formed single CP (−1, −1) (second row, extreme right in Figure 5), all solutions acquire a countable set of movable singular points.
Approaching the transition from D ++ to D 00 through D +0 , that is, the convergence and merging of the nodal and saddle singular points that are located on the left and right sides of the rectangle in Figure 6 occurs with decreasing and tending to 0 of discriminant D 2 (D 2 = 1, 0.64, 0.36, 0.16, 0.04 , first and second rows left to right in Figure 6, respectively) when the distance between the y-coordinates of the singular points tends to 0. At D 2 = 0, in D +0 , there is a merging (a) of two nodal and saddle singular points 1 and 3 as a result of which one point 1 is formed, a stable "one-sided" node, and (b) two nodal and saddle singular points 2 and 4 as a result of which one point 2 is formed, an unstable "one-sided" node (second row, extreme right in Figure 6). Additionally, a gap is formed at y = 0 and a discontinuous "fold" characteristic of D +0 .
The transition from D 0+ to D 00 , that is, the convergence and merging of the nodal singular points that are located on the same vertical line in Figure 5 occurs with decreasing and tending to 0 of discriminant D 2 (D 2 = 1, 0.64, 0.36, 0.16, 0.04, third and fourth rows left to right in Figure 5, respectively), so that the distance between the y-coordinates of the singular points tends to 0. At D 2 = 0, in D 00 , two nodal "one-sided" singular points 1 and 2 merge, resulting in one point 3, a stable "one-sided" node. Additionally, the already existing discontinuity at y = 0 is supplemented by a discontinuity at x = 0, and a "cross" is formed from the discontinuity lines. The transition from D +0 to D 00 , that is, the convergence and merging of nodal singular points that are located on the same horizontal line in Figure 6 occurs with decreasing and tending to 0 of the discriminant D 1 (D 1 = 1, 0.64, 0.36, 0.16, 0.04, third and fourth rows left to right in Figure 6, respectively), so that the distance between the x-coordinates of the singular points tends to 0. At D 1 = 0, in D 00 , two nodal "one-sided" singular points 1 and 2 merge, resulting in one point 3, a stable "one-sided" node. Additionally, the already existing discontinuity at y = 0 in D +0 is supplemented by a discontinuity at x = 0 and a "cross" is formed from the discontinuity lines.
The transition from D 00 to D 0− shown in Figure 7 occurs with an increase in modulus of negative values of discriminant D 2 (D 2 = −0.16, −1, −2.56, −5.76, , first and second rows left to right in Figure 7, respectively), so that as |D 2 | increases, the tangensoids accumulate near the singular point (second row extreme right in Figure 7). Additionally, with an increase in the number (condensation) of tangensoids, one can observe the appearance of two special lines, asymptotes (y 1 = −0.5, y 2 = 0.5, which are the real parts of the zeros of the trinomial), between which the tangensoids accumulate.
The transition from D 0− to D +− shown in Figure 7 occurs when the positive values of discriminant D 1 increase (D 1 = 0.09, 0.25, 0.49, 0.81 , third and fourth rows left to right in Figure 7, respectively), so as the value of D 1 increases, one can observe how two singular points 2 and 3 are formed from one point 1 and begin to move away from each other, while the tangensoids, in turn, accumulate in their vicinity.

Conclusions
A qualitative theory has been constructed of two-dimensional quadratic-polynomial integrable DSs. It has been shown that basic qualitative features of solutions are governed by the signs of discriminants of the trinomials entering the DS right-hand side. The latter has made it possible to formulate the discriminant criterion as a main axis of the qualitative theory.
Nine (all possible) DS solution classes have been identified and described possessing essentially different qualitative properties, as well as the DS critical and singular points and solutions.
Visualization has been performed, resulting in an atlas of phase portraits. The detailed analysis of the qualitative behavior has revealed all possible catastrophic scenarios that result from the transition between the solution classes.
Author Contributions: Conceptualization, methodology, resources, and data curation, A.S.; validation, formal analysis, investigation, writing-original draft preparation, and writing-review and editing, Y.S. All authors have read and agreed to the published version of the manuscript.