Exactly Solvable Quadratic Differential Equation Systems Through Generalized Inversion

We study the autonomous systems of quadratic differential equations of the form x˙i(t)=x(t)TAix(t)+viTx(t)\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${\dot{x}}_i(t)={\textbf{x}}(t)^T {\textbf{A}}_i {\textbf{x}}(t) + {\textbf{v}}_i^T {\textbf{x}}(t)$$\end{document} with x(t)=(x1(t),\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${\textbf{x}}(t) = (x_1(t),$$\end{document}x2(t),…,xi(t),⋯)\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$x_2(t), \ldots ,x_i(t),\dots )$$\end{document} which, in general, cannot be solved exactly. In the present paper, we introduce a subclass of analytically solvable quadratic systems, whose solution is realized through a multi-dimensional generalization of the inversion which transforms a quadratic system into a linear one. We provide a constructive algorithm which, on one hand, decides whether the system of differential equations is analytically solvable with the inversion transformation and, on the other hand, provides the solution. The presented results apply for arbitrary, finite number of variables.


Introduction
Linear, ordinary, autonomous systems of differential equations are celebrated for their analytical integrability which makes them excellent tools in almost all fields of science.This solvability breaks down when entering the realm of non-linear ordinary differential equation systems where no general solution method exists.In most cases, non-linear systems are handled with sophisticated numerical methods such as Runge-Kutta schemes [1].
Despite of the great accuracy of numerical methods, it is also desirable to find analytical techniques to solve non-linear systems.Beside the inherent beauty of exact solutions, they can also play the role of benchmark for numerical methods.In this work, we focus on firstorder, quadratic, autonomous differential equation systems (QDEs) which emerge in various research fields including models of population dynamics [2], fluid dynamics [3], control systems [4] and even quantum dynamics [5].Classification of two-variable quadratic systems have been studied extensively [6,7,8,9].
Email address: bacsi.adam@sze.hu( Ádám Bácsi) There exists no general method to find analytical solution of QDEs.Except for the special case, when the number of variables is just one.By denoting this variable with x(t), the quadratic differential equation is written as with a = 0 and v being constants.An additional constant term on the right-handside could be eliminated by shifting x(t).By introducing y(t) = x(t) −1 , the equation is transformed to the linear equation of and, hence, the quadratic differential equation ( 1) becomes exactly solvable.
In this paper, we investigate a class of QDEs with more than one variable which can be solved exactly by transforming the QDE to a linear differential equation system (LDE).The transformation is realized through a multi-variable extension of y(t) = x(t) −1 , similar to the generalized spherical inversions presented in Refs.[10,11].

Generalized inversion transformation
Let us consider some differentiable functions x 1 (t), x 2 (t), . . ., x n (t) whose dynamics is governed by the QDE ẋi where A i are n × n matrices which are chosen to be symmetric and v T i are row vectors of n components.The entries of A i and v i are independent from t and are assumed to be real without loss of generality 1 .Our goal is to determine if the variables x i (t) can be transformed with a multivariable variant of the inversion in such a way that the QDE (3) is transformed to an LDE.The multivariable inversion is defined as with some symmetric matrix B which does not depend on time.Note that we use the same B matrix for all i.This ensures that the transformation is easily inverted by The transformations (4) and (5) will henceforth be referred to as B-transformation and inverse B-transformation, respectively.The variables y i (t) are expected to obey the linear differential equations ẏi = m T i y + w i .
with time-independent row vectors m T i and the scalars w i .Since Eq. ( 6) is exactly solvable for arbitrary initial conditions, the QDE can also be solved by applying the inverse Btransformation (5).The goal of the present paper is to provide a procedure based on which one can decide whether a given QDE, defined through A i and v i , can be B-transformed to an LDE.
To start our study, let us calculate the time derivative of Eq. ( 5) ẋi = ẏi y T By − y i (y T By) ẏT By + y T B ẏ .
Substituting the LDE differential equations ( 6) and then transforming all y i variables back to x i by means of Eq. ( 4), we obtain with the matrix M built from the row vectors m T i and the vector w consisting of the constants w i .
By comparing Eq. ( 8) with the original QDE (3), one can note that Eq. ( 8) is not necessarily a quadratic differential equation.The last term has an x-dependent denominator and, hence, is neither quadratic nor linear in general.This term scales with x as a linear term but becomes an actual linear term only if holds with some scalar −λ.The minus sign is chosen for later convenience.The condition (9) means that B must be the eigenmatrix [12] of M with the eigenvalue of −λ.In this case, the differential equations ( 8) become where • denotes dyadic product and e T i is the unit row vector with 1 in the ith entry and zero otherwise.The equation is already a quadratic differential equation containing quadratic and linear terms in x.
Comparing the linear terms of Eqs. ( 10) and (3), we find where V is constructed from the row vectors v T i and I is the n×n identity matrix.Substituting (11) into the condition (9), we obtain indicating that the matrix B must be a symmetric eigenmatrix of V as well.
Beside the linear terms, the quadratic terms of Eq. ( 10) and (3) must also equal.Since both A i and the kernel in the quadratic terms in Eq. ( 10) are symmetric, they must equal for all i.Note that the second and third terms are matrices full of zeros except for the ith row and ith column, respectively.Therefore, omitting the ith row and column in the equation leads to the proportionality condition where Ãi i ( Bi ) is the (n − 1) × (n − 1) matrix which is obtained from A i (B) by skipping the ith row and column.For the ith column of Eq. (13), we have To summarize, the QDE as given in Eq. ( 3) can be transformed with a B-transformation to an LDE if we find a symmetric matrix B fulfilling Eq. ( 12) with some constant λ and, at the same time, obeying Eqs. ( 14) and (15) with some constants w i .The linear terms of the LDE can be calculated by expressing V from Eq. ( 11).

