# Symbolic solution of systems of polynomial differential equations via the
Cauchy–Riemann equation. Applications to kinetic differential equations

Kelvin Kiprono  Department of Analysis and Operations Research, Institute of
Mathematics, Budapest University of Technology and Economics, Műegyetem rkp.
3., H-1111 Budapest, Hungary  János Tóth  Department of Analysis and
Operations Research, Institute of Mathematics, Budapest University of
Technology and Economics, Műegyetem rkp. 3., H-1111 Budapest, Hungary
Chemical Kinetics Laboratory, Institute of Chemistry, ELTE Eötvös Loránd
University, Pázmány P. sétány 1/A, H-1117 Budapest

###### Abstract

The differential equations of chemical kinetics are systems of nonlinear
(polynomial) differential equations, therefore their solutions cannot usually
be found in symbolic form. Here we offer a method to solve classes of kinetic
differential equations based on the Cauchy–Riemann equations. It turns out
that the method can be used to symbolically solve some polynomial differential
equations that are not necessarily kinetic, as well.

MSC: 34A05, 80A30, 34A34, 92C45

###### Contents

  1. 1 Introduction
  2. 2 The method
  3. 3 Special cases
    1. 3.1 Linear systems, first-order reactions
    2. 3.2 Second-degree systems, second-order reactions
    3. 3.3 Third-degree systems, third-order reactions
    4. 3.4 Rth-degree systems
  4. 4 Application of multivariate complex functions
  5. 5 Discussion and Outlook
    1. 5.1 Further directions
      1. 5.1.1 Partial differential equations
      2. 5.1.2 Approximations
    2. 5.2 Final remarks
  6. 6 Appendix
    1. 6.1 Useful formulas
    2. 6.2 Addendum to the second-degree case
    3. 6.3 Addendum to an example of the second-degree case
    4. 6.4 Addendum to the third-degree case
    5. 6.5 An example of the four-variable case

##  1 Introduction

The solutions of the differential equations of chemical reaction kinetics
generally cannot be explicitly given symbolically, except in a few, very
special cases. Here we consider the solution of a kinetic differential
equation as explicit or symbolic, if the concentrations of the species are
given as a function of time, not necessarily via elementary functions, and
infinite series also including. We categorize solutions into this class also
if the inverse function of a concentration vs. time function is given. (See a
review on closed form solutions [1].) It is almost as useful as these if the
concentration of one species is given as the function of the concentration of
another. As is customary, we define a solution as symbolic if a finite
algorithm (typically consisting of fewer than 10 steps in practice) is
provided, allowing for the construction of the solution, as seen in the
derivation of the solution of cubic or similar equations, or simple
differential equations. However, it is important to note that this definition
may imply that solutions with symbolic parameters (or with "unlucky" numeric
values) cannot be practically calculated in a usable form.

Preliminary work in this field includes [29] presenting a large collection of
explicitly solvable cases. A major merit of that paper is that the author also
attached a chemical example with references to all the cases he was able to
solve. In some cases (as it is known from the theory of differential
equations) he was only able to provide the inverse of a concentration vs. time
function. Also, it happened sometimes that he could give the _selectivity
curves_ , i.e. the concentration of one species as the function of the
concentration of another. These are the projections of the trajectories in
mathematical terms.

Another large collection has been published in [26]. These authors have shown
how to solve the case when all the reaction steps have an order not larger
than one: all the reaction steps are of order one or zero. No wonder that one
of their tools is Laplace transformation. Although the solution of linear
differential equations that can be given is generally considered to be
explicit, it is only relatively explicit, so far as the eigenvalues and
eigenvectors are known—which is not the case if the number of variables is
larger than four.

Recently, Lente [19, 20, 21, 22] solved quite a few kinetic differential
equations: those describing reactions with two steps and not being of the
order larger than two. The explicit solutions contain special functions like
the hypergeometric function, exponential integral, etc., and also infinite
series. Again, in some cases, he could give the _selectivity curves_. The work
by Park [24] has similar, but more limited ambitions. Hernandez and Lubenia
[15] could provide a symbolic representation of the stationary states of
kinetic differential equations.

Kamke has gathered the classical collection of almost 2000 symbolically
solvable ordinary differential equations [18]. (He also published a volume on
partial differential equations: [17].) These days a minimal requirement for a
mathematical program package, like e.g. the Wolfram Language, is to be capable
of reproducing Kamke’s solutions and going beyond. A more recent and more
modern collection can be found on the Internet: [25]. Still, this is not the
end of the story because papers studying first integrals [23] or symmetries
[12, 13] also contribute to this area. Let us note that in pure mathematics,
finding explicit solutions to ordinary differential equations is a separate
area, papers belonging to this are classified under 34A05 of AMS MSC, see e.g.
the works by Calogero and his coworkers [2, 3, 4, 5]. As our work goes beyond
the boundaries of reaction kinetics, the present results are an addition to
this kind of Painlevé approach as well, as [2, p. 1] puts it.

It is quite common in the literature of reaction kinetics that a special,
individual equation (not a parameterized class) is solved, we are not
interested in these cases as these may count to millions.

In the present paper, we aim to find kinetic differential equations of the
Cauchy–Riemann–Erugin type (cf. [35, V. §2, $4^{o}$], and [35, Ex. 1091.,
1092.], furthermore [34, p. 116],[30]): two-variable equations with a right-
hand side the components of which can be considered to be the real and
imaginary part of a differentiable complex function. Then, we are going to
solve the kinetic differential equations belonging to this class. Next, we
shall give a few (second and third-order) complex chemical reactions whose
induced kinetic differential equation can be solved by the method. To find a
realization of a given kinetic differential equation is far from being a
trivial problem, see e.g. [6, 7, 9, 14, 16].

A restriction of our method is that the reaction rate coefficients should
fulfill some equalities, this is not quite uncommon either in reaction
kinetics or in the theory of differential equations. The method is extended to
four-variable equations. Although kinetic differential equations form a
special subclass of polynomial differential equations, our method will only be
easier to apply to them, the results of the method extend beyond this field.

If our equation is not of the Cauchy–Riemann–Erugin type but is close to that
then the solutions will also be close. However, the trajectories may suffer
quantitative or even qualitative distortions.

The structure of our paper is as follows. Section 2 presents the essence of
our method. Section 3 starts with a trivial example: the equation of first-
order reactions in two variables. The next part is about second and third-
order reactions, providing more useful results from the point of view of
kinetics because one has more free parameters in higher-order cases. This
section contains a nice compact form of the original real equations and the
complex equation as well (Theorem 1) that may be the result most interesting
from the mathematical point of view. Section 3 applies the multivariate
Cauchy–Riemann equation to extend the method to four (or, $2n$) variable
cases. Then, perspectives of approximations obtained by the method, partial
differential equations, and possible applications of quaternions are also
treated or at least mentioned in the Discussion and Outlook section. In all
cases, we start without the restriction to kinetic differential equations [10,
31], the results will only be specialized at the end. Most of the general and
special calculations are relegated to the Appendix.

##  2 The method

The steps of the method are as follows.

  1. 1.

Start from a two-variable differential equation with a polynomial right-hand
side

$\dot{x}=u\circ(x,y),\quad\dot{y}=v\circ(x,y),$ (1)

and apply the restrictions on the coefficients coming from the assumption that
the right-hand sides are the real and imaginary parts of a differentiable
complex function $f,$ i.e. they fulfill the Cauchy–Riemann conditions:

$\partial_{1}u=\partial_{2}v,\quad\partial_{2}u=-\partial_{1}v.$ (2)

  2. 2.

Rewrite the differential equations into a single equation for the complex
function $f.$

  3. 3.

Solve the equation for $z:=x+iy.$

  4. 4.

Determine the real and imaginary parts of the solution.

  5. 5.

Determine the general form of the induced kinetic differential equation of
two-species first-order reactions, taking care of the absence of negative
cross-effect, and specialize the above-obtained results to this special case.

  6. 6.

Finally, try to find "simple" realizations for the pair of real differential
equations: with as small many complexes and reaction steps as possible. Note
that this is not our final goal here.

Let us mention in passing that the relations (2) can be rephrased as

$f^{\prime}J=Jf^{\prime}\text{ with }J:=\begin{bmatrix}0&1\\\
-1&0\end{bmatrix}.$ (3)

For the multivariate case of two complex functions and two complex variables
such as the one below.

$\displaystyle f_{1}(z_{1},z_{2})$
$\displaystyle=u_{1}(x_{1},x_{2},y_{1},y_{2})+iv_{1}(x_{1},x_{2},y_{1},y_{2})$
$\displaystyle f_{2}(z_{1},z_{2})$
$\displaystyle=u_{2}(x_{1},x_{2},y_{1},y_{2})+iv_{2}(x_{1},x_{2},y_{1},y_{2})$

we have:

$f^{\prime}J=Jf^{\prime}\text{ with }J:=\begin{pmatrix}\mbox{\Large
0}&\vline&\begin{matrix}1&0&0&0\\\ 0&1&0&0\\\ 0&0&1&0\\\
0&0&0&1\end{matrix}\\\ \hline\cr\begin{matrix}-1&0&0&0\\\ 0&-1&0&0\\\
0&0&-1&0\\\ 0&0&0&-1\end{matrix}\par&\vline&\mbox{\Large 0}\\\ \end{pmatrix}$

In general, if we have $n$ multivariate complex functions with $n$ complex
variables, then its Jacobian matrix is given by the following block matrix:

$J:=\left[\begin{array}[]{c|c}\mbox{\Large 0}&\mathbf{I_{2n}}\\\
\hline\cr-\mathbf{I_{2n}}&\mbox{\Large 0}\end{array}\right]$ (4)

##  3 Special cases

We treat in some detail cases that are either relevant from the point of view
of reaction kinetics or are not too complicated allowing us to carry out the
calculations symbolically up to the end.

