# Liouvillian integrability of polynomial differential systems

Xiang Zhang  Department of Mathematics, and MOE-LSC, Shanghai Jiao Tong
University, Shanghai 200240, People’s Republic of China<EMAIL_ADDRESS>
###### Abstract.

M.F. Singer [Liouvillian first integrals of differential equations, Trans.
Amer. Math. Soc. 333 (1992), 673–688] proved the equivalence between
Liouvillian integrability and Darboux integrability for two dimensional
polynomial differential systems. In this paper we will extend Singer’s result
to any finite dimensional polynomial differential systems. We prove that if an
$n$–dimensional polynomial differential system has $n-1$ functionally
independent Darboux Jacobian multiplier then it has $n-1$ functionally
independent Liouvillian first integrals. Conversely if the system is
Liouvillian integrable then it has a Darboux Jacobian multiplier.

###### Key words and phrases:

Liouville integrability; Darboux integrability; Jacobian multiplier; Galois
group.
To appear in Transactions of the American Mathematical Society

###### 2010 Mathematics Subject Classification:

34A34, 37C10, 34C14, 37G05.

##  1\. Background and statement of the main results

The theory of integrability for differential systems is classic and it is
useful in the study of dynamics of differential system. Integrability has
different definitions in different fields. Here we mainly concern the
algebraic aspects of integrability for polynomial differential systems, which
involves analysis, algebraic geometry, the field extension and so on. For
further information on this subject, we refer readers to Daboux [7, 8],
Jouanolou [12], Prelle and Singer [21], Singer [23], Schlomiuk [22], Llibre
[14], Dumortier and Llibre et al [9], Christopher et al [3, 6] and Llibre and
Zhang [16, 18, 19].

Darboux theory of integrability was established by Darboux [7, 8] in 1878 for
polynomial differential systems of degree $n$ by using the invariant algebraic
curves (resp. surfaces or hypersurfaces) in dimension 2 (resp. 3 or $n>3$).
Jouanolou [12] in 1979 extended the Darboux’s theory to construct rational
first integrals with the help of algebraic geometry. An elementary proof of
Jouanolou’s result was provided respectively by Christopher and Llibre [5] in
2000 for two dimensional differential systems and by Llibre and Zhang [17] in
2010 for any finite dimensional differential systems. On further extensions to
Darboux theory of integrability, Christopher, Llibre and Pereira [6] in 2007
took into account not only the number of invariant algebraic curves but also
their multiplicities for two dimensional differential systems. Llibre and
Zhang [16] further extended Christopher et al’s result in [6] to any finite
dimensional differential systems, where there are some deep characterizations
on the number of exponential factors and the multiplier of invariant algebraic
hypersurfaces. Darboux theory of integrability has important applications in
the center–focus problem, dynamical analysis and so on, see for instance [4,
14, 22, 24] and the reference therein.

Darboux theory of integrability has a nice extension to Weierstrass
integrability, see e.g. [10], where they used Weierstrass polynomials to
replace the usual polynomials. By definition the former include the latter as
a special one. In [1] Blázquez–Sanz and Pantazi provided a new approach to
study the Darboux integrability of polynomial differential systems of degree
$m$, where they replaced the dimension of $\mathbb{C}_{m-1}[x]$ which the
cofactors of Darboux polynomials and exponential factors are located in by the
rank of a matrix associated to these cofactors. Here $\mathbb{C}_{m-1}[x]$
denotes the linear space formed by polynomials in $x\in\mathbb{C}^{n}$ of
degree no more than $m-1$. Recently Darboux theory of integrability was also
successfully extended to nonautonomous differential systems which are
polynomial ones in space variables with coefficients the smooth functions of
the time, see e.g. [15, 11], where they extended the notion of invariant
algebraic hypersurfaces in the phase space to polynomial invariant
hypersurfaces in the extended space including the time.

Prelle and Singer [21] in 1983 proved that if a polynomial differential system
has an elementary first integral then it has a first integral of a very simple
form. As a corollary of their results, one gets that if a planar polynomial
differential system has an elementary first integral, then it has an
integrating factor of the form $f_{1}^{m_{1}}\ldots f_{p}^{m_{p}}$ with
$f_{i}\in\mathbb{C}[x,y]$ and $m_{i}\in\mathbb{Z}$. This shows the equivalence
between the existence of elementary first integrals and the Darboux
integrating factors for planar polynomial differential systems.

Singer [23] in 1992 proved that a planar polynomial differential system has a
Liouvillian first integral if and only if it has an integrating factor of the
form

$R(x,y)=\exp\left(\int U(x,y)dx+V(x,y)dy\right),$

where $U(x,y),V(x,y)$ are rational functions in $x,y$. For a simple proof to
Singer’s result, see [3] and [4, Theorem 3.2].

In this paper we will extend Singer’s result to any finite dimensional
polynomial differential systems.

Consider polynomial differential systems

(1.1) $\dot{x}=P(x),\qquad x\in\mathbb{C}^{n},$

where $P(x)=(P_{1}(x),\ldots,P_{n}(x))$ are vector–valued polynomial
functions. We call $m:=\max\\{\mbox{deg}P_{1},\ldots,\mbox{deg}P_{n}\\}$ the
degree of polynomial differential systems (1.1). In what follows we also use

$\mathcal{X}_{P}=P_{1}(x)\frac{\partial}{\partial
x_{1}}+\ldots+P_{n}(x)\frac{\partial}{\partial x_{n}},$

to represent the vector field associated to system (1.1). For simplifying
notations, in what follows we denote ${\partial}/{\partial x_{i}}$ by
$\partial_{i}$ and $(\partial_{1},\ldots,\partial_{n})$ by $\partial$.

Denote by $\mathbb{C}[x]$ the ring of polynomials in $x$ with coefficients in
$\mathbb{C}$. A polynomial $f(x)\in\mathbb{C}[x]$ is called a Darboux
polynomial of $\mathcal{X}_{P}$ if there exists a $k(x)\in\mathbb{C}[x]$ such
that

$\mathcal{X}_{P}(f)(x)=k(x)\,f(x),\qquad x\in\mathbb{C}^{n}.$

The polynomial $k$ is called cofactor of $f$. A function of the form

$\exp\left(\frac{g}{h}\right)f_{1}^{l_{1}}\ldots f_{r}^{l_{r}},\quad\mbox{
with }g,\,h,\,f_{i}\in\mathbb{C}[x],\,\,\,l_{i}\in\mathbb{C},\,\,i=1,\ldots,r$

is called a Darboux function. For a Darboux function, we always require that
its factors $f_{i}$ are irreducible and relatively different, and $g,h$ are
relative coprime.

A Darboux first integral of (1.1) is a Darboux function and it is a first
integral of (1.1). Note that a first integral is not necessary to be defined
in the full space but in a full Lebesgue measure subset of $\mathbb{C}^{n}$.
System (1.1) is Darboux integrable if it has $n-1$ functionally independent
Darboux first integrals.

A smooth function $J(x)$ is a Jacobian multiplier of system (1.1) if

$\partial_{1}(JP_{1})+\ldots+\partial_{n}(JP_{n})=0.$

A Darboux Jacobian multiplier of system (1.1) is a Jacobian multiplier of the
system and it is a Darboux function. For planar polynomial differential
systems, a Jacobian multiplier is usually called an integrating factor. If a
planar polynomial differential system has a Darboux integrating factor, it is
also called Darboux integrable.

For stating our results we recall the definition of Liouvillian functions.

A differential field  $(K,\,\Delta)$ consists of the field $K$ and the set
$\Delta$ of commutative derivatives defined on $K$. In this paper all
mentioned fields have characteristic $0$.

A differential field extension of a differential field $(K,\,\Delta)$ is a
differential field $(L,\,\Delta^{\prime})$ with the properties that $K\subset
L$ and for $\forall\,\delta^{\prime}\in\Delta^{\prime}$ we have
$\left.\delta^{\prime}\right|_{K}\in\Delta$. Because of the relation between
the derivatives of differential field $(K,\,\Delta)$ and its field extension
$(L,\,\Delta^{\prime})$, we also use $\Delta$ to represent $\Delta^{\prime}$.
For simplifying notations we also use $L/K$ to denote the differential field
extension $(L,\Delta)$ of $(K,\Delta)$.

For a field extension $L/K$,

  * •

$\alpha\in L$ is called

    * –

an algebraic element of $K$, if there exists a polynomial with coefficients in
$K$ such that $F(\alpha)=0$.

    * –

a transcendental element of $K$, if $\alpha$ is not an algebraic element over
$K$.

  * •

If each element of $L$ is algebraic over $K$, we call  $L/K$ an algebraic
extension of field $K$.

  * •

$L$ can be considered as a vector space over $K$: the elements of $L$ are
treated as vectors, and elements of $K$ are treated as scalars, and the
summation of vectors is that of elements of field and the product of elements
of $L$ and $K$ is that of elements of field $L$.

    * –

The dimension of this vector space is called degree of this differential field
extension, denoted by $[L:K]$.

  * •

If $[L:K]\in\mathbb{N}$, we call $L/K$ finite field extension.

  * •

Let $S\subset L$,

    * –

$K(S)$ denotes the minimal subfield of $L$ including $K$ and $S$.

    * –

If $S$ contains only one element, we call  $K(S)$ the minimal field extension
of $K$.

A differential field extension $L/K$ is Liouvillian, if this differential
field extension can be written in the tower form

$K=K_{0}\subset K_{1}\subset\ldots\subset K_{r}=L,$

such that

  * $(a)$

$K_{i+1}$ is a finite algebraic extension of $K_{i}$, or

  * $(b)$

$K_{i+1}=K_{i}(t)$, where $t$ is a transcendental element of $K_{i}$
satisfying: for each $\delta\in\Delta$, $\dfrac{\delta t}{t}\in K_{i}$, or

  * $(c)$

$K_{i+1}=K_{i}(t)$, where $t$ is a transcendental element of $K_{i}$
satisfying: for $\delta\in\Delta$, $\delta t\in K_{i}$.

A Liouvillian first integral of (1.1) is a Liouvillian function and is a first
integral of (1.1). System (1.1) is Liouvillian integrable if it has $n-1$
functionally independent Liouvillian first integrals.

Now we can state our main results. The first one characterizes the existence
of Liouvillian first integrals via Darboux Jacobian multipliers.

######  Theorem 1.1.

If polynomial differential system (1.1) has $n-1$ functionally independent
Darboux Jacobian multipliers, then they have $n-1$ functionally independent
Liouvillian first integrals.

The next one shows that Liouvillian integrability implies the existence of
Darboux Jacobian multipliers.

######  Theorem 1.2.

If system (1.1) is Liouvillian integrable, i.e. it has $n-1$ functionally
independent Liouvillian first integrals, then the system has a Darboux
Jacobian multiplier.

In the rest of this paper we will prove our main results.

##  2\. Proof of the main results

###  2.1. Proof of Theorem 1.1

Let $J_{1}(x),\,\ldots,\,J_{n-1}(x)$ be $n-1$ functionally independent Darboux
Jacobian multipliers of system (1.1). Then we have

(2.1) $\mathcal{X}_{P}(J_{l})=-J_{l}\,\mbox{div}P,\qquad l=1,\ldots,n-1,$

where $\mbox{div}P=\partial_{1}P_{1}(x)+\ldots+\partial_{n}P_{n}(x)$ is the
divergence of the vector fields $P(x)$. Recall that $\partial_{i}P_{i}$
denotes the partial derivative of the function $P_{i}$ with respect to
$x_{i}$.

From the definition of Darboux functions and some direct calculations we get
that

$\frac{\partial_{i}J_{l}}{J_{l}}\in\mathbb{C}(x),\qquad
l\in\\{1,\ldots,n-1\\},\,\,\,i\in\\{1,\ldots,n\\}.$

Recall that $\mathbb{C}(x)$ is the field of rational functions in $x$. So it
follows from the condition $(b)$ of Liouvillian extension of field that
$J_{l}$ for $l=1,\ldots,n-1$ are Liouvillian functions. Furthermore, some easy
calculations show that

$\frac{J_{l}}{J_{k}},\qquad\mbox{for }\quad 1\leq l\neq k\leq n-1,$

are non–trivial Liouvillian first integrals of the vector field
$\mathcal{X}_{P}$, i.e. $\frac{J_{l}}{J_{k}}$ is not a constant and
$\mathcal{X}_{P}\left(\frac{J_{l}}{J_{k}}\right)\equiv 0$.

