# Construction of Lyapunov Functions Using Vector Field Decomposition

Yuanyuan Liu<EMAIL_ADDRESS>Northwestern Polytechnical
University, Xi’an 710072, China.

###### Abstract

In the present paper, a novel vector field decomposition based approach for
constructing Lyapunov functions is proposed. For a given dynamical system, if
the defining vector field admits a decomposition into two mutually orthogonal
vector fields, one of which is curl-free and the other is divergence-free,
then the potential function of the former can serve as a Lyapunov function
candidate, since its positive definiteness will reflect the stability of the
system. Moreover, under some additional conditions, its sublevel sets will
give the exact attraction domain of the system. A sufficient condition for the
existence of the proposed vector field decomposition is first obtained for
$2$-dimensional systems by solving a partial differential equation and then
generalized to $n$-dimensional systems. Furthermore, the proposed vector field
decomposition always exists for linear systems and can be obtained by solving
a specific algebraic Riccati equation.

###### keywords:

dynamical systems, Lyapunov functions, vector field decomposition,
differential forms, algebraic Riccati equation

††journal: arXiv

##  1 Introduction

As is well known, _Lyapunov functions_ (LFs) play a key role in the analysis
of stability of dynamical systems. Since the Russian mathematician A. M.
Lyapunov introduced the famous Lyapunov direct method in his seminal work
entitled _The General Problem of Stability of Motion_ [1], finding effective
methods for constructing LFs has become a central issue of modern control
theory.

The simplest and most commonly used methods may be the _variable gradient
method_ and _Krasovskii’s method_. The variable gradient method assumes a
certain form for the gradient of the unknown LF which usually depends on some
unknown parameters that need to be determined by the properties of a LF.
Although integrating the gradient may produce a LF, there is no unified
approach for the selection of the gradient. It relies heavily on experiences
and physical intuition. Krasovskii’s method is essentially a linearization
method, and hence the resulting LF usually can only be used in a neighborhood
(which may be very small) of the equilibrium point. Another kind of effective
method is the energy-based method, such as _Energy-Casimir method_ [2] and
_Hamiltonian function method_ [3, 4]. However, there exist some complex
dynamical systems to which none of the aforementioned methods can be applied
and one has to seek for numerical methods. Commonly used numerical methods
include the radial basis functions based method [5], the linear programming
based method [6], and the iterative based method [7]. Other kinds of numerical
methods can be found in [8, 9, 10, 11] and references therein. Numerical
methods are very useful in practice, though, they provide little insight to
the structure of dynamical systems.

Apart from stability, there is another issue that is of great importance in
the study of stability of dynamical systems, namely, the _domain of
attraction_ (DA). Since in practice, it may happen that a dynamical system is
stable but with very narrow stability margin, and becomes unstable once
disturbed even by slightest external disturbance. However, finding the exact
DA of a dynamical system is anything but easy. In fact, the only method
noticed by us that does give an accurate characterization of the DA is given
by Zubov [12]. But it has to solve a partial differential equation (known as
_Zubov’s equation_) subjecting to some additional conditions which is
difficult in general. Thus, one turns to estimate the DA instead. For example,
_Zubov’s Method_ [12] estimates the DA of a dynamical system by approximately
solving Zubov’s equation using the series expansion of the LF. Generalizations
of Zubov’s method to systems subject to perturbations or control inputs can be
found in [13, 14]. Traditionally, estimates of the DA are given by sublevel
sets of LFs. But such estimates are usually too conservative, justifying the
term _inner estimates_. Less conservative estimates can be found in [15] which
are given by positively invariant sets that are not necessarily sublevel sets
of a LF. Sometimes, however, _outer estimates_ containing the DA are more
useful in practice such as collision avoidance [16]. A method for finding
outer estimates of the DA can be found in [17], where the DA is computed by
solving an infinite-dimensional convex linear programming problem over the
space of measures. This method is then extended to obtain both inner and outer
estimates of the DA [16]. Computing both inner and outer estimates, compared
to just one or the other, enables assessing the tightness of the estimates
obtained and provides a valuable insight into achievable performance and/or
safety of a given constrained control system. For a detailed discussion on the
estimates of the DA, please refer to the book [18].

In this paper, we shall propose a novel vector field decomposition based
approach for constructing LFs. This approach stand a chance to be applied to a
large class of dynamical systems (depends on the existence of the desired
vector field decomposition). Also, it will provide some insights to the
structure of the motion of the system. In fact, if the desired vector field
decomposition does exist for a given dynamical system, then naturally there is
a corresponding LF candidate and the motion of the system can be decomposed
into two directions, one within the level surfaces of the LF candidate and the
other orthogonal to them. Moreover, under some additional conditions, the DA
of the dynamical system is fully characterized by a sublevel set of the LF. We
shall also give a sufficient condition for the existence of the proposed
vector filed decomposition for general nonlinear dynamical systems. This
condition is first obtained for the $2$-dimensional cases by solving a partial
differential equation and then generalized to the $n$-dimensional cases by
treating the defining vector field as a differential $1$-form. Finally, we
shall prove that for linear systems, the proposed vector field decomposition
always exists and can be found explicitly by solving a specific algebraic
Riccati equation.

The organization of the paper is as follows. The proposed vector field
decomposition based approach is presented in Section 2. A sufficient condition
for the existence of the proposed vector field decomposition for nonlinear
systems is given in Section 3. The existence of the proposed vector field
decomposition for linear systems is proved in Section 4. Finally, the
conclusions are given in Section 5.

##  2 A vector field decomposition based approach

Consider a general nonlinear dynamical system of the form

$\dot{x}=f(x),$ (1)

where $x\in\mathbb{R}^{n}$ is the system state vector and
$f:\mathbb{R}^{n}\to\mathbb{R}^{n}$ is a smooth map (the smoothness can be
relaxed but we assume it for simplicity). To see the motivation behind our
vector field decomposition based approach, observe that if the dynamical
system (1) admits a smooth LF, say $V$, then the map $f$, viewed as a vector
field (hereafter, called the defining vector field of the system (1)) on
$\mathbb{R}^{n}$, can be naturally decomposed into the sum of two vector
fields as

$f(x)=g(x)+h(x),\quad x\in\mathbb{R}^{n},$ (2)

where $-g(x)=\nabla V(x),x\in\mathbb{R}^{n}$ is the _gradient field_ (and
hence _curl-free_) of $V$. Moreover, we have, by the definition of LF,

$\left<\nabla V(x),f(x)\right>=-\left<g(x),g(x)+h(x)\right>\leq 0,\quad
x\in\mathbb{R}^{n},$ (3)