###  3.1 Linear systems, first-order reactions

Let us apply the method in a case when it is not needed, but transparently
shows the essential steps.

  1. 1.

The general form of two-variable linear differential equations is:

$\dot{x}=a+bx+cy,\quad\dot{y}=A+Bx+Cy$

with the real coefficients $a,b,c,A,B,C.$ Assuming that the right-hand sides
of the above differential equations are real and imaginary parts of the
differentiable complex function $f,$ i.e. they are satisfying the Cauchy-
Riemann conditions, we get $C=b,B=-c,$ or

$\dot{x}=a+bx+cy,\quad\dot{y}=A-cx+by.$

Without excluding relevant interesting cases we assume $b^{2}+c^{2}\neq 0.$ In
this very special (linear) case, one can solve the system of equations (1)
with the initial condition

$x(0)=x_{0},\quad y(0)=y_{0}$ (5)

as

$\displaystyle x(t)$
$\displaystyle=\frac{e^{bt}(\cos(ct)X+\sin(ct)Y)-ab+Ac}{b^{2}+c^{2}},$ (6)
$\displaystyle y(t)$ $\displaystyle=\frac{e^{bt}(\cos(ct)Y-\sin(ct)X)-ac-
Ab}{b^{2}+c^{2}},$ (7)

with

$X:=ab-Ac+x_{0}(b^{2}+c^{2}),\quad Y:=ac+Ab+y_{0}(b^{2}+c^{2}).$ (8)

We shall not be as lucky as this later.

  2. 2.

The equation for $z:=x+iy$ if $x$ and $y$ obeys (1) and (5) is as follows:

$\dot{z}=a+iA+(b-ic)z,z(0)=x_{0}+iy_{0}.$ (9)

  3. 3.

The solution of (9) can easily be found such as

$z(t)=e^{(b-ic)t}\left(z(0)+\frac{a+iA}{b-ic}\right)-\frac{a+iA}{b-ic},$ (10)

  4. 4.

that can be rewritten using Euler’s formula and Lemma 1 as

$\displaystyle x(t)$ $\displaystyle=e^{bt}\left(\cos(bt)(x_{0}+\frac{ab-
Ac}{b^{2}+c^{2}})+\sin(bt)(y_{0}+\frac{ac+Ab}{b^{2}+c^{2}})\right)-\frac{ab-
Ac}{b^{2}+c^{2}},$ (11) $\displaystyle y(t)$
$\displaystyle==e^{bt}\left(\sin(bt)(x_{0}+\frac{ab-
Ac}{b^{2}+c^{2}})+\cos(bt)(y_{0}+\frac{ac+Ab}{b^{2}+c^{2}})\right)-\frac{ac+Ab}{b^{2}+c^{2}}$
(12)

which is the same as given in (6)–(7).

  5. 5.

To have a kinetic (or, Hungarian) polynomial differential equation no negative
cross-effect ([14]) can be present, thus we have the restrictions $a,c,A,B\geq
0.$ These, together with the Cauchy–Riemann conditions imply that $c=0.$
Again, we are excluding here the trivial case $b=0.$ Specializing (6)–(7) to
this case gives

$\displaystyle x(t)$ $\displaystyle=(x_{0}+\frac{a}{b})e^{bt}-\frac{a}{b},$
(13) $\displaystyle y(t)$
$\displaystyle=(y_{0}+\frac{A}{b})e^{bt}-\frac{A}{b}.\quad(t\in\mathbb{R}).$
(14)

  6. 6.

Now we find a realization of the differential equation (of the one that is not
only of the Cauchy–Riemann type but also kinetic)

$\dot{x}=a+bx\quad\dot{y}=A+by$ (15)

by finding a reaction inducing it. A solution with a minimal number of
complexes and reaction steps is seen in Fig. 1 when $b<0.$

Figure 1: A reaction with a minimal number of complexes and reaction steps
inducing the kinetic differential equation (15) with $b<0$.

When $b>0,$ we have the reaction in Fig. 2

Figure 2: A reaction with a minimal number of complexes and reaction steps
inducing the kinetic differential equation (15) with $b>0$.

###  3.2 Second-degree systems, second-order reactions

Let us apply the method.

  1. 1.

The general form of two-variable differential equations with a second-degree
right-hand side is:

$\dot{x}=a+bx+cy+dx^{2}+exy+fy^{2},\quad\dot{y}=A+Bx+Cy+Dx^{2}+Exy+Fy^{2}.$

Assuming that components of the right-hand side of the above differential
equation form the real and imaginary parts of a differentiable complex
function $f$ satisfying the Cauchy–Riemann conditions, we have:

$b+2dx+ey=C+Ex+2Fy,\quad c+ex+2fy=-B-2Dx-Ey.$

Equating the coefficients on the right–hand side and those of the left
hand–side in each of the equations above we get:

$B=-c,C=b,D=-\frac{e}{2},E=2d=-2f,F=\frac{e}{2}.$

Since $E=2d$ and $E=-2f$ it follows that $f=-d,$ thus we can eliminate $f.$ We
obtain the following differential equation:

$\displaystyle\dot{x}=a+bx+cy+dx^{2}+exy-dy^{2},$ (16)
$\displaystyle\dot{y}=A-cx+by-\frac{e}{2}x^{2}+2dxy+\frac{e}{2}y^{2}.$ (17)

Now it does not seem possible to solve this system with classical methods.

  2. 2.

However, the equation for $z:=x+iy$ if $x$ and $y$ obeys (16)–(17) and (5) is
as follows:

$\dot{z}=a+iA+(b-ic)z+(d-i\frac{e}{2})z^{2},\quad z(0)=x_{0}+iy_{0}.$ (18)

  3. 3.

The solution of (18) can be found, see the Appendix, Subsection 6.2.

  4. 4.

The final result giving the solution of the real system of equations (16)–(17)
together with the initial condition (5) can be given, but it is so complicated
that we only give it in special cases.

  5. 5.

Because of the absence of a negative cross-effect, we have the restrictions
$a,A,\geq 0,d,e\leq 0$ and $c=0.$ Specializing the general solution to this
case gives the solution of the kinetic differential equation.

  6. 6.

Now we find a realization of the differential equation (of the one that is not
only of the Cauchy–Riemann type but also kinetic)

$\dot{x}=a+bx+dx^{2}+exy-
dy^{2},\quad\dot{y}=A+by-\frac{e}{2}x^{2}+2dxy+\frac{e}{2}y^{2}+dy.$ (19)

with $a,A\geq 0;d,e\leq 0,$ and $b<0,b=0,b>0.$

A relatively simple realization is the canonic realization, see Fig. 3.

Figure 3: The canonic realization of Eq. (19) for $b<0,b=0,b>0$ respectively.

Let us make a few remarks on the realization problem. The right-hand sides of
the equation (19) are sums of monomials with the exponents

$\begin{bmatrix}0\\\ 0\end{bmatrix},\begin{bmatrix}1\\\
0\end{bmatrix},\begin{bmatrix}0\\\ 1\end{bmatrix},\begin{bmatrix}2\\\
0\end{bmatrix},\begin{bmatrix}1\\\ 1\end{bmatrix},\begin{bmatrix}0\\\
2\end{bmatrix}.$

The complexes corresponding to these vectors, i.e.

$\ce{0},\ce{X},\ce{Y},\ce{2X},\ce{X+Y},\ce{2Y}$

will surely play a role as a reactant complex. As to the realizations, we may
have different goals:

  1. 1.

No more complexes be used.

  2. 2.

Let the number of reaction steps be as small as possible.

  3. 3.

Let the product complexes (if we are forced to add new ones) have integer
coefficients.

  4. 4.

Let the realization be reversible, weakly reversible, complex balanced,
detailed balanced, let it have as small a deficiency as possible.

There is abundant literature on the topics from [14] through [7] until recent
papers. Here we shall be content with the canonic realizations above, because
they are not too bad in the sense that they only use two additional complexes
beyond the absolutely necessary ones.

###  3.3 Third-degree systems, third-order reactions

As third-order reactions are not always excluded, it is worth studying this
case, as well. Furthermore, the number of free (adjustable) parameters is much
larger here, therefore the explicit solutions can be used for many reactions.

  1. 1.

The general form of two-variable differential equations with a third-degree
right-hand side is:

$\displaystyle\dot{x}$
$\displaystyle=a+bx+cy+dx^{2}+exy+fy^{2}+gx^{3}+hx^{2}y+jxy^{2}+ky^{3},$
$\displaystyle\dot{y}$
$\displaystyle=A+Bx+Cy+Dx^{2}+Exy+Fy^{2}+Gx^{3}+Hx^{2}y+Jxy^{2}+Ky^{3}.$

The Cauchy–Riemann conditions mean here:

$\displaystyle b+2dx+ey+3gx^{2}+2hxy+jy^{2}$
$\displaystyle=C+Ex+2Fy+Hx^{2}+2Jxy+3Ky^{2}$ (20) $\displaystyle
c+ex+2fy+hx^{2}+2jxy+3ky^{2}$ $\displaystyle=-(B+2Dx+Ey+3Gx^{2}+2Hxy+Jy^{2}).$
(21)

Comparing the coefficients on the left-hand side and right-hand side of both
equations we obtain the following results on coefficients:

$\displaystyle B=-c,C=b,D=-\frac{e}{2},E=2d,E=-2f,F=\frac{e}{2},$ (22)
$\displaystyle G=-\frac{h}{3},H=3g,H=-j,J=h,J=-3k,K=\frac{j}{3}.$ (23)

Since $E=2d$ and $E=-2f$ it follows that $f=-d,$ thus we can eliminate $f.$
Similarly, $H=3g$ and $H=-j$ imply $j=-3g;$ $J=h$ and $J=-3k$ imply
$k=-\frac{h}{3};$ therefore

