# Computation of Normal Forms for Systems
with Many Parameters

Tatjana Petek1,4 and Valery G. Romanovski1,2,3111Corresponding author, email:
<EMAIL_ADDRESS>
1Faculty of Electrical Engineering and Computer Science,
University of Maribor, Koroška cesta 46, SI-2000 Maribor, Slovenia
2Center for Applied Mathematics and Theoretical Physics,
Mladinska 3, SI-2000 Maribor, Slovenia
3Faculty of Natural Science and Mathematics,
University of Maribor, Koroška cesta 160, SI-2000 Maribor, Slovenia
4Institute of Mathematics, Physics and Mechanics,
Jadranska 19, SI-1000 Ljubljana, Slovenia

###### Abstract

There are two ways to compute Poincaré-Dulac normal forms of systems of ODEs.
Under the original approach used by Poincaré the normalizing transformation is
explicitly computed. On each step, the normalizing procedure requires the
substitution of a polynomial to a series. Under the other approach, a normal
form is computed using Lie transformations. In this case, the changes of
coordinates are performed as actions of certain infinitesimal generators. In
both cases, on each step the homological equation is solved in the vector
space of polynomial vector fields $\mathcal{V}^{n}_{j}$ where each component
of the vector field is a homogeneous polynomial of degree $j$. We present the
third way of computing normal forms of polynomial systems of ODEs where the
coefficients of all terms are parameters. Although we use Lie transforms the
homological equation is solved not in $\mathcal{V}^{n}_{j}$ but in the vector
space of polynomial vector fields where each component is a homogeneous
polynomial in the parameters of the system. It is shown that the space of the
parameters is a kind of dual space and, the computation of normal forms can be
performed in the space of parameters treated as the space of generalized
vector fields. The approach provides a simple way to parallelize the normal
form computations opening the way to compute normal forms up to higher order
than under previously known two approaches.

Mathematical Subject Classification 2010: 34C20, 37G05

Keywords: Poincaré-Dulac normal form, Lie transform, Lie algebra of vector
fields

##  1 Introduction

The normal form theory of ordinary differential equations is a powerful tool
for studying the qualitative behavior of nonlinear systems near equilibrium
points and local integrability. The central idea of the theory, which
originates from the works of Poincaré, Dulac and Lyapunov [13, 7, 9], is to
perform a change of coordinates or a succession of coordinate transformations,
to place the original system into a form most amenable to study.

There are many books devoted to various aspects of the normal form theory (see
e.g. [2, 4, 6, 8, 10, 11]). Some discussions on the recent progress can be
found in the survey papers [1, 18].

Three main questions in the theory of normal form are to analyze the structure
of normal form, compute it and analyze the convergence of the obtained series.
Our paper concerns only algorithms for computing of normal forms regardless of
their structure and convergence. Usually, computations of normal forms require
huge amount of algebraic operations, so the algorithms used in normal form
theory are complex and computationally expensive. The complexity of the
algorithms can make it challenging to compute normal forms for systems with a
large number of variables. The complexity of computations fast increases for
equations dependent on parameters. For these reasons, normal form computations
rarely can be performed by hand but are carried out using computer algebra
systems.

There are two main approaches to computations of normal form. The classical
one, used by Poincaré, Dulac and Lyapunov requires substitution series into
series. It appears that the most computationally efficient algorithms obtained
using this approach are given in [21, 19]. The other one is based on Lie
transformations and was developed starting from the works of Birkhoff [3],
Steinberg [17], Chen [5] and others. Using this way a normalization is
performed as a sequence of linear transformations.

Throughout this paper, fix $\mathbb{F}$ to be either the field of all real
numbers or the field of all complex numbers. Let $\mathcal{V}^{n}$ be the
space of vector fields $v:\mathbb{F}^{n}\to\mathbb{F}^{n}$ which coordinates
are power series in $x_{1},\dots,x_{n}$ (for brevity we will say ”in $x$”)
vanishing at the origin, which can be convergent or merely formal, and
$\mathcal{V}^{n}_{j}$ be the vector space of polynomial vector fields on
$\mathbb{F}^{n}$, that is, vector fields $v:\mathbb{F}^{n}\to\mathbb{F}^{n}$
such that each component $v_{i}$, $i=1,\dots,n$, of $v$ is a homogeneous
polynomial of degree $j$. Note that $\mathcal{V}^{n}$ and
$\mathcal{V}^{n}_{j}$ are vector spaces over $\mathbb{F}$ and
$\mathcal{V}^{n}=\oplus_{j=0}^{\infty}\mathcal{V}^{n}_{j}$.

Consider the differential equation

$\dot{x}=u(a,x),$ (1)

where $a=(a_{1},\dots,a_{\ell})$ is an $\ell$-tuple of parameters (regarded as
element of $\mathbb{F}^{\ell}$), $x=(x_{1},\dots,x_{n})$,
$u(a,x)\in\mathcal{V}^{n}$ and all terms of $u(a,x)$ depend polynomially on
parameters $a_{1},\dots,a_{\ell}$ (on $a$ for brevity).

The space $\mathcal{V}^{n}$ can also be viewed as a module over the ring of
power series in $x_{1},\dots,x_{n}$. It has natural grading by the degree of
polynomials, that is, equation (1) can be written as

$\dot{x}=\sum_{j=1}^{\infty}u_{j}(a,x),$ (2)

where $u_{j}(a,x)$ is a homogeneous polynomial of degree $j$ in $x$, that is,
$u_{j}(a,x)\in\mathcal{V}^{n}_{j}$. As a rule, the normalization of system (1)
is performed according to this grading – polynomials $u_{j}$ are changed into
a more amendable form for studying for $j=1,2,3,\dots$, consequently.

In this paper, we use another approach which is based on the grading of the
power series on the right hand side of (1) by the degree of polynomials viewed
as polynomials in the parameters of the system. Indeed, since by our
assumption each term of $u(a,x)$ depends polynomially on parameters $a$, each
component of the vector field $u(a,x)$ is a power series in
$a_{1},\dots,a_{\ell},x_{1},\dots,x_{n}$, so we can write (1) in the form

$\dot{x}=\sum_{s=1}^{\infty}\bar{u}_{s}(a,x),$ (3)

where $\bar{u}_{s}$ is a homogeneous polynomial of degree $s$ in the
parameters $a_{1},\dots,a_{\ell}$ of the system.

We will work with formal vector fields using the grading of the formal series
by polynomials which are homogeneous with respect to the parameters $a$ of
$u(a,x)$ (generally, they are _not homogeneous_ with respect to the variable
$x$). This idea looks unusual at first glance, but it makes perfect sense
since when a vector field depends on many parameters, then it becomes more
important to have a tool to handle the monomials with respect to the
parameters rather than with respect to the unknown functions $x_{i}$.

However, the module of formal vector fields with components from
$\mathbb{F}[[a,x]]$ is too large. We define a submodule which is graded using
the structure of the monomials on the right-hand side of (3). As a result, we
will obtain an algorithm for computing normal forms which has two advantages
with respect to the traditional ones:

  * –

It allows performing the parallelization of normal forms computations. Namely,
for equation (1), the terms of powers series in a normal form have the form
$p(a)x^{r}$, where $x^{r}$ is a resonant monomial and $p(a)=\sum_{\alpha\in
Supp(p)}p_{\alpha}a^{\alpha}$ is a polynomial in $a$. Using traditional
methods, one can compute only the whole term $p(a)x^{r}$. Using our technique,
it is possible to compute any term $p_{\alpha}a^{\alpha}$ of $p(a)$ without
computing the whole polynomial $p(a)$. Since the computation of normal form
using computer algebra systems is extremely time and memory consuming, this
parallelization should allow computing normal forms up to higher order than
using the traditional approaches.

  * –

The obtained algorithm is more convenient for computations by hand than the
ones known before because it does not require the computation of derivatives
and collection of similar terms, since only arithmetic operations with numbers
are used.

As we will see, the space of parameters of the differential equations is a
kind of dual space and all computations can be performed in this dual space.

The paper is organized as follows. In the next section, we present some basics
of the normal forms. In Section 3 a module which we use for our normal form
computations is defined and our normalization algorithm is presented. In
Section 4 we introduce the notion of generalized vector fields, reformulate
the algorithm in terms of such vector fields and give an example of the normal
form computation.

##  2 Preliminaries

Consider the autonomous differential equation (we can also say $n$-dimensional
autonomous system)

$\dot{x}=f(x),\hskip 8.53581ptx\in\mathbb{F}^{n},$ (4)

where

$f(x)=Ax+\sum_{j=2}^{\infty}f_{j}(x),$

$f_{j}(x)$ is an $n$-dimensional vector valued homogeneous polynomial of
degree $j$, that is, $f_{j}(x)\in\mathcal{V}^{n}_{j}$ and
$f(x)\in\mathcal{V}^{n}$. In this paper, we assume that $A$ is a _diagonal
matrix_ with the eigenvalues $\lambda_{1},\dots,\lambda_{n}$,
$x=(x_{1},\dots,x_{n})^{T}.$

Let

$\lambda=(\lambda_{1},\dots,\lambda_{n})^{T}\in\mathbb{F}^{n}$

be the column vector of eigenvalues of $A$. Set
$\mathbb{N}_{0}=\mathbb{N}\cup\\{0\\}$. For
$\alpha=(\alpha_{1},\dots,\alpha_{n})\in\mathbb{N}_{0}^{n}$ we denote
$\langle\lambda,\alpha\rangle=\sum_{i=1}^{n}\alpha_{i}\lambda_{i}$,
$|\alpha|=\alpha_{1}+\dots+\alpha_{n}$ and
$x^{\alpha}=x_{1}^{\alpha_{1}}x_{2}^{\alpha_{2}}\cdots x_{n}^{\alpha_{n}}$.
Throughout this paper $x$ is a column vector and $x^{\alpha}$ is a monomial.

It is well known that any formal invertible change of coordinates of the form

$x=y+\tilde{h}(y)=y+\sum_{j=2}^{\infty}\tilde{h}_{j}(y),$ (5)