where $\left<\cdot,\cdot\right>$ denotes the usual Euclidean inner product.
Intuitively, equation (3) means that the sum vector field $f=g+h$ always point
‘inside’ the level surfaces of $V$, i.e., the trajectories of (1) always move
from higher level surfaces to lower ones, despite the vector field $h$ may
point ‘inside’, ‘outside’, or tangent to the level surfaces. It is noteworthy,
however, that in the critical case where the vector field $h$ tangent to the
level surfaces of $V$, we have

$\left<g(x),h(x)\right>=0,\quad x\in\mathbb{R}.$ (4)

and (3) always holds since $-\left<g(x),g(x)\right>$ is negative semidefinite.

The geometric meaning of (4) is that the vector field $g$ is pointwise
orthogonal to $h$ in $\mathbb{R}^{n}$. Therefore, we would like to decompose
the vector field $f$ into a sum of two vector fields as $f=g+h$ such that $g$
is curl-free and orthogonal to $h$ everywhere in $\mathbb{R}^{n}$, i.e.,
equation (4) holds. Moreover, in order to rule out the trivial solution
$g(x)\equiv 0$, we shall also require $h$ to be _divergence-free_ (which
implies that the vector fields $f$ and $g$ have the same divergence, a
property that is not satisfied by the trivial solution $g(x)\equiv 0$ in
general).

Suppose the defining vector field $f$ of system (1) has the proposed
decomposition, then there exists a unique smooth function
$V:\mathbb{R}^{n}\rightarrow\mathbb{R}$ satisfying $\nabla V(x)=-g(x)$ and
$V(0)=0$ such that $V$ is decreasing along the trajectories of system (1),
since

$\frac{\mathrm{d}}{\mathrm{d}t}V(x)=\nabla
V(x)f(x)=\left<-g(x),g(x)+h(x)\right>=-\left\|g(x)\right\|^{2}\leq
0,\quad\forall x\in\mathbb{R}^{n}.$

Thus, if $V(x)$ happens to be positive definite, then it is a LF of the
system. Actually, we have

######  Theorem 1.

Consider the dynamical system (1). Suppose the defining vector field $f$
admits a decomposition $f=g+h$ in a neighborhood $U\subset\mathbb{R}^{n}$ of
$0$ such that $g$ is curl-free, $h$ is divergence-free and
$\left<h(x),g(x)\right>=0$ for all $x\in U$. Let $V:U\rightarrow\mathbb{R}$ be
the smooth function satisfying $\nabla V(x)=-g(x)$ and $V(0)=0$.

  1. 1.

If $V(x)$ is positive definite on $U\backslash\left\\{0\right\\}$, then the
equilibrium point $x=0$ is Lyapunov stable.

  2. 2.

If in addition, $g(x)\neq 0$ for all $x\in U\backslash\\{0\\}$, then the
equilibrium point $x=0$ is asymptotically stable.

  3. 3.

Finally, if there exists a scalar $\varepsilon>0$ such that $g(x)\neq 0$ for
all $x\in\mathcal{B}_{\varepsilon}(0)\backslash\\{0\\}$, and for every
sufficiently small $\delta>0$ there exists $x_{0}\in\mathbb{R}^{n}$ such that
$\|x_{0}\|<\delta$ and $V(x_{0})<0$, then the equilibrium point $x=0$ is
unstable.

###### Proof.

1) and 2) follow directly from Lyapunov’s stability theorem (see [2], Theorem
3.1), while 3) follows directly from Lyapunov’s instability theorem (see [2]
Theorem 3.12). ∎

It seems that the proposed scheme not only complicates the searching for a
scalar function $V$ satisfying Lyapunov’s stability theorem, but also strongly
restricts the set of systems and LFs the scheme works for. This issue is clear
for the case of linear systems studied in Section 4 where the method reduces
to solve a restricted case of the algebraic Lyapunov equation. Nevertheless,
it has the virtue that it gives an accurate characterization of the DA.
Indeed, we have

######  Theorem 2.

Assume the defining vector field $f$ of (1) admits a decomposition $f=g+h$ as
in Theorem 1 and $V:\mathbb{R}^{n}\rightarrow\mathbb{R}$ be the smooth
function satisfying $\nabla V(x)=-g(x)$ and $V(0)=0$. Let $c\in\mathbb{R}$ be
positive and $\Omega_{c}$ the path component of $\\{x:V(x)<c\\}$ containing
the equilibrium point $x=0$. If the closure $\bar{\Omega}_{c}\subset U$,

$\begin{split}g(x)&=0,\quad x\in\partial\Omega_{c},\\\ g(x)&\neq 0,\quad
x\in\Omega_{c}\backslash\\{0\\},\\\ V(x)&>0,\quad
x\in\Omega_{c}\backslash\\{0\\},\end{split}$

and $V$ is radially unbounded, then $\Omega_{c}$ is the DA of system (1).

###### Proof.

First of all, the equilibrium point $x=0$ is asymptotically stable according
to Theorem 1. Let $\mathcal{D}_{a}$ be the DA of system (1). Since $V(x)$ is
positive definite on $\Omega_{c}\backslash\\{0\\}$, there exists
$\varepsilon>0$ such that
$\Omega_{\varepsilon}\coloneqq\\{x\in\Omega_{c}:V(x)<\varepsilon\\}\subseteq\mathcal{D}_{a}$.
To prove that $\Omega_{c}=\mathcal{D}_{a}$, we first show that
$\Omega_{c}\subseteq\mathcal{D}_{a}$. To this end, it suffices to show that
for each $\varepsilon\leq c^{\prime}<c$ and
$x_{0}\in\bar{\Omega}_{c^{\prime}}\backslash\Omega_{\varepsilon}$, the
trajectory $x(t)$ issued from $x_{0}$, i.e., that satisfies $x(0)=x_{0}$, will
enter $\Omega_{\varepsilon}$ in a finite time $t$. By the radially
unboundedness of $V$, the set
$\bar{\Omega}_{c^{\prime}}\backslash\Omega_{\varepsilon}$ is compact. Thus,
there exists $\delta>0$ such that $\|g(x)\|^{2}>\delta$ for all
$x\in\bar{\Omega}_{c^{\prime}}\backslash\Omega_{\varepsilon}$, since $g(x)\neq
0$ on this set. Consequently, $x(t)$ will enter $\Omega_{\varepsilon}$ at
$t=(c^{\prime}-\varepsilon)/\delta$ since we have

$\frac{\mathrm{d}}{\mathrm{d}t}V(x)=\left<\nabla
V(x),f(x)\right>=-\|g(x)\|^{2}<-\delta$

and $V(x)\leq c^{\prime}$ in
$\bar{\Omega}_{c^{\prime}}\backslash\Omega_{\varepsilon}$.