$\displaystyle f$ $\displaystyle=-d,j=-3g,k=-\frac{h}{3};$ (24) $\displaystyle
B$ $\displaystyle=-c,C=b,D=-\frac{e}{2},E=2d,F=\frac{e}{2},$ (25)
$\displaystyle G$ $\displaystyle=-\frac{h}{3},H=3g,J=h,K=-g.$ (26)

We obtain the following differential equation:

$\displaystyle\dot{x}$ $\displaystyle=a+bx+cy+dx^{2}+exy-
dy^{2}+gx^{3}+hx^{2}y-3gxy^{2}-\frac{h}{3}y^{3},$ (27) $\displaystyle\dot{y}$
$\displaystyle=A-cx+by-\frac{e}{2}x^{2}+2dxy+\frac{e}{2}y^{2}-\frac{h}{3}x^{3}+3gx^{2}y+hxy^{2}-gy^{3}.$
(28)

Again, it does not seem possible to solve this system with classical methods.

  2. 2.

However, the equation for $z:=x+iy$ if $x$ and $y$ obeys (28) and (5) is as
follows:

$\dot{z}=(h-i\frac{g}{3})z^{3}+(d-i\frac{e}{2})z^{2}+(b-ic)z+(a+iA),\quad
z(0)=x_{0}+iy_{0}.$ (29)

  3. 3.

The implicit solution of (29) can be found such as

$\displaystyle(z(t)-z_{1})(z(t)-z_{2})(z(t)-z_{3})$ (30)
$\displaystyle=(z(0)-z_{1})(z(0)-z_{2})(z(0)-z_{3})\exp{(z_{1}-z_{2})(z_{1}-z_{3})(z_{2}-z_{3})(a-i\frac{b}{3})t},$
(31)

where $z_{1},z_{2}$ and $z_{3}$ are the roots of the polynomial