with $\tilde{h}_{j}\in\mathcal{V}^{n}_{j}$, $j=2,3,\dots,$ brings (4) to a
system of a similar form,

$\dot{y}=Ay+\tilde{g}(y),$ (6)

where

$\tilde{g}(y)=\sum_{j=2}^{\infty}\tilde{g}_{j}(y),\qquad\tilde{g}_{j}(y)\in\mathcal{V}^{n}_{j},\qquad
j=2,3,\dots.$

Let $e_{k}$ be the $n$–dimensional column vector with its $k$-th component
equal to $1$ and the others all equal to zero. It is said that the monomial
$x^{\alpha}$ ($|\alpha|>1$) in $e_{k}^{T}f(x)$, $k=1,\dots,n$ is a _resonant_
monomial if $k$ and $\alpha$ satisfy

$\langle\lambda,\alpha\rangle-\lambda_{k}=0.$ (7)

A term of $f(x)$ is called _resonant_ if it is a resonant monomial multiplied
by a scalar from $\mathbb{F}$ (or a parameter) and _nonresonant_ otherwise.
Similarly, the monomial $y^{\alpha}$ in $e_{k}^{T}\tilde{g}(y)$ or in
$e_{k}^{T}\tilde{h}(y)$ is resonant if (7) holds.

######  Definition 1.

We say that system (6) is in the _Poincaré–Dulac normal form_ (or, simply in
the normal form) if $\tilde{g}(y)$ contains only resonant terms.

By Poincaré–Dulac theorem any system (4) can be brought to a normal form by a
substitution (5) [13, 7]. Transformation (5) which brings system (4) to a
Poincaré–Dulac normal form is called a normalization or a normalizing
transformation. In general, the normalization of system (6) is not unique.
However, the normalization containing only nonresonant terms is unique.

For $v(x)\in\mathcal{V}^{n}$ we denote by $Dv(x)$ the $n\times n$ matrix of
partial derivatives of $v(x)$. The linear operator
$\pounds:\mathcal{V}^{n}\to\mathcal{V}^{n}$ defined by

$(\pounds v)(x)=Dv(x)Ax-Av(x)$ (8)

is called the _homological operator_. It is easy to see that
$\pounds(\mathcal{V}^{n}_{j})\subseteq\mathcal{V}^{n}_{j}$, so let

$\pounds_{j}:\mathcal{V}^{n}_{j}\to\mathcal{V}^{n}_{j}$

be the restriction of $\pounds$ to $\mathcal{V}^{n}_{j}$. It is known (see
e.g. [2, 10]) that the $\pounds_{j}$ is semisimple with the eigenvalues
$\beta_{\alpha,i}=\langle\alpha,\lambda\rangle-\lambda_{i}$ and (basis)
eigenvectors $p_{i,\alpha}=e_{i}x^{\alpha}$, where $|\alpha|=j$ and
$i=1,\dots,n$. The kernel of $\pounds_{j}$ is then spanned by all monomials
$e_{i}x^{\alpha}$, where $i$ and $\alpha$ satisfy $\beta_{\alpha,i}=0$. Then
we can decompose $\mathcal{V}^{n}_{j}$ as

$\mathcal{V}^{n}_{j}={\rm im}\pounds_{j}\oplus\ker\pounds_{j}.$

According to this decomposition, we can write every
$v_{j}\in\mathcal{V}^{n}_{j}$ uniquely as

$v_{j}=\pounds_{j}\tilde{h}_{j}+\tilde{g}_{j}$ (9)

for some $\tilde{h}_{j}\in\mathcal{V}^{n}_{j}$, i.e.
$\pounds\tilde{h}_{j}\in{\rm im}\pounds_{j}$ and
$\tilde{g}_{j}\in\ker\pounds_{j}$ is the (natural) projection of $v_{j}$ onto
$\ker\pounds_{j}$. In other words, in $\tilde{g}_{j}$ can only appear terms
$e_{i}x^{\alpha}$ with $\langle\alpha,\lambda\rangle-\lambda_{i}=0$, thus
being resonant terms.

Assume that equation (4) has already been normalized to order $j-1$, so (6) is
a normal form of (4) up to order $j-1$. Then to perform the normalization to
order $j$, one solves the equation

$\pounds_{j}\tilde{h}_{j}=f_{j}-\tilde{g}_{j},$ (10)

known as the homological equation, for $\tilde{h}_{j}$ and $\tilde{g}_{j}$
(the solution is unique, providing we choose all resonant terms of
$\tilde{h}_{j}$ to be zero). With this $\tilde{g}_{j}$ (6) is in the normal
form to order $j$, and to perform the next step of the normalization, one
needs to do in (4) the transformation

$x\to x+\tilde{h}_{j}(x)$

and re-expand the series $f(x)$ on the right hand of (4) up to order $j+1$.

The computational procedure described above is simple, however, it requires
substitutions of polynomials into series and, therefore, it is very laborious
from the computational point of view.

Another setting for computing the normal form is using the Lie brackets. For
two vector fields $v(x)$ and $w(x)$ in $\mathcal{V}^{n}$

$[w,v]=D(v)w-D(w)v$

is the Lie bracket. For a given $w(x)\in\mathcal{V}^{n}$, the operator

$\text{ad}\,w(x)=[w(x),\cdot]\,:\,{\mathcal{V}^{n}}\to{\mathcal{V}^{n}}$

is the so-called adjoint map. In particular, $\pounds=\text{ad\,}Ax$ and
equation (6) is in the normal form if and only if

$(\text{ad\,}Ax)\,\tilde{g}_{k}=0,{\ \rm for\ all\ }k=2,3,\dots.$

Let

$\frac{dx}{ds}=b(x),\qquad b(x)\in{\mathcal{V}^{n}},\qquad b(0)=0,$ (11)

and let $\psi_{s}(x)$ be the flow of (11) satisfying $\psi_{0}(x)=x$. By the
Theorem on Lie series for vector fields (see [10, Theorem 4.4.8] or [12, p.
202]) after the transformation

$x\to\psi_{s}(x)$ (12)

the right-hand side of equation (4) is changed to

$\psi_{s}^{\prime}(x)^{-1}f(\psi_{s}(x))=f(x)+s(\text{ad\,}b(x))f(x)+\frac{s^{2}}{2}(\text{ad\,}b(x))^{2}f(x)+\dots.$
(13)

Thus, if equation (4) is in the normal form form up to order $j-1$, to perform
the next step of the normalization, one solves the homological equation (10)
for $\tilde{h}_{j}$ and $\tilde{g}_{j}$ and performs in (4) the transformation
(12) with $s=1$ and $\psi_{s}(x)$ being the flow of

$\frac{dx}{ds}=\tilde{h}_{j}(x).$

Denoting $\psi_{1}(x)$ by $\exp(\tilde{h}_{j})$ we can write the normalization
to order $m$ as

$x=H_{m}(y)\circ\cdots\circ H_{2}(y),$

where $H_{k}=\exp({\tilde{h}_{k}})$, $k=2,\dots,m$.

An algorithmic description of this normalizing procedure can be found, for
instance, in [16, 20].

##  3 A new algorithm for computing of normal forms

###  3.1 A grading of the formal power series module

Assume that the terms of the function $f(x)$ in (4) depends polynomially on
parameters $a_{1},\dots,a_{\ell}$ and let $a=(a_{1},\dots,a_{\ell})$, that is,
the terms of $f(x)=f(a,x)$ are of the form

$f^{(i)}_{(\mu_{1},\dots,\mu_{\ell},\beta_{1},\dots,\beta_{n})}a_{1}^{\mu_{1}}\cdots
a_{\ell}^{\mu_{\ell}}x_{1}^{\beta_{1}}\cdots
x_{n}^{\beta_{n}}\qquad(i=1,\dots,n).$ (14)

Then the coefficients of the resonant monomials in the normal form are
polynomials in parameters $a_{1},\dots a_{\ell}$.

As an example, we consider the following quite arbitrarily chosen differential
equation

$\displaystyle\dot{x}_{1}=$
$\displaystyle\phantom{-}x_{1}+a^{(1)}_{10}x_{1}^{2}+a^{(1)}_{01}x_{1}x_{2}+a^{(1)}_{-13}x_{2}^{3}=x_{1}(1+a^{(1)}_{10}x_{1}+a^{(1)}_{01}x_{2}+a^{(1)}_{-13}x_{1}^{-1}x_{2}^{3}),$
(15) $\displaystyle\dot{x}_{2}=$ $\displaystyle-
x_{2}+a^{(2)}_{10}x_{1}x_{2}+a^{(2)}_{01}x_{2}^{2}+a^{(2)}_{02}x_{2}^{3}=x_{2}(-1+a^{(2)}_{10}x_{1}+a^{(2)}_{01}x_{2}+a^{(2)}_{02}x_{2}^{2}).$

Let us remark that there is a deep reason for the given notation of parameters
which will be very important in the sequel.

Computing the normal form of (15) up to order 5 we obtain

$\dot{x}={\rm diag}(1,-1)x+g_{2}(x)+g_{3}(x)+g_{4}(x)+g_{5}(x),$ (16)

where $x=(x_{1},x_{2})^{T}$, $g_{2}(x)=g_{4}(x)=0$ and