Next, we prove $\mathcal{D}_{a}\subseteq\Omega_{c}$. This can be achieved by
showing that for each $x_{0}\in U\backslash\Omega_{c}$ the trajectory $x(t)$
issued from $x_{0}$ cannot across the boundary $\partial\Omega_{c}$ of
$\Omega_{c}$ in any finite time. Suppose on the contrary there does exist a
finite $t>0$ such that $x(t)$ crosses the boundary $\partial\Omega_{c}$ into
$\Omega_{c}$, then the vector field $f$ intersects $\partial\Omega_{c}$
transversally at $x(t)$. Since $f=g+h$ and $g(x)=0$ at $x(t)$, we conclude
that the vector field $h$ intersects $\partial\Omega_{c}$ transversally at
$x(t)$. But then for a sufficiently small $\varepsilon>0$, the vector field
$h$ will intersects the lever surface $\\{x\in U:V(x)=c^{\prime}\\}$
transversally at $x(t+\varepsilon)$, where $c^{\prime}=V(x(t+\varepsilon))<c$.
This violates the orthogonality of vector fields $g$ and $h$, a contradiction.
∎

######  Example 1.

Consider the following nonlinear dynamical system adopted from [19],

$\begin{split}\dot{x}_{1}=-x_{1}[1+(x_{1}^{2}+x_{2}^{2})-(x_{1}^{2}+x_{2}^{2})^{2}]+x_{2},\\\
\dot{x}_{2}=-x_{2}[1+(x_{1}^{2}+x_{2}^{2})-(x_{1}^{2}+x_{2}^{2})^{2}]-x_{1}.\end{split}$

The defining vector field $f$ of it has the desired decomposition $f=g+h$ with

$g=-[1+(x_{1}^{2}+x_{2}^{2})-(x_{1}^{2}+x_{2}^{2})^{2}](x_{1},x_{2})^{T}$

curl-free and $h=(x_{2},-x_{1})^{T}$ divergence-free. Thus, by Theorem 1, the
function

$V(x)=\frac{1}{2}(x_{1}^{2}+x_{2}^{2})\left[1+\frac{1}{2}(x_{1}^{2}+x_{2}^{2})-\frac{1}{3}(x_{1}^{2}+x_{2}^{2})^{2}\right]$

is a LF of the system and the equilibrium point $x=0$ is asymptotically
stable. Let
$U=\\{x\in\mathbb{R}^{2}:x_{1}^{2}+x_{2}^{2}<\frac{3+\sqrt{57}}{4}\\}$ and
$\Omega=\\{x\in\mathbb{R}^{2}:x_{1}^{2}+x_{2}^{2}<\frac{1+\sqrt{5}}{2}\\}$.
Then $V$ is positive definite on $U\backslash\\{0\\}$ and $g(x)=0$ for
$x\in\partial\Omega$. Moreover, $\Omega$ is the path component of a sublevel
set of $V$ that contains $x=0$, and the conditions of Theorem 2 are satisfied.
Therefore, $\Omega$ is the DA of the system. This can be seen via the phase
portrait of the system.

At this point, a natural question is whether the decomposition required by
Theorem 1 and 2 exists. It is a well known result in physics (especially in
electromagnetics) that one can always decompose a vector field into the sum of
a divergence-free and a curl-free vector field, and such decomposition is not
unique. But with the additional orthogonality condition (4), the question
turns out to be much harder. Nevertheless, we shall give a sufficient
condition for the existence of the proposed decomposition in the next section.
Furthermore, a complete answer for linear systems will be given in Section 4.

##  3 A sufficient condition for nonlinear systems

In this section, we give a sufficient condition for the existence of the
proposed vector field decomposition for general nonlinear systems. To this
end, we first address the simpler $2$-dimensional cases and then generalize
the result to $n$-dimensional cases.

###  3.1 $2$-dimensional cases

Given a smooth vector field $f$ on $\mathbb{R}^{2}$, suppose it can be
decomposed into the sum of two smooth vector fields as $f=g+h$, where $g$ is
curl-free, $h$ is divergence-free, and $\left<g(x),h(x)\right>=0$ for all
$x\in\mathbb{R}^{2}$. Then, by the irrotationality of $g$, i.e.,

$\frac{\partial g_{1}}{\partial x^{2}}=\frac{\partial g_{2}}{\partial x^{1}},$

there exists (see the Poincaré Lemma in [20]) a scaler function
$\Phi:\mathbb{R}^{2}\to\mathbb{R}$ such that

$g_{1}=\frac{\partial\Phi}{\partial x^{1}}\quad\text{and}\quad
g_{2}=\frac{\partial\Phi}{\partial x^{2}},$ (5)

where $g_{i}$ and $x^{i},i=1,2$ denote the $i$-th component of $g$ and $x$,
respectively. Similarly, there exists a scalar function
$\Psi:\mathbb{R}^{2}\rightarrow\mathbb{R}$ such that

$h_{1}=-\frac{\partial\Psi}{\partial x^{2}}\quad\text{and}\quad
h_{2}=\frac{\partial\Psi}{\partial x^{1}},$ (6)

since $h$ is divergence-free. Plugging (5) and (6) into $\left<g,h\right>=0$
we find that

$-\frac{\partial\Phi}{\partial x^{1}}\frac{\partial\Psi}{\partial
x^{2}}+\frac{\partial\Phi}{\partial x^{2}}\frac{\partial\Psi}{\partial
x^{1}}=0,\quad x\in\mathbb{R}^{2}.$ (7)

If the vector field $g$ is nondegenerate on an open set
$U\subset\mathbb{R}^{2}$, then considered as a partial differential equation
of $\Psi$, the classical differential equation theory tells us that the
general solution of equation (7) can be found by integrating the corresponding
characteristic system of ordinary differential equation (see [21, 22] for
details), i.e.,

$\frac{dx^{1}}{{\partial\Phi}/{\partial
x^{2}}}=\frac{dx^{2}}{-{\partial\Phi}/{\partial x^{1}}},\quad x\in U.$ (8)

Equation (8) can be solved explicitly and the solution is

$\Phi(x)=c,\quad x\in U,$

where $c$ is the constant of integration. And any solution of (7) can be
written as a function of $c$, i.e., $\Psi(x)=\alpha(c)=\alpha(\Phi(x))$ for a
smooth function $\alpha:\mathbb{R}\rightarrow\mathbb{R}$. Therefore, we have

$\begin{split}f_{1}&=\frac{\partial\Phi}{\partial
x^{1}}-\alpha^{\prime}(\Phi)\frac{\partial\Phi}{\partial x^{2}},\\\
f_{2}&=\frac{\partial\Phi}{\partial
x^{2}}+\alpha^{\prime}(\Phi)\frac{\partial\Phi}{\partial x^{1}},\\\
\end{split}$ (9)