$(z\mapsto(h-i\frac{g}{3})z^{3}+(d-i\frac{e}{2})z^{2}+(b-ic)z+(a+iA),$

and $z(0):=x_{0}+iy_{0}.$ To get an explicit solution, the application of the
Cardano formula is needed, and the result would not be really useful. Even the
specialization to the kinetic case is not simpler. To show the steps of the
process, we do all the calculations in Subsection 6.4 of the Appendix for a
few cases.

  4. 4.

Now we find a realization of the differential equation of the one that is not
only of the Cauchy–Riemann type but also kinetic. The above equation is
kinetic if $c=0;a,A\geq 0;d,h,e\leq 0.$ However, there are no restrictions to
the signs of $b$ and $g.$

$\displaystyle\dot{x}$ $\displaystyle=a+bx+dx^{2}+exy-
dy^{2}+gx^{3}+hx^{2}y-3gxy^{2}-\frac{h}{3}y^{3},$ (32) $\displaystyle\dot{y}$
$\displaystyle=A+by-\frac{e}{2}x^{2}+2dxy+\frac{e}{2}y^{2}-\frac{h}{3}x^{3}+3gx^{2}y+hxy^{2}-gy^{3}.$
(33)

We shall treat the case $b<0$ and $g<0$ in detail and leave the three
remaining cases $b\leq 0,g\geq 0$, and $b\geq 0,g\leq 0$, and $b\geq 0,g\geq
0$, and also the cases when some of the constants are zero, to the reader. A
relatively simple realization is the canonic realization: Note that beyond the
complexes corresponding to the exponents of the monomials on the right-hand
sides we needed four more complexes: 2X + 2Y,X + 3 Y, 3 X + Y, and 4Y. Fig. 4
shows the FHJ graph for the canonic realization of the Eq. (33).

Figure 4: Canonic realization in the third order case when $b<0,g<0$.

Finally, we mention that it is mathematically possible to apply our method in
the fourth-order case, but it is not worth the price.

###  3.4 Rth-degree systems

The general case is interesting from the point of view of differential
equations, certainly not from the point of view of formal reaction kinetics.
Some of the earlier results can be summarized and generalized in a nice
statement.

######  Theorem 1.

Let $R\in\mathbb{N}_{0},$ and suppose that the polynomials $u,v$ defined as

$\displaystyle u(x,y)$
$\displaystyle:=\sum_{r=0}^{R}\sum_{s=0}^{r}a_{r}^{s}x^{r-s}y^{s},\quad
v(x,y):=\sum_{r=0}^{R}\sum_{s=0}^{r}A_{r}^{s}x^{r-s}y^{s}$
$\displaystyle(a_{r}^{s},A_{r}^{s}\in\mathbb{R},s=0,1,2,\dots,r;r=1,2,\dots,R)$

satisfy the Cauchy–Riemann equations (2). Then, the system of differential
equations

$\dot{x}=u\circ(x,y),\quad\dot{y}=v\circ(x,y)$ (34)

can be simplified to

$\displaystyle\dot{x}$
$\displaystyle=a^{0}_{0}+\sum_{r=1}^{R}\left(a^{0}_{r}\mathfrak{Re}[(x+iy)^{r}]+\frac{a^{1}_{r}}{r}\mathfrak{Im}[(x+iy)^{r}]\right),$
(35) $\displaystyle\dot{y}$
$\displaystyle=A^{0}_{0}+\sum_{r=1}^{R}\left(a^{0}_{r}\mathfrak{Im}[(x+iy)^{r}]-\frac{a^{1}_{r}}{r}\mathfrak{Re}[(x+iy)^{r}]\right)$
(36)

Furthermore, for the complex-valued function $z:=x+iy$ one has

$\dot{z}=a^{0}_{0}+iA^{0}_{0}+\sum_{r=1}^{R}(a_{r}^{0}-i\frac{a_{r}^{1}}{r})z^{r}.$

######  Remark 1.

Note that the number of independent parameters left is $2R+1,$ showing that
the Cauchy–Riemann equations mean strong restrictions for the relations
between the parameters. Furthermore, the coefficients of the terms with
different degrees are independent, therefore it is enough to do the
calculations with a single term.

Proof.

  1. 1.

A direct comparison of the coefficients in the Cauchy–Riemann equations gives

$\displaystyle A^{s}_{r}$ $\displaystyle=\frac{r-s+1}{s}a^{s-1}_{r}\quad$
$\displaystyle(s=1,2,3,\dots,r;r=1,2,3,\dots,R),$ (37) $\displaystyle
A^{s}_{r}$ $\displaystyle=-\frac{s+1}{r-s}a^{s+1}_{r}\quad$
$\displaystyle(s=0,1,2,\dots,r-1;r=1,2,3,\dots,R).$ (38)

Eqs.(37) and (38) imply

$\displaystyle
a^{s}_{r}=-\frac{(r-s+2)(r-s+1)}{s(s-1)}a^{s-2}_{r}\quad(s=2,3,\dots,r;r=2,3,\dots,R).$
(39)

The recursive application of (39) (with $i:=\sqrt{-1}$) leads to

$\displaystyle a^{s}_{r}$ $\displaystyle=i^{s}\binom{r}{s}a^{0}_{r},\quad$
$\displaystyle(s=0,2,4,6,\dots;r=2,3,\dots,R);$ (40) $\displaystyle a^{s}_{r}$
$\displaystyle=i^{s-1}\frac{1}{r}\binom{r}{s}a^{1}_{r},\quad$
$\displaystyle(s=1,3,5,\dots;r=2,3,\dots,R).$ (41)

Therefore the $r\mathrm{th}$ term of the first equation of the system (34) can
be written in the following ways:

$\displaystyle
i^{0}\binom{r}{0}a^{0}_{r}x^{r}+i^{0}\frac{1}{r}\binom{r}{1}a^{1}_{r}x^{r-1}y+i^{2}\binom{r}{2}a^{0}_{r}x^{r-2}y^{2}+i^{2}\frac{1}{r}\binom{r}{3}a^{1}_{r}x^{r-3}y^{3}$
(42)
$\displaystyle+i^{4}\binom{r}{4}a^{0}_{r}x^{r-4}y^{4}+i^{4}\frac{1}{r}\binom{r}{5}a^{1}_{r}x^{r-5}y^{5}+\dots$
(43)
$\displaystyle=a^{0}_{r}\left(i^{0}\binom{r}{0}x^{r}+i^{2}\binom{r}{2}x^{r-2}y^{2}+i^{4}\binom{r}{4}x^{r-4}y^{4}+\dots\right)$
(44)
$\displaystyle-i\frac{a^{1}_{r}}{r}\left(i^{1}\binom{r}{1}x^{r-1}y+i^{3}\binom{r}{3}x^{r-3}y^{3}+i^{5}\binom{r}{5}x^{r-5}y^{5}+\dots\right)$
(45)
$\displaystyle=a^{0}_{r}\mathfrak{Re}[(x+iy)^{r}]+\frac{a^{1}_{r}}{r}\mathfrak{Im}[(x+iy)^{r}].$
(46)

  2. 2.

Similarly, the recursive application of (37) leads to

$\displaystyle A^{s}_{r}$
$\displaystyle=i^{s-2}\frac{1}{r}\binom{r}{s}a^{1}_{r},\quad(s=2,4,6,\dots;r=2,3,\dots,R;)$
(47) $\displaystyle A^{s}_{r}$
$\displaystyle=i^{s-1}\binom{r}{s}a^{0}_{r},\quad(s=1,3,5,\dots;r=2,3,\dots,R.)$
(48)

Therefore the $r\mathrm{th}$ term of the second equation of the system (34)
can be written in the following ways:

$\displaystyle-\frac{1}{r}\binom{r}{0}a^{1}_{r}x^{r}+\binom{r}{1}a^{0}_{r}x^{r-1}y+\frac{1}{r}\binom{r}{2}a^{1}_{r}x^{r-2}y^{2}-\binom{r}{3}a^{0}_{r}x^{r-3}y^{3}$
(49)
$\displaystyle-\frac{1}{r}\binom{r}{4}a^{1}_{r}x^{r-4}y^{4}+\binom{r}{5}a^{0}_{r}x^{r-5}y^{5}+\dots$
(50)
$\displaystyle=a^{0}_{r}\left(i^{0}\binom{r}{1}x^{r-1}y+i^{2}\binom{r}{3}x^{r-3}y^{3}+i^{4}\binom{r}{5}x^{r-5}y^{5}+\dots\right)$
(51)
$\displaystyle-\frac{a^{1}_{r}}{r}\left(i^{0}\binom{r}{0}x^{r}+i^{2}\binom{r}{2}x^{r-2}y^{2}+i^{4}\binom{r}{4}x^{r-4}y^{4}+\dots\right)$
(52)
$\displaystyle=a^{0}_{r}\mathfrak{Im}[(x+iy)^{r}]-\frac{a^{1}_{r}}{r}\mathfrak{Re}[(x+iy)^{r}].$
(53)

  3. 3.

Finally, the equation for $z:=x+iy$ if $x$ and $y$ obeys (2) and (5) is as
follows:

$\displaystyle\dot{z}$
$\displaystyle=a^{0}_{r}\mathfrak{Re}[(x+iy)^{r}]-i\frac{a^{1}_{r}}{r}\mathfrak{Im}[(x+iy)^{r}]+i(a^{0}_{r}\mathfrak{Im}[(x+iy)^{r}]+\frac{a^{1}_{r}}{r}\mathfrak{Re}[(x+iy)^{r}])$
(54) $\displaystyle=\sum_{r=0}^{R}(a_{r}^{0}-\frac{a_{r}^{1}}{r})z^{r},\quad
z(0)=x_{0}+iy_{0}.$ (55)

##  4 Application of multivariate complex functions

Here we extend our method to treat differential equations in multiple
variables [27, 32]. Suppose, $u$ and $v$ are the real and imaginary parts of
the multivariate complex function $f.$ Then, the Cauchy–Riemann equations are
as follows.

$\frac{\partial u}{\partial x_{j}}=\frac{\partial v}{\partial
y_{j}},\quad\frac{\partial u}{\partial y_{j}}=-\frac{\partial v}{\partial
x_{j}}.$ (56)

Consider more specifically the case when we have two complex functions and two
complex variables such as

$\displaystyle f_{1}(z_{1},z_{2})$
$\displaystyle=u_{1}(x_{1},x_{2},y_{1},y_{2})+iv_{1}(x_{1},x_{2},y_{1},y_{2}),$
$\displaystyle f_{2}(z_{1},z_{2})$
$\displaystyle=u_{2}(x_{1},x_{2},y_{1},y_{2})+iv_{2}(x_{1},x_{2},y_{1},y_{2}),$

then the equalities are:

$\displaystyle\frac{\partial u_{1}}{\partial x_{1}}=\frac{\partial
v_{1}}{\partial y_{1}},\quad\frac{\partial u_{1}}{\partial
y_{1}}=-\frac{\partial v_{1}}{\partial x_{1}},\quad\frac{\partial
u_{1}}{\partial x_{2}}=\frac{\partial v_{1}}{\partial
y_{2}},\quad\frac{\partial u_{1}}{\partial y_{2}}=-\frac{\partial
v_{1}}{\partial x_{2}},$ (57) $\displaystyle\frac{\partial u_{2}}{\partial
x_{1}}=\frac{\partial v_{2}}{\partial y_{1}},\quad\frac{\partial
u_{2}}{\partial y_{1}}=-\frac{\partial v_{2}}{\partial
x_{1}},\quad\frac{\partial u_{2}}{\partial x_{2}}=\frac{\partial
v_{2}}{\partial y_{2}},\quad\frac{\partial u_{2}}{\partial
y_{2}}=-\frac{\partial v_{2}}{\partial x_{2}},$ (59)

and we would have a system of differential equations in real quantities:

$\displaystyle\dot{x}_{1}=u_{1}(x_{1},x_{2},y_{1},y_{2}),\quad$
$\displaystyle\dot{y}_{1}=v_{1}(x_{1},x_{2},y_{1},y_{2}),$ (60)
$\displaystyle\dot{x}_{2}=u_{2}(x_{1},x_{2},y_{1},y_{2}),\quad$
$\displaystyle\dot{y}_{2}=v_{2}(x_{1},x_{2},y_{1},y_{2}).$ (61)

If we were able to reduce the number of equations to two (in complex
variables), it would be a small step forward. Let us see the linear case
first, to get some experience. The system of real equations is:

$\displaystyle\dot{x}_{1}$ $\displaystyle=ax_{1}+bx_{2}+cy_{1}+dy_{2},\quad$
$\displaystyle\dot{y}_{1}$ $\displaystyle=ex_{1}+fx_{2}+gy_{1}+hy_{2},$
$\displaystyle\dot{x}_{2}$ $\displaystyle=Ax_{1}+Bx_{2}+Cy_{1}+Dy_{2},\quad$
$\displaystyle\dot{y}_{2}$ $\displaystyle=Ex_{1}+Fx_{2}+Gy_{1}+Hy_{2},$

and the Cauchy–Riemann equations imply

$a=g,c=-e,b=h,d=-f,\quad A=G,C=-E,B=H,D=-F.$

The complex equations are:

$\dot{z}_{1}=(a-ci)z_{1}+(b-di)z_{2},\quad\dot{z}_{2}=(A-Ci)z_{1}+(B-Di)z_{2}.$

Promising. But not from the point of reaction kinetics. Why? To have a kinetic
equation one should have

$c=d=e=f=C=D=E=F=0,$

thus the system falls apart:

$\displaystyle\dot{x}_{1}$ $\displaystyle=ax_{1}+bx_{2},$
$\displaystyle\dot{x}_{2}$ $\displaystyle=Ax_{1}+Bx_{2},$
$\displaystyle\dot{y}_{1}$ $\displaystyle=gy_{1}+hy_{2},$
$\displaystyle\dot{y}_{2}$ $\displaystyle=Gy_{1}+Hy_{2},$

In the quadratic case, we start from

$\displaystyle\dot{x}_{1}$
$\displaystyle=j_{1}x_{1}^{2}+j_{2}x_{2}^{2}+j_{3}y_{1}^{2}+j_{4}y_{2}^{2}$
$\displaystyle+j_{5}x_{1}x_{2}+j_{6}y_{1}y_{2}+j_{7}x_{1}y_{1}+j_{8}x_{1}y_{2}+j_{9}x_{2}y_{1}+j_{10}x_{2}y_{2}$
$\displaystyle+j_{11}x_{1}+j_{12}x_{2}+j_{13}y_{1}+j_{14}y_{2}+j_{15},$ (62)
$\displaystyle\dot{y}_{1}$
$\displaystyle=k_{1}x_{1}^{2}+k_{2}x_{2}^{2}+k_{3}y_{1}^{2}+k_{4}y_{2}^{2}$
$\displaystyle+k_{5}x_{1}x_{2}+k_{6}y_{1}y_{2}+k_{7}x_{1}y_{1}+k_{8}x_{1}y_{2}+k_{9}x_{2}y_{1}+k_{10}x_{2}y_{2}$
$\displaystyle+k_{11}x_{1}+k_{12}x_{2}+k_{13}y_{1}+k_{14}y_{2}+k_{15},$ (63)
$\displaystyle\dot{x}_{2}$
$\displaystyle=J_{1}x_{1}^{2}+J_{2}x_{2}^{2}+J_{3}y_{1}^{2}+J_{4}y_{2}^{2}$
$\displaystyle+J_{5}x_{1}x_{2}+J_{6}y_{1}y_{2}+J_{7}x_{1}y_{1}+J_{8}x_{1}y_{2}+J_{9}x_{2}y_{1}+J_{10}x_{2}y_{2}$
$\displaystyle+J_{11}x_{1}+J_{12}x_{2}+J_{13}y_{1}+J_{14}y_{2}+J_{15},$ (64)
$\displaystyle\dot{y}_{2}$
$\displaystyle=K_{1}x_{1}^{2}+K_{2}x_{2}^{2}+K_{3}y_{1}^{2}+K_{4}y_{2}^{2}$
$\displaystyle+K_{5}x_{1}x_{2}+K_{6}y_{1}y_{2}+K_{7}x_{1}y_{1}+K_{8}x_{1}y_{2}+K_{9}x_{2}y_{1}+K_{10}x_{2}y_{2}$
$\displaystyle\quad+K_{11}x_{1}+K_{12}x_{2}+K_{13}y_{1}+K_{14}y_{2}+K_{15}.$
(65)

and the Cauchy–Riemann equations imply

$\displaystyle k_{1}=-k_{3}=-\frac{j_{7}}{2},$ $\displaystyle
k_{2}=-k_{4}=-\frac{j_{10}}{2},$ $\displaystyle-k_{5}=k_{6}=j_{8}=j_{9},$
$\displaystyle k_{7}=2j_{1}=-2j_{3},$ $\displaystyle
k_{8}=k_{9}=j_{5}=-j_{6},$ $\displaystyle k_{10}=2j_{2}=-2j_{4},$
$\displaystyle k_{11}=-j_{13},$ $\displaystyle k_{12}=-j_{14},$ $\displaystyle
k_{13}=j_{11},$ $\displaystyle k_{14}=j_{12},$ $\displaystyle
K_{1}=-K_{3}=-\frac{J_{7}}{2},$ $\displaystyle
K_{2}=-K_{4}=-\frac{J_{10}}{2},$ $\displaystyle K_{5}=-K_{6}=-J_{8}=-J_{9},$
$\displaystyle K_{7}=2J_{1}=-2J_{3},$ $\displaystyle
K_{8}=K_{9}=J_{5}=-J_{6},$ $\displaystyle K_{10}=2J_{2}=-2J_{4},$
$\displaystyle K_{11}=-J_{13},$ $\displaystyle K_{12}=-J_{14},$ $\displaystyle
K_{13}=J_{11},$ $\displaystyle K_{14}=J_{12}.$

To use as few parameters as possible: altogether 2*11 (instead of 2*30), and
these are

$j_{1},j_{2},j_{5},j_{7},j_{8},j_{10},j_{11},j_{12},j_{13},j_{14},j_{15}$

and

$J_{1},J_{2},J_{5},J_{7},J_{8},J_{10},J_{11},J_{12},J_{13},J_{14},J_{15},$

we shall use these relations.

$\displaystyle k_{1}=-\frac{j_{7}}{2},$ $\displaystyle
k_{2}=-\frac{j_{10}}{2},$ $\displaystyle k_{3}=\frac{j_{7}}{2},$
$\displaystyle k_{4}=\frac{j_{10}}{2},$ $\displaystyle k_{5}=-j_{8},$
$\displaystyle k_{6}=j_{8},$ $\displaystyle k_{7}=2j_{1},$ $\displaystyle
k_{8}=j_{5},$ $\displaystyle k_{9}=j_{5},$ $\displaystyle k_{10}=2j_{2}$
$\displaystyle k_{11}=-j_{13},$ $\displaystyle k_{12}=-j_{14},$ $\displaystyle
k_{13}=j_{11},$ $\displaystyle k_{14}=j_{12},$ $\displaystyle j_{3}=-j_{1},$
$\displaystyle j_{4}=-j_{2}$ $\displaystyle j_{6}=-j_{5}$ $\displaystyle
j_{9}=j_{8},$ $\displaystyle K_{1}=-\frac{J_{7}}{2},$ $\displaystyle
K_{2}=-\frac{J_{10}}{2},$ $\displaystyle K_{3}=\frac{J_{7}}{2},$
$\displaystyle K_{4}=\frac{J_{10}}{2},$ $\displaystyle K_{5}=-J_{8},$
$\displaystyle K_{6}=J_{8},$ $\displaystyle K_{7}=2J_{1},$ $\displaystyle
K_{8}=J_{5}$ $\displaystyle K_{9}=J_{5},$ $\displaystyle K_{10}=2J_{2},$
$\displaystyle K_{11}=-J_{13},$ $\displaystyle K_{12}=-J_{14},$ $\displaystyle
K_{13}=J_{11},$ $\displaystyle K_{14}=J_{12},$ $\displaystyle J_{3}=-J_{1},$
$\displaystyle J_{4}=-J_{2}$ $\displaystyle J_{6}=-J_{5}$ $\displaystyle
J_{9}=J_{8}.$

Thus, our system obeys the Cauchy–Riemann conditions are:

$\displaystyle\dot{x}_{1}$
$\displaystyle=j_{1}x_{1}^{2}+j_{2}x_{2}^{2}-j_{1}y_{1}^{2}-j_{2}y_{2}^{2}$
$\displaystyle+j_{5}x_{1}x_{2}-j_{5}y_{1}y_{2}+j_{7}x_{1}y_{1}+j_{8}x_{1}y_{2}+j_{8}x_{2}y_{1}+j_{10}x_{2}y_{2}$
$\displaystyle+j_{11}x_{1}+j_{12}x_{2}+j_{13}y_{1}+j_{14}y_{2}+j_{15},$ (66)
$\displaystyle\dot{y}_{1}$
$\displaystyle=-\frac{j_{7}}{2}x_{1}^{2}-\frac{j_{10}}{2}x_{2}^{2}+\frac{j_{7}}{2}y_{1}^{2}+\frac{j_{10}}{2}y_{2}^{2}$
$\displaystyle-
j_{8}x_{1}x_{2}+j_{8}y_{1}y_{2}+2j_{1}x_{1}y_{1}+j_{5}x_{1}y_{2}+j_{5}x_{2}y_{1}+2j_{2}x_{2}y_{2}$
$\displaystyle-j_{13}x_{1}-j_{14}x_{2}+j_{11}y_{1}+j_{12}y_{2}+k_{15},$ (67)
$\displaystyle\dot{x}_{2}$
$\displaystyle=J_{1}x_{1}^{2}+J_{2}x_{2}^{2}-J_{1}y_{1}^{2}-J_{2}y_{2}^{2}$
$\displaystyle+J_{5}x_{1}x_{2}-J_{5}y_{1}y_{2}+J_{7}x_{1}y_{1}+J_{8}x_{1}y_{2}+J_{8}x_{2}y_{1}+J_{10}x_{2}y_{2}$
$\displaystyle+J_{11}x_{1}+J_{12}x_{2}+J_{13}y_{1}+J_{14}y_{2}+J_{15},$ (68)
$\displaystyle\dot{y}_{2}$
$\displaystyle=-\frac{J_{7}}{2}x_{1}^{2}-\frac{J_{10}}{2}x_{2}^{2}+\frac{J_{7}}{2}y_{1}^{2}+\frac{J_{10}}{2}y_{2}^{2}$
$\displaystyle-
J_{8}x_{1}x_{2}+J_{8}y_{1}y_{2}+2J_{1}x_{1}y_{1}+J_{5}x_{1}y_{2}+J_{5}x_{2}y_{1}+2J_{2}x_{2}y_{2}$
$\displaystyle-J_{13}x_{1}-J_{14}x_{2}+J_{11}y_{1}+J_{12}y_{2}+K_{15}.$ (69)

Therefore we have simpler complex equations as follows:

$\displaystyle\dot{z}_{1}$
$\displaystyle=(j_{1}-i\frac{j_{7}}{2})z_{1}^{2}+(j_{2}-i\frac{j_{10}}{2})z_{2}^{2}+(j_{5}-ij_{8})z_{1}z_{2}+(j_{11}-ij_{13})z_{1}$
$\displaystyle\quad+(j_{12}-ij_{14})z_{2}+(j_{15}+ik_{15}),$ (70)
$\displaystyle\dot{z}_{2}$
$\displaystyle=(J_{1}-i\frac{J_{7}}{2})z_{1}^{2}+(J_{2}-i\frac{J_{10}}{2})z_{2}^{2}+(J_{5}-iJ_{8})z_{1}z_{2}+(J_{11}-iJ_{13})z_{1}$
$\displaystyle\quad+(J_{12}-iJ_{14})z_{2}+(J_{15}+iK_{15}).$ (71)

Assuming also that our system is kinetic, we have

$\displaystyle j_{1},j_{7}\leq
0;j_{2}=j_{5}=j_{8}=j_{10}=j_{13}=j_{14}=0;j_{12},j_{15},k_{15}\geq 0$ (72)
$\displaystyle J_{2},J_{10}\leq
0;J_{1}=J_{5}=J_{7}=J_{8}=J_{13}=J_{14}=0;J_{11},J_{15},K_{15}\geq 0$ (73)
$\displaystyle\dot{x}_{1}$
$\displaystyle=j_{1}x_{1}^{2}-j_{1}y_{1}^{2}+j_{7}x_{1}y_{1}+j_{11}x_{1}+j_{12}x_{2}+j_{15},$
(74) $\displaystyle\dot{y}_{1}$
$\displaystyle=-\frac{j_{7}}{2}x_{1}^{2}+\frac{j_{7}}{2}y_{1}^{2}+2j_{1}x_{1}y_{1}+j_{11}y_{1}+j_{12}y_{2}+k_{15},$
(75) $\displaystyle\dot{x}_{2}$
$\displaystyle=J_{2}x_{2}^{2}-J_{2}y_{2}^{2}+J_{10}x_{2}y_{2}+J_{11}x_{1}+J_{12}x_{2}+J_{15},$
(76) $\displaystyle\dot{y}_{2}$
$\displaystyle=-\frac{J_{10}}{2}x_{2}^{2}+\frac{J_{10}}{2}y_{2}^{2}+2J_{2}x_{2}y_{2}+J_{11}y_{1}+J_{12}y_{2}+K_{15}.$
(77)

The complex equations will also be simpler:

$\displaystyle\dot{z}_{1}$
$\displaystyle=(j_{1}-i\frac{j_{7}}{2})z_{1}^{2}+j_{11}z_{1}+j_{12}z_{2}+(j_{15}+ik_{15}),$
(78) $\displaystyle\dot{z}_{2}$
$\displaystyle=(J_{1}-i\frac{J_{7}}{2})z_{2}^{2}+J_{11}z_{1}+J_{12}z_{2}+(J_{15}+iK_{15}).$
(79)

Here, we have reduced the four–variable real system into a two–variable
complex system. In general, given a $2n$–variable real system, where $n$ is a
positive integer, we can reduce it to a $n$-variable complex system using the
technique above.

If all the right-hand sides of the reduced system are homogeneous polynomials
of the same degree then a further (final) step of reduction is also possible.
We show how to start the calculations with a simple example. Suppose that we
have arrived at two equations with homogeneous second-degree polynomials:

$\dot{z}_{1}=az_{1}^{2}+bz_{1}z_{2}+cz_{2}^{2}\quad\dot{z}_{2}=Az_{1}^{2}+Bz_{1}z_{2}+Cz_{2}^{2}.$
(80)

Assuming that there is a function $Z$ such that $z_{2}=Z\circ z_{1},$ we have
$\dot{z_{2}}=Z^{\prime}\circ z_{1}\dot{z_{1}},$ and

$Z^{\prime}\circ z_{1}=\frac{\dot{z_{2}}}{\dot{z_{1}}}=\frac{A+B\frac{Z\circ
z_{1}}{z_{1}}+C(\frac{Z\circ z_{1}}{z_{1}})^{2}}{a+b\frac{Z\circ
z_{1}}{z_{1}}+c(\frac{Z\circ z_{1}}{z_{1}})^{2}}$ (81)

Let $U(\zeta):=\frac{Z(\zeta)}{\zeta},$ then $Z(\zeta)=\zeta U(\zeta),$
therefore ${Z^{\prime}}(\zeta)=\zeta{U^{\prime}}(\zeta)+U(\zeta),$ thus

$\zeta{U^{\prime}}(\zeta)+U(\zeta)=\frac{A+BU(\zeta)+C(U(\zeta))^{2}}{a+bU(\zeta)+c(U(\zeta))^{2}},$
(82)

or

$\zeta{U^{\prime}}(\zeta)=\frac{A+BU(\zeta)+C(U(\zeta))^{2}}{a+bU(\zeta)+c(U(\zeta))^{2}}-U(\zeta),$
(83)

and

${U^{\prime}}(\zeta)\frac{a+bU(\zeta)+c(U(\zeta))^{2}}{A+(B-a)U(\zeta)+(C-b)U(\zeta)^{2}-cU(\zeta)^{3}}=\frac{1}{\zeta}$
(84)

which is a separable differential equation for the function $U.$ Solving via
the Cardano formula one gets a relationship $\Phi(z_{1}(t),z_{2}(t))=0,$
(which is a first integral of the equation (80)) that does not seem to be
easily rearranged for the real and imaginary parts of the two functions
$z_{1}$ and $z_{2}.$

##  5 Discussion and Outlook

###  5.1 Further directions

####  5.1.1 Partial differential equations

Consider the system of reaction-diffusion equations

$\displaystyle\frac{\partial x(t,\mathbf{r})}{\partial
t}=u(x(t,\mathbf{r}),y(t,\mathbf{r}))+D_{x}\Delta x(t,\mathbf{r})$ (85)
$\displaystyle\frac{\partial y(t,\mathbf{r})}{\partial
t}=v(x(t,\mathbf{r}),y(t,\mathbf{r}))+D_{y}\Delta y(t,\mathbf{r}),$ (86)

Assuming (2), and $D_{x}=D_{y},$ there exists $f:=u+iv,$ furthermore one one
gets a single, complex reaction-diffusion equation for the function $z:=x+iy.$

$\frac{\partial z(t,\mathbf{r})}{\partial t}=f(z(t,\mathbf{r}))+D\Delta
z(t,\mathbf{r}).$ (87)

Note that the restriction that the two diffusion constants have the same value
$D$ is not too stringent: in aqueous solutions, this is approximately true for
many materials.

####  5.1.2 Approximations

This example will also serve as an example of the second-degree case. Here we
carry out all the calculations for a second-degree case with specified
(numerical) coefficients. One could substitute the concrete values of the
parameters into the general formula, however—as in general when solving
differential equations—it may be better to reproduce the derivation of the
results in the concrete case.

Fig. 5 shows the FHJ graph of the reaction (88).

Figure 5: FHJ graph of the reaction (88)

$\ce{2X->[1]X+Y<->[2][1]0},\quad\ce{2Y->[1]X+Y}$ (88)

The induced kinetic differential equation of the reaction (88) with the
initial concentrations $x(0)=2$ and $y(0)=1$ is

$\displaystyle\dot{x}=1-x^{2}-2xy+y^{2},\quad\dot{y}=1+x^{2}-2xy-y^{2},$ (89)
$\displaystyle x(0)=2,\quad y(0)=1$

leading to the complex differential equation

$\dot{z}=1+i+(-1+i)z^{2},\quad z(0)=2+i.$ (90)

The solution of Eq. (90) is

$\displaystyle
z(t)=\frac{(1+i)(-2ie^{2\sqrt{2}t}+3(1+\sqrt{2})e^{4\sqrt{2}t}+3-3\sqrt{2})}{-8e^{2\sqrt{2}t}+3(2+\sqrt{2})e^{4\sqrt{2}t}+6-3\sqrt{2}}\quad(t\in\mathbb{R}).$

Therefore

$\displaystyle x(t)$
$\displaystyle=\frac{2e^{2\sqrt{2}t}+3(1+\sqrt{2})e^{4\sqrt{2}t}+3-3\sqrt{2}}{-8e^{2\sqrt{2}t}+3(2+\sqrt{2})e^{4\sqrt{2}t}+6-3\sqrt{2}}\quad(t\in\mathbb{R}),$
(91) $\displaystyle y(t)$
$\displaystyle=\frac{-2e^{2\sqrt{2}t}+3(1+\sqrt{2})e^{4\sqrt{2}t}+3-3\sqrt{2}}{-8e^{2\sqrt{2}t}+3(2+\sqrt{2})e^{4\sqrt{2}t}+6-3\sqrt{2}}\quad(t\in\mathbb{R}).$
(92)

Note that the denominator in the above equations is never zero 2.

Now consider the equation

$\displaystyle\dot{x}=1-(1+\varepsilon)x^{2}-2xy+y^{2},\quad\dot{y}=1+x^{2}-2xy-y^{2},$
(93) $\displaystyle x(0)=2,\quad y(0)=1$

While Eq. (93) is an equation fulfilling the Cauchy–Riemann condition for
$\varepsilon=0,$ in general, it is not. Note that it is also a kinetic
differential equation for any real $\varepsilon$. As one has the explicit
solution of the unperturbed equation one might hope that it is a good
approximation of the solution to the perturbed equation.

Figure 6: Differences between the first coordinates of the solution of (93)
with $\varepsilon=1$ and $\varepsilon=0$.

Figure 7: Differences between the second coordinates of the solution of (93)
with $\varepsilon=1$ and $\varepsilon=0$.

Figures 7 and 7 show the differences between the first, respectively, second
coordinates of the solutions of (93) with $\varepsilon=1$ and $\varepsilon=0.$
Note that for $\varepsilon=0$ the solutions can be determined symbolically, as
it is a Cauchy–Riemann–Erugin equation.

###  5.2 Final remarks

The closest mathematical paper to ours has been written by Calogero and
Payandeh [3]. Now, we will systematically compare the two works. Their main
result can be summarized in the following statement.

######  Theorem 2.

Consider the equations Eq. (1b)–Eq. (1c) of [3]

$\displaystyle\dot{x_{1}}$
$\displaystyle=c_{11}x_{1}^{2}+c_{12}x_{1}x_{2}+c_{13}x_{2}^{2}+c_{14}x_{1}+c_{15}x_{2}+c_{16},$
(94) $\displaystyle\dot{x_{2}}$
$\displaystyle=c_{21}x_{1}^{2}+c_{22}x_{1}x_{2}+c_{23}x_{2}^{2}+c_{24}x_{1}+c_{25}x_{2}+c_{26},$
(95)

and assume that the equalities

$\displaystyle 4c_{13}c_{21}-c_{12}c_{22}=0$ (96) $\displaystyle
2(-c_{12}+2c_{23})c_{21}+(2c_{11}-c_{22})c_{22}=0$ (97) $\displaystyle
c_{24}(2c_{11}-c_{22})+2c_{21}(c_{25}-c_{14})=0$ (98) $\displaystyle
c_{12}c_{24}-2c_{15}c_{21}=0$ (99)

hold. Furthermore, assume that none of the denominators in the steps
(39a)–(44f) in [3], describing the algorithm, are different from zero. Then,
the explicit symbolic solution of (94)–(95) can be constructed by applying the
steps of the algorithm (39a)–(44f) in [3].

If (94)–(95) is a Cauchy–Riemann–Erugin system, then the coefficients fulfill
(96)–(99).

To show that this method implies our method under some conditions, we try to
obtain the above constraints by assuming that the right-hand side of (94)–(95)
satisfy Cauchy–Riemann equations:

$\displaystyle 2c_{11}-c_{22}=0,\quad 2c_{23}-c_{12}=0,\quad
c_{14}-c_{25}=0,\quad 2c_{21}=-c_{12},\quad 2c_{13}=-c_{22},\quad
c_{15}=-c_{24}.$

Simple substitutions show the statement.

It is important to note that this does not imply that all the
Cauchy–Riemann–Erugin systems can be solved by the method of Calogero and
Payandeh; for example, the linear systems cannot, because they do not satisfy
the conditions about the denominators.

In the following cases, we demonstrate that the far-reaching method presented
by Calogero and Payandeh may also lead to a dead end, awaiting further
simplification methods.

In the example

$\displaystyle\dot{x}_{1}$ $\displaystyle=1+2x+2x^{2}+y+2xy+2y^{2},$ (100)
$\displaystyle\dot{y}$ $\displaystyle=1+x+x^{2}+2y+4xy+y^{2}$ (101)

the Cauchy–Riemann conditions are not fulfilled, but the Calogero–Payandeh
conditions are, including the fact that none of the corresponding denominators
are different from zero,therefore this equation can be solved by their method,
and cannot be solved by ours. What is more, there are no "hidden" $i$ numbers
under the square roots:

$\displaystyle x(t)$
$\displaystyle=\frac{57(3e^{t}-2)+\sqrt{15}(2e^{t}-25)\sin(\sqrt{15}t)-3(62e^{t}-53)\cos(\sqrt{15}t)}{3(e^{t}-2)(7\sqrt{15}\sin(\sqrt{15}t)+71\cos(\sqrt{15}t)-76)},$
(102) $\displaystyle y(t)$
$\displaystyle=\frac{-57(e^{t}+2)+\sqrt{15}(23e^{t}-25)\sin(\sqrt{15}t)+3(9e^{t}+53)\cos(\sqrt{15}t)}{3(e^{t}-2)(7\sqrt{15}\sin(\sqrt{15}t)+71\cos(\sqrt{15}t)-76)}.$
(103)

The denominators in [3] are usually more or less complicated expressions
(polynomials) of the coefficients, except the denominator in Eq. (41) of [3]
that seems to be very hard to analyze.

In summary, there is a significant overlap in the scopes of the two methods
concerning two-variable second-degree systems. However, we have also provided
some ideas for cases outside this class of systems. If you truly need the
symbolic solution of an equation of the form (94)–(95) you should try
everything offered by the two papers. Finally, one should not overlook the
computational difficulties. For instance, in the case of (89), we were able to
arrive at an explicit solution (91)–(92), but we struggled in vain to
demonstrate its identity with the result obtained by the method of [3],
despite attempting all the simplifying functions available in the Wolfram
language.

We have seen that in the linear case, our method requires lengthier
calculations than the traditional approach. However, in the second-degree
case, our method is effective as opposed to the traditional one.

In [5] the authors can solve equations of specific form under restrictions on
the parameters. In one of their example, they find (explicitly) a periodic
solution, and to arrive at this they use an equation for a complex-valued
function and expand it into real and imaginary parts. [33] go in the opposite
direction: they use a five-variable real equation to prove the boundedness of
the complex Lorenz system.

Our method can also be used to approximate solutions of (kinetic or other)
differential equations. If our polynomial equation is close to a
Cauchy–Riemann–Erugin system, then one can solve that equation, and use the
Peano inequality to estimate the difference between the obtained solution and
the solution we are interested in. If we are in the position to make the
difference between the right-hand sides then we can make the difference of the
solutions arbitrarily small on any finite interval.

There is some hope that some special four-variable systems can be solved
similarly by applying quaternionic calculus, e.g. [8, Theorem 3.1] might be
useful for this purpose. See also [11, 28].

##  6 Appendix

###  6.1 Useful formulas

Before anything else, we present two formulas we are recurrently using.

######  Lemma 1.

  1. Suppose $\xi$ and $\eta\neq 0$ are real numbers. Then the square roots of $\xi+i\eta$ are as follows:

$\displaystyle\pm\left(\frac{\sqrt{\sqrt{\xi^{2}+\eta^{2}}+\xi}}{\sqrt{2}}+i\frac{\eta}{\sqrt{2}\sqrt{\sqrt{\xi^{2}+\eta^{2}}+\xi}}\right)$
(104) $\displaystyle=$
$\displaystyle\pm\left(\frac{\eta}{\sqrt{2}\sqrt{\sqrt{\xi^{2}+\eta^{2}}-\xi}}+i\frac{\sqrt{\sqrt{\xi^{2}+\eta^{2}}-\xi}}{\sqrt{2}}\right)$
(105)

###  6.2 Addendum to the second-degree case

Let us start with the solution of (18)

$z(t)=z_{1}+\frac{z_{2}-z_{1}}{1-\frac{z(0)-z_{2}}{z(0)-z_{1}}e^{(d-i\frac{e}{2})(z_{2}-z_{1})t}}=z_{2}+\frac{z_{1}-z_{2}}{1-\frac{z(0)-z_{1}}{z(0)-z_{2}}e^{(d-i\frac{e}{2})(z_{1}-z_{2})t}}$
(106)

where $z_{1}$ and $z_{2}$ are the roots of the polynomial
$z\mapsto(a+iA)+(b-ic)z+(d-i\frac{e}{2})z^{2},$ and $z(0):=x_{0}+iy_{0}.$ To
get an explicit solution, further—although elementary—lengthy calculations are
needed. Let us apply Euler’s formula.

$\displaystyle z(t)=$ (107) $\displaystyle
z_{1}+\frac{z_{2}-z_{1}}{1-\frac{z(0)-z_{2}}{z(0)-z_{1}}e^{(d-i\frac{e}{2})(z_{2}-z_{1})t}}=$
(108) $\displaystyle
x_{1}+iy_{1}+\frac{x_{2}-x_{1}+i(y_{2}-y_{1})}{1-\frac{x_{0}-x_{2}+i(y_{0}-y_{2})}{x_{0}-x_{1}+i(y_{0}-y_{1})}e^{(d(x_{2}-x_{1})+\frac{e}{2}(y_{2}-y_{1}))t}\left(\cos(\Omega
t)+i\sin(\Omega t)\right)}$ (109)

with $\Omega:=-\frac{e}{2}(x_{2}-x_{1})+d(y_{2}-y_{1}),$ and with
$x_{1}+iy_{1}:=z_{1},x_{2}+iy_{2}:=z_{2}.$

######  Proposition 1.

Assume that the real numbers $a,b,c,d,e,A$ are such that $(d-i\frac{e}{2})\neq
0.$ Then, the roots of the polynomial

$z\mapsto a+iA+(b-ic)z+(d-i\frac{e}{2})z^{2}$ (110)

are as follows

$\displaystyle z_{1,2}$
$\displaystyle=\frac{-b+ic\pm(q+ir)}{2d-ie}=\frac{-b\pm q+i(c\pm r)}{2d-ie}$
(111) $\displaystyle=\left(\frac{\pm 2dq-2bd-ec\mp
er}{4d^{2}+e^{2}}+i\frac{2cd\pm 2dr\pm eq-be}{4d^{2}+e^{2}}\right)$ (112)

with

$\displaystyle q$
$\displaystyle:=\sqrt{\frac{1}{2}(b^{2}-c^{2}-4ad+2eA+\sqrt{(b^{2}-c^{2}-4ad+2eA)^{2}+4(-\frac{ae}{2}-2bC+dA)^{2}})},$
(113) $\displaystyle r$
$\displaystyle:=\frac{-\sqrt{2}(-\frac{ae}{2}-2bC+dA)}{\sqrt{b^{2}-c^{2}-4ad+2eA+\sqrt{(b^{2}-c^{2}-4ad+2eA)^{2}+4(-\frac{ae}{2}-2bC+dA)^{2}})}}.$
(114)