$\displaystyle g_{3}(x)$
$\displaystyle=((-a^{(1)}_{10}a^{(1)}_{01}+a^{(1)}_{01}a^{(2)}_{10})x_{1}^{2}x_{2},(-a^{(1)}_{01}a^{(2)}_{10}+a^{(2)}_{01}a^{(2)}_{10})x_{1}x_{2}^{2})^{T}$
$\displaystyle=(x_{1}(-a^{(1)}_{01}a^{(1)}_{10}+a^{(1)}_{01}a^{(2)}_{10})x_{1}x_{2},\
x_{2}(-a^{(1)}_{01}a^{(2)}_{10}+a^{(2)}_{01}a^{(2)}_{10})x_{1}x_{2})^{T},$
$\displaystyle g_{5}(x)$
$\displaystyle=(((a^{(1)}_{01})^{2}a^{(1)}_{10}a^{(2)}_{10}+a^{(1)}_{01}a^{(1)}_{10}a^{(2)}_{01}a^{(2)}_{10}-2(a^{(1)}_{01})^{2}(a^{(2)}_{10})^{2})x_{1}^{3}x_{2}^{2},$
$\displaystyle\phantom{-}(-a^{(1)}_{01}a^{(1)}_{10}a^{(2)}_{01}a^{(2)}_{10}+a^{(1)}_{10}a^{(2)}_{02}a^{(2)}_{10}+2(a^{(1)}_{01})^{2}(a^{(2)}_{10})^{2}-a^{(1)}_{01}a^{(2)}_{01}(a^{(2)}_{10})^{2}+2a^{(2)}_{02}(a^{(2)}_{10})^{2})x_{1}^{2}x_{2}^{3})^{T}$
$\displaystyle=(x_{1}((a^{(1)}_{01})^{2}a^{(1)}_{10}a^{(2)}_{10}+a^{(1)}_{01}a^{(1)}_{10}a^{(2)}_{01}a^{(2)}_{10}-2(a^{(1)}_{01})^{2}(a^{(2)}_{10})^{2})x_{1}^{2}x_{2}^{2},$
$\displaystyle\phantom{-}x_{2}(-a^{(1)}_{01}a^{(1)}_{10}a^{(2)}_{01}a^{(2)}_{10}+a^{(1)}_{10}a^{(2)}_{02}a^{(2)}_{10}+2(a^{(1)}_{01})^{2}(a^{(2)}_{10})^{2}-a^{(1)}_{01}a^{(2)}_{01}(a^{(2)}_{10})^{2}+2a^{(2)}_{02}(a^{(2)}_{10})^{2})x_{1}^{2}x_{2}^{2})^{T}.$

We see from the example that the expression for the normal form is ”dominated”
by parameters in the sense that there are 12 monomials involved parameters,
whereas for the phase variables we have just four monomials $x_{1}^{2}x_{2}$,
$x_{1}x_{2}^{2}$, $x_{1}^{3}x_{2}^{2}$, $x_{1}^{2}x_{2}^{3}$.

Consider polynomial equations (4) written in the form

$\dot{x}_{k}=\lambda_{k}x_{k}+x_{k}\sum_{{\bar{\imath}}\in\Omega_{k}}a^{(k)}_{\bar{\imath}}x^{\bar{\imath}}\qquad(k=1,\dots,n),$
(17)

where $\Omega_{k}$, $k=1,\dots,n$, is a fixed ordered set of multi-indices
($n$-tuples) $\bar{\imath}=(i_{1},\dots,i_{n})$, whose $k$-th entry is from
$\mathbb{N}_{-1}=\\{-1\\}\cup\mathbb{N}_{0}$ and all other entries are from
$\mathbb{N}_{0}$. The symbols $a_{\bar{\imath}}^{(k)}$ are parameters of the
family (17). Let $\Omega=\Omega_{1}\cup\Omega_{2}\cup\dots\cup\Omega_{n}$. If
$|\Omega_{k}|$ is the number of elements of $\Omega_{k}$, then
$\ell=|\Omega_{1}|+\dots+|\Omega_{n}|$ is the number of parameters
$a^{(k)}_{\bar{\imath}}$ in (17). We will write (17) in a short form as the
equation

$\dot{x}=Ax+F(a,x),$ (18)

where

$F(a,x)=(x_{1}\sum_{{\bar{\imath}}\in\Omega_{1}}a^{(1)}_{\bar{\imath}}x^{\bar{\imath}},\dots,x_{n}\sum_{{\bar{\imath}}\in\Omega_{n}}a^{(n)}_{\bar{\imath}}x^{\bar{\imath}})^{T}$

and $A\in\mathbb{F}^{n,n}$ is a diagonal matrix.

Let $\tilde{L}$ be the $\ell\times n$ matrix which rows are all $n$-tuples
$\bar{\imath}$ from $\Omega_{1},\Omega_{2},\ldots,\Omega_{n}$ successively,
(that is, the rows of $\tilde{L}$ are the subscripts $\bar{\imath}$ of the
parameters $a^{(k)}_{\bar{\imath}}$ of (17) written as row vectors). For
$\nu=(\nu_{1},\dots,\nu_{\ell})\in\mathbb{N}^{\ell}_{0}$ we denote by $L$ the
map

$L(\nu)=\nu\tilde{L}.$ (19)