where $\alpha^{\prime}$ denotes the derivative of $\alpha$ with respect to
$\Phi$. Clearly, if a vector field $f$ can be written in the form of (9), then
it has the desired decomposition.

On the other hand, if $h$ is nondegenerate on an open set $U$, then by the
same argument we will find that $\Phi=\beta(\Psi)$ on $U$ for a smooth
function $\beta:\mathbb{R}\to\mathbb{R}$ and consequently, we have

$\begin{split}f_{1}&=\beta^{\prime}(\Psi)\frac{\partial\Psi}{\partial
x^{1}}-\frac{\partial\Psi}{\partial x^{2}},\\\
f_{2}&=\beta^{\prime}(\Psi)\frac{\partial\Psi}{\partial
x^{2}}+\frac{\partial\Psi}{\partial x^{1}}.\end{split}$ (10)

Combining equations (9) and (10) yields the following sufficient condition for
the existence of the proposed vector field decomposition.

######  Theorem 3.

Let $f=(f_{1},f_{2})^{T}$ be a smooth vector field on $\mathbb{R}^{2}$. If
there exist smooth functions $\Theta:\mathbb{R}^{2}\rightarrow\mathbb{R}$ and
$\alpha,\beta:\mathbb{R}\rightarrow\mathbb{R}$ such that

$\begin{split}f_{1}&=\beta(\Theta)\frac{\partial\Theta}{\partial
x^{1}}-\alpha(\Theta)\frac{\partial\Theta}{\partial x^{2}},\\\
f_{2}&=\beta(\Theta)\frac{\partial\Theta}{\partial
x^{2}}+\alpha(\Theta)\frac{\partial\Theta}{\partial x^{1}}\end{split}$ (11)

then $f$ admits the proposed decomposition $f=g+h$ with
$g=\beta(\Theta)(\frac{\partial\Theta}{\partial
x^{1}},\frac{\partial\Theta}{\partial x^{2}})^{T}$ and
$h=\alpha(\Theta)(-\frac{\partial\Theta}{\partial
x^{2}},\frac{\partial\Theta}{\partial x^{1}})^{T}$.

For instance, the defining vector field of the nonlinear system given in
Example 2 satisfies the conditions of Theorem 3, and hence admits the proposed
decomposition as we have seen before.

######  Remark 1.

It is clear that any smooth vector field $f$ on $\mathbb{R}^{2}$ can be
written in the form

$\begin{split}f_{1}&=\beta(x)\frac{\partial\Theta}{\partial
x^{1}}-\alpha(x)\frac{\partial\Theta}{\partial x^{2}},\\\
f_{2}&=\beta(x)\frac{\partial\Theta}{\partial
x^{2}}+\alpha(x)\frac{\partial\Theta}{\partial x^{1}},\\\ \end{split}$ (12)

for some smooth functions $\Theta:\mathbb{R}^{2}\to\mathbb{R}$ and
$\alpha,\beta:\mathbb{R}^{2}\to\mathbb{R}$. What Theorem 3 tells us is that if
the function $\Theta$ can be chosen in a manner such that $\alpha$ and $\beta$
are constant on each level set of $\Theta$, then the existence of the proposed
decomposition is guaranteed. The question is, under this restriction, whether
equation (12) is still sufficient to represent all vector fields on
$\mathbb{R}^{2}$.

###  3.2 $n$-dimensional cases

In order to generalize Theorem 3 to $n$-dimensional cases, we have to
introduce some more notations. For high dimensional cases, it is more
convenient to treat a vector filed $f$ on $\mathbb{R}^{n}$ as a _differential
$1$-form_, i.e.,

$f=\sum_{i=1}^{n}f_{i}\mathrm{d}x^{i}.$

Under this new notation, the curl of $f$ can be written as

$\mathrm{curl}f=\mathrm{d}f=\sum_{i=1}^{n}\sum_{j=1}^{n}\frac{\partial
f_{i}}{\partial x^{j}}\mathrm{d}x^{j}\wedge\mathrm{d}x^{i}=\sum_{1\leq i<j\leq
n}\left(\frac{\partial f_{i}}{\partial x^{j}}-\frac{\partial f_{j}}{\partial
x^{i}}\right)\mathrm{d}x^{j}\wedge\mathrm{d}x^{i},$ (13)

where $\mathrm{d}$ is the _exterior derivative operator_. If $f$ is curl-free,
i.e., $\mathrm{d}f=0$, then there exist a $0$-form (i.e., a scaler function)
$\Phi$ such that $f=\mathrm{d}\Phi$.

In addition, let $\mathrm{d}^{\ast}$ be the $L^{2}$-adjoint of the exterior
derivative $\mathrm{d}$, then the divergence of $f$ can be written as

$\begin{split}\mathrm{div}f&=-\mathrm{d}^{\ast}f=\ast\mathrm{d}\ast f\\\
&=\ast\mathrm{d}\sum_{i=1}^{n}(-1)^{i-1}f_{i}\mathrm{d}x^{1}\wedge\cdots\wedge\widehat{\mathrm{d}x^{i}}\wedge\cdots\wedge\mathrm{d}x^{n}\\\
&=\ast\sum_{i=1}^{n}\frac{\partial f_{i}}{\partial
x^{i}}\mathrm{d}x^{1}\wedge\cdots\wedge\mathrm{d}x^{n}\\\
&=\sum_{i=1}^{n}\frac{\partial f_{i}}{\partial x^{i}},\end{split}$ (14)

where $\ast$ is the _Hodge star operator_ (the precise definition of
$\mathrm{d}^{\ast}$ and $\ast$ as well as their relationship can be found in
[23]). If $f$ is divergence free, i.e., $\mathrm{d}^{\ast}f=0$, then there
exists a $2$-form

$\psi=\sum_{1\leq i<j\leq n}\Psi_{ij}\mathrm{d}x^{i}\wedge\mathrm{d}x^{j}$

such that $f=\mathrm{d}^{\ast}\psi$. Now we claim:

######  Theorem 4.

Let $f=\sum_{i=1}^{n}f_{i}\mathrm{d}x^{i}$ be a vector field on
$\mathbb{R}^{n}$, if there exist smooth functions
$\Theta:\mathbb{R}^{n}\rightarrow\mathbb{R}$,
$\alpha_{ij}:\mathbb{R}\rightarrow\mathbb{R}$, and
$\beta:\mathbb{R}\to\mathbb{R}$ such that

$f=\mathrm{d}\Phi+\mathrm{d}^{\ast}\psi$

hold with

$\Phi=\beta(\Theta)\quad\text{and}\quad\Psi_{ij}=\alpha_{ij}(\Theta),$ (15)

then $f$ admits the proposed decomposition $f=g+h$ with $g=\mathrm{d}\Phi$ and
$h=\mathrm{d}^{\ast}\psi$.