We claim that

$\frac{J_{1}}{J_{n-1}},\,\ldots,\,\frac{J_{n-2}}{J_{n-1}},$

are functionally independent. Indeed, assume that

$c_{1}\partial\left(\frac{J_{1}}{J_{n-1}}\right)+\ldots+c_{n-2}\partial\left(\frac{J_{n-2}}{J_{n-1}}\right)=0.$

Since

$\partial\left(\frac{J_{l}}{J_{n-1}}\right)=\frac{J_{n-1}\,\partial
J_{l}-J_{l}\partial J_{n-1}}{J_{n-1}^{2}},$

we have

$\displaystyle c_{1}J_{n-1}\partial J_{1}+\ldots+c_{n-2}J_{n-1}\partial
J_{n-2}$
$\displaystyle\qquad\qquad\,\,-(c_{1}J_{1}+\ldots+c_{n-2}J_{n-2})\partial
J_{n-1}=0.$

So by the functional independence of $J_{1},\ldots,J_{n-1}$ we must have

$c_{1}J_{n-1}=\ldots=c_{n-2}J_{n-1}=c_{1}J_{1}+\ldots+c_{n-2}J_{n-2}=0,$

in a full Lebesgure measure subset of $\mathbb{C}^{n}$. Consequently

$c_{1}=\ldots=c_{n-2}=0,$

in a full Lebesgure measure subset of $\mathbb{C}^{n}$. This proves the claim.

Using the last claim, we assume without loss of generality that

$\displaystyle y_{i}$ $\displaystyle=$
$\displaystyle\frac{J_{i}}{J_{n-1}},\qquad i=1,\ldots,n-2,$ $\displaystyle
y_{n-1}$ $\displaystyle=$ $\displaystyle x_{n-1},$ $\displaystyle y_{n}$
$\displaystyle=$ $\displaystyle x_{n},$

are invertible, at least in some full Lebesgue measure subset $\Omega$ of
$\mathbb{C}^{n}$. Denote by $y=G(x)$ this last transformation. Then under it
the differential system (1.1) is equivalent to

$\displaystyle\dot{y}_{i}$ $\displaystyle=$ $\displaystyle 0,\quad\qquad
i=1,\ldots,n-2,$ (2.2) $\displaystyle\dot{y}_{n-1}$ $\displaystyle=$
$\displaystyle P_{n-1}\circ G^{-1}(y),$ $\displaystyle\dot{y}_{n}$
$\displaystyle=$ $\displaystyle P_{n}\circ G^{-1}(y).$

Clearly system (2.1) has the first integrals $I_{i}(y)=y_{i}$,
$i=1,\ldots,n-2$. In addition, we can prove that system (2.1) has the Jacobian
multiplier

$M(y)=J_{n-1}\circ G^{-1}(y)D_{y}G^{-1}(y),$

where $D_{y}G^{-1}(y)$ denotes the Jacobian matrix of $G^{-1}$ with respect to
$y$. This shows that the two dimensional differential system

$\displaystyle\dot{y}_{n-1}$ $\displaystyle=$ $\displaystyle P_{n-1}\circ
G^{-1}(I_{1},\ldots,I_{n-2},y_{n-1},y_{n})=:g_{n-1}(y_{n-1},y_{n}),$
$\displaystyle\dot{y}_{n}$ $\displaystyle=$ $\displaystyle P_{n}\circ
G^{-1}(I_{1},\ldots,I_{n-2},y_{n-1},y_{n})=:g_{n}(y_{n-1},y_{n}),$

has the integrating factor

$V(y_{n-1},y_{n})=J_{n-1}\circ
G^{-1}(y)D_{y}G^{-1}(I_{1},\ldots,I_{n-2},y_{n-1},y_{n}),$

where we take $I_{1},\ldots,I_{n-2}$ as constants. Hence this last two
dimensional differential system has the first integral.

$I_{n-1}(y_{n-1},y_{n})=\int Vg_{n}dy_{n-1}-Vg_{n-1}dy_{n}.$

Obviously, $I_{n-1}$ is functionally independent of $I_{1},\dots,I_{n-2}$,
because the latter are independent of $y_{n-1}$ and $y_{n}$.

Next we prove that $I_{n-1}$ is a Liouvillian function. First we prove that
$G^{-1}(x)$ is a Liouvillian function. Indeed, by

$G^{-1}\circ G(x)=x,$

we have

$D_{y}G^{-1}(G(x))D_{x}G(x)=E.$

where $E$ is the $n$–dimensional unit matrix. Since $G$ is Liouvillian, and so
is $(D_{x}G(x))^{-1}$. Hence we have

$D_{y}G^{-1}(G(x))=(D_{x}G(x))^{-1},$

is a Liouvillian function. This shows that $G^{-1}(y)$ is a Liouvillian
function. Furthermore it follows from the above construction that
$g_{n-1},g_{n}$ and $V$ are Liouvillian functions. This proves that $I_{n-1}$
is a Liouvillian function.

Applying the transformation $y=G(x)$ to $I_{1}(y),\ldots,I_{n-1}(y)$, we get
$n-1$ functionally independent Liouvillian first integrals

$H_{1}(x):=I_{1}\circ G(x),\quad\ldots,\quad H_{n-1}(x)=I_{n-1}\circ G(x),$

of differential system (1.1). We complete the proof of the theorem.

###  2.2. Proof of Theorem 1.2

For proving Theorem 1.2 and readers’s convenience, we recall some notions.

Given a field $K$,

  * •

A  separating field of a polynomial $p(x)$ over $K$ is a minimal field
extension of $K$ such that $p(x)$ can be decomposed into product of linear
factors over this field extension, i.e. $p(x)=\prod(x-a_{i})$, $a_{i}\in L$,
$L/K$ is the minimal field extension such that this decomposition can happen.

  * •

We say that an algebraic field extension $L/K$ of $K$ is normal, if $L$ is a
separating field of polynomials in $K[x]$.

  * •

The normal closure of an algebraic field extension $L/K$ is a field extension
$\overline{L}$ of $L$ such that $\overline{L}/K$ is normal, and $\overline{L}$
is the minimal field extension satisfying this property.

  * •

Field automorphism over field $K$ is a bijective map $\varphi:\,\,K\rightarrow
K$ which keeps the algebraic properties of $K$, i.e. $\varphi$ satisfies that
$\varphi(0_{K})=0_{K}$, $\varphi(1_{K})=1_{K}$,
$\varphi(a+b)=\varphi(a)+\varphi(b)$ and $\varphi(ab)=\varphi(a)\varphi(b)$.

  * •

The set of all field automorphisms over field $K$ fixing elements of a
subfield $K^{\prime}\subset K$ forms a group under the composition of maps.
This group is called Galois group.

  * •

The order of a group is the number of elements of a group $G$, denoted by
$|G|$.

Now we can prove Theorem 1.2. It will follows from the following lemmas.

######  Lemma 2.1.

If system (1.1) is Liouvillian integrable, then it has a Jacobian multiplier
of the form

$J=\exp\left(\int U_{1}dx_{1}+\ldots+U_{n}dx_{n}\right),$

with $U_{i}\in\mathbb{C}(x)$, $i=1,\ldots,n$, and

$\partial_{i}U_{j}=\partial_{j}U_{i},\qquad 1\leq j<i\leq n.$

###### Proof.

Assume that system (1.1) has the functionally independent Liouvillian first
integrals $H_{1},\ldots,H_{n-1}$, which are defined in a full Lebesgue measure
subset $\Omega$ of $\mathbb{C}^{n}$. By definition of first integrals we have

(2.3) $\mathcal{X}_{P}(H_{i})(x)\equiv 0,\quad x\in\Omega,\qquad
i=1,\ldots,n-1.$

From independence of $H_{1},\ldots,H_{n-1}$ we can assume without loss of
generality that

$\Gamma:=\det\left(\partial_{1}\mathcal{H},\cdots,\partial_{n-1}\mathcal{H}\right)\neq
0,\qquad x\in\Omega,$

with

$\mathcal{H}:=(H_{1},\ldots,H_{n-1})^{T},$

where $T$ denotes the transpose of a matrix, and

$\partial_{i}\mathcal{H}:=(\partial_{i}H_{1},\ldots,\partial_{i}H_{n-1})^{T},\qquad
i=1,\ldots,n.$

Set for $i=1,\ldots,n-1$

$\Gamma_{i}:=\det\left(\partial_{1}\mathcal{H},\cdots,\partial_{i-1}\mathcal{H},\partial_{n}\mathcal{H},\partial_{i+1}\mathcal{H},\cdots,\partial_{n-1}\mathcal{H}\right).$

Clearly $\Gamma$ and $\Gamma_{i}$, $i=1,\ldots,n$, are Liouvillian functions.

By the Cramer’s rule we get from (2.3) that

$P_{i}(x)=-\frac{\Gamma_{i}}{\Gamma}P_{n}(x),\quad i=1,\ldots,n-1.$

Hence we have

$\Gamma(P_{1}(x),\ldots,P_{n-1})=-(\Gamma_{1},\ldots,\Gamma_{n-1})P_{n}(x).$

Since $P_{1},\ldots,P_{n-1},P_{n}$ are relative coprime, and $\Gamma$ and
$\Gamma_{i}$ are Liouvillian functions, so there exists a Liouvillian function
$h(x)$ such that

$h(x)\Gamma=P_{n}(x).$

Consequently we have

$P_{i}(x)=-h(x)\Gamma_{i},\quad i=1,\ldots,n-1.$

Set

(2.4) $A_{i}:=\frac{\partial_{i}h}{h},\quad i=1,\ldots,n;\qquad
A:=(A_{1},\ldots,A_{n}).$

Then $A_{i}$ is Liouvillian for each $i\in\\{1,\ldots,n\\}$.

We claim that

(2.5) $\displaystyle\partial_{i}A_{j}$ $\displaystyle=$
$\displaystyle\partial_{j}A_{i},\qquad 1\leq j<i\leq n,$ (2.6)
$\displaystyle\langle A,P\rangle$ $\displaystyle:=$ $\displaystyle
A_{1}P_{1}+\ldots+A_{n}P_{n}=\mbox{div}P.$

The equality (2.5) can be proved by direct calculations via the fact

$\partial_{j}\partial_{i}h=\partial_{i}\partial_{j}h\qquad\mbox{ for all
}\,\,\,i,j\in\\{i,\ldots,n\\}.$

For proving (2.6), some computations show that

$\displaystyle\mbox{div}P$ $\displaystyle=$
$\displaystyle\partial_{1}(-h\Gamma_{1})+\ldots+\partial_{n-1}(-h\Gamma_{n-1})+\partial_{n}(h\Gamma)$
$\displaystyle=$ $\displaystyle
A_{1}P_{1}+\ldots+A_{n}P_{n}+h(\partial_{n}\Gamma-\partial_{1}\Gamma_{1}-\ldots-\partial_{n-1}\Gamma_{n-1}).$

Next we only need to prove that

$\partial_{n}\Gamma-\partial_{1}\Gamma_{1}-\ldots-\partial_{n-1}\Gamma_{n-1}=0.$

Since

$\displaystyle\partial_{n}\Gamma$ $\displaystyle=$
$\displaystyle\sum\limits_{i=1}\limits^{n-1}\det\left(\partial_{1}\mathcal{H},\ldots,\partial_{i-1}\mathcal{H},\partial_{n}\partial_{i}\mathcal{H},\partial_{i+1}\mathcal{H},\ldots,\partial_{n-1}\mathcal{H}\right),$
$\displaystyle\partial_{i}\Gamma_{i}$ $\displaystyle=$
$\displaystyle\det\left(\partial_{1}\mathcal{H},\ldots,\partial_{i-1}\mathcal{H},\partial_{i}\partial_{n}\mathcal{H},\partial_{i+1}\mathcal{H},\ldots,\partial_{n-1}\mathcal{H}\right)$
$\displaystyle+\sum\limits_{j=1,j\neq
i}\limits^{n-1}\det\left(\partial_{1}\mathcal{H},\ldots,\partial_{j-1}\mathcal{H},\partial_{i}\partial_{j}\mathcal{H},\partial_{j+1}\mathcal{H},\ldots,\partial_{i-1}\mathcal{H},\right.$
$\displaystyle\qquad\qquad\qquad\qquad\qquad\left.\partial_{n}\mathcal{H},\partial_{i+1}\mathcal{H},\ldots,\partial_{n-1}\mathcal{H}\right),$

we have