Then $L(\nu)$ is the row vector $L(\nu)=(L_{1}(\nu),\dots,L_{n}(\nu)$. Let us
add that $L$ is additive map on the monoid $\mathbb{N}_{0}^{\ell}$, hence
$L(0)=0$ and $L(j\mu)=jL(\mu)$ for any $j\in\mathbb{N}_{0}$ and
$\mu\in\mathbb{N}_{0}^{\ell}$.

Denote by ${a}$ the ordered (according to the order in $\Omega_{k}$,
$k=1,2,\dots,n$) $\ell$-tuple of parameters of equation (18),

${a}=(a^{(1)}_{\bar{\imath}^{(1)}},a^{(1)}_{\bar{\imath}^{(2)}},\dots,a^{(n)}_{\bar{\imath}^{(\ell)}})$
(20)

and by $\mathbb{F}[{a}]$ the ring of polynomials in variables
$a^{(1)}_{\bar{\imath}^{(1)}},\dots,a^{(n)}_{\bar{\imath}^{(\ell)}}$ over
$\mathbb{F}$. Any monomial in parameters of system (18) has the form

$a^{\nu}=(a^{(1)}_{\bar{\imath}^{(1)}})^{\nu_{1}}(a^{(1)}_{\bar{\imath}^{(2)}})^{\nu_{2}}\cdots(a^{(n)}_{\bar{\imath}^{(\ell)}})^{\nu_{\ell}},$
(21)

so, $a^{\nu}$ is the abbreviation of the expression on the right hand side of
(21). We also define the L-preimage of $\Omega$ to be

$\omega=L^{-1}(\Omega)$ (22)

being non-empty since $L(e_{1})=\bar{\imath}^{(1)}\in\Omega$.

In the case of system (15) we have $\Omega_{1}=((1,0,),(0,1),(-1,3))$,
$\Omega_{2}=((0,2),(1,0,),(0,1))$ and $\ell=6$. Then

$\tilde{L}=\left(\begin{matrix}1&0&-1&0&1&0\\\
0&1&3&2&0&1\end{matrix}\right)^{T}$

and
$a=(a^{(1)}_{10},a^{(1)}_{01},a^{(1)}_{-13},a^{(2)}_{02},a^{(2)}_{10},a^{(2)}_{01})$.
For any
$(\nu_{1},\nu_{2},\nu_{3},\nu_{4},\nu_{5},\nu_{6})\in\mathbb{N}_{0}^{6}$,

$a^{\nu}=(a^{(1)}_{10})^{\nu_{1}}(a^{(1)}_{01})^{\nu_{2}}(a^{(1)}_{-13})^{\nu_{3}}(a^{(2)}_{02})^{\nu_{4}}(a^{(2)}_{10})^{\nu_{5}}(a^{(2)}_{01})^{\nu_{6}}$

so,

$L(\nu)=(\nu_{1}-\nu_{3}+\nu_{5},\nu_{2}+3\nu_{3}+2\nu_{4}+\nu_{6}).$

We see that the range of $L$, in this example, is equal to
$\mathrm{im}L=\mathbb{Z}\times\mathbb{N}_{0}$.

In the general case, $\mathrm{im}L$ is an additive submonoid of
$\mathbb{Z}^{n}$. Moreover, by our assumptions in each column of $\tilde{L}$
at most one $-1$ can appear, and $|L(\mu)|\geq 1$ for every
$\mu\in\mathbb{N}_{0}^{n}$. Therefore, $L(\mu)=0\in\mathbb{N}_{-1}^{n}$
implies that $\mu=0$.

######  Definition 2.

For $m\in\mathrm{im}L\subseteq\mathbb{Z}^{n}$, a (Laurent) polynomial $p(a)$,
$p=\sum_{\nu\in\text{Supp}(p)}p^{(\nu)}a^{\nu}$, is an $m$-_polynomial_ if for
every $\nu\in\text{Supp}(p)\subset\mathbb{N}^{\ell}_{0}$, $L(\nu)=m$. The zero
polynomial is said to be an $m$-polynomial for every $m\in\mathbb{Z}^{n}$.

Remark. Observe that by our definition, in (15) every parameter $a^{(j)}_{rs}$
(that is, $a^{(1)}_{10},a^{(1)}_{01},...$ ) is an $(r,s)$-polynomial.
Moreover, the monomials $a^{(j)}_{rs}$ have the following property: if
$a^{(j)}_{rs}$ is the $q$-th element of (20), then $a^{(j)}_{rs}=a^{e_{q}}$,
$L(e_{q})=(r,s)$ and $a^{(j)}_{rs}x_{1}^{r}x_{2}^{s}=a^{e_{q}}x^{L(e_{q})}$
($e_{q}$ is the $q$-th unit basis vector of $\mathbb{F}^{\ell}$). A very
important observation is that all monomials $g_{3}$ and $g_{5}$ of (16), after
$x_{1}$ and $x_{2}$ being factored out, appear in the form
$a^{\nu}x^{L(\nu)}$. It is one of the cornerstones of our investigation.

As it has just been used, $e_{q}\in\mathbb{N}_{0}^{\ell}$ will be considered
as a row vector with 1 in the $q$-th position and zeroes elsewhere, and a
column vector (being defined before) if being the standard basis vector in
$\mathbb{F}^{n}$.

For a given $m\in\mathbb{Z}^{n}$ let $R_{m}$ be the subset of $\mathbb{F}[a]$
consisting of all $m$-polynomials. Denote by $R$ the direct sum of $R_{m}$,

$R=\oplus_{m\in\mathbb{Z}^{n}}R_{m}.$

Since

$R_{m_{1}}R_{m_{2}}\subseteq R_{m_{1}+m_{2}},$

$R$ is a graded ring. Clearly, by the remark preceding Definiton 2 we have
$R_{0}=\mathbb{F}$. Moreover, $R_{m}$ as well as $R$ are vector spaces over
$\mathbb{F}$ for the usual addition and multiplication by numbers from
$\mathbb{F}$.

Let

$\mathcal{M}=\\{(m_{1},\dots,m_{n})\in\mathbb{N}_{-1}^{n}:|m|\geq
0,m\in\mathrm{im}L,m_{j}=-1\text{ for at most one }j\\}.$

Let $M_{j}$ be the space of formal vector fields of the form

$\dot{x}_{j}=x_{j}\sum_{\begin{subarray}{c}m\in\mathcal{M},\\\ m_{i}\geq 0{\
\rm if\ }i\neq j\end{subarray}}p_{j}^{(m)}(a)x^{m}\qquad(j=1,\dots,n)$

where $p_{j}^{(m)}(a)\in R_{m}$ for all $j=1,\dots,n$.

Further, let

$M=M_{1}\times\dots\times M_{n}.$

Denoting $p^{(m)}(a)=\sum_{j=1}^{n}e_{j}p_{j}^{(m)}(a)$, for every member of
$M$ there exists a finite or infinite set $\mathcal{L}\subset\mathcal{M}$ such
that it can be written as

$\sum_{m\in{\mathcal{L}}}x\odot p^{(m)}(a)x^{m}$ (23)

with $\odot$ denoting the Hadamard multiplication. It can be easily seen that
$M$ is an additive group and moreover, $M$ is a module over the ring
$(\mathbb{F}^{n},+,\odot)$ with the multiplicative unit
$\bar{1}=\sum_{i=1}^{n}e_{i}\in\mathbb{F}^{n}$. Obviously, $M$ is also a
vector space over $\mathbb{F}$.

We note that the grading introduced in this section was implicitly used in
[14] in order to derive an efficient algorithm for computing focus (Lyapunov)
quantities of planar polynomial vector fields.

###  3.2 Lie brackets in $M$

We consider elements of $M$ as (formal) vector fields, so we can compute Lie
brackets of elements of $M$.

######  Lemma 1.

Any element ${\bf\Theta}=(\Theta^{(1)},\dots,\Theta^{(n)})\in M$ can be
written in the form

${\bf\Theta}=\sum_{\mu\in\omega}(\theta_{\mu}\odot x)a^{\mu}x^{L(\mu)},$ (24)

where $\theta_{\mu}=(\theta_{\mu,1},\dots,\theta_{\mu,n})$, $\omega$ is a
finite or infinite subset of $\mathbb{N}^{\ell}_{0}$ such that if
$\mu\in\omega$ then $L(\mu)\in{\mathcal{M}}$. Additionally, if
$L_{j}(\mu)=-1$, then $\theta_{\mu}=\theta_{\mu,j}e_{j}$.

###### Proof.

Every ${\bf\Theta}\in M$ can be represented in the form (23). We rewrite the
polynomials $p^{(m)}(a)$, applying that each component $p^{(m)}_{j}(a)$ is an
$m$-polynomial, as

$p^{(m)}(a)=\sum_{j=1}^{n}p_{j}^{(m)}(a)e_{j}=\sum_{j=1}^{n}\sum_{\mu:L(\mu)=m}\theta_{\mu,j}e_{j}a^{\mu}=\sum_{\mu:L(\mu)=m}\theta_{\mu}a^{\mu}$

where $\theta_{\mu}=\sum_{j=1}^{n}\theta_{\mu,j}e_{j}$. It follows that

${\bf\Theta}=\sum_{m\in\mathcal{L}}x\odot\sum_{\mu:L(\mu)=m}\theta_{\mu}a^{\mu}x^{m}=\sum_{\mu:L(\mu)\in\mathcal{L}}(\theta_{\mu}\odot
x)a^{\mu}x^{L(\mu)}$

where $\omega=\\{\mu:L(\mu)\in\mathcal{L}\subseteq\mathcal{M}\\}$.

The last assertion follows by the definition of the set $\mathcal{M}$. ∎

As every ${\bf\Theta}\in M$ is a (formal) vector field, we can compute Lie
brackets of such vector fields:

$[{\bf\Theta},{\bf\Phi}]:=(D{\bf\Phi}){\bf\Theta}-(D{\bf\Theta}){\bf\Phi}$

for any ${\bf\Theta},{\bf\Phi}\in M$.

######  Lemma 2.

If ${\bf\Theta}=(\theta_{\mu}\odot x)a^{\mu}x^{L(\mu)}$ and
${\bf\Phi}=(\phi_{\nu}\odot x)a^{\nu}x^{L(\nu)}$, where
$\mu,\nu\in\mathbb{N}_{0}^{\ell}$,
$\theta_{\mu}=(\theta_{1},\dots,\theta_{n})^{T}$ and
$\phi_{\nu}=(\phi_{1},\dots,\phi_{n})^{T}$, then

$[{\bf\Theta},{\bf\Phi}]=(\alpha_{\mu+\nu}\odot x)a^{\mu+\nu}x^{L(\mu+\nu)}\in
M,$ (25)

with

$\alpha_{\mu+\nu}=\langle L(\nu),\theta_{\mu}\rangle{\phi_{\nu}}-\langle
L(\mu),{\phi_{\nu}}\rangle\theta_{\mu}.$ (26)

###### Proof.

We recall that for $f:\mathbb{F}^{n}\to\mathbb{F}^{n}$ the derivative of the
Hadamard product reads

$D(x\odot f(x))={\rm diag}(f)+{\rm diag}(x)Df(x).$ (27)

Let $x^{-\bar{1}}=(x_{1}^{-1},\dots,x_{n}^{-1})^{T}$. Using (27) we get

$\displaystyle D(x\odot\phi_{\nu}\,x^{L(\nu)})$
$\displaystyle=\mathrm{diag}(\phi_{\nu})x^{L(\nu)}+\mathrm{diag}(x)\phi_{\nu}(L(\nu)^{T}\odot
x^{-\bar{1}})^{T}x^{L(\nu)}$
$\displaystyle=\mathrm{diag}(\phi_{\nu})x^{L(\nu)}+(\phi_{\nu}\odot
x)(L(\nu)^{T}\odot x^{-\bar{1}})^{T}x^{L(\nu)}$

Taking into account (27), we have

$(D{\bf\Phi}){\bf\Theta}=(D((\phi_{\nu}\odot
x)a^{\nu}x^{L(\nu)})){\bf\Theta}=(D(x\odot\phi_{\nu}x^{L(\nu)})a^{\nu})){\bf\Theta}=\\\
(a^{\nu}({\rm diag}(\phi_{\nu})x^{L(\nu)}+(\phi_{\nu}\odot x)(L(\nu)^{T}\odot
x^{-\bar{1}})^{T}x^{L(\nu)})){\bf\Theta}=\\\
a^{\nu+\mu}(\phi_{\nu}\odot\theta_{\mu}\odot
x)x^{L(\nu+\mu)}+a^{\nu+\mu}(\phi_{\nu}\odot x)(L(\nu)^{T}\odot
x^{-\bar{1}})^{T}(x\odot\theta_{\mu})x^{L(\nu+\mu)}=\\\
a^{\nu+\mu}x^{L(\nu+\mu)}(\phi_{\nu}\odot\theta_{\mu}\odot x+\langle
L(\nu),\theta_{\mu}\rangle\phi_{\nu}\odot x)$

and similarly,

$(D{\bf\Theta}){\bf\Phi}=a^{\nu+\mu}x^{L(\nu+\mu)}(\phi_{\nu}\odot\theta_{\mu}\odot
x+\langle L(\mu),\phi_{\nu}\rangle\theta_{\mu}\odot x).$ (28)

Now, (25) takes place.

To verify the last assertion that $[{\bf\Theta},{\bf\Phi}]\in M$ one has to
recall the definition of the set $\mathcal{M}$. The vector field
$[{\bf\Theta},{\bf\Phi}]$ would have not been in $M$ when $L(\mu+\nu)$ was not
in $\mathcal{M}$. We have to consider two possibilities. Either
$L_{i}(\mu)=L_{i}(\nu)=-1$ for some $i$ or, for some $i\neq j$ we have
$L_{i}(\mu)=-1$, $L_{i}(\nu)=0$, $L_{j}(\mu)=0$ and $L_{j}(\nu)=-1$.

In the first case, $\theta_{\mu}=\theta_{i}e_{i}$, $\phi_{\nu}=\phi_{i}e_{i}$
and so,

$\alpha_{\mu+\nu}=L_{i}(\nu)\theta_{i}\phi_{i}e_{i}-L_{i}(\mu)\phi_{i}\theta_{i}e_{i}=0.$

In the second case, we have $\theta_{\mu}=\theta_{i}e_{i}$,
$\phi_{\nu}=\phi_{j}e_{j}$ and again, $\alpha_{\mu+\nu}=0$. ∎

The following statement is a direct corollary of Lemma 2.

######  Corollary 1.

If ${\bf\Theta},{\bf\Phi}\in M$, then $[{\bf\Theta},{\bf\Phi}]\in M$.
Moreover, if

${\bf\Phi}=\sum_{\nu\in\omega_{1}}(\phi_{\nu}\odot x)a^{\nu}x^{L(\nu)},\ \
{\bf\Theta}=\sum_{\mu\in\omega_{2}}(\theta_{\mu}\odot x)a^{\mu}x^{L(\mu)},$

then

$[{\bf\Theta},{\bf\Phi}]=\sum_{\mu\in\omega}\sum_{\nu\in\omega_{1}}\alpha_{\mu+\nu}\
a^{\mu+\nu}x^{L(\mu+\nu)}.$ (29)

with $\alpha_{\mu+\nu}$ being defined by (26).

Let $\mathcal{U}_{s}^{\ell}$, $s\geq 0$, be the space of polynomial vector
fields of the form

${\bf\Theta}_{s}=\sum_{|\mu|=s}(\theta(\mu)\odot x)a^{\mu}x^{L(\mu)}.$ (30)