Algorithm
Let us provide a detailed description of the procedure based on the findings of the previous section.The starting point of the algorithm is Eq. ( 3) determined by the matrices A i and vectors v i .
1.As a first step, one has to build up the matrix V from the row vectors v T i and find all symmetric eigenmatrices based on V T B + BV = λB.Note that the eigenmatrices are definite only up to an overall multiplication factor similarly to the case of usual eigenvectors.
To obtain the eigenmatrices, one may calculate the vector-eigenvalues s j of V T and the corresponding right-handside eigenvectors r j , i.e., V T r j = s j r j .The sum of the geometric multiplicity of all eigenvalues is denoted by N .The matrix-eigenvalues of V are given by s 1 + s 1 , s 1 + s 2 , . . ., s 1 + s n , s 2 + s 2 , s 2 + s 3 , . . .,s N + s N and the eigenmatrix corresponding to s j + s m is P jm = r j • r T m + r m • r T j /2.Note that the number of symmetric eigenmatrices is N (N + 1)/2 which is utmost n(n + 1)/2 when n eigenvectors are found.
In case of degenerate matrix eigenvalues λ, all linear combinations within the subspace are potential B matrices.For an example, see Sec. 5.

2.
For each potential eigenmatrix B, the proportionality condition (14) must be checked for all i with some constants w i , which may also be zero.
If Eq. ( 14) cannot be fulfilled with any potential B matrix, then the QDE is not solvable by B-transforming to an LDE.If Eq. ( 14) is satisfied with a B matrix and some w i constants, one may proceed to Step 3.

3.
For the potential pairs of B and w surviving Step 2, one has to check Eq. ( 15) for all i.
If it holds true for all i, then the QDE (3) can be transformed to an LDE of the form of Eq. ( 6) where the constants w i are the coefficients of proportionality found in Eq. ( 14) in Step 2 and the vectors m T i are the row vectors of M = V − λI with λ the eigenvalue corresponding to the eigenmatrix B.
4. The solution of the QDE can be obtained by solving first the linear equation (6) as where y 0 = x 0 /(x T 0 Bx 0 ) is the initial condition with x 0 = x(t = 0) and y p (t) = e Mt − I M −1 w if M is invertible.If M is not invertible, the function y p (t) has linear time-dependence in the nullspace of M. Using the inverse B transformation, the solution of the QDE is obtained as Before presenting examples, let us make a remark on the region where the B-transformation is undefined, i.e., where x T Bx = 0.This region can be any quadratic hypersurface depending on the specific structure of the B matrix.The region also involves the nullspace of B, i.e., the points where Bx = 0.In Sec. 4, the hypersurface contains two straight lines while in Sec. 5, it consists of two cone-shaped surfaces.This raises the question how the system behaves if the initial condition x 0 is an element of the region, i.e. x T 0 Bx 0 = 0. First, one can prove that if the initial condition of the QDE is on the hypersurface, the dynamics is constrained to the hypersurface for the whole time-evolution.To justify this, we define b(t) = x(t) T Bx(t) whose dynamics is derived from the quadratic (3) as ḃ = b λ − 2w T Bx .If the initial condition is on the hypersurface, i.e., b(0) = 0, then the differential equation solves to b(t) = 0.
Second, it can also be shown that inside the hypersurface, the quadratic differential equation is analytically solvable and the solution is obtained simply by taking x T 0 Bx 0 = 0 in Eq. ( 17) which reads as Note that (18) does not obviously solve the quadratic differential equation because it was derived from Eq. ( 17) which was computed through the B transformation.However, by substituting (18) directly into (3) and taking advantage of x T 0 Bx 0 = 0 and the properties of B, one can prove that the QDE is solved indeed by (18).Hence, if a QDE is exactly solvable by a B-transformation, the integrability is preserved also in the region where the generalized inversion is undefined.