Proof.

  1. 1.

To go from (111) to (112) apply the rule for ratios.

  2. 2.

Next, let us apply (104) with
$\xi:={b^{2}-c^{2}-4ad+2eA},\quad\eta:={-2(-\frac{ae}{2}-2bC+dA)}$ to get
(112).

Using the explicit formula for the roots one can further expand the solution
to (18), but if one makes all the necessary substitutions with $q$ and $r$ as
defined above, the results will be useless. However, we shall show examples
that the procedure can be finished in all special cases.

###  6.3 Addendum to an example of the second-degree case

######  Proposition 2.

The denominator in the equations (90) is never zero.

Proof. Let $y:=e^{2\sqrt{2}t}$, then we have the following expression in the
denominator.

$3(2+\sqrt{2})y^{2}-8y+(6-3\sqrt{2})$ (115)

The zeros of (115) above are:

$y_{1}=\frac{4-i\sqrt{2}}{6+3\sqrt{2}},\quad
y_{2}=\frac{4+i\sqrt{2}}{6+3\sqrt{2}}$

But $y:=e^{2\sqrt{2}t}$ is always a real number for all $t\in\mathbf{R}$ and
cannot be equal to the above complex roots. Therefore, $x(t)$ and $y(t)$ are
always defined for all real numbers $t\in\mathbf{R}$.