Then, obviously, $M=\oplus_{s=0}^{\infty}\mathcal{U}_{s}^{\ell}$.

For any $s\in\mathbb{N}_{0}$, let

$U_{s}=\\{(\bar{1}\odot x)a^{\mu}x^{L(\mu)}:|\mu|=s\\}.$ (31)

Since any vector field ${\bf\Theta}_{s}$ of the form (30) can be written as

${\bf\Theta}_{s}=\sum_{|\mu|=s}\theta(\mu)\odot((\bar{1}\odot
x)a^{\mu}x^{L(\mu)}),$

we can treat $\mathcal{U}_{s}^{\ell}$ as a free module generated by $U_{s}$
over the ring $(\mathbb{F}^{n},+,\odot)$ with the operation of multiplication
by the elements of the ring being the Hadamard product. The dimension of
$\mathcal{U}_{s}^{\ell}$ as the module is the same as the dimension of the
space of homogeneous polynomials of degree $s$ in $\ell$ variables, that is,
it is ${\ell+s-1\choose s}$.

######  Definition 3.

We say that ${\bf\Theta}\in M$ is of level $s$ if
${\bf\Theta}\in\mathcal{U}_{s}^{\ell}$ and ${\bf\Theta}\in M$ is of the level
at least $s$ if each term

$\Theta=(\theta\odot x)a^{\mu}x^{L(\mu)}$

of ${\bf\Theta}$ is in some of $\mathcal{U}_{s+j}^{\ell}$, where
$j\in\mathbb{N}_{0}$.

By Lemma 2 if ${\bf\Theta}\in\mathcal{U}_{s}^{\ell}$,
${\bf\Phi}\in\mathcal{U}_{t}^{\ell}$ and $\text{ad}\ {\bf\Theta}$ is the
adjoint operator acting on $\bf\Phi$ by

$(\text{ad\, }{\bf\Theta})\,{\bf\Phi}=[{\bf\Theta},{\bf\Phi}],$

then $(\text{ad\, }{\bf\Theta})^{i}{\bf\Phi}$ is an element of
$\mathcal{U}_{is+t}^{\ell}$, that is,

$(\text{ad\,
}{\bf\Theta})^{i}:\mathcal{U}_{t}^{\ell}\normalcolor\to\mathcal{U}_{is+t}^{\ell}.$
(32)

Remark. The map $(\text{ad\, }{\bf\Theta})^{i}$ ”lifts” the space
$\mathcal{U}_{t}^{\ell}$ to the space $\mathcal{U}_{is+t}^{\ell}.$ For this
reason we speak about levels in the normalization process.

We denote the set of $\mu$’s appearing in (30) by $\sigma(s)$,

$\sigma(s):=\\{\mu\in\mathbb{N}^{\ell}_{0}\ :|\mu|=s\\}.$

###  3.3 The normal form algorithm

Any equation of the form (18) can be written as a differential equation of the
form

$\dot{x}=a_{0}(x)+a_{1}(x)+a_{2}(x)+\dots=a(x),$ (33)

where $a_{0}(x)=\lambda\odot x,$ $a_{s}(x)\in\mathcal{U}_{s}^{\ell}$ for
$s\geq 1$,

$a_{s}(x)=\sum_{\mu\in\sigma(s)}(\alpha_{\mu}\odot x)a^{\mu}x^{L(\mu)},$

where $\alpha_{\mu}\in\mathbb{F}^{n}.$ Since
$a_{0}(x)\in\mathcal{U}_{0}^{\ell}$ and the nonlinear part $a_{1}(x)$ of (18)
is from $\mathcal{U}_{1}^{\ell}$, equation (18) is of the form (33).

For instance, system (15) is written as equation (33) with

$a_{0}(x)={1\choose-1}\odot{x_{1}\choose x_{2}},$ $\displaystyle a_{1}(x)=$
$\displaystyle{1\choose 0}\odot{x_{1}\choose x_{2}}a^{(1)}_{10}x_{1}+{1\choose
0}\odot{x_{1}\choose x_{2}}a^{(1)}_{01}x_{2}+{1\choose 0}\odot{x_{1}\choose
x_{2}}a^{(1)}_{-13}x_{1}^{-1}x_{2}^{3}+$ $\displaystyle{0\choose
1}\odot{x_{1}\choose x_{2}}a^{(2)}_{10}x_{1}+{0\choose 1}\odot{x_{1}\choose
x_{2}}a^{(2)}_{01}x_{2}+{0\choose 1}\odot{x_{1}\choose
x_{2}}a^{(2)}_{02}x_{2}^{2}.$

Note that (17) is just a special case of (33).

Below we will present an algorithm which normalizes equation (33) up to a
certain level.

######  Lemma 3.

If

$\langle L(\mu),\lambda\rangle=0,$ (34)

then all entries of the vector field

$(\alpha_{\mu}\odot x)a^{\mu}x^{L(\mu)}$ (35)

are resonant terms.

###### Proof.

The $i$-th entry of (35) is

$\alpha^{(\mu)}_{i}x_{i}x^{L(\mu)}a^{\mu}=\alpha^{(\mu)}_{i}x_{1}^{L_{1}(\mu)}\cdots
x_{i-1}^{L_{i-1}(\mu)}x_{i}^{L_{i}(\mu)+1}x_{i+1}^{L_{i+1}(\mu)}\cdots
x_{n}^{L_{n}(\mu)}a^{\mu}.$

Since

$\langle{L_{1}(\mu)},\dots,{L_{i-1}(\mu)},{L_{i}(\mu)+1},{L_{i+1}(\mu)},\dots,{L_{n}(\mu)},\lambda\rangle-\lambda_{i}=\langle
L(\mu),\lambda\rangle,$

the statement of the lemma takes place. ∎

Lemma 3 justifies the following definition.

######  Definition 4.

It is said that a term of the form (35) of the right-hand side of (33) is
_resonant_ if (34) holds.

Under known approaches a normalization is performed up to terms of some degree
of $x$. However, it is possible to perform the normalization up to some degree
of the polynomials in the parameters of the system.

######  Definition 5.

We say that equation (33) is in the normal form up to level $s$ if all non-
resonant terms in $a_{1}(x),\dots,a_{s}(x)$ are equal to zero.

Clearly, if equation (33) is in the normal form for all levels
$s\in\mathbb{N}$, then it is in the normal form in the sense of Definition 1.

By definition we say that the operator $\pounds^{a}:M\to M$ which acts on
elements of the form (24) by

$\pounds^{a}(\sum_{\mu\in\omega}(\theta_{\mu}\odot
x)a^{\mu}x^{L(\mu)})=\sum_{\mu\in\omega}\langle
L(\mu),\lambda\rangle(\theta_{\mu}\odot x)a^{\mu}x^{L(\mu)})$ (36)

is the homological operator  of equation (33) (the superscript $a$ indicates
that the operator is related to the equation (33) with the right-hand sides
being an element of the module $M$, so it is different from the homological
operator of (4) defined by (8). It is not difficult to check that
$\pounds^{a}$ is a homomorphism of the module $M$ over $\mathbb{F}^{n}$.
Clearly being additive, we also have that for any $\phi\in\mathbb{F}^{n}$

$\displaystyle\pounds^{a}(\phi\odot\sum_{\mu\in\omega}(\theta{(\mu)}\odot
x)a^{\mu}x^{L(\mu)})$ $\displaystyle=\sum_{\mu\in\omega}\langle
L(\mu),\lambda\rangle(\phi\odot\theta{(\mu)}\odot x)a^{\mu}x^{L(\mu)}$
$\displaystyle=\phi\odot\pounds^{a}(\sum_{\mu\in\omega}(\theta{(\mu)}\odot
x)a^{\mu}x^{L(\mu)})).$

The restriction of $\pounds^{a}$ on $\mathcal{U}_{s}^{\ell}$ is denoted by
$\pounds^{a}_{s}$. Obviously,
$\pounds^{a}_{s}:\mathcal{U}_{s}^{\ell}\to\mathcal{U}_{s}^{\ell}$. From (36)
it is clear that the set $U_{s}$ defined by (31) is the set of basis
eigenvectors of $\pounds^{a}_{s}$ and we can split $\mathcal{U}_{s}^{\ell}$ as

$\mathcal{U}_{s}^{\ell}={\rm im}\pounds_{s}^{a}\oplus\ker\pounds_{s}^{a}.$

Assume that equation (33) is in the normal form up to level $s-1$, $s\geq 1$,
that is, for terms of the form (35) appearing in (33) if $|\mu|\leq s-1$ and
$\langle L(\mu),\lambda\rangle\neq 0$, then $\alpha{(\mu)}=0$. Then the
homological equation

$\pounds^{a}_{s}(h_{s})=a_{s}-g_{s}$

can be solved for $h_{s}$ and $g_{s}$ as follows:

$h_{s}(x)=\sum_{\begin{subarray}{c}\mu:\mu\in\sigma(s),\\\ \langle
L(\mu),\lambda\rangle\neq 0\end{subarray}}\frac{1}{\langle
L(\mu),\lambda\rangle}(\alpha_{\mu}\odot x)a^{\mu}x^{L(\mu)},$ (37)
$g_{s}(x)=\sum_{\begin{subarray}{c}\mu:\mu\in\sigma(s),\\\ \langle
L(\mu),\lambda\rangle=0\end{subarray}}(\alpha_{\mu}\odot x)a^{\mu}x^{L(\mu)}.$
(38)

######  Theorem 1.

Assume that equation (33) is in the normal form up to level $s-1$, $s\geq 1$,
and let

$H_{s}(x)=\exp(h_{s}(x)),$ (39)

where $h_{s}$ is defined by (37). Then performing the substitution
$y=H_{s}(x)$ and changing $y$ to $x$ we obtain from (33) an equation, which
the right hand side is from $M$ and is in the normal form up to level $s$.

###### Proof.

According to the theorem on Lie series for vector fields after transformation
(39) we obtain from (18) the vector field