Example, two-dimensional quadratic system
In this section, we consider the quadratic equations of from which we read out The first step is to obtain the eigenmatrices of V.The right-handside eigenvectors of V T are given by with the eigenvalues of s 1 = −2 and s 2 = 1.The potential B matrices are given by which are non-degenerate eigenmatrices.The second step is to check the proportionality condition of Eq. ( 14) and determine the coefficients w i .Proportionality check, Eq. ( 14) The table shows that all eigenmatrices obey the proportionality condition (14).Note that for a two-dimensional QDE, the proportionality check always simplifies to comparison of 1 × 1 matrices, which is in most cases trivially fulfilled.For higher dimensions, however, (14) might mean a much stricter condition.For an example, see Sec. 5.
In the example, we continue with the third step by checking Eq. (15) for each i.In the table below, rhs stands for the right-handside of Eq. (15) evaluated with B and w.
The investigation shows that a 1 1 and a 2 2 are reproduced only by P 12 .The result of the algorithm is that the QDE can be transformed to an LDE by B = P 12 .The linear system is obtained as which can be solved analytically for arbitrary initial conditions.The analytical solution of the quadratic system can be given based on Eq. ( 17).The phase portrait is shown in Fig. 1.The two dashed lines comprise the region where the B transformation is undefined, i.e., where x T Bx = (x T r 1 ) • (x T r 2 ) = 0.  19).The green dots indicate the fix points of the system.The dashed lines comprise the region where the B transformation is undefined.

Three-dimensional example
In this section, the algorithm is demonstrated in an example with three variables x 1 (t), x 2 (t) and x 3 (t).The quadratic differential equations are given by ẋ1 from which we can read out and The eigenvectors of V T are simply obtained as r i = e i with the eigenvalues of s 1 = 5, s 2 = 2 and s 3 = −1.Hence, the eigenmatrices of V matrices are as follows.
Note that the eigenmatrices P 22 and P 13 are degenerate.Therefore, their linear combination is also a potential B matrix with an arbitrary value of b.
We continue the algorithm with Step 2, the proportionality check.which are exactly solvable for arbitrary initial conditions.The analytical solution can be given based on Eq. ( 17).Note that the B transformation is undefined in the region where x T Bx = x 1 x 3 + x 2 2 = 0.This quadratic equation determines two cone-shaped surfaces touching each other at the origin as shown in Fig. 2.

Conclusion
We studied the system of quadratic differential equations of the form of Eq. (3).Although these differential equation systems cannot be solved in general, we have presented a

− 1 . 2 Figure 1 :
Figure1: Phase portrait of the quadratic system, Eq. (19).The green dots indicate the fix points of the system.The dashed lines comprise the region where the B transformation is undefined.

if b = 1
It has been found that the quadratic system transformed to a linear system with the matrix and the resulting linear differential equations are given as ẏ1 = y 1 + 7 ẏ2 = −2y 2 + 2 ẏ3 = −5y 3 − 1 (28)

Figure 2 :
Figure 2: The hypersurface determined by x T Bx = 0 consists of two cone-shaped surfaces.The hypersurface contains the two special directions determined by ℓ1 (blue) and ℓ3 (red).