###### Proof.

Apparently, $g$ is divergence-free and $f$ is curl-free, since
$\mathrm{d}\circ\mathrm{d}=0$ and $\mathrm{d}^{\ast}\circ\mathrm{d}^{\ast}=0$.
Moreover, we have

$g=\mathrm{d}\Phi=\sum_{i=1}^{n}\frac{\partial\Phi}{\partial
x^{i}}\mathrm{d}x^{i}$

and

$\begin{split}h=&\mathrm{d}^{\ast}\psi=(-1)^{n(2+1)+1}\ast\mathrm{d}\ast\sum_{1\leq
i<j\leq n}\Psi_{ij}\mathrm{d}x^{i}\wedge\mathrm{d}x^{j}\\\
=&(-1)^{n+1}\ast\mathrm{d}\sum_{1\leq i<j\leq
n}(-1)^{i+j-1}\Psi_{ij}\mathrm{d}x^{1}\wedge\cdots\wedge\widehat{\mathrm{d}x^{i}}\wedge\cdots\wedge\widehat{\mathrm{d}x^{j}}\wedge\cdots\wedge\mathrm{d}x^{n}\\\
=&(-1)^{n+1}\ast\sum_{1\leq i<j\leq
n}(-1)^{j}\frac{\partial\Psi_{ij}}{\partial
x^{i}}\mathrm{d}x^{1}\wedge\cdots\wedge\widehat{\mathrm{d}x^{j}}\wedge\cdots\wedge\mathrm{d}x^{n}\\\
&+(-1)^{n+1}\ast\sum_{1\leq i<j\leq
n}(-1)^{i-1}\frac{\partial\Psi_{ij}}{\partial
x^{j}}\mathrm{d}x^{1}\wedge\cdots\wedge\widehat{\mathrm{d}x^{i}}\wedge\cdots\wedge\mathrm{d}x^{n}\\\
=&(-1)^{n+1}\sum_{1\leq i<j\leq n}(-1)^{n}\frac{\partial\Psi_{ij}}{\partial
x^{i}}\mathrm{d}x^{j}-\sum_{1\leq i<j\leq
n}(-1)^{n}\frac{\partial\Psi_{ij}}{\partial x^{j}}\mathrm{d}x^{i}\\\
=&\sum_{1\leq i<j\leq n}\left(\frac{\partial\Psi_{ij}}{\partial
x^{j}}\mathrm{d}x^{i}-\frac{\partial\Psi_{ij}}{\partial
x^{i}}\mathrm{d}x^{j}\right).\end{split}$

Therefore,

$\begin{split}\left<g,h\right>&=\left<\sum_{k=1}^{n}\frac{\partial\Phi}{\partial
x^{k}}\mathrm{d}x^{k},\sum_{1\leq i<j\leq
n}\left(\frac{\partial\Psi_{ij}}{\partial
x^{j}}\mathrm{d}x^{i}-\frac{\partial\Psi_{ij}}{\partial
x^{i}}\mathrm{d}x^{j}\right)\right>\\\ &=\sum_{k=1}^{n}\sum_{1\leq i<j\leq
n}\frac{\partial\Phi}{\partial x^{k}}\left(\frac{\partial\Psi_{ij}}{\partial
x^{j}}\delta^{ki}-\frac{\partial\Psi_{ij}}{\partial
x^{i}}\delta^{kj}\right)\\\ &=\sum_{1\leq i<j\leq
n}\left(\frac{\partial\Psi_{ij}}{\partial x^{j}}\frac{\partial\Phi}{\partial
x^{i}}-\frac{\partial\Psi_{ij}}{\partial x^{i}}\frac{\partial\Phi}{\partial
x^{j}}\right)\\\ &=\sum_{1\leq i<j\leq
n}\frac{\mathrm{d}\alpha_{ij}}{\mathrm{d}\Theta}\frac{\mathrm{d}\beta}{\mathrm{d}\Theta}\left(\frac{\partial\Theta}{\partial
x^{j}}\frac{\partial\Theta}{\partial x^{i}}-\frac{\partial\Theta}{\partial
x^{i}}\frac{\partial\Theta}{\partial x^{j}}\right)\\\ &=0,\end{split}$

wehre in the second last equation we have used the condition (15). Thus,
$f=g+h$ is the desired vector field decomposition. ∎

######  Remark 2.

The remarkable _Hodge decomposition theorem_ (see Corollary 3.4.2 in [23])
says that any vector field can be decomposed into the sum of a curl-free, a
divergence-free, and a harmonic vector field, such that they are orthogonal to
each other with respect to the $L^{2}$ inner product. One can hardly overlook
the intimate connection between the proposed decomposition and the Hodge
decomposition.

##  4 Existence for linear systems

To prove the existence for linear systems, we shall first recall some basics
about algebraic Riccati equation. The contents are mainly based on the book
[24] of I. Gohberg, P. Lancaster and L. Rodman. A general _continuous
algebraic Riccati equation_ (CARE) is a quadratic matrix equation of the form

$XDX+XA+BX-C=0,$ (16)

where $D,A,B,C$ are real or complex matrices of sizes $n\times m,n\times
n,m\times m$, and $m\times n$, respectively. Note that the coefficients of
CARE (16) fit into a square matrix of size $(n+m)\times(n+m)$:

$T=\begin{pmatrix}A&D\\\ C&-B\\\ \end{pmatrix}.$

The solutions of (16) can be fully characterized by the Jordan chains of $T$.

######  Definition 1 ([24]).

Let $T$ be an $n\times n$ matrix and $\lambda$ an eigenvalue of $T$, if there
exist $\eta_{i}\in\mathbb{C}^{n}$ satisfying $\eta_{i}\neq 0,i=1,\cdots,l$,
such that

$T\eta_{1}=\lambda\eta_{1},\ \ \text{and }\ \
T\eta_{i}=\lambda\eta_{i}+\eta_{i-1},\quad i=2,\cdots,l,$

then the ordered $l$-tuple $(\eta_{1},\cdots,\eta_{l})$ is called a _Jordan
chain_ of $T$ corresponding to the eigenvalue $\lambda$.

The following constructive theorem then gives the solutions of CARE (16).

######  Theorem 5 ([24]).

Equation (16) has a solution $X\in\mathbb{C}^{m\times n}$ if and only if there
exist vectors $v_{1},\cdots,v_{n}\in\mathbb{C}^{m+n}$ that form a set of
Jordan chains for $T$ and, if

$v_{j}=\begin{pmatrix}y_{j}\\\ z_{j}\\\ \end{pmatrix},\quad j=1,2,\cdots,$

where $y_{j}\in\mathbb{C}^{n}$, then $y_{1},\cdots,y_{n}$ form a basis for
$\mathbb{C}^{n}$. Furthermore, if