$a(x)+(\text{ad\ }h_{s})a(x)+\sum_{i=2}^{\infty}\frac{1}{i!}(\text{ad\
}h_{s})^{i}a(x).$ (40)

By (32) the last summand is of the level at least $s+1$ and for the first two,
we have

$a(x)+(\text{ad\ }h_{s})a(x)=a(x)+[h_{s}(x),a(x)]=\\\
a_{0}(x)+a_{1}(x)+\dots+a_{s-1}(x)+a_{s}(x)+[h_{s}(x),a_{0}(x)]+\dots,$

where the dots stand for the terms of the level at least $s+1$.

By (25), (37) and (38) $a_{s}(x)+[h_{s}(x),a_{0}(x)]=g_{s}(x)$. Since by our
assumption (33) is in the normal form up to level $s-1$, the expression (40)
is in the normal form up to level $s$. ∎

As a direct outcome of Theorem 1 we have the following statement.

######  Corollary 2.

There are polynomial maps $H_{1}(x),\dots,H_{s}(x)$, such that equation (18)
is transformed to an equation which is in the normal form up to level $s$ by
the transformation $y=H_{s}\circ\cdots\circ H_{1}$.

From the results presented above we have the following algorithm for the
normalization of equation (33) up to level $m$.

Algorithm A.

Set $a_{0}(x):=Ax$, ${\bf\Xi}_{k}(x):=0$ for $k=1,2,\dots,m.$

For $s=1,\dots,m$ do the following:

  1. (i)

Define $h_{s}$ and $g_{s}$ by (37) and (38);

  2. (ii)

Compute

${\bf\Xi}=\sum_{k=0}^{s-1}\sum_{i=1}^{\lfloor\frac{m-k}{s}\rfloor}\frac{1}{i!}(\text{ad\
}h_{s}(x))^{i}a_{k}(x)$

and write ${\bf\Xi}=\sum_{i=s+1}^{m}{\bf\Xi}_{i},$ where
${\bf\Xi}_{i}\in\mathcal{U}^{\ell}_{i}$;

  3. (iii)

Let $a_{s}=g_{s}$, $a_{i}={\bf\Xi}_{i}$ for $i=s+1,\dots,m$.

We see that the algorithm is different from Algorithm 1, in particular, the
calculations are performed in different spaces. In the next section, we
reformulate the algorithm in a simpler setting and give an example of
computing the normal form of equation (15).

######  Proposition 1.

If system (33) is in normal form up to level $s$, then it is in the normal
form up to order at least $s+1$.

###### Proof.

For each multiindex $\bar{\imath}\in\Omega$ it holds $|\bar{\imath}|\geq 1$.
Therefore, by (19), if a term of the form $(\theta\odot x)a^{\mu}x^{L(\mu)}$
is of level $s$ then $|L(\mu)|\geq s$. Thus, all resonant terms in the normal
form of level $s$ or higher are of order $s+1$ or higher. ∎

##  4 Generalized vector fields

We treat $\mathbb{N}_{0}^{\ell}$ as an ordered set (ordered, for instance,
using the degree lexicographic order).

######  Definition 6.

Let $\alpha$ be a map defined on some subset $\omega$ of
$\mathbb{N}_{0}^{\ell}$

$\alpha:\omega\subset\mathbb{N}_{0}^{\ell}\to\mathbb{F}^{n},$

that is, $\alpha$ assigns to every $\nu\in\omega$ an $n$-tuple

$\alpha_{\nu}=(\alpha_{1}{(\nu)},\dots,\alpha_{n}{(\nu)}).$

We say that an $n$-tuple of formal power series

$\hat{\alpha}=\sum_{\nu\in\omega}\alpha_{\nu}a^{\nu},$ (41)

where $\omega=\text{Supp}(\hat{\alpha})$, is a _generalized formal vector
field_.

Note that expression (41) looks similar to a usual vector field depending on
variables
$a^{(1)}_{\bar{\imath}^{(1)}}),a^{(1)}_{\bar{\imath}^{(2)}},\dots,a^{(n)}_{\bar{\imath}^{(\ell_{n})}}$
defined on $\mathbb{F}^{\ell}$, however, it is not a vector field in these
variables in the usual sense, because the usual vector field is defined by
assigning to a vector from $\mathbb{F}^{\ell}$ a vector of the same dimension,
but if a series (41) converges it assigns to a point from $\mathbb{F}^{\ell}$
a vector from $\mathbb{F}^{n}$.

Since the objects defined above are not vector fields in the usual sense, we
call them generalized vector fields.

We denote the set of all formal vector fields defined by (41) by
$\mathcal{A}$. It is not difficult to check that $\mathcal{A}$ is a module
over the ring $(\mathbb{F}^{n},+,\odot)$. The zero vector in $\mathcal{A}$ is
a series (41) where $\alpha_{\nu}=0\in\mathbb{F}^{n}$ for all
$\nu\in\mathbb{N}_{0}^{\ell}$. For any $k\in\mathbb{N}_{0}$, let
$\mathcal{A}_{k}$ be the subset of all elements of $\mathcal{A}$ of the form

$\sum_{\mu:|\mu|=k}\alpha_{\mu}a^{\mu}.$

Then $\mathcal{A}_{k}$ is a module over the ring $(\mathbb{F}^{n},+,\odot)$
and $\mathcal{A}$ is a direct sum of $\mathcal{A}_{k}$, $k=0,1,\dots$.

Recall that we consider $M$ as the direct sum of modules
$\mathcal{U}^{\ell}_{s}$ over $\mathbb{F}^{n}$, $s=0,1,2,\dots$, and define a
module homomorphism

$\mathfrak{T}:\mathcal{A}\to M$
$\mathfrak{T}\left(\sum_{\mu\in\omega}\theta_{\mu}a^{\mu}\right)=\sum_{\mu\in\omega}(\theta_{\mu}\odot
x)a^{\mu}x^{L(\mu)}.$ (42)

Clearly, $\mathfrak{T}$ is an isomorphism.

We define the Lie bracket of
$\hat{\theta}=\sum_{\mu\in\omega}\theta_{\mu}a^{\mu}$ and
$\hat{\phi}=\sum_{\nu\in\omega_{1}}\phi_{\nu}a^{\nu}$ by

$[\hat{\theta},\hat{\phi}]=\mathfrak{T}^{-1}([\mathfrak{T}(\hat{\theta}),\mathfrak{T}(\hat{\phi})]).$

Then, by (29) and (42)

$[\hat{\theta},\hat{\phi}]=\sum_{\mu\in\omega}\sum_{\nu\in\omega_{1}}\left(\langle
L(\nu),\theta_{\mu}\rangle\phi_{\nu}-\langle
L(\mu),\phi_{\nu}\rangle\theta_{\mu}\right)a^{\mu+\nu}.$ (43)

Since

$\mathfrak{T}\left([\hat{\psi},[\hat{\theta},\hat{\phi}]]\right)=[\mathfrak{T}(\hat{\psi}),\mathfrak{T}([\hat{\theta},\hat{\phi}])]=[\mathfrak{T}(\hat{\psi}),[\mathfrak{T}(\hat{\theta}),\mathfrak{T}(\hat{\phi})])],$

it is easily seen that for the Lie bracket in $\mathcal{A}$ defined by (43)
the Jacobi identity holds, so $\mathcal{A}$, being also a vector space over
$\mathbb{F}$, is a Lie algebra and so $\mathfrak{T}$ defines a Lie algebra
isomorphism.

We could have defined the Lie bracket in $\mathcal{A}$ in a different way. The
$n\times n$ matrix
$\mathfrak{D}(\theta_{\mu}a^{\mu}):=\theta_{\mu}L(\mu)a^{\mu}$, associated to
the monomial $\theta_{\mu}a^{\mu}$, can be extended to a
$\mathbb{F}^{n}$-module homomorphism (or linear operator) from
$\mathbb{F}^{n}$ to $\mathbb{F}^{n}$ inducing the Lie bracket. Indeed, the
setting

$[\hat{\theta},\hat{\phi}]=\mathfrak{D}(\phi_{\nu}a^{\nu})\theta_{\mu}a^{\mu}-\mathfrak{D}(\theta_{\mu}a^{\mu})\phi_{\nu}a^{\nu}$

produces the same as (43) which can be observed from

$\displaystyle\mathfrak{D}(\phi_{\nu}a^{\nu})\theta_{\mu}a^{\mu}-\mathfrak{D}(\theta_{\mu}a^{\mu})\phi_{\nu}a^{\nu}$
$\displaystyle=\phi_{\nu}L(\nu)\theta_{\mu}a^{\mu+\nu}-\theta_{\mu}L(\mu)\phi_{\nu}a^{\mu+\nu}$
$\displaystyle=\left(\langle L(\nu),\theta_{\mu}\rangle\phi_{\nu}-\langle
L(\mu)\phi_{\nu}\rangle\theta_{\mu}\right)a^{\mu+\nu}.$

One may expect that this operator is the derivative $D(\hat{\Theta})$
evaluated at $x=\bar{1}$ but it is not the case. However, the connection
exists, the second term of (28) evaluated at $x=\bar{1}$, coincides with
$\mathfrak{D}(\theta_{\mu}a^{\mu})$. So, we may think of $\mathfrak{D}$ as a
derivative.

We will compute normal forms in terms of such generalized vector fields in the
following way.

Let $\hat{\alpha}$ be the image of the right hand side of (33) under the
isomorphism $\mathfrak{T}^{-1}$, so

$\mathfrak{T}^{-1}(a(x))=\sum_{k=0}^{\infty}\hat{\alpha}_{k},$

where

$\hat{\alpha}_{0}=\lambda,\quad\hat{\alpha}_{k}=\sum_{\mu\in\sigma(k)}\alpha_{\mu}a^{\mu}\quad\text{for}\
k\geq 1.$ (44)

######  Definition 7.