###  6.4 Addendum to the third-degree case

FHJ graphs

Let us apply the proposed procedure in some typical cases.

  1. 1.

The real initial value problem

$\displaystyle\dot{x}=-8+12x-6x^{2}+6y^{2}+x^{3}-3xy^{2},\quad\dot{y}=12y-12xy+3x^{2}y-y^{3},$
(116) $\displaystyle x(0)=2,\quad y(0)=1$

leads to the complex differential equation

$\dot{z}=(z-2)^{3},\quad z(0)=2+i$

having the solution

$z(t)=2+\frac{i}{\sqrt{1+2t}},\quad(t\in]-\frac{1}{2},+\infty[),$ (117)

or

$x(t)=2,\quad y(t)=\frac{1}{\sqrt{1+2t}},\quad(t\in]-\frac{1}{2},+\infty[).$
(118)

Substituting this result into (116) verifies the calculations.

  2. 2.

One can proceed similarly if the starting initial value problem is

$\displaystyle\dot{x}=-3x+6xy+x^{3}-3xy^{2},\quad\dot{y}=1-3y-3x^{2}+3y^{2}+3x^{2}y-y^{3},$
(119) $\displaystyle x(0)=1,\quad y(0)=2.$

The complex differential equation is now:

$\dot{z}=(z-i)^{3},\quad z(0)=1+2i$