$Y=(y_{1},\cdots,y_{n}),\quad Z=(z_{1},\cdots,z_{n}),$

then every solution of (16) has the form $X=ZY^{-1}$ for some set of Jordan
chains $v_{1},\cdots,v_{n}$ for $T$ such that $Y$ is nonsingular.

Now we can prove the existence of the proposed vector field decomposition for
linear systems. Note that for linear systems we have $f(x)=Fx$ for a certain
$n\times n$ matrix $F\in\mathbb{R}^{n\times n}$. Let $g(x)=Gx$ for some
$G\in\mathbb{R}^{n\times n}$, then by the irrotationality of $g$ and equation
(13), we find that $G_{ij}=G_{ji}$, i.e., $G$ is a symmetric matrix. On the
other hand, according to (14) we find that

$\mathrm{tr}G=\mathrm{tr}F,$ (17)

since $h=f-g$ is divergence free. Moreover, the orthogonality condition (4)
becomes to $x^{T}G^{T}(F-G)x=0$ which implies that $G^{T}(F-G)$ is a skew-
symmetric matrix, that is $G^{T}(F-G)=-(F-G)^{T}G$, or equivalently

$GF+F^{T}G=2G^{2},$ (18)

where we have utilized the symmetry of matrix $G$. So we are going to find a
symmetric solution of (18) subjecting to constraint (17).

Note that equation (18) is a special form of the algebraic Riccati equation
(16). Indeed, if we take $m=n$, $D=2I,X=G,A=-F,B=-F^{T}$, and $C=0$, then
equation (16) becomes to equation (18). Therefore, all the solutions of (18)
can be found by Theorem 5. Among these solutions, we have to prove that there
exists a symmetric one such that constraint (17) is fulfilled. To begin with,
let us look at a simple illustrative example.

######  Example 2.

Let $D=2I_{2\times 2},C=0_{2\times 2},A=\begin{pmatrix}1&1\\\ 0&0\\\
\end{pmatrix}$, and $B=A^{T}$, then $T=PJP^{-1}$ with

$J=\begin{pmatrix}0&1&0&0\\\ 0&0&0&0\\\ 0&0&-1&0\\\ 0&0&0&1\\\
\end{pmatrix},\quad\text{and}\quad P=\begin{pmatrix}2&-2&2&0\\\ -2&0&0&2\\\
0&0&0&1\\\ 0&-1&0&1\\\ \end{pmatrix},$

be the Jordan normal form and the corresponding transformation matrix of $T$,
respectively. There are $4$ sets of Jordan chains for $T$ satisfying the
conditions of Theorem 5, i.e.,
$(P_{1},P_{2}),(P_{1},P_{3}),(P_{1},P_{4}),(P_{3},P_{4})$, where $P_{i}$
denotes the $i$-th column of $P$. Hence, there are four solutions for equation
(16), i.e.,

$X_{1}=\begin{pmatrix}0&0\\\ \frac{1}{2}&\frac{1}{2}\\\ \end{pmatrix},\quad
X_{2}=\begin{pmatrix}0&0\\\ 0&0\\\ \end{pmatrix},\quad
X_{3}=\begin{pmatrix}\frac{1}{2}&\frac{1}{2}\\\ \frac{1}{2}&\frac{1}{2}\\\
\end{pmatrix},\quad\text{and}\quad X_{4}=\begin{pmatrix}0&\frac{1}{2}\\\
0&\frac{1}{2}\\\ \end{pmatrix}.$

Observe that solution $X_{3}$ is the one with the same trace as $A$, and it
corresponds to the eigenvalues $0$ and $1$ of $T$ which are exactly the
eigenvalues of $A$.

Indeed, the observation is justified by the following lemma.

######  Lemma 6.

Let $X$ be a solution of equation (16) and
$v_{1},\cdots,v_{n}\in\mathbb{C}^{m+n}$ the corresponding set of Jordan chains
for $T$ satisfying the conditions of Theorem 5. If
$\lambda_{1},\cdots,\lambda_{n}$ are the eigenvalues of $T$ corresponding to
$v_{1},\cdots,v_{n}$, then we have

$\mathrm{tr}\left(A+DX\right)=\sum_{i=1}^{n}{\lambda_{i}}.$

###### Proof.

Since $X$ is a solution of (16), we have, according to Proposition 14.4.1 in
[24], that the equation

$T\begin{pmatrix}I\\\ X\\\ \end{pmatrix}=\begin{pmatrix}I\\\ X\\\
\end{pmatrix}Q$

holds with $Q=A+DX$. Let $J$ be the Jordan normal form of $T$ and $P$ the
corresponding transformation matrix, i.e. $T=PJP^{-1}$. Then, multiplying both
sides of the above equation by $P^{-1}$ on the left yields

$JP^{-1}\begin{pmatrix}I\\\ X\\\ \end{pmatrix}=P^{-1}\begin{pmatrix}I\\\ X\\\
\end{pmatrix}Q.$

From Theorem 5, we know that

$\begin{pmatrix}I\\\ X\\\ \end{pmatrix}=vY^{-1},$

where $v=\left(v_{1},\cdots,v_{n}\right)$ and $Y$ is the same as defined in
Theorem 5. Therefore,

$JP^{-1}v=P^{-1}v\tilde{Q}$ (19)

with $\tilde{Q}=Y^{-1}QY$. Since $v_{1},\cdots,v_{n}$ form a set of Jordan
chains for $T$, there exist $1\leq i_{1},\cdots,i_{n}\leq m+n$ such that
$v_{j}=P_{i_{j}},1\leq j\leq n$, where $P_{i_{j}}$ is the $i_{j}$-th column of
matrix $P$, and $i_{j}\neq i_{k}$ whenever $j\neq k$. Then we have
$P^{-1}v=\left(e_{i_{1}},\cdots,e_{i_{n}}\right)$, where
$e_{1},\cdots,e_{m+n}$ are the standard bases of $\mathbb{R}^{m+n}$.
Consequently, equation (19) is equivalent to $\tilde{J}=\tilde{Q}$, where
$\tilde{J}=\left(J_{i_{k}i_{l}}\right),1\leq k,l,\leq n$ is the matrix
consists of entries of $J$ at $i_{k}$-th row and $i_{l}$-th column. It is
clear that $\tilde{J}$ is itself in Jordan normal form since
$v_{1},\cdots,v_{n}$ form a set of Jordan chains of $T$. Moreover, the
diagonal entries of $\tilde{J}$ are exactly the eigenvalues
$\lambda_{1},\cdots,\lambda_{n}$, hence
$\mathrm{tr}\left(A+DX\right)=\mathrm{tr}\tilde{Q}=\mathrm{tr}\tilde{J}=\sum_{i=1}^{n}{\lambda_{i}}$,
as expected. ∎