$\displaystyle\partial_{1}\Gamma_{1}+\ldots+\partial_{n-1}\Gamma_{n-1}-\partial_{n}\Gamma$
$\displaystyle=$
$\displaystyle\sum\limits_{i=1}\limits^{n-1}\sum\limits_{j=1,j\neq
i}\limits^{n-1}\det\left(\partial_{1}\mathcal{H},\ldots,\partial_{j-1}\mathcal{H},\partial_{i}\partial_{j}\mathcal{H},\partial_{j+1}\mathcal{H},\ldots,\partial_{i-1}\mathcal{H},\right.$
$\displaystyle\qquad\qquad\qquad\qquad\qquad\quad\left.\partial_{n}\mathcal{H},\partial_{i+1}\mathcal{H},\ldots,\partial_{n-1}\mathcal{H}\right)$
$\displaystyle=$ $\displaystyle 0,$

where in the last equality we have used the fact that

$\displaystyle\det\left(\partial_{1}\mathcal{H},\ldots,\partial_{j-1}\mathcal{H},\partial_{i}\partial_{j}\mathcal{H},\partial_{j+1}\mathcal{H},\ldots,\partial_{i-1}\mathcal{H},\right.$
$\displaystyle\qquad\qquad\qquad\qquad\qquad\quad\quad\left.\partial_{n}\mathcal{H},\partial_{i+1}\mathcal{H},\ldots,\partial_{n-1}\mathcal{H}\right)$
$\displaystyle\qquad+\det\left(\partial_{1}\mathcal{H},\ldots,\partial_{j-1}\mathcal{H},\partial_{n}\mathcal{H},\partial_{j+1}\mathcal{H},\ldots,\partial_{i-1}\mathcal{H},\right.$
$\displaystyle\qquad\qquad\qquad\qquad\qquad\quad\quad\left.\partial_{j}\partial_{i}\mathcal{H},\partial_{i+1}\mathcal{H},\ldots,\partial_{n-1}\mathcal{H}\right)=0.$

This proves the claim.

Set

$U_{i}=-A_{i},\qquad i=1,\ldots,n,$

with $A_{i}$ defined in (2.4). By (2.5) and the Stokes’s theorem (see e.g. [2,
p.3]), we get that

(2.7) $J=\exp\left(\int U_{1}dx_{1}+\ldots+U_{n}dx_{n}\right),$

is well defined in any connected subsets where $U_{1},\ldots,U_{n}$ are
defined. Furthermore we can check that $J$ in (2.7) is a Jacobian multiplier
of system (1.1) if and only if (2.6) holds. So, in what follows we only need
to prove that there exist rational functions $B_{i}$ for $i=1,\ldots,n$
instead of $A_{i}$ such that the equalities (2.5) and (2.6) hold.

According to the Liouvillian extension of field in the tower form, all the
$A_{i}$ belong to some tower for $i=1,\ldots,n$. We distinguish three
different cases according to the definition of the tower.

$(a)$ $A_{i}\in K_{l+1}$, $i=1,\ldots,n$, and $K_{l+1}$ is an algebraic
extension of the field $K_{l}$. We will prove that there exist $B_{i}\in
K_{l}$ instead of $A_{i}$ for $i=1,\ldots,n$ such that (2.5) and (2.6) hold.

Let $\overline{K}_{l+1}$ be the normal closure of $K_{l+1}$, and $\mathcal{G}$
be the Galois group formed by the automorphisms of $\overline{K}_{l+1}$ fixing
$K_{l}$. Then it follows from a result of Artin (see Lang [13, Theorem 1.1])
that $\mathcal{G}$ is of finite order, and denote by $N=|\mathcal{G}|$ the
order of the group. Note that $N\leq[K_{l+1}:K_{l}]$, the degree of the
algebraic extension of the field.

Since $P\in\mathbb{C}(x)^{n}$ and $\mathbb{C}(x)\subset K_{l}$, we get from
(2.5) and (2.6) that

(2.8)
$\displaystyle\begin{array}[]{rcl}\partial_{i}g(A_{j})&=&\partial_{j}g(A_{i}),\qquad\forall
g\,\in\mathcal{G},\\\
\left\langle\sum\limits_{g\in\mathcal{G}}g(A),P\right\rangle&=&\sum\limits_{g\in\mathcal{G}}g(A_{1})P_{1}+\ldots+\sum\limits_{g\in\mathcal{G}}g(A_{n})P_{n}=N\mbox{div}P,\end{array}$

where in the second equality we have used the fact that $g\in\mathcal{G}$
fixes $K_{l}$. Set

$B_{i}=\frac{1}{N}\sum\limits_{g\in\mathcal{G}}g(A_{i}),\qquad i=1,\ldots,n.$

Then $B_{i}\in K_{l}$ for $i=1,\ldots,n$, because all $B_{i}$ are fixed under
the action of all elements of the Galois group $\mathcal{G}$. Furthermore, we
get from (2.8) that (2.5) and (2.6) hold for $B_{i}$ instead of $A_{i}$ for
$i=1,\ldots,n$.

$(b)$ Assume $K_{l+1}=K_{l}(t)$ with $t$ a transcendental element over $K_{l}$
and $\partial_{i}t/t\in K_{l}$ for $i=1,\ldots,n$.

Since $A_{i}\in K_{l}(t)$ for $i=1,\ldots,n$, we can assume without loss of
generality that

$A_{i}=a_{i}(t)\in K_{l}(t),\qquad i=1,\ldots,n.$

Expanding $a_{j}(t)$ into Laurent series in $t$ gives

(2.9)
$a_{j}(t)=a_{0}^{(j)}+\sum\limits_{s\in\mathbb{Z}\setminus\\{0\\}}a_{s}^{(j)}t^{s},\qquad
j=1,\ldots,n,$

with $a_{s}^{(j)}\in K_{l}$ for $j=1,\dots,n$ and all $s$. Then we have

(2.10)
$\partial_{i}A_{j}=\partial_{i}a_{0}^{(j)}+\sum\limits_{s\in\mathbb{Z}\setminus\\{0\\}}\left(\partial_{i}a_{s}^{(j)}+sa_{s}^{(j)}p_{i}\right)t^{s},$

where $p_{i}\in K_{l}$ satisfying $\partial_{i}t/t=p_{i}\in K_{l}$ for
$i=1,\ldots,n$. Set

$B_{i}=a_{0}^{(i)},\qquad i=1,\ldots,n,$

we have $B_{i}\in K_{l}$. Substituting (2.9) and (2.10) into (2.5) and (2.6),
and equating the coefficients of $t^{0}$, we get that $B_{i}$ for
$i=1,\ldots,n$ satisfy (2.5) and (2.6) instead of $A_{i}$.

$(c)$ Assume that $K_{l+1}=K_{l}(t)$ with $t$ a transcendental element over
$K_{l}$ and $\partial_{i}t\in K_{l}$ for all $i\in\\{1,\ldots,n\\}$.

Similar to $(b)$ we set

$A_{j}=a_{j}(t)\in K_{l}(t),\qquad j=1,\ldots,n.$

Now the Laurent expansion in $t$ does not work, we choose the Laurent
expansion in $1/t$ of $a_{j}(t)$. Since $A_{i}\in K_{l}(t)$, we get from its
construction (2.4), i.e. $A_{i}=\partial_{i}h/h$, that the degree of numerator
in $t$ of $a_{j}(t)$ is less than or equal to the degree of its denominator.
Write $a_{j}(t)$ as

$\displaystyle a_{j}(t)$ $\displaystyle=$
$\displaystyle\frac{a_{j0}+a_{j1}t+\ldots+a_{jk}t^{k}}{b_{j0}+b_{j1}t+\ldots+b_{jl}t^{l}}$
$\displaystyle=$
$\displaystyle\frac{a_{j0}t^{-l}+a_{j1}t^{-(l-1)}+\ldots+a_{jk}t^{-(l-k)}}{b_{j0}t^{-l}+b_{j1}t^{-(l-1)}+\ldots+b_{jl}}.$

Since $l\geq k$, so the Laurent expansion in $t^{-1}$ of $a_{j}(t)$ has the
form

(2.11) $a_{j}(t)=\sum\limits_{s=-\infty}\limits^{0}a_{s}^{(j)}t^{s},\qquad
j=1,\ldots,n,$

with $a_{s}^{(j)}\in K_{l}$ for $j=1,\dots,n$ and all $s$. Direct calculation
shows that

(2.12)
$\partial_{i}A_{j}=\sum\limits_{s=-\infty}\limits^{0}\left(\partial_{i}a_{s-1}^{(j)}+sa_{s}^{(j)}q_{i}\right)t^{s-1}+\partial_{i}a_{0}^{(j)},$

where $q_{i}\in K_{l}$ satisfying $q_{i}=\partial_{i}t\in K_{l}$ for
$i=1,\ldots,n$.

Set

$B_{i}=a_{0}^{(i)},\qquad i=1,\ldots,n.$

Then $B_{i}\in K_{l}$. Using the expansions (2.11) and (2.12) we get from
(2.5) and (2.6) that

$\displaystyle\partial_{i}B_{j}$ $\displaystyle=$
$\displaystyle\partial_{j}B_{i},\qquad 1\leq j<i\leq n,$
$\displaystyle\mbox{div}P$ $\displaystyle=$ $\displaystyle
B_{1}P_{1}+\ldots+B_{n}P_{n}.$

Of course if all $B_{i}=0$, then $\mbox{div}P=0$, and so $J=1$ is a Jacobian
multiplier.

Summarizing the cases $(a)$, $(b)$ and $(c)$, and combining the definition of
Liouvillian functions in the tower form, by induction we get that there exist
$U_{1},\ldots,U_{n}\in K_{0}=\mathbb{C}(x)$ for which (2.5) and (2.6) hold
instead of $A_{1},\ldots,A_{n}$. We complete the proof of the lemma. ∎

The next result shows that the existence of Jacobian multipliers of the form
given in Lemma 2.1 implies the existence of Darboux Jacobian multipliers.

######  Lemma 2.2.

If polynomial differential system (1.1) has a Jacobian multiplier

$J=\exp\left(\int U_{1}dx_{1}+\ldots+U_{n}dx_{n}\right),$

with $U_{i}\in\mathbb{C}(x)$, and $\partial_{j}U_{i}=\partial_{i}U_{j}$ for
$1\leq i,j\leq n$, then it has a Darboux Jacobian multiplier

$\exp\left(\frac{g}{h}\right)\prod\limits_{i}f_{i}^{l_{i}},$

where $g,h,f_{i}\in\mathbb{C}[x,y]$, $l_{i}\in\mathbb{C}$.

###### Proof.

Since $U_{1},\ldots,U_{n}\in\mathbb{C}(x)$, we treat their numerators and
denominators as polynomials in $x_{1}$ with coefficients in
$\mathbb{C}[x_{2},\ldots,x_{n}]$. Let $K$ be the minimal normal algebraic
field extension of $\mathbb{C}(x_{2},\ldots,x_{n})$ such that it is the
separating field of the numerators and denominators of $U_{1},\ldots,U_{n}$.

By the properties on normal algebraic field extension, the rational functions
$U_{1},\ldots,U_{n}$ over $K$ can be expanded in

(2.13)
$U_{k}(x)=\sum\limits_{i=1}\limits^{r}\sum\limits_{j=1}\limits^{m}\frac{\alpha_{ij}^{(k)}}{(x_{1}-\beta_{i})^{j}}+\sum\limits_{i=0}\limits^{p}\xi_{i}^{(k)}x_{1}^{i},\qquad
k=1,\ldots,n,$

where $\alpha_{ij}^{(k)},\beta_{i},\xi_{i}^{(k)}\in K$, and parts of them can
be zero. Direct calculations show that for $l\in\\{2,\ldots,n\\}$

$\displaystyle\partial_{l}U_{1}$ $\displaystyle=$
$\displaystyle\sum\limits_{i=1}\limits^{r}\sum\limits_{j=1}\limits^{m}\left(\frac{\partial_{l}\alpha_{ij}^{(1)}}{(x_{1}-\beta_{i})^{j}}+\frac{j\alpha_{ij}^{(1)}\partial_{l}\beta_{i}}{(x_{1}-\beta_{i})^{j+1}}\right)+\sum\limits_{i=0}\limits^{p}\partial_{l}\xi_{i}^{(1)}x^{i},$
$\displaystyle\partial_{1}U_{l}$ $\displaystyle=$
$\displaystyle\sum\limits_{i=1}\limits^{r}\sum\limits_{j=1}\limits^{m}\frac{-j\alpha_{ij}^{(l)}}{(x_{1}-\beta_{i})^{j+1}}+\sum\limits_{i=0}\limits^{p}i\xi_{i}^{(l)}x_{1}^{i-1}.$