It is said that the generalized vector field
$\hat{\alpha}=\sum_{k=0}^{\infty}\hat{\alpha}_{k}$ (where $\hat{\alpha}_{k}$
is of the form (44)) is in the normal form up to level $s$ if the coefficients
of all non-resonant terms in $\hat{\alpha}_{1},\dots,\hat{\alpha}_{s}$ are
equal to zero. The generalized vector field $\hat{\alpha}_{0}$ will be always
considered to be in the normal form up to level $0$.

Assume that $\hat{\alpha}$ is in the normal form up to level $s-1$, $s\geq 1$.
Let

$\displaystyle\hat{\eta}_{s}$
$\displaystyle=\sum_{\begin{subarray}{c}\mu:\mu\in\sigma(s),\\\ \langle
L(\mu),\lambda\rangle\neq 0\end{subarray}}\frac{1}{\langle
L(\mu),\lambda\rangle}\alpha_{\mu}a^{\mu},$ (45)
$\displaystyle\hat{\zeta}_{s}$
$\displaystyle=\sum_{\begin{subarray}{c}\mu:\mu\in\sigma(s),\\\ \langle
L(\mu),\lambda\rangle=0\end{subarray}}\alpha_{\mu}a^{\mu}.$

that is $\hat{\eta}_{s}=\mathfrak{T}^{-1}(h_{s}(x)),$ and
$\hat{\zeta}_{s}=\mathfrak{T}^{-1}(g_{s}(x)),$ where $h_{s}(x)$ and $g_{s}(x)$
are defined by (37) and (38), respectively. Since $\mathfrak{T}$ is a Lie
algebra isomorphism we can reformulate Algorithm A as follows.

Algorithm B.

Set $\hat{\alpha}_{0}:=\lambda$, $\hat{\xi}_{k}:=0$ for $k=1,2,\dots,m$.

For $s=1,\dots,m$ do the following:

  1. (i)

Define $\hat{\eta}_{s}$ and $\hat{\zeta}_{s}$ by (45);

  2. (ii)

Compute

$\hat{\xi}=\sum_{k=0}^{s-1}\sum_{i=1}^{\lfloor\frac{m-k}{s}\rfloor}\frac{1}{i!}(\text{ad\
}\hat{\eta}_{s})^{i}\hat{\alpha}_{k}$

(where $\text{ad\,}\hat{\eta}_{s}:=[\hat{\eta}_{s},\cdot]$ is the adjoint
operator acting on $\mathcal{A}$) and represent $\hat{\xi}$ in the form
$\hat{\xi}=\sum_{i=s}^{m}\hat{\xi}_{i}$, where
$\hat{\xi}_{i}\in\widehat{U}_{i}^{\ell}$;

  3. (iii)

Let $\hat{\alpha}_{s}=\hat{\zeta}_{s}$,
$\hat{\alpha}_{s+1}=\hat{\xi}_{s+1},\dots,\hat{\alpha}_{m}=\hat{\xi}_{m}$.

The obtained vector field $\hat{\alpha}$ is in the normal form up to level
$m$.

By Proposition 1 if system (18) is in normal form up to level $s$ then it is
in the normal form up to order at least $s+1$. Obviously, when $\hat{\alpha}$
is in the normal form up to level $m$, a normal form of (18) can be built up
from $\hat{\alpha}$ using the following procedure:

Set $g_{k}(x)=(0,0,\dots,0)^{T}$ for $k=1,\dots,m$. For
$\mu\in\cup_{k=1}^{m}{\sigma(k)}$ do the following:
if $\alpha_{\mu}\neq 0,\ |L(\mu)|=k$, then
$g_{k}(x)=g_{k}(x)+(\alpha_{\mu}\odot x)a^{\mu}x^{L(\mu)}$.

Clearly, the obtained equation

$\dot{x}=Ax+\sum_{k=1}^{m}g_{k}(x)$

is the normalization of (18) up to order $m$.

Example. As an example of applying the algorithm we consider the computation
of normal form of system (15). For $\mu\in\mathbb{N}_{0}^{\ell}$ we will use
the abbreviation $[\mu]=[\mu_{1},\dots,\mu_{\ell}]:=a^{\mu}$.

By Proposition 1 in order to compute the normal form of (15) up to order 5 it
is sufficient to compute the normal form of (15) up to level 4.

At the level $0$ the set $\sigma(0)$ consists of only one vector,
$(0,0,0,0,0,0)$ with

$\hat{\alpha}_{0}=\alpha_{(0,0,0,0,0,0)}={1\choose-1}[0,0,0,0,0,0].$

Passing to the level 1, $\sigma(1)$ is the set of vectors

$e_{1}^{T},\dots,e_{6}^{T},$ (46)

which form the standard basis of $\mathbb{Z}^{6}$. The vector field
$\hat{\alpha}_{1}$ is obtained by using the nonlinear terms of (15):

$\displaystyle\hat{\alpha}_{1}$ $\displaystyle={1\choose
0}[e_{1}^{T}]+{1\choose 0}[e_{2}^{T}]+{1\choose 0}[e_{3}^{T}]+{0\choose
1}[e_{4}^{T}]+{0\choose 1}[e_{5}^{T}]+{0\choose 1}[e_{6}^{T}]$
$\displaystyle={1\choose 0}a^{(1)}_{10}+{1\choose 0}a^{(1)}_{01}+{1\choose
0}a^{(1)}_{-13}+{0\choose 1}a^{(2)}_{02}+{0\choose 1}a^{(2)}_{10}+{0\choose
1}a^{(1)}_{01}.$

Then by (i) of Algorithm B for $s=1$ we have $\zeta_{1}=\bar{0}$ and

$\displaystyle\hat{\eta}_{1}$ $\displaystyle={1\choose 0}[e_{1}^{T}]-{1\choose
0}[e_{2}^{T}]-\frac{1}{4}{1\choose 0}[e_{3}^{T}]-\frac{1}{2}{0\choose
1}[e_{4}^{T}]+{0\choose 1}[e_{5}^{T}]-{0\choose 1}[e_{6}^{T}]$
$\displaystyle={1\choose 0}a^{(1)}_{10}-{1\choose
0}a^{(1)}_{01}-\frac{1}{4}{1\choose 0}a^{(1)}_{-13}-\frac{1}{2}{0\choose
1}a^{(2)}_{02}+{0\choose 1}a^{(2)}_{10}-{0\choose 1}a^{(1)}_{01}.$

When we know the level $s$, that is the set $\sigma(s)$, the next level, the
set $\sigma(s+1)$, is obtained by adding to the elements of $\sigma(s)$, one
of vectors (46).

According to the Algorithm B we have to set

$\hat{\alpha}_{2}=\hat{\alpha}_{3}=\hat{\alpha}_{4}=(0,0,\dots,0)^{T}.$ (47)

Next, we compute $\hat{\xi}_{1}+\dots+\hat{\xi}_{4}$ according to (ii); that
is, we compute the sum

$\displaystyle\phantom{\hat{\alpha}_{0}+(\text{ad\,}\hat{\eta}_{1})\hat{\alpha}_{0}+}\frac{1}{2}(\text{ad\,}\hat{\eta}_{1})^{2}\hat{\alpha}_{0}+\frac{1}{3!}(\text{ad\,}\hat{\eta}_{1})^{3}\hat{\alpha}_{0}+\frac{1}{4!}(\text{ad\,}\hat{\eta}_{1})^{4}\hat{\alpha}_{0}+$
(48)
$\displaystyle\phantom{\hat{\alpha}_{1}+}(\text{ad\,}\hat{\eta}_{1})\hat{\alpha}_{1}+\frac{1}{2}(\text{ad\,}\hat{\eta}_{1})^{2}\hat{\alpha}_{1}+\frac{1}{3!}(\text{ad\,}\hat{\eta}_{1})^{3}\hat{\alpha}_{1}+$
$\displaystyle\hat{\alpha}_{2}+(\text{ad\,}\hat{\eta}_{1})\hat{\alpha}_{2}+\frac{1}{2}(\text{ad\,}\hat{\eta}_{1})^{2}\hat{\alpha}_{2}+$
$\displaystyle\hat{\alpha}_{3}+(\text{ad\,}\hat{\eta}_{1})\hat{\alpha}_{3}+$
$\displaystyle\hat{\alpha}_{4}.$

Then, for level 1, we have

$(\text{ad\,}\hat{\eta}_{1})\hat{\alpha}_{0}+\hat{\alpha}_{1}=\bar{0}.$ (49)

One can observe that the terms on the diagonals of (48) are at the same
levels, that is, they are from $\mathcal{U}^{6}_{2},\mathcal{U}^{6}_{3}$ and
$\mathcal{U}^{6}_{4}$, respectively. Thus, for level 2, we obtain

$\hat{\alpha}_{2}+(\text{ad\,}\hat{\eta}_{1})\hat{\alpha}_{1}+\frac{1}{2}(\text{ad\,}\hat{\eta}_{1})^{2}\hat{\alpha}_{0}.$
(50)

By (47) $\hat{\alpha}_{2}=0$ and for the second term of the sum given above

$\displaystyle(\text{ad\,}\hat{\eta}_{1})\hat{\alpha}_{1}=$
$\displaystyle{-2\choose 0}[1,1,0,0,0,0]+{-5/2\choose
0}[1,0,1,0,0,0]+{3/4\choose 0}[0,1,1,0,0,0]+$ $\displaystyle{1/2\choose
0}[0,1,0,1,0,0]+{2\choose-2}[0,1,0,0,1,0]+{-3/4\choose 0}[0,0,1,1,0,0]+$
$\displaystyle{15/4\choose-5/4}[0,0,1,0,1,0]+{-9/4\choose
0}[0,0,1,0,0,1]+{0\choose 3}[0,0,0,1,1,0]+$
$\displaystyle{0\choose-1/2}[0,0,0,1,0,1]+{0\choose 2}[0,0,0,0,1,1].$

From (49) we observe that
$(\text{ad\,}\hat{\eta}_{1})\,\hat{\alpha}_{0}=-\hat{\alpha}_{1}$ and using
this in (50) gives