######  Remark 3.

In the special case of (18), the eigenvalues of $F$ are also eigenvalues of
$T$. If corresponding to them there exists a solution $G$, then the lemma says
that $\mathrm{tr}(-F+2G)=\mathrm{tr}F$ which is exactly (17).

Now, we prove the existence of the proposed vector field decomposition for a
special class of linear systems where $F$ satisfies some additional
conditions.

######  Theorem 7.

Let $F\in\mathbb{R}^{n\times n}$ be a real matrix and
$\lambda_{1},\cdots,\lambda_{n}$ the eigenvalues of it. If
$\lambda_{i}+\lambda_{j}\neq 0$ for any $1\leq i,j\leq n$, then equation (18)
always has a real symmetric solution $G$ such that
$\mathrm{tr}G=\mathrm{tr}F$.

###### Proof.

Let $-F=PJP^{-1}$ and $F^{T}=UKU^{-1}$, where $J,K$ are Jordan normal forms
and $P,U$ the corresponding transformation matrices of $-F$ and $F^{T}$,
respectively. Since $\lambda_{i}+\lambda_{j}\neq 0$ for $1\leq i,j\leq n$,
$-F$ and $F^{T}$ have no eigenvalue in common. Thus, there is an invertible
matrix $W$ such that $T=WLW^{-1}$ with

$T=\begin{pmatrix}-F&2I\\\ 0&F^{T}\\\ \end{pmatrix}\quad\text{and}\quad
L=\begin{pmatrix}K&0\\\ 0&J\end{pmatrix}$

a Jordan normal form of $T$. Indeed, we can choose

$W=\begin{pmatrix}W_{11}&W_{12}\\\ W_{21}&W_{22}\end{pmatrix}$

with $W_{11}=XU,W_{12}=P,W_{21}=U,W_{22}=0$, and $X$ the unique (See [25, 26]
for the existence and uniqueness.) solution of the _Sylvester equation_

$XF^{T}+FX=2I.$ (20)

We claim that the solution $X$ is invertible. Indeed, $X$ is symmetric since
apparently $X^{T}$ is also a solution of (20). Thus, (20) becomes to

$(FX)^{T}+FX=2I,$

from which we conclude that $FX$ and consequently $X$ is invertible. But then
$W_{11}$ is also invertible and the condition of Theorem 5 is satisfied with
$Y=W_{11}$ and $Z=W_{21}=U$. Therefore $G=ZY^{-1}=X^{-1}$ is a solution of
(18). Note that $\mathrm{tr}G=\mathrm{tr}F$ (See Remark 3) and $G$ is
symmetric since $X$ is. And the reality of $G$ is guaranteed by Theorem
14.10.2 in [24]. So the proof is complete. ∎

Based on this result, we show that the proposed vector field decomposition
exists for arbitrary linear system.

######  Theorem 8.

For any real matrix $F\in\mathbb{R}^{n\times n}$, equation (18) always has a
real symmetric solution $G$ such that $\mathrm{tr}F=\mathrm{tr}G$.

###### Proof.

Suppose (18) has a real symmetric solution $G$ such that
$\mathrm{tr}F=\mathrm{tr}G$. Then there exists an orthogonal matrix $Q$ such
that $G=Q\Lambda Q^{T}$, where $\Lambda$ is a diagonal matrix such that
$\mathrm{tr}\Lambda=\mathrm{tr}F$. Plugging this into (18) and multiplying
both sides on the left by $Q^{T}$ and right by $Q$ yield

$\Lambda Q^{T}FQ+Q^{T}F^{T}Q\Lambda=2\Lambda^{2},$

which implies that

$\Lambda Q^{T}FQ=\Lambda^{2}+B$ (21)

for a certain skew-symmetric matrix $B$.

Conversely, if there exist an orthogonal matrix $Q$, a skew-symmetric matrix
$B$, and a diagonal matrix $\Lambda$ satisfying
$\mathrm{tr}\Lambda=\mathrm{tr}F$ such that equation (21) holds, then (18) is
satisfied with $G=Q\Lambda Q^{T}$. Thus, to prove the theorem, it suffices to
show that for any matrix $F\in\mathbb{R}^{n\times n}$, such matrices $Q,B$ and
$\Lambda$ do exist. To this end, we assume the set of eigenvalues of $F$ is
$\sigma=\\{\lambda_{1},\cdots,\lambda_{n}\\}$ and construct two subsets
$\sigma_{1}$ and $\sigma_{2}$ of it iteratively as follows.

  1. 1.

Let $\sigma_{1}=\sigma$ and $\sigma_{2}=\emptyset$.

  2. 2.

If we can find a pair of eigenvalues $\lambda_{i}$ and $\lambda_{j}$ ($i,j$
may be equal) in $\sigma_{1}$ such that $\lambda_{i}+\lambda_{j}=0$, then go
to step 3. Else, go to step 4.

  3. 3.

Remove $\lambda_{i}$ and $\lambda_{j}$ from $\sigma_{1}$ and put them into
$\sigma_{2}$, return to Step 2.

  4. 4.

End!

Apparently, the process ends after a finite number of steps and yields two
sets $\sigma_{1},\sigma_{2}$ with the following properties:

  1. 1.

If $\lambda\in\sigma_{1}$ (or $\in\sigma_{2}$), then
$\bar{\lambda}\in\sigma_{1}$ (or $\in\sigma_{2}$, respectively) too, where
$\bar{\lambda}$ denotes the complex conjugate of $\lambda$.

  2. 2.

$\lambda_{i}+\lambda_{j}\neq 0$ ($i,j$ may be equal) for any two eigenvalues
$\lambda_{i},\lambda_{j}\in\sigma_{1}$.

  3. 3.

$\sum_{\lambda\in\sigma_{2}}{\lambda}=0$.

  4. 4.

$\sigma=\sigma_{1}\cup\sigma_{2}$.

According to the well-known real Schur decomposition theorem (See Theorem
7.4.1 in [27]), there exists an orthogonal matrix $U\in\mathbb{R}^{n\times n}$
such that

$U^{T}FU=\begin{pmatrix}R_{11}&0&\cdots&0\\\ R_{21}&R_{22}&\cdots&0\\\
\vdots&\vdots&\ddots&\vdots\\\ R_{m1}&R_{m2}&\cdots&R_{mm}\\\
\end{pmatrix}\eqqcolon R\in\mathbb{R}^{n\times n},$

where each $R_{ii}$ is a $1\times 1$ matrix or a $2\times 2$ matrix having
complex conjugate eigenvalues. In view of properties 1)-4) of $\sigma_{1}$ and
$\sigma_{2}$, we can, without loss of generality, always assume that $R$ can
be divided into blocks as