having the solution

$\displaystyle z(t)$
$\displaystyle=i+\frac{2}{\sqrt{1-8t}},\quad(t\in]-\infty,\frac{1}{8}[),$

with

$x(t)=\frac{2}{\sqrt{1-8t}},\quad(t\in]-\infty,\frac{1}{8}[),y(t)=1,\quad(t\in\mathbb{R})$
(120)

As can be seen from the above examples, if the right-hand side of the complex
equation has a triple root, the solutions of the real system can be expressed
symbolically. One can never have this nice form with kinetic differential
equations.

######  Proposition 3.

The right-hand side $(z-a-ib)^{3};\quad a,b\in\mathbb{R};a^{2}+b^{2}\neq 0$
always comes from a real system containing negative cross-effect(s).

Proof. Since

$\displaystyle(z-a-ib)^{3}=$
$\displaystyle-a^{3}+3ab^{2}+3(a^{2}-b^{2})x-6aby-3ax^{2}+6bxy+3ay^{2}+x^{3}-3xy^{2}$
$\displaystyle+$ $\displaystyle
i(-3a^{2}b+b^{3}+6abx+3(a^{2}-b^{2})y-3bx^{2}-6axy+3by^{2}+3x^{2}y-y^{3})$

$6ab\leq 0$ and $6ab\geq 0$ implies that $6ab=0.$ If $a=0,$ then $b^{3}\geq 0$
and $-3b\geq 0$ should hold, thus $b=0$ would be implied. If $b=0,$ then
$-a^{3}\geq 0$ and $3a\geq 0$ should hold, thus $a=0$ would be implied.

######  Remark 2.

In other cases, one may either have terrible formulas based on the Cardano
formula, or one may be content with an implicit solution for the real
components. E.g. the system