$\displaystyle\frac{1}{2}(\text{ad\ }\eta_{1})^{2}\alpha_{0}+(\text{ad\
}\eta_{1})\alpha_{1}$ $\displaystyle={-1\choose 0}[1,1,0,0,0,0]+{-5/4\choose
0}[1,0,1,0,0,0]$ $\displaystyle+{3/8\choose 0}[0,1,1,0,0,0]+{1/4\choose
0}[0,1,0,1,0,0]$ $\displaystyle+{1\choose-1}[0,1,0,0,1,0]+{-3/8\choose
0}[0,0,1,1,0,0]$ $\displaystyle+{15/8\choose-5/8}[0,0,1,0,1,0]+{-9/8\choose
0}[0,0,1,0,0,1]$ $\displaystyle+{0\choose
3/2}[0,0,0,1,1,0]+{0\choose-1/4}[0,0,0,1,0,1]$ $\displaystyle+{0\choose
1}[0,0,0,0,1,1].$

Continuing the computations following Algorithm B we obtain

$\displaystyle\hat{\alpha}_{1}$ $\displaystyle=0;$
$\displaystyle\hat{\alpha}_{2}$
$\displaystyle=\left(-a^{(1)}_{01}a^{(1)}_{10}+a^{(1)}_{01}a^{(2)}_{10},-a^{(1)}_{01}a^{(2)}_{10}+a^{(2)}_{01}a^{(2)}_{10}\right)^{T};$
$\displaystyle\hat{\alpha}_{3}$
$\displaystyle=\left(0,a^{(1)}_{10}a^{(2)}_{02}a^{(2)}_{10}+2a^{(2)}_{02}(a^{(2)}_{10})^{2}\right)^{T};$
$\displaystyle\hat{\alpha}_{4}$
$\displaystyle=\left((a^{(1)}_{01})^{2}a^{(1)}_{10}a^{(2)}_{10}+a^{(1)}_{01}a^{(1)}_{10}a^{(2)}_{01}a^{(2)}_{10}-2(a^{(1)}_{01})^{2}(a^{(2)}_{10})^{2}\right.,$
$\displaystyle\left.-a^{(1)}_{01}a^{(1)}_{10}a^{(2)}_{01}a^{(2)}_{10}+2(a^{(1)}_{01})^{2}(a^{(2)}_{10})^{2}-a^{(2)}_{01}(a^{(2)}_{10})^{2}\right)^{T}.$

so, the normal form up to level 4 is

$\hat{\alpha}=\hat{\alpha}_{0}+\hat{\alpha}_{1}+\hat{\alpha}_{2}+\hat{\alpha}_{3}+\hat{\alpha}_{4}.$
(51)

Remark. From Algorithm A it is obvious that in the case of equation (18),
where $F(a,x)\in\mathcal{V}^{n}_{j}$ the normal form produced by Algorithm A
is the same as the usual Poincaré-Dulac normal form (produced e.g. algorithms
of [10, 16, 20]). Normal form (51) agrees with (16). However, computing up to
the seventh order we see that Algorithm B and the algorithms of [10, 16, 20]
produce normal forms which differ in terms of order 7.

A great advantage of Algorithm B is that it allows parallel computations of
terms of normal forms, that is, it is possible to compute one term of the
normal form at once based on the following statement.

######  Proposition 2.

If for some $\kappa\in\mathbb{N}_{0}^{\ell}$ we have $\langle
L(\kappa),\lambda\rangle=0$, then the coefficient $\alpha_{\kappa}$ in the
normal form $\hat{\alpha}$ is computed using only the terms
$\alpha_{\mu}a^{\mu}$ such that $a^{\mu}$ divides $a^{\kappa}$, i.e.
$\mu_{j}\leq\kappa_{j}$ for all $j=1,2,\dots,\ell$.

###### Proof.

The correctness of the statements follows from (43) and Algorithm B. ∎

Thus, in order to compute a coefficient $\alpha_{\kappa}$ of the normal form,
we first look for the set of $\mu$’s involved in the computation of
$\alpha_{\kappa}$ by means of Algorithm B. Denote this set $\omega_{\kappa}$.
The set $\omega_{\kappa}$ can be found using, for instance, the following
procedure:

Let $|\kappa|=s$. Set $p=1$, $\tau_{\kappa}(s)=\\{\kappa\\}.$

While $p<s$ do
set $\tau_{\kappa}(s-p)=\emptyset$;
for $\mu=(\mu_{1},\dots,\mu_{\ell}):$
for $i=1,\dots,\ell$: if $\mu_{i}-i\geq 0$ then
$\tau_{\kappa}(s-p)=\tau_{\kappa}(k)\cup\\{\mu\\}$;
set $p=p+1$.

The output of the procedure is the sequence of the sets $\tau_{\kappa}(i)$,
$i=1,\dots,s-1$, where $\tau_{\kappa}(i)$ is a subset of elements of level
$i$. Then

$\tau_{\kappa}=\cup_{i=1}^{s-1}\tau_{\kappa}(i)$

is the subset of $\mathbb{N}_{0}^{\ell}$ needed in the computation of
$\alpha_{\kappa}$ and, in order to compute the $\alpha_{\kappa}$, one just
uses Algorithm B, where the Lie brackets are computed with $\omega_{1}$ and
$\omega_{2}$ in (43) being subsets of $\tau_{\kappa}$.

Thus, based on Proposition 2 and the procedure described above, one can easily
parallelize normal form computations.

## Acknowledgments

The first author is supported by the Slovenian Research Agency (core research
program P1-0288) and the second author is supported by the Slovenian Research
Agency (core research program P1-0306).

We acknowledge COST (European Cooperation in Science and Technology) actions
CA15140 (Improving Applicability of Nature-Inspired Optimisation by Joining
Theory and Practice (ImAppNIO)) and IC1406 (High-Performance Modelling and
Simulation for Big Data Applications (cHiPSet)) and Tomas Bata University in
Zlín for accessing Wolfram Mathematica.

## References

  * [1] A. Algaba, E. Freire, E. Gamero, Characterizing and computing normal forms using Lie transforms, Dyn. Contin. Discrete Impuls. Syst. Ser. A: Math. Anal. 8 (2001) 449–475.
  * [2] Y. N. Bibikov, Local Theory of Nonlinear Analytic Ordinary Differential Equations. Lecture Notes in Math. Vol. 702. Springer–Verlag, New York, London, 1979.
  * [3] G. D. Birkhoff, Dynamical systems, Amer. Math. Soc. Colloq. Publ. 9, Amer. Math. Soc., Providence, R. I., 1927.
  * [4] A. D. Brjuno, _A Local Method of Nonlinear Analysis for Differential Equations_. Moscow: Nauka, 1979; _Local Methods in Nonlinear Differential Equations._ Translated from the Russian by William Hovingh and Courtney S. Coleman. Springer-Verlag, Berlin, 1989\.
  * [5] K.T. Chen, Equivalence and decomposition of vector fields about an elementary critical point, Amer. J. Math., 85 (1963) 693–722, https://doi.org/10.2307/2373115.
  * [6] S-N. Chow, C. Li, D. Wang, Normal forms and bifurcations of planar vector fields. Cambridge University Press, Cambridge, 1994, https://doi.org/10.1017/CBO9780511665639.
  * [7] H. Dulac, Recherches sur les points singuliers des equations différentielles. J. Ecole Polytechn. Ser. II, 9 (1904) 1–25.
  * [8] M. Han, P. Yu, Normal forms, Melnikov functions and bifurcations of limit cycles. Springer, London, 2012, https://doi.org/10.1007/978-1-4471-2918-9.
  * [9] A. M. Liapunov, _Stability of Motion_. With a contribution by V. A. Pliss. Translated by F. Abramovici and M. Shimshoni. Academic Press, New York, 1966.
  * [10] J. Murdock, Normal forms and unfolding for local dynamical systems. Springer Monographs in Mathematics, Springer-Verlag, New York, 2003, https://doi.org/10.1007/b97515.
  * [11] A. H. Nayfeh, Method of Normal Form. Wiley, New York, 1993.
  * [12] P. J. Olver, Applications of Lie Groups to Differential Equations, Springer-Verlag, New York-Berlin-Heidelberg, 1993, https://doi.org/10.1007/978-1-4684-0274-2.
  * [13] H. Poincaré. Mémoire sur les courbes définies par une équation différentielle, _J. Math. Pures et Appl._ (Sér. 3) 7 (1881) 375–422; (Sér. 3) 8 (1882) 251–296; (Sér. 4) 1 (1885) 167–244; (Sér. 4) 2 (1886) 151–217.
  * [14] V. G. Romanovskii, Calculation of Lyapunov numbers in the case of two pure imaginary roots. Differential Equations, 29 (1993), 782–784.
  * [15] V. G. Romanovski and D. S. Shafer. _The Center and Cyclicity Problems: A Computational Algebra Approach_. Birkhäuser, Boston, 2009, https://doi.org/10.1007/978-0-8176-4727-8.
  * [16] J. Sanders, Computational Aspects of Lie Group Representations and Related Topics, Proceedings of the 1990 Computational Algebra Seminar in CWI Tract 84, 129 (1991) 129–142.
  * [17] S. Sternberg, Infinite Lie groups and formal aspects of dynamical systems. J. Math. Mech., 10 No. 3 (1961) 451–474.
  * [18] L. Stolovitch, Progress in normal form theory. Nonlinearity 22 No. 7, (2009) R77–99, https://doi.org/:10.1088/0951-7715/22/7/R01.
  * [19] Y. Tian, P. Yu, An explicit recursive formula for computing the normal forms associated with semisimple cases, Commun. Nonlinear Sci. Numer. Simul. 19(7) (2014) 2294–2308, https://doi.org/10.1016/j.cnsns.2013.11.019.
  * [20] S. Walcher, On transformations into normal form, J. Math. Anal. Appl. 180 (1993) 617–632, https://doi.org/10.1006/jmaa.1993.1420.
  * [21] P. Yu, A simple and efficient method for computing center manifold and normal forms associated with semisimple cases. Dyn. Contin. Discrete Impuls. Syst. Ser. B: Appl Algorithm 10 (2003) 273–86.