Using the assumption $\partial_{1}U_{l}=\partial_{l}U_{1}$, and comparing the
coefficients of $(x_{1}-\beta_{i})^{-j}$ and $x^{i}$, we get that for
$l\in\\{2,\ldots,n\\}$

(2.14)
$\partial_{l}\alpha_{i,j+1}^{(1)}+j\alpha_{ij}^{(1)}\partial_{l}\beta_{i}+j\alpha_{ij}^{(l)}=0,\qquad\partial_{l}\xi_{i}^{(1)}=(i+1)\xi_{i+1}^{(l)}.$

The first equality with $j=0$ of (2.14) shows that
$\alpha_{i1}^{(1)}\in\mathbb{C}$, because $\alpha_{ij}^{(k)}$ are functions in
$x_{2},\ldots,x_{n}$ as prescribed.

Set

$\displaystyle\Phi(x)$ $\displaystyle=$
$\displaystyle\sum\limits_{i=1}\limits^{r}\alpha_{i1}^{(1)}\log(x_{1}-\beta_{i})+\sum\limits_{i=1}\limits^{r}\sum\limits_{j=2}\limits^{m}\frac{-1}{j-1}\frac{\alpha_{ij}^{(1)}}{(x_{1}-\beta_{i})^{j-1}}$
$\displaystyle+\sum\limits_{i=1}\limits^{p}\frac{\xi_{i}^{(1)}}{i+1}x_{1}^{i+1}+\int\xi_{0}^{(2)}dx_{2}+\ldots+\xi_{0}^{(n)}dx_{n},$

where the integration represents any primitive function of
$\xi_{0}^{(2)}dx_{2}+\ldots+\xi_{0}^{(n)}dx_{n}$. Direct calculations show
that $\partial_{1}\Phi=U_{1}$. For $l>1$ since

$\displaystyle\partial_{l}\Phi(x)$ $\displaystyle=$
$\displaystyle-\sum\limits_{i}\frac{\alpha_{i1}^{(1)}\partial_{l}\beta_{i}}{x_{1}-\beta_{i}}+\sum\limits_{i,j}\frac{-1}{j-1}\left(\frac{\partial_{l}\alpha_{ij}^{(1)}}{(x_{1}-\beta_{i})^{j-1}}+\frac{(j-1)\alpha_{ij}^{(1)}\partial_{l}\beta_{i}}{(x_{1}-\beta_{i})^{j}}\right)$
$\displaystyle+\sum\limits_{i}\frac{\partial_{l}\xi_{i}^{(1)}}{i+1}x^{i+1}+\xi_{0}^{(l)}.$

Using the equalities (2.14) and by some calculations, we get that

$\partial_{l}\Phi=U_{l},\qquad l=2,\ldots,n.$

These show that

$\Phi(x)=\int U_{1}dx_{1}+\ldots+U_{n}dx_{n},$

with possible a constant difference.

Denote by $\mathcal{G}$ the group of automorphisms over $K$ which keep
$\mathbb{C}(y)$, where $y=(x_{2},\ldots,x_{n})$. Since $K$ is the minimal
normal algebraic field extension of $\mathbb{C}(y)$, we get from the
properties of field extensions that $\mathcal{G}$ is a finite group. Denote by
$N=|\mathcal{G}|$, the order of $\mathcal{G}$. Set

$\Psi=\frac{1}{N}\sum\limits_{\sigma\in\mathcal{G}}\sigma(\Phi).$

Since $\sigma\in\mathcal{G}$ is an automorphism over the algebraic field
extension $K$ of $\mathbb{C}(y)$, it follows that

$\displaystyle\sigma\left(\alpha_{i1}\log(x_{1}-\beta_{i})\right)$
$\displaystyle=$ $\displaystyle\alpha_{i1}\log(x_{1}-\sigma(\beta_{i})),$
$\displaystyle\sigma\left(\int\gamma_{0}^{(2)}dx_{2}+\ldots+\gamma_{0}^{(n)}dx_{n}\right)$
$\displaystyle=$
$\displaystyle\int\sigma(\gamma_{0}^{(2)})dx_{2}+\ldots+\sigma(\gamma_{0}^{(n)})dx_{n},$

where the second equality may have a constant difference.

Since $\sigma(\partial_{l}\Phi)=\sigma(U_{l})$, we have

(2.15) $\partial_{l}\sigma(\Phi)=\sigma(U_{l})=U_{l},\qquad l=1,\ldots,n,$

where in the second equality we have used the facts that the numerators and
denominators of $U_{l}$’$s$ can be written in the polynomials of $x_{1}$ with
coefficients in $\mathbb{C}(y)$ and $\sigma$ keeps $\mathbb{C}(y)$.

The equalities (2.15) show that

$\partial_{l}\Psi=U_{l},\qquad l=1,\ldots,n.$

Moreover we have

$\Psi(x)=\sum\limits_{i=1}\limits^{r_{0}}c_{i}\log R_{i}(x)+R(x)+\int
S_{2}(y)dx_{2}+\ldots+S_{n}(y)dx_{n},$

where $c_{i}\in\mathbb{C}$, $R_{i},\,R\in\mathbb{C}(x)$ and
$S_{i}\in\mathbb{C}(y)$. Recall that $y=(x_{2},\ldots,x_{n})$. By the
expansions of $U_{k}$’$s$ in (2.13) and $\partial_{l}U_{k}=\partial_{k}U_{l}$,
it follows that

$\partial_{l}\xi_{0}^{(k)}(y)=\partial_{k}\xi_{0}^{(l)}(y),\qquad 2\leq
k,l\leq n.$

So we have

$\sum\limits_{\sigma\in\mathcal{G}}\sigma(\xi_{0}^{(k)})\in\mathbb{C}(y),\qquad
k=2,\ldots,n,$

and for $2\leq k,l\leq n$

$\partial_{l}\left(\sum\limits_{\sigma\in\mathcal{G}}\sigma(\xi_{0}^{(k)})\right)=\sum\limits_{\sigma\in\mathcal{G}}\sigma(\partial_{l}\xi_{0}^{(k)})=\sum\limits_{\sigma\in\mathcal{G}}\sigma(\partial_{k}\xi_{0}^{(l)})=\partial_{k}\left(\sum\limits_{\sigma\in\mathcal{G}}\sigma(\xi_{0}^{(l)})\right).$

These show that

$\partial_{x_{i}}S_{j}(y)=\partial_{x_{j}}S_{i}(y),\quad 2\leq i,j\leq n.$

Now for the integration $\int S_{2}(y)dx_{2}+\ldots+S_{n}(y)dx_{n}$ with
$y=(x_{2},\ldots,x_{n})$ we are in the same conditions as those of integration
$\int U_{1}(y)dx_{1}+\ldots+U_{n}(y)dx_{n}$, so working in a similar way as
that in the above proof we get that there exists a function $\Psi_{1}(y)$ such
that

$\partial_{l}\Psi_{1}(y)=S_{l}(y),\qquad l=2,\ldots,n,$

and

$\Psi_{1}(y)=\sum\limits_{i=1}\limits^{r_{1}}d_{i}\log T_{i}(y)+T(y)+\int
W_{3}(z)dx_{3}+\ldots+W_{n}(z)dx_{n},$

where $d_{i}\in\mathbb{C}$, $T_{i},\,T\in\mathbb{C}(y)$, and
$W_{i}\in\mathbb{C}(z)$ with $z=(x_{3},\ldots,x_{n})$ satisfy

$\partial_{x_{i}}W_{j}(z)=\partial_{x_{j}}W_{i}(z),\qquad 3\leq i,\,j\leq n.$

By induction we can prove that

$\displaystyle\Psi(x)$ $\displaystyle=$
$\displaystyle\sum\limits_{i=1}\limits^{r_{0}}c_{i}\log R_{i}(x)+R(x)$
$\displaystyle+\sum\limits_{j=2}^{n}\left(\sum\limits_{i=1}\limits^{r_{j}}c_{i}^{(j)}\log
R_{i}^{(j)}(x_{j},\ldots,x_{n})+R^{(j)}(x_{j},\ldots,x_{n})\right),$

where $c_{i},c_{i}^{(j)}\in\mathbb{C}$, and
$R_{i}^{(j)},R^{(j)}\in\mathbb{C}(x_{j},\ldots,x_{n})$. Recall that
$\partial_{l}\Psi(x)=U_{l}(x)$ for $l=1,\ldots,n$. Furthermore we have

$\displaystyle\exp\left(\Psi\right)$ $\displaystyle=$
$\displaystyle\exp\left(R(x)+\sum\limits_{j=2}^{n}R^{(j)}(x_{j},\ldots,x_{n})\right)$
$\displaystyle\times\prod\limits_{i=1}\limits^{r_{0}}\left(R_{i}(x)\right)^{c_{i}}\prod\limits_{j=2}^{n}\prod\limits_{i=1}^{r_{j}}\left(R_{i}^{(j)}(x_{j},\ldots,x_{n})\right)^{c_{i}^{(j)}}.$

Since $R,\,R^{(j)},\,R_{i},\,R_{i}^{(j)}\in\mathbb{C}(x)$, it follows that
$\exp(\Psi(x))$ is a Darboux function, and consequently is a Darboux Jacobian
multiplier. This proves Lemma 2.2. ∎

Summarizing Lemmas 2.1 and 2.2, we complete the proof of Theorem 1.2.

Acknowledgements. The author is partially supported by NNSF of China grant
11271252, RFDP of Higher Education of China grant 20110073110054, and
FP7-PEOPLE-2012-IRSES-316338 of Europe.

## References

  * [1] D. Blázquez-Sanz and Ch. Pantazi, A note on the Darboux theory of integrability of non-autonomous polynomial differential systems, Nonlinearity 25 (2012), 2615–2624.
  * [2] R. Bott and L.W. Tu, Differential Forms in Algebraic Topology, Springer–Verlag, New York, 1982.
  * [3] C. Christopher, Liouvillian first integrals of second order polynomial differential equations, Electron. J. Differential Equations 1999, no. 49, 1–7.
  * [4] C. Christopher and C. Li, Limit Cycles of Differential Equations, Birkhäuser, Basel, 2007.
  * [5] C. Christopher and J. Llibre, Integrability via invariant algebraic curves for planar polynomial differential systems, Ann. Diff. Eqns. 16 (2000), 5–19.
  * [6] C. Christopher, J. Llibre and J.V. Pereira, Multiplicity of invariant algebraic curves in polynomial vector fields, Pacific J. Math. 229 (2007), 63–117.
  * [7] G. Darboux, Mémoire sur les équations différentielles algébriques du premier ordre et du premier degré (Mélanges), Bull. Sci. Math. 2 (1878), 60–96; 123–144; 151–200.
  * [8] G. Darboux, De l’emploi des solutions particulières algébriques dans l’intégration des systèmes d’équations différentielles algébriques, C. R. Math. Acad. Sci. Paris 86 (1878), 1012–1014.
  * [9] F. Dumortier, J. Llibre and J.C. Artés, Qualitative theory of planar differential systems, UniversiText, Springer–Verlag, New York, 2006.
  * [10] J. Giné and M. Grau, Weierstrass integrability of differential equations, Appl. Math. Lett. 23 (2010), 523–526.
  * [11] J. Giné, M. Grau and J. Llibre, On the extensions of the Darboux theory of integrability, Nonlinearity 26 (2013), 2221–2229.
  * [12] J.P. Jouanolou, Equations de Pfaff algébriques, Lecture Notes in Math. 708, Springer–Verlag, New York/Berlin, 1979.
  * [13] S. Lang, Algebra , Graduate Texts in Mathematics 211 (third ed.), Springer-Verlag, New York, 2002.
  * [14] J. Llibre, Integrability of polynomial differential systems, in Handbook of differential equations, Elsevier, Amsterdam, 2004, pp.437–532.
  * [15] J. Llibre and Ch. Pantazi, Darboux theory of integrability for a class of nonautonomous vector fields, J. Math. Phys. 50 (2009), 102705, 19 pp.
  * [16] J. Llibre and X. Zhang, Darboux Theory of Integrability in $\mathbb{C}^{n}$ taking into account the multiplicity, J. Differential Equations 246 (2009), 541–551.
  * [17] J. Llibre and X. Zhang, Rational first integrals in the Darboux theory of integrability in $\mathbb{C}^{n}$, Bull. Sci. Math. 134 (2010), 189–195.
  * [18] J. Llibre and X. Zhang, On the Darboux integrability of polynomial differential systems, Qual. Theory Dyn. Syst. 11 (2012), 129–144.
  * [19] Y. Pan and X. Zhang, Algebraic aspects of integrability for polynomial differential systems, J. Appl. Anal. Comp. 3 (2013), 51–69.
  * [20] H. Poincaré, Sur l’intégration des équations différentielles du premier ordre et du premier degré I and II, Rendiconti del Circolo Matematico di Palermo 5 (1891), 161–191; 11 (1897), 193–239.
  * [21] M.J. Prelle and M.F. Singer, Elementary first integrals of differential equations, Trans. Amer. Math. Soc. 279 (1983), 215–229.
  * [22] D. Schlomiuk, Algebraic particular integrals, integrability and the problem of the center, Trans. Amer. Math. Soc. 338 (1993), 799–841.
  * [23] M.F. Singer, Liouvillian first integrals of differential equations, Trans. Amer. Math. Soc. 333 (1992), 673–688.
  * [24] X. Zhang, Global structure of quaternion polynomial differential equations, Commun. Math. Phys. 303 (2011), 301–316.