$\displaystyle\dot{x}$ $\displaystyle=1-3y-2x^{2}+4xy+2y^{2}+x^{3}-3xy^{2},$
(121) $\displaystyle\dot{y}$
$\displaystyle=-1+3x-2x^{2}-4xy+2y^{2}+3x^{2}y-y^{3}$ (122)

leads to the complex differential equation:

$\dot{z}=(z-1)(z-i)(z-i-1),\quad z(0)=x_{0}+iy_{0},$ (123)

the solution of which can either be obtained using the Cardano formula, but
the result is too complicated to be useful, or it can be left in the form to
produce an implicit solution, first for $z(t),$ then for $x(t),y(t).$

###  6.5 An example of the four-variable case

Let us start with a system of equations with a homogeneous right-hand side.

$\displaystyle\dot{x}_{1}$
$\displaystyle=ax_{1}^{2}+bx_{2}^{2}+cy_{1}^{2}+dy_{2}^{2}+ex_{1}x_{2}+fy_{1}y_{2}+gx_{1}y_{1}+hx_{1}y_{2}+ix_{2}y_{1}+jx_{2}y_{2},$
(124) $\displaystyle\dot{y}_{1}$
$\displaystyle=kx_{1}^{2}+lx_{2}^{2}+my_{1}^{2}+ny_{2}^{2}+px_{1}x_{2}+qy_{1}y_{2}+rx_{1}y_{1}+sx_{1}y_{2}+tx_{2}y_{1}+ux_{2}y_{2},$
(125) $\displaystyle\dot{x}_{2}$
$\displaystyle=Ax_{1}^{2}+Bx_{2}^{2}+Cy_{1}^{2}+Dy_{2}^{2}+Ex_{1}x_{2}+Fy_{1}y_{2}+Gx_{1}y_{1}+Hx_{1}y_{2}+Ix_{2}y_{1}+Jx_{2}y_{2},$
(126) $\displaystyle\dot{y}_{2}$
$\displaystyle=Kx_{1}^{2}+Lx_{2}^{2}+My_{1}^{2}+Ny_{2}^{2}+Px_{1}x_{2}+Qy_{1}y_{2}+Rx_{1}y_{1}+Sx_{1}y_{2}+Tx_{2}y_{1}+Ux_{2}y_{2}.$
(127)

and the Cauchy–Riemann equations imply

$\displaystyle
a=-c=\frac{r}{2},e=t=s=-f,h=q=i=-p,\frac{g}{2}=m=-k,b=-d=\frac{u}{2},n=-l=\frac{j}{2},$
(128) $\displaystyle
A=-C=\frac{R}{2},E=T=S=-F,H=Q=I=-P,\frac{G}{2}=M=-K,B=-D=\frac{U}{2},N=-L=\frac{J}{2}.$
(129)

therefore we have simpler complex equations as follows:

$\displaystyle\dot{z}_{1}$
$\displaystyle=(a-i\frac{g}{2})z_{1}^{2}+(b-i\frac{j}{2})z_{2}^{2}+(e-ih)z_{1}z_{2}$
(130) $\displaystyle\dot{z}_{2}$
$\displaystyle=(A-i\frac{G}{2})z_{1}^{2}+(B-i\frac{J}{2})z_{2}^{2}+(E-iH)z_{1}z_{2}$
(131)

Now we can again use the method described in 4.

## References

  * [1] J. M. Borwein and R. E. Crandall.  Closed forms: What they are and why we care.  Notices of the AMS, 60(1):50–65, 2013.
  * [2] F. Calogero, R. Conte, and F. Leyvraz.  New algebraically solvable systems of two autonomous first-order ordinary differential equations with purely quadratic right-hand sides.  Journal of Mathematical Physics, 61(10):102704, 2020.
  * [3] F. Calogero and F. Payandeh.  Solution of the system of two coupled first-order ODEs with second-degree polynomial right-hand sides.  Mathematical Physics, Analysis and Geometry, 24(3):29, 2021.
  * [4] F. Calogero and F. Payandeh.  Solvable systems of two coupled first-order ODEs with homogeneous cubic polynomial right-hand sides.  Journal of Mathematical Physics, 62(1), 2021.
  * [5] F. Calogero and F. Payandeh.  Explicitly solvable systems of first-order ordinary differential equations with homogeneous right-hand sides, and their periodic variants.  Theoretical and Mathematical Physics, 213(1):1317–1330, 2022.
  * [6] Gh. Craciun, J. Jin, and P. Y. Yu.  An algorithm for finding weakly reversible deficiency zero realizations of polynomial dynamical systems.  SIAM Journal on Applied Mathematics, 83(4):1717–1737, 2023.
  * [7] Gh. Craciun, M. D. Johnston, G. Szederkényi, E. Tonello, J. Tóth, and P. Y. Yu.  Realizations of kinetic differential equations.  Mathematical Biosciences and Engineering, 17:862–892, 2019.
  * [8] C. A. Deavours.  The quaternion calculus.  The American Mathematical Monthly, 80:995–1008, 1973.
  * [9] A. Deshpande.  Source-only realizations, weakly reversible deficiency one networks, and dynamical equivalence.  SIAM Journal on Applied Dynamical Systems, 22(2):1502–1521, 2023\.
  * [10] M. Feinberg.  Foundations of Chemical Reaction Network Theory.  Springer, 2019.
  * [11] G. Gentili, C. Stoppato, and D. C. Struppa.  Regular functions of a quaternionic variable.  Springer Nature, 2022.
  * [12] J. Giné and V. G. Romanovski.  Integrability conditions for Lotka–Volterra planar complex quintic systems.  Nonlinear Analysis: Real World Applications, 11(3):2100–2105, 2010\.
  * [13] J. Giné and C. Valls.  Analytic integrability of cubic–linear planar polynomial differential systems.  Journal of Differential Equations, 260(2):1690–1716, 2016.
  * [14] V. Hárs and J. Tóth.  On the inverse problem of reaction kinetics.  Qualitative theory of differential equations, 30:363–379, 1981\.
  * [15] B. S. Hernandez and P. V. N. Lubenia.  Deriving positive steady states of chemical reaction networks for classes of non-mass action kinetics.  arXiv preprint arXiv:2304.11594, 2023.
  * [16] M. D. Johnston, D. Siegel, and G. Szederkényi.  Computing weakly reversible linearly conjugate chemical reaction networks with minimal deficiency.  Mathematical Biosciences, 241(1):88–98, 2013.
  * [17] E. Kamke.  Differentialgleichungen, Volume II: Lösungsmethoden und Lösungen.  American Mathematical Soc., 1974.
  * [18] E. Kamke.  Differentialgleichungen: Lösungsmethoden und Lösungen.  Springer-Verlag, 2013.
  * [19] G. Lente.  Analytical solutions for the rate equations of irreversible two-step consecutive processes with second order later steps.  J. Math. Chem., 53:1759–1771, 2015.
  * [20] G. Lente.  Deterministic Kinetics in Chemistry and Systems Biology: The Dynamics of Complex Reaction Networks.  Springer, 2015.
  * [21] G. Lente.  Kinetics of irreversible consecutive processes with first order second steps: analytical solutions.  J. Math. Chem., 53:1172–1183, 2015.
  * [22] G. Lente.  Analytical solutions for the rate equations of irreversible two-step consecutive processes with mixed second order later steps.  J. Math. Chem., 55(3):832–848, 2017.
  * [23] I. Nagy and J. Tóth.  Quadratic first integrals of kinetic differential equations.  J. Math. Chem., 52:93–114, 2014.  https://doi.org/10.1007/s10910-013-0247-3.
  * [24] Tae Jun Park.  Kinetics of reversible consecutive reactions.  Bulletin of the Korean Chemical Society, 34(1):243–245, 2013.
  * [25] A. D. Polyanin.  Exact solutions of mathematical equations, 2004–2017.  https://eqworld.ipmnet.ru/en/solutions.htm.
  * [26] N. M. Rodiguin and E. N. Rodiguina.  Consecutive Chemical Reactions: Mathematical Analysis and Development (English edn.).  D. Van Nostrand Co., New York, 1964.
  * [27] L. I. Ronkin.  Elements of the Theory of Multivariate Analytic Functions.  Naukova Dumka, Kyiv, 1977.
  * [28] A. Sudbery.  Quaternionic analysis.  Mathematical Proceedings of the Cambridge Philosophical Society, 99:169–186, 1986.
  * [29] Z. G. Szabó.  Kinetic characterization of complex reaction systems.  In Ch. Bamford and C. F. H. Tipper, editors, Comprehensive Chemical Kinetics, volume 2, chapter 1, pages 1–80. Elsevier, Amsterdam, Oxford, New York, 1969.
  * [30] J. Tóth and P. Érdi.  Kinetic symmetries: Some hints.  In M. Moreau and P. Turq, editors, Chemical Reactivity in Liquids: Fundamental Aspects, pages 517–522. Springer, 1988.
  * [31] J. Tóth, A. L. Nagy, and D. Papp.  Reaction Kinetics: Exercises, Programs and Theorems. Mathematica for Deterministic and Stochastic Kinetics.  Springer Nature, Berlin, Heidelberg, New York, 2018.
  * [32] W. Wirtinger.  Zur formalen Theorie der Funktionen von mehreren komplexen Veränderlichen.  Mathematische Annalen, 97:357–375, 1927.
  * [33] F. Zhang and G. Zhang.  Boundedness solutions of the complex Lorenz chaotic system.  Applied Mathematics and Computation, 243:12–23, 2014.
  * [34] Ю. С. Богданов and Ю. Б. Сыроид.  Дифференциальные уравнения.  Вышэйшая Школа, Минск, 1983.
  * [35] Н. М. Матвеев.  Сборник задач и упражнений по обыкновенным дифференциальным уравнениям.  Вышэйшая Школа, Минск, 1977.