$R=\begin{pmatrix}\mathcal{R}_{11}&0\\\ \mathcal{R}_{21}&\mathcal{R}_{22}\\\
\end{pmatrix},$

where $\mathcal{R}_{11}\in\mathbb{R}^{p\times p}$ is a matrix with
$\sigma_{1}$ as its set of eigenvalues (which implies that
$\mathrm{tr}\mathcal{R}_{11}=\mathrm{tr}F$ in view of property 3)) and
$\mathcal{R}_{22}\in\mathbb{R}^{(n-p)\times(n-p)}$ a matrix with $\sigma_{2}$
as its set of eigenvalues. Since $\mathcal{R}_{11}$ satisfies the condition of
Theorem 7, there exist, according to the discussion above, an orthogonal
matrix $\tilde{Q}$, a diagonal matrix $\tilde{\Lambda}$, and a skew-matrix
$\tilde{B}$ such that
$\tilde{\Lambda}\tilde{Q}^{T}\mathcal{R}_{11}\tilde{Q}=\tilde{\Lambda}^{2}+\tilde{B}$.
Then, a simple calculation shows that (21) is satisfied with

$\Lambda=\begin{pmatrix}\tilde{\Lambda}&0\\\ 0&0\\\ \end{pmatrix},\quad
B=\begin{pmatrix}\tilde{B}&0\\\ 0&0\\\ \end{pmatrix},\quad\text{and}\quad
Q=U\begin{pmatrix}\tilde{Q}&0\\\ 0&I\\\ \end{pmatrix},$

which completes the proof. ∎

##  5 Conclusions

In this paper, a novel approach for constructing Lyapunov functions was
proposed. The advantage of this approach is that under some additional
conditions, the resulting Lyapunov function can give an accurate
characterization of the attraction domain of the given dynamical system. The
approach make use of a specific decomposition of the defining vector field of
the system. For nonlinear systems, a sufficient condition for the existence of
such decompositions was given. For linear systems, it was proved that such
decompositions always exist and can be found explicitly by solving a specific
algebraic Riccati equation.

## References

  * [1] A. M. Lyapnov, The General Problem of the Stability of Motion, Kharkov Mathematical Society, Kharkov, Russia, 1892.
  * [2] M. Haddad, V. Chellaboina, Nonlinear Dynamical Systems and Control, Princeton University Press, William Street, Princeton, New Jersey, 2007.
  * [3] Y. Wang, C. Li, D. Cheng, Generalized Hamiltonian realization of time-invariant nonlinear systems, Automatica 39 (8) (2003) 1437–1443.
  * [4] B. Maschke, R. Ortega, A. J. Van Der Schaft, Energy-based Lyapunov functions for forced Hamiltonian systems with dissipation, IEEE Transactions on automatic control 45 (8) (2000) 1498–1502.
  * [5] P. Giesl, Construction of Global Lyapunov Functions Using Radial Basis Functions, Springer-Verlag, Berlin, Heidelberg, 2007.
  * [6] P. Giesl, S. Hafstein, Construction of Lyapunov functions for nonlinear planar systems by linear programming, Journal of Mathematical Analysis and Applications 388 (1) (2012) 463–479.
  * [7] C. Argáez, P. Giesl, S. F. Hafstein, Iterative construction of complete Lyapunov functions., in: SIMULTECH, 2018, pp. 211–222.
  * [8] K. Matsue, T. Hiwaki, N. Yamamoto, On the construction of Lyapunov functions with computer assistance, Journal of Computational and Applied Mathematics 319 (2017) 385–412.
  * [9] S. F. Hafstein, An algorithm for constructing Lyapunov functions., Electronic Journal of Differential Equations 2007.
  * [10] M. Anghel, F. Milano, A. Papachristodoulou, Algorithmic construction of Lyapunov functions for power system stability analysis, IEEE Transactions on Circuits and Systems I: Regular Papers 60 (9) (2013) 2533–2546.
  * [11] P. Giesl, N. Mohammed, Verification estimates for the construction of Lyapunov functions using meshfree collocation, Discrete & Continuous Dynamical Systems-B 24 (9) (2019) 4955.
  * [12] V. I. Zubov, Methods of AM Lyapunov and their application, P. Noordhoff, 1964\.
  * [13] F. Camilli, L. Grüne, F. Wirth, A generalization of Zubov’s method to perturbed systems, SIAM Journal on Control and Optimization 40 (2) (2001) 496–515.
  * [14] F. Camilli, L. Grüne, F. Wirth, Calculating the domain of attraction: Zubov’s method and extensions, in: Proc. of the Stability and Control Processes Conference, 2005, pp. 27–36.
  * [15] G. Valmorbida, J. Anderson, Region of attraction estimation using invariant sets and rational Lyapunov functions, Automatica 75 (2017) 37–45.
  * [16] M. Korda, D. Henrion, C. N. Jones, Controller design and region of attraction estimation for nonlinear dynamical systems, IFAC Proceedings Volumes 47 (3) (2014) 2310–2316.
  * [17] D. Henrion, M. Korda, Convex computation of the region of attraction of polynomial control systems, IEEE Transactions on Automatic Control 59 (2) (2013) 297–312.
  * [18] G. Chesi, Domain of attraction: analysis and control via SOS programming, Vol. 415, Springer Science & Business Media, 2011.
  * [19] H. K. Khalil, Nonlinear Systems, Prentice Hall, Upper Saddle River, 2002.
  * [20] J. M. Lee, Introduction to smooth manifolds, Springer Science & Business Media, 2012.
  * [21] Y. N. Grigoriev, V. F. Kovalev, S. V. Meleshko, N. H. Ibragimov, Symmetries of integro-differential equations: with applications in mechanics and plasma physics, Springer, 2010.
  * [22] P. Olver, Applications of Lie Groups to Differential Equations, Springer, 1986\.
  * [23] J. Jost, Riemannian Geometry and Geometric Analysis, Springer-Verlag, Heidelberg, Dordrecht, London, New York, 2011.
  * [24] I. Gohberg, P. Lancaster, L. Rodman, Indefinite Linear Algebra and Applications, Birkhäuser Verlag, Basel, Boston, Berlin, 2005.
  * [25] D. G. Luenberger, Invertible solutions to the operator equation $TA-BT=C$, Proceedings of the American Mathematical Society 16 (6) (1965) 1226–1229.
  * [26] V. E. Shestopal, Solution of the matrix equation $AX-XB=C$, Mathematical Notes of the Academy of Sciences of the Ussr 19 (3) (1976) 275–276.
  * [27] G. H. Golub, C. F. Van Loan, Matrix computations. 4th ed., The Johns Hopkins University Press, 2013.