OMAE 2020 the ASME 2020 39th International Conference on Ocean, Offshore and
Arctic Engineering 28-3 June-July 2020 Fort Lauderdale USA OMAE2020-18556

# Deep Convolutional Recurrent Autoencoders for Flow Field Prediction

Sandeep R. Bukka Email address of<EMAIL_ADDRESS>Allan Ross
Magee  Keppel-NUS Corp Lab
National University of Singapore
Department of Civil and Environmental Engineering
Singapore  Rajeev K. Jaiman Currently at University of British Columbia
Department of Mechanical Engineering
National University of Singapore
Singapore

###### Abstract

In this paper, an end-to-end nonlinear model reduction methodology is
presented based on the convolutional recurrent autoencoder networks. The
methodology is developed in the context of overall data-driven reduced order
model framework proposed in the paper. The basic idea behind the methodology
is to obtain the low dimensional representations via convolutional neural
networks and evolve these low dimensional features via recurrent neural
networks in time domain. The high dimensional representations are constructed
from the evolved low dimensional features via transpose convolutional neural
networks. With an unsupervised training strategy, the model serves as an end
to end tool which can evolve the flow state of the nonlinear dynamical system.
The convolutional recurrent autoencoder network model is applied on the
problem of flow past bluff bodies for the first time. To demonstrate the
effectiveness of the methodology, two canonical problems namely the flow past
plain cylinder and the flow past side-by-side cylinders are explored in this
paper. Pressure and velocity fields of the unsteady flow are predicted in
future via the convolutional recurrent autoencoder model. The performance of
the model is satisfactory for both the problems. Specifically, the multiscale
nature and the gap flow dynamics of the side-by-side cylinders are captured by
the proposed data-driven model reduction methodology. The error metrics, the
normalized squared error and the normalized reconstruction error are
considered for the assessment of the data-driven framework.

Keywords: Data-driven modeling, reduced-order model, recurrent neural
networks, convolutional neural networks, autoencoders

## INTRODUCTION

In the past decade, the hardware and techniques in experiments and
computations have become much more advanced and resulted in the generation of
massive amounts of data. It is definitely a welcoming note to have more data,
but the gigantic scale of modern datasets poses serious challenges in terms of
analysis. The demand for data driven methods rises from these challenges and
the primary motivation of such methods is to come up with efficient models
which can characterize the data in meaningful ways. One other way to put the
analysis via data-driven methods is learning from experience. In nature,
various animals ranging from fish, cheetah, birds and other aerial, aquatic
and terrestrial lifeforms are able to achieve remarkable levels of efficiency
in terms of speed, manoeuvring and control by manipulating the fluid
environment around them. The critical point here is to question the learning
process of these organisms.

It is very clear that they do not have any prior knowledge of calculus or the
Navier-Stokes equations. Humans also achieved remarkable engineering feats of
fluid manipulation before the knowledge of current equation based analysis
came into existence. The feasible explanation for such learning is that it
came from intuition and experience [1]. It is also an undeniable fact that the
current approach of physics-based engineering is an epitome of human
achievement which revolutionized the performance of the engineering systems.
However, there are serious hurdles associated with physics-based engineering
such as high-dimensionality and non-linearity, which defy closed-form
solutions and limit real-time optimization and control efforts. With the
recent surge in computational capabilities and advancements in machine
learning, deep learning, artificial intelligence, learning from experience has
taken the spotlight.

The marine/offshore industry is going through a paradigm shift with regard to
automation and autonomy. While automation generally implies to the activity
performed devoid of human intervention, autonomy is the ability of a system to
achieve goals while operating independently of external control (i.e., self-
sufficient and self-directed system). Major industries such as automobile,
aerospace, health care, and logistics operations are already seeing the
benefits of having smart autonomous systems. The marine/offshore industry is
not an exception. Since the presence of technologies related to smart
autonomous systems (SAS) is burgeoning at a rapid pace, the mechanics
community needs to integrate their tools and techniques with these SAS
technologies so that the systems can be made more efficient, reliable and
adaptable to marine/offshore applications. The construction of an efficient
ROM and its integration with the sensory data is necessary to achieve a smart
autonomous system [2].

An overall representation of the proposed data-driven framework is presented
in Figure 1. The framework encompasses the previously reported POD-RNN hybrid
model in [2]. The main idea of this hybrid model comprises the following two
steps:

  * •

1\. The identification of a low-dimensional manifold $\bm{\Phi}$ embedded in
$\mathbb{R}^{N}$ on which most of the data is supported. This yields, in some
sense, an optimal low-dimensional representations $\textbf{A}=f(\textbf{S})$
of the data S in terms of the intrinsic coordinates on $\bm{\Phi}$, and

  * •

2\. The identification of a dynamic model which effectively evolves the low
dimensional representation A on the manifold $\bm{\Phi}$.

While the POD based model reduction is optimal and can generate physically
interpretable modes it is not good enough for practical problems. In such
cases, the number of dominant modes will increase and their computation can
pose serious challenges. Up to 27 modes were required to capture the dominant
energy of the side-by-side cylinder whereas only 4 modes were required for the
flow past an isolated (single) cylinder. Similarly it was shown in [3], that
123 modes were required for the flow past a square cylinder at $Re=22000$. As
a natural consequence, the question of completely bypassing POD in building
reduced order model arises and the alternative technique based on
convolutional neural networks (CNN) will be explored in this paper for
extracting the low dimensional features.

In addition to the POD-RNN, and also inspired from the overall framework, is a
complete data-driven one, where both the low-dimensional representation of the
state variable and its temporal evolution are learned via machine learning
techniques. Such approach has been explored in [4] in which an autoencoder is
used to learn a low-dimensional representation of the high-dimensional state,

$\textbf{A}=\mathcal{E}(\textbf{S}),$ (1)

where $\textbf{S}\in\mathbb{R}^{N}$ high-dimensional state of the system,
$\textbf{A}\in\mathbb{R}^{N_{A}},N_{A}<N$, $\mathcal{E}$ is a nonlinear
enocoder network and a linear recurrent model is used to evolve the low
dimensional features

$\textbf{A}^{n+1}=\textbf{K}\textbf{A}^{n},$ (2)

where $\textbf{K}\in\mathbb{R}^{N_{A}\times N_{A}}$. This approach was first
introduced in the context of learning a dictionary of functions used in
extended dynamic mode decomposition to approximate the Koopman operator of a
nonlinear system [5]. Based on the framework introduced in [6] [4] [7], [8]
introduced a deep convolutional autoencoder architecture which provides
certain advantages in identifying low-dimensional representations of the input
data. Figure 2 depicts the schematic of the model reduction via convolutional
recurrent autoencoder network. Benchmark problems of lid-driven cavity,
viscous Burgers equation were explored via the above approach in [8]. However,
such problems are far away of the practical scenarios. There is a need to
address the development of efficient models for more complex problems such as
flow past side-by-side cylinders, flexible offshore riser, and ocean wave
predictions. Such problems contain more complex physics with multiscale
features and various physical phenomena such as vortex-induced vibration,
vortex shedding, gap flow, flip flopping, bifurcation and turbulence. In this
paper, we extend the capabilities of the convolutional recurrent autoencoder
model by applying it for the first time on the more complex problems of fluid-
structure interaction. [8] contains more details on the model reduction via
convolutional recurrent networks and they are also presented here for the
completeness of the present paper.

Figure 1: A generic sketch of data-driven reduced order framework using
POD/CNN and RNN techniques Figure 2: Schematic of Convolutional Recurrent
Autoencoder network

## NUMERICAL METHODOLOGY

This section deals with the formulation used for high-dimensional full-order
model followed by the implementation of deep convolutional recurrent
autoencoder model for a coupled fluid-structure problem. To begin with, the
formulation of the coupled variational fluid-structure solver is outlined in
the following section.

### Full-order model formulation

A numerical scheme implementing Petrov-Galerkin finite-element and semi-
discrete time stepping are adopted in the present work [9, 10]. The
incompressible Navier-Stokes equations are used in the Arbitrary Lagrangian-
Eulerian (ALE) reference frame and formulated in the following form:

$\displaystyle\rho^{\mathrm{f}}\left(\frac{\partial\textbf{u}^{\mathrm{f}}}{\partial
t}\bigg{\rvert}_{\widehat{x}}+\left(\textbf{u}^{\mathrm{f}}-\textbf{w}\right)\cdot\bm{\nabla}\textbf{u}^{\mathrm{f}}\right)$
$\displaystyle=\bm{\nabla}\cdot\bm{\sigma}^{\mathrm{f}}+\textbf{b}^{\mathrm{f}}\quad\text{on}\quad\Omega^{f}(t),$
(3) $\displaystyle\bm{\nabla}\cdot\textbf{u}^{\mathrm{f}}$
$\displaystyle=0\quad\text{on}\quad\Omega^{\mathrm{f}}(t),$ (4)

where $\textbf{u}^{\mathrm{f}}=\textbf{u}^{\mathrm{f}}(\textbf{x},t)$ and
$\textbf{w}=\textbf{w}(\textbf{x},t)$ are the fluid and mesh velocities,
respectively. In Eq. (3) the partial time derivative with respect to the ALE
referential coordinate $\widehat{x}$ is constant. Here
$\textbf{b}^{\mathrm{f}}$ represents the body force per unit mass and
$\bm{\sigma}^{\mathrm{f}}$ is the Cauchy stress tensor for a Newtonian fluid
which is defined as

$\bm{\sigma}^{\mathrm{f}}=-\textit{p}\textbf{I}+\mu^{\mathrm{f}}\left(\bm{\nabla}\textbf{u}^{\mathrm{f}}+(\bm{\nabla}\textbf{u}^{\mathrm{f}})^{T}\right),$
(5)

where p, $\mu^{\mathrm{f}}$ and I are the hydrodynamic pressure, the dynamic
viscosity of the fluid, and the identity tensor, respectively. A rigid-body
structure submerged in the fluid experiences unsteady fluid forces and
consequently may undergo flow-induced vibrations if the body is mounted
elastically. To simulate the translational motion of a rigid body about its
center of mass, the equation along the Cartesian axes is given by:

$\textbf{m}\cdot\frac{\partial\textbf{u}^{\mathrm{s}}}{\partial
t}+\textbf{c}\cdot\textbf{u}^{\mathrm{s}}+\textbf{k}\cdot\left(\bm{\phi}^{\mathrm{s}}(\textbf{z}_{0},t)-\textbf{z}_{0}\right)=\textbf{F}^{\mathrm{s}}+\textbf{b}^{\mathrm{s}}\>\text{on}\>\Omega^{\mathrm{s}},$
(6)

where m, c, k, $\textbf{F}^{\mathrm{s}}$ and $\textbf{b}^{\mathrm{s}}$ are the
mass, damping coefficient, and stiffness coefficient vectors for the
translational motions, fluid traction, and body forces on the rigid body,
respectively. Here $\Omega^{\mathrm{s}}$ represents the domain occupied by the
rigid body and $\textbf{u}^{\mathrm{s}}(t)$ represents the velocity of the
immersed rigid body. The fluid and the structural equations are coupled by the
continuity of velocity and traction along the fluid-structure interface.

The new position of the rigid body is updated through a position vector
$\bm{\phi}^{\mathrm{s}}$ , which maps the initial position $\textbf{z}_{0}$ of
the rigid body to its new position at time t. Let $\gamma$ be the Lagrangian
point on $\Gamma$ and its corresponding mapping position vector to the new
position after the motion of the rigid body is $\phi(\gamma,t)$ at time $t$.
Since the position and flow field around the moving rigid body is updated
continuously, the no-slip and traction continuity conditions should be
satisfied on the fluid-body interface $\Gamma$

$\displaystyle\textbf{u}^{\mathrm{f}}(\bm{\phi}^{\mathrm{s}}(\textbf{z}_{0},t),t)$
$\displaystyle=\textbf{u}^{\mathrm{s}}(\textbf{z}_{0},t),$ (7)
$\displaystyle\int_{\phi(\gamma,t)}\bm{\sigma}^{\mathrm{f}}(\textbf{x},t)\cdot\textbf{n}\text{d}\Gamma+\int_{\gamma}\textbf{F}^{\mathrm{s}}\text{d}\Gamma$
$\displaystyle=0\quad\forall\>\gamma\>\in\>\Gamma,$ (8)

where n is the outer normal to the fluid-body interface. The characterization
of the moving fluid-body interface is constructed by means of the ALE
technique [11] The movement of the internal ALE nodes is constructed by
solving a continuum hyperelastic model for the fluid mesh such that the mesh
quality does not deteriorate as the displacement of the body increases. For
the spatial and material mapping problem, we use classical Neo-Hookean
material properties for the ALE variational formulation [12] which does not
entail any additional user-defined remeshing parameter. The weak variational
form of Eq. (3) is discretized in space using $\mathbb{P}_{2}/\mathbb{P}_{1}$
isoparametric finite elements for the fluid velocity and pressure. The second-
order backward-differencing scheme is used for the time discretization of the
Navier-Stokes system [13]. A partitioned staggered scheme is considered for
the full-order simulations of fluid-structure interaction [14]. The above
coupled variational formulation completes the presentation of full-order model
for high-fidelity simulation. The employed in-house FSI solver has been
extensively validated in [15, 16].

## Model reduction via convolutional recurrent autoencoders

### Overview

Fully-connected autoencoders are used in many applications for dimensionality
reduction. The high dimensionality of DNS-level input data can act as a major
bottleneck for applying fully-connected autoencoders. The input data should be
flattened into an 1D array when we consider autoencoders. This might result in
elimination of local spatial relations between the data. In other words,
dominant flow features, often found in FSI simulations might be ignored or
might not be captured efficiently. Therefore, direct application of fully-
connected autoencoders on DNS-level input data is generally prohibited. The
current model seeks to extract dominant flow features present in many fluid
flow and fluid-structure interaction data by employing convolutional neural
networks.

In particular, rather than applying a fully-connected autoencoder directly to
complex, high-dimensional simulation or experimental data, we apply it to a
vectorized feature map of a much lower-dimension obtained from a deep
convolutional network acting directly on the high-dimensional data. Wrapping
the fully-connected autoencoder with a convolutional neural network has two
significant advantages [8]:

  * •

1\. The local approach of each convolutional layer helps to exploit local
correlations in field values. Thus, much the same way finite-difference
stencils can capture local gradients, each filter $\textbf{K}^{f}$ in a filter
bank computes local low-level features from a small subset of the input.

  * •

2\. The shared nature of each filter bank both allows identification of
similar features throughout the input domain and reduction of the overall
number of trainable parameters compared to a fully-connected layer with the
same input size.

Figure 3: A sketch of 12 layer convolutional recurrent autoencoder model

### Proposed architecture

As shown in figure 3, we construct a 12 layer convolutional autoencoder model.
The first 4 layers together forms the convolution aspect of the model. The two
dimensional input $S\in\mathbb{R}^{N_{x}\times N_{y}}$, with $N_{x}=N_{y}=64$
is first passed through this convolutional block. All the layers in the
convolution block used a filter bank $\textbf{K}^{f}\in\mathbb{R}^{4\times 4}$
with the number of filters in the first layer $f$ increasing from 4 in the
first layer to 32 in the fourth layer. A stride of 2 and a padding of 1 is
used across all the layers.

The process of transpose convolution is carried out in the last four layers of
the network. The tranpose convolution is not exactly the reverse of regular
convolutions but it up-samples the low dimensional representations across four
layers gradually to reconstruct the high dimensional state. The output size of
the transpose convolutional layer is exactly the inverse of regular
convolutional layer and is given by $N_{o}=(N_{i}-1)*s+f-2*P$, where $s$ is
the stride, $f$ is the filter size, $P$ is the padding, $N_{i}$ is the input
size and $N_{o}$ is the output size. Table 1 outlines the architecture of the
convolutional encoder and decoder subgraphs. In this work, we consider the
sigmoid activation function $\sigma(s)=\frac{1}{1+\exp^{-s}}$ for each layer
of the autoencoder.

Layer | filter size | filters | stride | padding
---|---|---|---|---
1 | $4\times 4$ | 4 | 2 | 1
2 | $4\times 4$ | 8 | 2 | 1
3 | $4\times 4$ | 16 | 2 | 1
4 | $4\times 4$ | 32 | 2 | 1
9 | $4\times 4$ | 16 | 2 | 1
10 | $4\times 4$ | 8 | 2 | 1
11 | $4\times 4$ | 4 | 2 | 1
12 | $4\times 4$ | 1 | 2 | 1
Table 1: Filter sizes and strides of the layers in convolutional autoencoder
network. Layers (1 to 4) correspond to encoder part and layers (9 to 12)
correspond to decoder part

A 4 layer fully connected regular autoencoder with a standard feed forward
neural network is also used in the full network. The first two layers performs
encoding by taking the vectorized form of 32 feature maps from the final layer
of convolution block and returns a one dimensional vector
$\textbf{A}\in\mathbb{R}^{N_{A}}$. The last two layers decode the evolved
feature vector into high dimensional vector which will be reshaped as feature
maps and passed into the transpose convolution block. Figure 3 illustrates the
whole architecture with the changes in the size of the data as it flows
throughout the entire model.

The key innovation in using convolutional autoencoders in model reduction is
that it allows for nonlinear autoencoders and thus nonlinear model reduction
can be applied to a large input data in a way that exploits structures
inherent in many physical systems.

### Evolution of features

Now we proceed to the modeling the evolution of low-dimensional features A in
a computationally efficient manner. In terms of of analysis viewpoint, it is
beneficial to identify the the linear dynamics of A. However, we consider a
general case of learning feature dynamics in a nonlinear setting.

Consider a set of observations $\\{\textbf{s}^{n}\\}_{n=0}^{m}$,
$\textbf{s}^{n}\in\mathbb{R}^{N}$ obtained from a CFD simulation or through
experimental sampling. An optimal low-dimensional representation
$\textbf{A}^{n}\in\mathbb{R}^{N_{A}}$ where $N_{A}\ll N$ is obtained for each
observations. These low-dimensional representations are obtained via
convolutional autoencoder described in the previous section. The model for
temporal evolution of this low dimensional features is constructed in a
complete data-driven fashion.

The LSTM networks are employed to model the evolution of A. The details on the
LSTM networks are outlined in the previous chapter. In a typical machine
translation task, the size of the hidden states and number of layers are large
when compared to the size of the feature vectors in the current work. Also, we
are not interested in evolving the full high dimensional system. A single
layer LSTM network is found to be sufficient for the current problem of
evolving the feature vectors $\textbf{A}^{n}$. It is important to note that
the the evolution of A does not require information from the full state S,
thereby avoiding a costly reconstruction at every step.

Initializing with a known low-dimensional representation $\textbf{A}^{0}$, one
obtains a prediction for the following steps by iteratively applying

$\hat{\textbf{A}}^{n+1}=\bm{\mathcal{H}}(\hat{\textbf{A}}^{n})\quad
n=1,2,3,\dots$ (9)

where $\hat{\textbf{A}}^{1}=\bm{\mathcal{H}}(\textbf{A}^{0})$ and
$\bm{\mathcal{H}}$ represents the action of a standard LSTM cell. A graphical
representation of this model is depicted in Figure 2.

## Unsupervised training strategy

The unsupervised training approach is the key component of the current model.
The weights of both the convolutional autoencoder network and LSTM recurrent
neural network are adjusted in a joint fashion. The main challenge here is to
prevent overfitting of the CNN and RNN portion of the model. The construction
of the training dataset as well as the training and evaluation algorithms are
presented in this section.

Consider a dataset $\\{\textbf{s}^{1},\dots,\textbf{s}^{m}\\}$, where
$\textbf{s}\in\mathbb{R}^{N_{x}\times N_{y}}$ is a 2D snapshot of some
dynamical system (e.g., a pressure field defined on a 2D grid). This dataset
is broken up into a set of $\textbf{N}_{s}$ finite time training sequences
$\\{\textbf{S}^{1},\dots,\textbf{S}^{N_{s}}\\}$ where each training sequence
$\textbf{S}^{i}\in\mathbb{R}^{N_{x}\times N_{y}\times N_{t}}$ consists of
$N_{t}$ snapshots.

The fluctuations around the temporal mean are considered as followed generally
in the data science community for improving the training. For example, we
define

$\textbf{s}^{{}^{\prime}n}=\textbf{s}^{n}-\bar{\textbf{s}},$ (10)

where $\bar{\textbf{s}}=\frac{1}{m}\sum_{n=1}^{m}\textbf{s}^{n}$ is the
temporal average over the entire dataset and $\textbf{s}^{{}^{\prime}}$ are
the fluctuations around this mean. Next the scaling of the above data is
conducted as follows

$\textbf{s}^{{}^{\prime}n}_{s}=\frac{\textbf{s}^{{}^{\prime}n}-\textbf{s}_{min}^{{}^{\prime}}}{\textbf{s}_{max}^{{}^{\prime}}-\textbf{s}_{min}^{{}^{\prime}}}$
(11)

where each $\textbf{s}^{{}^{\prime}n}_{s}\in[0,1]^{N_{x}\times N_{y}}$. The
scaling is carried out to ensure that the data lies in the interval of
$(0,1)$. This requirement stems from the fact that a sigmoid activation
function is used in all the layers of convolutional autoenocoder. The training
dataset with the above modifications has the following form.

$\mathcal{S}=\\{\textbf{S}^{{}^{\prime}1}_{s},\dots,\textbf{S}^{{}^{\prime}N_{s}}_{s}\\}\in[0,1]^{N_{x}\times
N_{y}\times N_{t}\times N_{s}}$ (12)

where each training sample
$\textbf{S}^{{}^{\prime}i}_{s}=[\textbf{S}^{{}^{\prime}1}_{s,i},\dots,\textbf{S}_{s,i}^{N_{t}}]$
is a matrix consisting of the feature-scaled fluctuations.

The process of training entails two stages primarily. The convolutional
autoencoder takes an $N_{b}$\- sized batch of the training data
$\mathcal{S}^{b}\subset\mathcal{S}$, where
$\mathcal{S}^{b}\in[0,1]^{N_{x}\times N_{y}\times N_{t}\times N_{b}}$,and
outputs both the current $N_{b}$-sized batch of low-dimensional
representations of the training sequence

$\mathcal{A}^{b}=\\{\textbf{A}^{1},\dots,\textbf{A}^{N_{b}}\\}\in\mathbb{R}^{N_{A}\times
N_{t}\times N_{b}},$ (13)

where
$\textbf{A}^{i}=[\textbf{A}^{1}_{i},\dots,\textbf{A}_{i}^{N_{t}}]\in\mathbb{R}^{N_{A}\times
N_{t}}$ and a reconstruction $\hat{\mathcal{S}}^{b}$ of the original input
training batch in the first stage. This is represented as purple arrow in
figure 4. The second stage is indicated by the blue arrows in the same figure
4, where the first feature vector of each sequence is used to initialize and
iteratively update Equation 9 to get a reconstruction $\hat{\mathcal{A}}^{b}$
of the low-dimensional representations of the training batch Equation 13.

Figure 4: An illustration of training process of convolutional recurrent
autoencoder model

A loss function is constructed that equally weights the error in the full-
state reconstruction and the evolution of the low-dimensional representations.
The target of the training is to find the model parameters $\theta$ such that
for any sequence
$\textbf{S}_{s}^{{}^{\prime}}=[\textbf{S}_{s}^{{}^{\prime}1},\dots,\textbf{S}_{s}^{{}^{\prime}N_{t}}]$
and its corresponding low-dimensional representation
$\textbf{A}=[\textbf{A}^{1},\dots,\textbf{A}^{N_{t}}]$, minimizes the
following error between the truth and the predictions as follows:

$\displaystyle\mathcal{J}(\theta)$
$\displaystyle=[\mathcal{L}(\hat{\textbf{S}}_{s}^{{}^{\prime}},\textbf{S}_{s}^{{}^{\prime}},\hat{\textbf{A}},\textbf{A})]$
(14)
$\displaystyle=\left[\frac{\alpha}{N_{t}}\operatornamewithlimits{\sum}_{n=1}^{N_{t}}\frac{\|\textbf{S}_{s}^{{}^{\prime}n}-\hat{\textbf{S}}_{s}^{{}^{\prime}n}\|_{2}^{2}}{\|\textbf{S}_{s}^{{}^{\prime}n}\|_{2}^{2}}+\frac{\beta}{N_{t}}\frac{\|\textbf{A}^{n}-\hat{\textbf{A}}^{n}\|_{2}^{2}}{\|\textbf{A}^{n}\|_{2}^{2}}\right]$

where $\alpha=\beta=0.5$. The values of $\alpha$ and $\beta$ are chosen in
order to give equal weightage to the errors in full state reconstruction and
feature prediction respectively. The predictions are denoted by the symbol
$\hat{}$. In practice, the error is approximated by averaging
$\mathcal{L}(\hat{\textbf{S}}_{s}^{{}^{\prime}},\textbf{S}_{s}^{{}^{\prime}},\hat{\textbf{A}},\textbf{A})$
over all samples in a training batch during each backward pass. The
convolutional autoencoder performs a regular forward pass and also constructs
the low dimensional representations simultaneously at every training step.
These low dimensional representations are used to train the RNN.

The ADAM optimizer [17] is used for updating the weights during the training
process. It is a modified version of the traditional stochastic gradient
descent that computes adaptive learning rates for different parameters using
estimates of first and second moments of the gradients. Algorithm 1 outlines
the offline training of the convolutional recurrent autoencoder in more
detail. Figure 4 illustrates the overall schematic of the training process for
easy comprehension. This model was built and trained using the open-source
deep learning library TensorFlow.

Figure 5: An illustration of prediction process of convolutional recurrent
autoencoder model

The online prediction is straightforward and a schematic of its process is
depicted in figure 5 for the ease of understanding. Firstly, a low-dimensional
representation $\textbf{A}^{0}\in\mathbb{R}^{N_{A}}$ is constructed using the
encoder network for a given initial condition
$\textbf{S}^{0}\in[0,1]^{N_{x}\times N_{y}}$ and a set of trained parameters
$\theta^{*}$. Equation 9 is applied iteratively for $N_{t}$ steps with
$\textbf{A}^{0}\in\mathbb{R}^{N_{A}}$ as the initial solution to get
predictions of low dimensional representations. Finally, the full-dimensional
state is reconstructed $\hat{\textbf{S}}^{n}$ from the low dimensional
representation $\hat{\textbf{A}}^{n}$ at every time step or at any specific
instance. The online prediction procedure is outlined in Algorithm 2.

Input : Training dataset $\mathcal{X}\in[0,1]^{N_{x}\times N_{y}\times
N_{t}\times N_{s}}$, number of train-steps $N_{train}$, batch size $N_{b}$

Result : Trained model parameters $\theta$

Randomly initialize $\theta$;

for  _$i\in\\{1,\dots,N_{train}\\}$_ do

Randomly sample batch from training data: $\mathcal{X}^{b}\subset\mathcal{X}$;

Flatten batch-mode: $\mathcal{X}^{b_{AE}}\leftarrow flatten(\mathcal{X}^{b})$
s.t. $\mathcal{X}^{b_{AE}}\in[0,1]^{N_{x}\times N_{y}\times(N_{t}.N_{b})}.$;

Encoder forward pass: $\tilde{\mathcal{H}}^{b}\leftarrow
f_{enc}(\mathcal{X}^{b_{AE}})$ where
$\tilde{\mathcal{H}}^{b}\in\mathbb{R}^{N_{h}\times(N_{t}.N_{b})}$;

Decoder forward pass: $\hat{\mathcal{X}}^{b_{AE}}\leftarrow
f_{dec}(\tilde{\mathcal{H}}^{b})$;

Reshape low-dimensional features: $\mathcal{H}^{b}\in\mathbb{R}^{N_{h}\times
N_{t}\times N_{b}}\leftarrow$ reshape $(\tilde{\mathcal{H}}^{b})$;

Initialize RNN subgraph loop: $\hat{\textbf{h}}_{i}^{2}\leftarrow
f_{LSTM}(\textbf{h}_{i}^{1})$ for
$i\in\\{1,\dots,N_{b}\\}$,$\textbf{h}_{i}^{1}\subset\mathcal{H}^{b}$;;

for  _$n\in\\{2,\dots,N_{t-1}\\}$_ do

$\hat{\textbf{h}}_{i}^{n+1}\leftarrow f_{LSTM}(\hat{\textbf{h}}_{i}^{n})$ for
$i\in\\{1,\dots,N_{b}\\}$,$\hat{\textbf{h}}_{i}^{n}\subset\hat{\mathcal{H}}^{b}$;

end for

Using $\mathcal{X}^{b},\hat{\mathcal{X}}^{b},\mathcal{H}^{b}$ and
$\hat{\mathcal{H}}^{b}$ calculate approximate gradient $\hat{\textbf{g}}$ of
equation;

Update parameters: $\bm{\theta}\leftarrow ADAM(\hat{\textbf{g}})$

end for

Algorithm 1 Convolutional Recurrent Autoencoder training algorithm [8]

Input : Initial condition $\textbf{x}^{0}\in[0,1]^{N_{x}\times N_{y}}$, number
of prediction steps $N_{t}$.

Result : Model prediction
$\hat{\textbf{X}}=[\hat{\textbf{x}}^{1},\dots,\hat{\textbf{x}}^{N_{t}}]\in[0,1]^{N_{x}\times
N_{y}\times N_{t}}$

Load trained parameters $\bm{\theta}^{*}$;

Endcoder forward pass: $\textbf{h}^{0}\leftarrow f_{enc}(\textbf{x}^{0})$;

Initialize RNN subgraph loop: $\hat{\textbf{h}}^{1}\leftarrow
f_{LSTM}(\textbf{h}^{0})$;

for  _$n\in\\{1,\dots,N_{t-1}\\}$_ do

$\hat{\textbf{h}}^{n+1}\leftarrow f_{LSTM}(\hat{\textbf{h}})^{n}$

end for

Decoder forward pass: $\hat{\textbf{X}}\leftarrow f_{dec}(\hat{\textbf{H}})$,
where
$\hat{\textbf{H}}=[\hat{\textbf{h}}^{1},\dots,\hat{\textbf{h}}^{N_{t}}]$;

Algorithm 2 Convolutional Recurrent Autoencoder prediction algorithm [8]

## Flow past an isolated cylinder

Figure 6: Schematic diagram of the problem setup: flow past a plain cylinder

(a)  (b)

Figure 7: (a) Unstructured mesh used in full order model simulation, (b)
Square uniform grid used as input function in the convolutional recurrent
autoencoder model. Figure 8: Normalized squared error $E_{f}$ for the velocity
field in X-direction ($U$): flow past a plain cylinder

(a)  (b)  (c)

Figure 9: Comparison of truth and predicted fields along with normalized
reconstruction error $E_{n}$ at (a) $t=1050s$, (b) $t=1125s$, (c) $t=1200s$
for velocity field in X-direction ($U$): flow past plain cylinder  Figure 10:
Normalized squared error $E_{f}$ for the Pressure field ($P$): flow past a
plain cylinder

The end to end nonlinear model reduction tool based on convolutional recurrent
autoencoder netoworks is applied first on the canonical problem of the flow
past a plain cylinder. The schematic of the problem set up is the same as the
one depicted in Fig. 6. The simulations were carried out via finite element
Navier-Stokes solver described above to generate the training and test data.

The Reynolds number of the problem is taken as $Re=100$. The full order
simulation was carried out for $1250s$ with a time step of $0.05s$. A total of
$5000$ snapshots of the simulation were collected for every $0.25s$. First we
collect the velocity field in X-direction ($U$) for the analysis. The first
and foremost step before proceeding with the training of the model is to map
the unstructured data coming from full order simulations onto a square uniform
grid of size $64\times 64$. It is essential that the input to the
convolutional recurrent autoencoder to be a square matrix of uniform grid. We
use the size $64\times 64$ throughout all our numerical examples for the sake
of simplicity. We map the data from the unstructured grid to the uniform grid
via $MATLAB$ function $``griddata".$ Figure 7 depicts both the unstructured
grid used in the full order simulations and the uniform grid which is used an
input to the convolutional recurrent autoencoder model. The presence of the
rigid body is ensured by masking the region covered by the cylinder with a
mandatory function which zeroes the predictions made by the model in the
points inside the cylinder. The procedure described in the previous sections
is followed to generate the feature scaled data for training.

$\mathcal{U}=\\{\textbf{U}^{{}^{\prime}1}_{s},\dots,\textbf{U}^{{}^{\prime}N_{s}}_{s}\\}\in[0,1]^{N_{x}\times
N_{y}\times N_{t}\times N_{s}}$ (15)

where each training sample
$\textbf{U}^{{}^{\prime}i}_{s}=[\textbf{u}^{{}^{\prime}1}_{s,i},\dots,\textbf{u}_{s,i}^{{}^{\prime}N_{t}}]$
with $N_{x},\;N_{y}=64$ and $N_{t}=25$, $N_{s}=150$. Of the total $5000$
samples, $3000$ are used for training whereas $1000$ are used for validation
and the final $1000$ are used for testing.

Four convolutional recurrent autoencoder models were trained using the above
dataset, with low-dimensional representations of sizes
$N_{h}=8,\;16,\;32,\;64$. All the four models are trained on single Intel
E5-2690v3 (2.60GHz, 12 cores) CPU node for $N_{train}=500000$ iterations. The
training took approximately $48$ hours. The online prediction is carried for
1000 steps with the trained model. Figure 8 shows the normalized squared error
for the velocity field in X-direction $(U)$.

$E_{f}=\frac{\|\textbf{u}^{{}^{\prime}n}_{s}-\hat{\textbf{u}}^{{}^{\prime}n}_{s}\|^{2}_{2}}{\|\textbf{u}_{s}^{{}^{\prime}n}\|^{2}_{2}}$
(16)

It is interesting to note that there is an order of magnitude difference in
normalized squared error $E_{f}$ for different sizes of low-dimensional
representations. The highest is with $N_{h}=8$, followed by $N_{h}=64$,
$N_{h}=16$ and $N_{h}=32$. The error is diverging rapidly for the case of
$N_{h}=16$, whereas it is diverging slowly for the case of $N_{h}=64$. We have
chosen $N_{h}=32$ as the final trained model to depict the flow field
predictions in detail.

Figure 9 depicts the comparison of truth and predicted values for the velocity
field in $X$-direction ($U$). The normalized reconstruction error $E_{n}$ is
constructed by taking the absolute value of differences between and truth and
predictions and normalizing it with $L_{2}$ norm of the truth data and is
given by

$E_{n}=\frac{|\textbf{u}_{s}^{{}^{\prime}n}-\hat{\textbf{u}}_{s}^{{}^{\prime}n}|}{\|\textbf{u}_{s}^{{}^{\prime}n}\|_{2}}$
(17)

The predictions from the convolutional recurrent autoencoder model are quite
satisfactory. The alternating vortices are clearly captured in the predictions
as shown in figure 9. The predictions are also quite smooth without any
oscillations. The reconstruction error $E_{n}$ is quite low in the order of
$10^{-4}$ for the predictions.

The same study is repeated with the pressure field. All the other details
i.e., the size of training, validation and testing samples are the same as the
case with the velocity field. Four convolutional recurrent autoencoder models
are evaluated for the pressure field predictions by changing the sizes of the
low-dimensional representations. Figure 10 depicts the normalized squared
error $E_{f}$ for the pressure field predictions from convolutional recurrent
autoencoder model. The normalized squared error $E_{f}$ for this case is
calculated by replacing the variable u with p. The errors from all the four
models are the in the same order of magnitude. There is no evidence for
divergence in the error but the magnitude of errors are higher when compared
to velocity field predictions. Without any loss of generality, we have chosen
$N_{h}=32$ for the depiction of full flow field predictions.

Figure 11 depicts the comparison of truth and predicted values for the
pressure field $P$ along with the normalized reconstruction error $E_{n}$. It
is very clear that the predictions are devoid of any spurious oscillations.
The normalized reconstruction error $E_{n}$ is in the order of $10^{-3}$ and
is higher when compared to velocity field predictions.

(a)  (b)  (c)

Figure 11: Comparison of truth and predicted fields along with normalized
reconstruction error $E_{n}$ at (a) $t=1050s$, (b) $t=1125s$, (c) $t=1200s$
for pressure field ($P$): flow past plain cylinder

There appears to be some common points in the predictions of both velocity
field $U$ and the pressure field $P$ for the problem of flow past an isolated
single cylinder. The predictions are quite smooth and close to the truth data.
In both the cases, the near wake region behind the cylinder appears to have
majority of the reconstruction error. The nonlinearities associated with the
near wake can be the reason for the errors.

It should be noted that we are using closed recurrent neural network in the
convolutional recurrent autoencoder model i.e., the prediction from the
previous step is used as input to the next prediction. The main advantage of
using such a model is that, there is no requirement of true data to be fed
into the network in a time delayed fashion. One can predict for a finite time
horizon with a single true sample which is used for the initialization of the
network. We have also used closed loop recurrent neural network in the POD-RNN
for the problem of flow past plain cylinder and the results were satisfactory.
In that sense, one can argue that the convolutional recurrent autoencoder
model does not provide any significant improvements over the POD-RNN model for
this problem. However, the same inference is expected since the problem at
hand is the simplest of fluid flow problems and any standard prediction tool
should be able to perform at equal capabilities.

## Side-by-Side Cylinders Arrangement

The same closed loop recurrent neural network in the POD-RNN model failed in
the case of flow past a side-by-side cylinder arrangement. We have implemented
the encoder-decoder networks which takes a certain length of true data to
predict a finite time horizon [2]. The flow past a side-by-side cylinder
arrangement is known to be a more complex problem with a bi-stable solution,
flip-flopping regime and gap flow [15]. In the next section, the same
convolutional recurrent autoencoder model with closed loop recurrent neural
network is employed on the problem of flow past side-by-side cylinders.

Figure 12: Schematic diagram of the problem setup: the side-by-side cylinder
arrangement Figure 13: Normalized squared error $E_{f}$ for the velocity field
in X-direction ($U$): flow past side-by-side cylinders

The canonical problem of flow past side-by-side cylinders can act as a
representative case for many offshore multibody systems. The dynamics of the
problem is still a popular topic of research in the fluid mechanics community.
Therefore this problem has been chosen to test the capabilities of our data-
driven framework. The schematic of the problem set up is given in Fig. 12
where all the information about the domain boundaries and the boundary
conditions are clearly outlined.

The Reynolds number of the problem is considered as $Re=100$. Although it is
in the regime of low Reynolds number, the dynamics associated with this
problem are far more complex when compared to the flow past a single cylinder.
The full order simulation was carried out for $600s$ with a time step of
$0.02s$. A total of $3000$ snapshots of the simulation were collected for
every $0.2s$. Similar to the previous problem, we map the unstructured data
coming from full order simulations onto a square uniform grid of size
$64\times 64$. Both the velocity field $U$ and the pressure field $P$ are
considered for the training and prediction. The training data is generated
with $N_{s}=100$ samples, each with a size of $N_{t}=20$, which accounts for
first $2000$ samples. The next $600$ samples are used for the validation and
the final $400$ samples are used for testing the model.

Similar to the previous section, four models with different sizes of low
dimensional representation were trained. All the four models are trained on
single Intel E5-2690v3 (2.60GHz, 12 cores) CPU node for $N_{train}=1000000$
iterations. The training took approximately $48$ hours. The online prediction
is carried for 400 steps with the trained model. Figure 13 shows the
normalized squared error $E_{f}$ for the velocity field in $X$-direction $(U)$
and Fig. 14 shows the same error for pressure field $P$

In contrast to the case with an isolated plain cylinder, the velocity error
here for all the four models appears to be higher than the pressure error.
Also for the predictions of both $U$ and $P$, the errors are in the same order
of magnitude for all the four low dimensional representations. We have chosen
$N_{h}=64$ as the final trained model to depict the flow field predictions in
detail since it has a better performance when compared to other models.

Figure 14: Normalized squared error $E_{f}$ for pressure field ($P$): flow
past side-by-side cylinders

Figure 15 shows the comparison of truth and predicted values for the velocity
field in X-direction ($U$) along with the normalized reconstruction error
$E_{n}$ and Fig.16 illustrates the same comparison for the pressure field.
Overall predictions from the convolutional recurrent autoencoder model are
quite satisfactory. It is remarkable to note that the convolutional recurrent
autoencoder model is able to capture the flip-flopping phenomenon in both
velocity $U$ and pressure $P$ predictions. The predictions are also quite
smooth without any oscillations. The reconstruction error $E_{n}$ is quite low
in the order of $10^{-3}$ for both the predictions. Similar to the isolated
cylinder, the reconstruction errors are concentrated in the near wake region
where there is strong presence of nonlinearity.

(a)  (b)  (c)

Figure 15: Comparison of truth and predicted fields along with normalized
reconstruction error $E_{n}$ at (a) $t=540s$, (b) $t=560s$, (c) $t=580s$ for
velocity field in X-direction ($U$): flow past side-by-side cylinders

(a)  (b)  (c)

Figure 16: Comparison of truth and predicted fields along with normalized
reconstruction error $E_{n}$ at (a) $t=540s$, (b) $t=560s$, (c) $t=580s$ for
pressure field ($P$): flow past side-by-side cylinders

It should be noted that the POD-RNN model with the closed loop recurrent
neural network completely failed to predict the flow past an arrangement of
side-by-side cylinders. Whereas, the convolutional recurrent autoencoder model
was able to predict the flow fields up to 400 time steps satisfactorily
without any divergence in the error. It is also able to capture the randomness
in the flow field with respect to dominant flow features. In the case of POD-
RNN model, the spatial and temporal parts of the problem are dealt
independently which might work for simple problems of the flow past a plain
isolated cylinder. Such linear decomposition procedures completely fail at
highly non-linear problems such as side-by-side cylinder. The improvement
compared to the POD-RNN model can be attributed to the low dimensional
features obtained by convolutional neural networks of convolutional recurrent
autoencoder model. They are very powerful in identifying dominant local flow
features [18]. Also, the complete end to end architecture of the network which
enables us to integrate the encoding, evolution and decoding in a complete
non-linear fashion is also the primary reason for this model to work for this
problem. This is a very promising result and motivates us to take forward this
concept of convolutional recurrent autoencoder models for more complicated
problems of fluid mechanics and fluid-structure interaction.

### CONCLUSION

An overall framework for the data-driven reduced order model is presented in
this work. The current work can be considered as a extension of our previous
work on POD-RNN [2]. The shortcomings of the POD-RNN model are discussed and
the idea of completely bypassing the POD for obtaining the low dimensional
features is presented via convolutional recurrent autoencoder networks. The
proposed convolutional recurrent autoencoder networks serves as an end to end
nonlinear model reduction tool for unsteady flow predictions. To demonstrate
the effectiveness, the proposed model was applied on the problem of the flow
past an isolated cylinder and more complicated problem of the flow past side-
by-side cylinders. Pressure and velocity fields were predicted ahead in time
for a finite time horizon and the model performed satisfactorily in agreement
with the error metrics. We are currently working to extend this framework to
more complicated cases of free-surface wave prediction, turbulent flows and
fluid-structure interaction of realistic geometries.

## ACKNOWLEDGEMENT

This work done at the Keppel-NUS Corporate Laboratory was supported by the
National Research Foundation, Keppel Offshore and Marine Technology Centre
(KOMtech) and National University of Singapore (NUS). The conclusions put
forward to reflect the views of the authors alone and not necessarily those of
the institutions within the Corporate Laboratory. The computational work for
this article was (fully/partially) performed on resources of the National
Supercomputing Centre, Singapore (https://www.nscc.sg).

## References

  * [1] Brunton, S. L., Noack, B. R., and Koumoutsakos, P., 2019.  “Machine learning for fluid mechanics”.  arXiv:1905.11075.
  * [2] Reddy, S. B., Magee, A. R., Jaiman, R. K., Liu, J., Xu, W., Choudhary, A., and Hussain, A., 2019.  “Reduced order model for unsteady fluid flows via recurrent neural networks”.  In ASME 2019 38th International Conference on Ocean, Offshore and Arctic Engineering, American Society of Mechanical Engineers Digital Collection.
  * [3] Miyanawala, T. P., and Jaiman, R. K., 2019.  “Decomposition of wake dynamics in fluid–structure interaction via low-dimensional models”.  Journal of Fluid Mechanics,  867, pp. 723–764.
  * [4] Otto, S. E., and Rowley, C. W., 2019.  “Linearly recurrent autoencoder networks for learning dynamics”.  SIAM Journal on Applied Dynamical Systems,  18(1), pp. 558–593.
  * [5] Li, Q., Dietrich, F., Bollt, E. M., and Kevrekidis, I. G., 2017.  “Extended dynamic mode decomposition with dictionary learning: A data-driven adaptive spectral decomposition of the koopman operator”.  Chaos: An Interdisciplinary Journal of Nonlinear Science,  27(10), p. 103111.
  * [6] Kani, J. N., and Elsheikh, A. H., 2017.  “Dr-rnn: A deep residual recurrent neural network for model reduction”.  arXiv preprint arXiv:1709.00939.
  * [7] Wang, Z., Xiao, D., Fang, F., Govindan, R., Pain, C. C., and Guo, Y., 2018.  “Model identification of reduced order fluid dynamics systems using deep learning”.  International Journal for Numerical Methods in Fluids,  86(4), pp. 255–268.
  * [8] Gonzalez, F. J., and Balajewicz, M., 2018.  “Deep convolutional recurrent autoencoders for learning low-dimensional feature dynamics of fluid systems”.  arXiv preprint arXiv:1808.01346.
  * [9] Jaiman, R., Guan, M., and Miyanawala, T., 2016.  “Partitioned iterative and dynamic subgrid-scale methods for freely vibrating square-section structures at subcritical reynolds number”.  Computers & Fluids,  133, pp. 68–89.
  * [10] Jaiman, R., Pillalamarri, N., and Guan, M., 2016.  “A stable second-order partitioned iterative scheme for freely vibrating low-mass bluff bodies in a uniform flow”.  Computer Methods in Applied Mechanics and Engineering,  301, pp. 187–215.
  * [11] Donea, J., 1983.  “Arbitrary Lagrangian–Eulerian finite element methods”.  Computational Methods for Transient Analysis, pp. 474–516.
  * [12] Kuhl, E., Askes, H., and Steinmann, P., 2004.  “An ale formulation based on spatial and material settings of continuum mechanics. part 1: Generic hyperelastic formulation”.  Computer Methods in Applied Mechanics and Engineering,  193(39), pp. 4207–4222.
  * [13] Liu, J., Jaiman, R. K., and Gurugubelli, P. S., 2014.  “A stable second-order scheme for fluid–structure interaction with strong added-mass effects”.  Journal of Computational Physics,  270, pp. 687–710.
  * [14] Jaiman, R., Geubelle, P., Loth, E., and Jiao, X., 2011.  “Transient fluid–structure interaction with non-matching spatial and temporal discretizations”.  Computers & Fluids,  50(1), pp. 120–135.
  * [15] Liu, B., and Jaiman, R. K., 2016.  “Interaction dynamics of gap flow with vortex-induced vibration in side-by-side cylinder arrangement”.  Physics of Fluids,  28(12), p. 127103.
  * [16] Mysa, R. C., Kaboudian, A., and Jaiman, R. K., 2016.  “On the origin of wake-induced vibration in two tandem circular cylinders at low reynolds number”.  Journal of Fluids and Structures,  61, pp. 76–98.
  * [17] Kingma, D. P., and Ba, J., 2014.  “Adam: A method for stochastic optimization”.  arXiv preprint arXiv:1412.6980.
  * [18] Miyanawala, T., and Jaiman, R., 2018.  “Decomposition of wake dynamics in fluid-structure interaction”.  arXiv preprint.