Construction of Lyapunov functions using Helmholtz-Hodge decomposition

TL;DR

This paper introduces a novel method for constructing Lyapunov functions using Helmholtz-Hodge decomposition, demonstrating its effectiveness and limitations in analyzing stability of vector fields.

Contribution

It presents a new approach to construct Lyapunov functions via HHD and explores the properties of orthogonal HHD vector fields, extending gradient field concepts.

Findings

Potential functions from HHD can serve as Lyapunov functions under stability conditions.

Orthogonal HHD vector fields generalize gradient vector fields with similar stability properties.

The method's limitations are examined through planar vector field analysis.

Abstract

The Helmholtz-Hodge decomposition (HHD) is applied to the construction of Lyapunov functions. It is shown that if a stability condition is satisfied, such a decomposition can be chosen so that its potential function is a Lyapunov function. In connection with the Lyapunov function, vector fields with strictly orthogonal HHD are analyzed. It is shown that they are a generalization of gradient vector fields and have similar properties. Finally, to examine the limitations of the proposed method, planar vector fields are analyzed.

Equations

$$\dot{\mathbf{x}}=\mathbf{F}(\mathbf{x}),$$

$$\dot{f}(\mathbf{x}(t)):=\frac{d}{dt}f\left(\mathbf{x}(t)\right)=\sum _i\frac{\partial f}{\partial x_i}\frac{d{x}_i(t)}{dt},$$

$$\dot{f}(\mathbf{x})=\nabla f(\mathbf{x})\cdot \mathbf{F}(\mathbf{x}).$$

$$\mathbf{F}=-\nabla V+\mathbf{u},$$

$$\omega (\mathbf{x})\subset \{\mathbf{y}\in \bar{D}|\nabla V(\mathbf{y})=0\}.$$

$$\alpha (\mathbf{x})\subset \{\mathbf{y}\in \bar{D}|\nabla V(\mathbf{y})=0\}.$$

$$\mathbf{F}=-\nabla \tilde{V}+\tilde{\mathbf{u}}.$$

$$\nabla \cdot \mathbf{F}$$

$$\mathrm{Hess}V$$

$$\frac{\partial \dot{V}}{\partial x_i}$$

$$\frac{{\partial }^2\dot{V}}{\partial x_i\partial x_j}=-2\sum _{l=1}^n\frac{{\partial }^{2}V}{\partial x_l\partial x_i}\frac{{\partial }^{2}V}{\partial x_l\partial x_j}+\sum _{l=1}^n\frac{{\partial }^{2}V}{\partial x_l\partial x_j}\frac{\partial u_l}{\partial x_i}+\sum _{l=1}^n\frac{{\partial }^{2}V}{\partial x_l\partial x_i}\frac{\partial u_l}{\partial x_j}.$$

$$\mathrm{Hess}\dot{V}=-2\left(\mathrm{Hess}V\right)\left(\mathrm{Hess}V\right)+{}^t\left(D{\mathbf{u}}_0\right)\left(\mathrm{Hess}V\right)+\left(\mathrm{Hess}V\right)D{\mathbf{u}}_0.$$

$$⟨\left({}^{t}D{\mathbf{F}}_0\left(\mathrm{Hess}V\right)+\left(\mathrm{Hess}V\right)D{\mathbf{F}}_0\right)\mathbf{x},\mathbf{x}⟩=\frac{2\rho }{n}⟨D{\mathbf{F}}_0\mathbf{x},\mathbf{x}⟩<0$$

$$⟨(\mathrm{Hess}\dot{V})\mathbf{x},\mathbf{x}⟩$$

$$=-2||(\mathrm{Hess}V)\mathbf{x}|{|}^2+2⟨D{\mathbf{u}}_0\mathbf{x},\left(\mathrm{Hess}V\right)\mathbf{x}⟩,$$

$$⟨(\mathrm{Hess}\dot{V})\mathbf{x},\mathbf{x}⟩$$

$$\le \left({\lambda }_{\mathbf{u}}^2-{\mu }_{\mathbf{F}}^2-{\mu }_V^2\right)||\mathbf{x}|{|}^2.$$

$${\bf F}(x,y)=\left(\begin{array}[]{cc}-a&b\\ c&-b\end{array}\right)\left(\begin{array}[]{c}x\\ y\end{array}\right).$$

$$V$$

$$\mathbf{u}$$

$${\bf F}=\left(\begin{array}[]{cc}-a&b\\ c&-b\end{array}\right){\bf x}-\nabla p+\left(\begin{array}[]{c}q\\ r\end{array}\right),$$

$$\frac{d}{dt}\left(\begin{array}[]{c}x\\ y\end{array}\right)=\left(\begin{array}[]{c}x\left(-1+4x^{2}+\frac{1}{4}y^{2}\right)+\frac{1}{8}y^{3}\\ y\left(-1+\frac{5}{2}x^{2}+\frac{3}{8}y^{2}\right)-6x^{3}\end{array}\right).$$

$$\Delta V=2-\frac{29}{2}{x}^2-\frac{11}{8}{y}^2,$$

$$V_1(x,y)=V_0(x,y)-\frac{29}{24}{x}^4-\frac{11}{96}{y}^4.$$

$$\frac{d}{dt}\left(\begin{array}[]{c}x\\ y\end{array}\right)=\left(\begin{array}[]{c}\mu x-\omega y-\left(x^{2}+y^{2}\right)x\\ \omega x+\mu y-\left(x^{2}+y^{2}\right)y\end{array}\right).$$

$$V=-\frac{\mu }{2}(x^2+y^2)+\frac{1}{4}(x^2+y^2{)}^2,$$

$${\bf u}=\left(\begin{array}[]{c}-\omega y\\ \omega x\end{array}\right).$$

$$\omega (\mathbf{x})\subset \{\mathbf{y}\in \bar{D}|\nabla V(\mathbf{y})=0\}.$$

$$\alpha (\mathbf{x})\subset \{\mathbf{y}\in \bar{D}|\nabla V(\mathbf{y})=0\}.$$

$$V^{\ast }\left(\mathbf{x},T\right):=\sup \{V(\mathbf{y})\mid \mathbf{y}\in \overline{{\cup }_{T\le t}\left\{{\varphi }^t(\mathbf{x})\right\}}\}.$$

$$V^{\ast }\left(\mathbf{x},T\right)=V\left({\varphi }^T(\mathbf{x})\right).$$

$$V(\mathbf{y})=\underset{T\to \infty }{\lim }{V}^{\ast }\left(\mathbf{x},T\right).$$

Full text

The first author is supported by Grant-in-Aid for JSPS Fellows (17J03931).

Tomoharu Suda

Graduate School of Human and Environmental Studies, Kyoto University

Yoshida-nihonmatsu-cho, Sakyo-ku, Kyoto, 606-8501, Japan

(Communicated by Bernold Fiedler)

1. Introduction

The Helmholtz–Hodge decomposition (HHD) is a decomposition of vector fields whereby they are expressed as the sum of a gradient vector field and a divergence-free vector field. It is applied to the study of Navier–Stokes equations (projection method) and to the detection and visualization of singularities of vector fields. For a review of the literature on HHD, the reader may refer to [1]. There are several methods for obtaining HHD. The most naive is by solving the Poisson equation, as its solution yields such a decomposition [8]. Thus, an HHD exists, provided that the Poisson equation has a solution. However, a certain boundary condition, whose choice is not obvious, must be prescribed. This condition is often imposed so that the two vector fields in the decomposition are orthogonal in the $L^{2}$ sense, and it is discussed in [5]. Other techniques to obtain an HHD involve the construction of basis functions or Green’s function [2, 7] and thus avoid the complications related to the choice of boundary conditions.

Although there are examples of HHD applied to the detection of singularities of vector fields [13, 16], few studies have been concerned with HHD from the viewpoint of dynamical systems. In the planar case, Demongeot, Glade, and Forest proposed a scheme for the analysis of planar vector fields with a limit cycle, in which they constructed a polynomial vector field so that it has a limit cycle resembling that of the original vector field [3, 4]. Mendes and Duarte approximated $n$-dimensional vector fields by divergence-free vector fields [6]. Even though these methods are different, their aim is to determine a function that best captures the behavior of the given vector field using HHD.

To describe vector fields using a function, a natural option is to consider a Lyapunov function. A Lyapunov function of a vector field is a smooth function that takes minimum value at an equilibrium point (more generally, on an invariant set), and decreases along solution curves in a neighborhood of the equilibrium point. The potential function of gradient vector fields is an example of Lyapunov function. Even though its existence is known for asymptotically stable equilibrium points of vector fields, the actual construction is not obvious except in the case of relatively simple vector fields. Several methods have been proposed for numerically constructing Lyapunov functions, and a comprehensive review can be found in [9]. They vary depending on the nature of the vector fields. In particular, methods based on linear programming are often applied to linear systems; however they may not be used for general vector fields. For an example of this construction, the reader is referred to [12]. If the vector fields are nonlinear, other techniques may be applied. Lyapunov functions can be constructed by considering a partial differential equation whose solution yields a Lyapunov function and then approximating this solution. A detailed discussion on such a method can be found in [10]. This approach can handle a variety of vector fields; however, it cannot readily be associated with the resulting function in terms of dynamics.

From the viewpoint of dynamical systems, it appears natural to choose an HHD whose gradient vector field is generated by a Lyapunov function. Indeed, the method of Demongeot, Glade, and Forest can be interpreted as an attempt to choose an HHD whose potential function is a Lyapunov function of a limit cycle [4]. Thus, it appears possible that a “good” choice of HHD yields a Lyapunov function. Furthermore, a Lyapunov function thereby obtained is naturally interpreted as the dissipative part of the original vector field. Hence, it is likely to offer insight into its dynamics.

In this study, the construction of Lyapunov functions using HHD is considered. It is shown that, under certain stability conditions, this is possible, and a method is provided for formulating it as an optimization problem. Moreover, strictly orthogonal HHDs are considered. This notion of orthogonality differs from the $L^{2}$-orthogonality commonly imposed in the literature meaning that the two vector fields obtained by the decomposition are orthogonal everywhere. It is proved that vector fields with strictly orthogonal HHD are similar to gradient vector fields and their behavior is completely determined in terms of potential functions that are Lyapunov functions of invariant sets. Finally, the limitations of the proposed method are examined in the case of planar vector fields.

Main Results. The terminology used in this study is first provided. It is assumed that all vector fields are $C^{2}$, unless indicated otherwise. Vectors will be denoted by bold fonts. Autonomous differential equations are considered and will be given in the form

[TABLE]

where ${\bf F}({\bf x})$ is a smooth vector field defined on $\bar{\Omega}$, where $\Omega\subset{\mathbb{R}}^{n}$ is an open domain which contains the origin. In what follows, unless otherwise stated, it is assumed that the domain $\Omega$ is bounded and has $C^{2}$-boundary.

Definition 1.1** (Orbital derivative).**

For a $C^{1}$-function $f$ on a domain $\Omega\subset{\mathbb{R}}^{n}$, its orbital derivative $\dot{f}$ along solutions of the equation (1) is a function defined by

[TABLE]

where ${\bf x}(t)$ is a solution of the equation (1). Since the vector field is autonomous, the orbital derivative $\dot{f}$ is given explicitly as a function of ${\bf x}$, namely,

[TABLE]

Lyapunov functions are defined as follows, requiring that they take a minimum value at an equilibrium point.

Definition 1.2** (Lyapunov functions).**

Let ${\bf x}_{*}$ be an equilibrium point and $U\subset{\mathbb{R}}^{n}$ a neighborhood of ${\bf x}_{*}$. A $C^{1}$-function $L:U\rightarrow{\mathbb{R}}$ is said to be a (local) Lyapunov function if the following hold:

1. (1)

For all ${\bf x}\in U-\{{\bf x}_{*}\}$, we have $L({\bf x})>L({\bf x}_{*}).$ 2. (2)

On the set $U-\{{\bf x}_{*}\}$, we have ${\dot{L}}\leq 0.$

In this study, the Helmholtz–Hodge decomposition is defined as follows.

Definition 1.3** (Helmholtz–Hodge decomposition).**

For a vector field ${\bf F}$ on a domain $\Omega\subset{\mathbb{R}}^{n}$, its Helmholtz–Hodge decomposition (HHD) is a decomposition of the form

[TABLE]

where $V:\Omega\rightarrow\mathbb{R}$ is a $C^{3}$ function and ${\bf u}$ is a vector field on $\Omega$ with the property $\nabla\cdot{\bf u}=0$. $V$ is called a potential function.

From the definition, it can be seen that the addition of a harmonic function $h$ to the potential function $V$ yields another decomposition, ${\bf F}=-\nabla\left(V+h\right)+\left({\bf u}+\nabla h\right)$. Therefore HHD is not unique in general. The existence of HHD is established if there is a solution of the Poisson equation obtained by taking the divergence of the equation in the definition. Such a solution exists if the boundary of the domain under consideration is, for example, sufficiently smooth. More details may be found in [8]. By the assumption on the smoothness of $\partial\Omega$, at least one HHD exists for the vector field considered on $\Omega$.

In the analysis of potential functions, the following concept is used.

Definition 1.4** (Singular values of a matrix).**

Let $A$ be a square real matrix of order $n$. Then, the matrix $\,{}^{t}AA$ is nonnegative definite; therefore, it has positive eigenvalues. The square root of an eigenvalue of $\,{}^{t}AA$ is called a singular value of $A$.

The singular values of a matrix are denoted by Greek letters, for example, $\lambda$ or $\mu$. It is easily verified that a matrix and its transpose have the same set of singular values. More details may be found in [15].

The main results of this study will now be stated. It is first shown that under a certain stability condition, it is possible to obtain an HHD with its potential function assuming a minimum at an equilibrium point.

Main Theorem 1**.**

Let ${\bf F}:\bar{\Omega}\to{\mathbb{R}}^{n}$ be a vector field with an equilibrium point at the origin. If $\nabla\cdot{\bf F}({\bf 0})<0$, then it has an HHD with the potential function $V$ having a minimum at the origin.

To obtain a Lyapunov function, a function whose orbital derivative decreases monotonically should be determined. This may be achieved by controlling the Jacobian matrix of the divergence-free vector field ${\bf u}$ and the Hessian of the potential function $V$. The next theorem provides a sufficient condition whereby the potential function is a Lyapunov function in a form suitable for optimization. The Jacobian matrix of ${\bf u}$ and ${\bf F}$ at the origin is denoted by $D{\bf u}_{\bf 0}$ and $D{\bf F}_{\bf 0}$, respectively.

Main Theorem 2**.**

Let ${\bf F}:\bar{\Omega}\to{\mathbb{R}}^{n}$ be a vector field with an equilibrium point at the origin and ${\bf F}=-\nabla V+{\bf u}$ be an HHD, assume $V$ has a minimum at the origin, and let the largest singular value of $D{\bf u}_{\bf 0}$ be $\lambda_{\bf u}$, the smallest singular value of $D{\bf F}_{\bf 0}$ be $\mu_{\bf F}$, and the smallest eigenvalue of ${\rm Hess}\,V$ be $\mu_{V}$. If $\lambda_{\bf u}^{2}-\mu_{V}^{2}<\mu_{\bf F}^{2}$, then $V$ is a (local) Lyapunov function of the origin.

A simple example of the use of Main Theorem 2 will be given in Example 1.

In connection with Lyapunov functions, vector fields with strictly orthogonal HHD are an example where the analysis based on potential function succeeds.

Definition 1.5**.**

HHD ${\bf F}=-\nabla V+{\bf u}$ is said to be strictly orthogonal on $D\subset{\mathbb{R}}^{n}$ if ${\bf u}({\bf x})\cdot\nabla V({\bf x})=0$ for all ${\bf x}\in D$.

Although vector fields with strictly orthogonal HHD exhibit behavior different from that of gradient vector fields, they are in fact a generalization of them, as the next theorem states.

Main Theorem 3**.**

Let $D\subset{\mathbb{R}}^{n}$ be a bounded domain. If a vector field ${\bf F}:{\mathbb{R}}^{n}\to{\mathbb{R}}^{n}$ has a strictly orthogonal HHD on $\bar{D}$ with potential function $V$, then the following hold:

1. (1)

If $D$ is forward invariant, then for all ${\bf x}\in D$,

[TABLE] 2. (2)

If $D$ is backward invariant, then for all ${\bf x}\in D$,

[TABLE]

This paper is organized as follows. In Section 2, it is shown that it is possible to construct Lyapunov functions by properly choosing an HHD. Moreover, Main Theorems 1 and 2 are proved. In Section 3, the properties of strictly orthogonal HHDs are considered. Main Theorem 3 is proved, and its compatibility with a boundary condition widely used in literature is further studied. In Section 4, to discuss the limitations of the analysis based on HHD, planar vector fields are analyzed using Fourier series expansions.

2. Construction of Lyapunov functions using HHD

In this section, the construction of Lyapunov functions near an equilibrium point is considered. The emphasis here is not on solving actual problems, but on showing that it is possible to construct Lyapunov functions using HHD. Without loss of generality, it may be assumed that the equilibrium point is the origin.

The aim here is to choose the potential functions of an HHD so that they can be used as Lyapunov functions. To achieve this, it is first ensured that the potential function assumes a minimum value at the equilibrium point.

Theorem 2.1** (Main Theorem 1).**

Let ${\bf F}:\bar{\Omega}\to{\mathbb{R}}^{n}$ be a vector field with an equilibrium point at the origin. If $\nabla\cdot{\bf F}({\bf 0})<0$, then it has an HHD such that the potential function $V$ has a minimum at the origin.

Proof.

An HHD of ${\bf F}$ is first chosen, that is,

[TABLE]

This is possible because an HHD exists in $\Omega$. If $\nabla\tilde{V}({\bf 0})\neq{\bf 0}$, then $\tilde{V}({\bf x})$ is replaced by $\tilde{V}({\bf x})-\nabla\tilde{V}({\bf 0})\cdot{\bf x}$, which also yields an HHD. It should be noted that $\nabla\tilde{V}({\bf 0})$ is merely a constant vector. Therefore, it may always be assumed that $\nabla\tilde{V}({\bf 0})={\bf 0}$.

The Hessian of $V$ is denoted by ${\rm Hess}\,V$. As

[TABLE]

then ${\rm tr}\,{\rm Hess}\,\tilde{V}$ is positive at the origin. Therefore, the matrix $\frac{\rho}{n}I$ is positive definite, where $\rho$ is the value of ${\rm tr}\,{\rm Hess}\,\tilde{V}$ at the origin. Let the symmetric matrix $A$ be defined by $-{\rm Hess}\,\tilde{V}+\frac{\rho}{n}I$. Then, ${\rm tr}\,A=0$, and $\,{}^{t}{\bf x}A{\bf x}\,$ is a harmonic function.

Let $V$ be defined by $\tilde{V}+\,^{t}{\bf x}A{\bf x}$. This function also yields an HHD ${\bf F}=-\nabla V+{\bf u}$, because $\,{}^{t}{\bf x}A{\bf x}$ is a harmonic function. As $\nabla V({\bf 0})={\bf 0}$ and the Hessian of $V$ at the origin is

[TABLE]

we conclude that $V$ takes a minimum value at the origin, owing to diagonality and positive definiteness. ∎

Remark 1**.**

The condition in the hypothesis of Main Theorem 1 is related to the Jacobian matrix of the vector field at the origin via the identity $\nabla\cdot{\bf F}({\bf 0})={\rm tr}D{\bf F}_{\bf 0}$. If this hypothesis fails, there are two possibilities: $\nabla\cdot{\bf F}({\bf 0})>0$ or $\nabla\cdot{\bf F}({\bf 0})=0$. In the former case, the origin cannot be a stable equilibrium point, and there is no Lyapunov function for it. In the latter case, the origin is a saddle, a center, or a degenerate equilibrium.

Even though the condition that $V$ takes a minimum value at the origin does not imply that it is a Lyapunov function, we have the following result.

Proposition 1**.**

Let ${\bf F}:\bar{\Omega}\to{\mathbb{R}}^{n}$ be a vector field with an equilibrium point at the origin, and let ${\bf F}=-\nabla V+{\bf u}$ be an HHD. If $V$ assumes a minimum at the origin, then $\dot{V}({\bf 0})=0$ and the origin is a critical point of $\dot{V}$. The Hessian of $\dot{V}$ is given by $\,{}^{t}\left(D{\bf F}_{\bf 0}\right)\left({\rm Hess}\,V\right)+\left({\rm Hess}\,V\right)D{\bf F}_{\bf 0}$.

Proof.

By the definition of the orbital derivative, we have $\dot{V}=\nabla V\cdot{\bf F}=-\left|\nabla V\right|^{2}+\nabla V\cdot{\bf u}$. As $V$ assumes a minimum at the origin, $\dot{V}({\bf 0})$ is zero.

The first-order derivatives of $\dot{V}$ are

[TABLE]

The first and last terms vanish at the origin because $V$ assumes a minimum. As ${\bf u}({\bf 0})={\bf F}({\bf 0})+\nabla V({\bf 0})=0$, the second term also vanishes. Therefore, the origin is a critical point of $\dot{V}$.

The second-order derivatives of $\dot{V}$ at the origin are

[TABLE]

Therefore,

[TABLE]

As $D{\bf F}_{\bf 0}=-\left({\rm Hess}\,V\right)+D{\bf u}_{\bf 0}$, this is equivalent to $\,{}^{t}\left(D{\bf F}_{\bf 0}\right)\left({\rm Hess}\,V\right)+\left({\rm Hess}\,V\right)D{\bf F}_{\bf 0}$. ∎

To summarize, an HHD of a vector field ${\bf F}$ with an equilibrium point at the origin yields a Lyapunov function if its potential function $V$ assumes a minimum at the origin and the matrix $\,{}^{t}D{\bf F}_{\bf 0}\left({\rm Hess}\,V\right)+\left({\rm Hess}\,V\right)D{\bf F}_{\bf 0}$ is negative definite. If these hypotheses are satisfied, $\dot{V}$ takes a maximum value of [math] at the origin; therefore, $\dot{V}\leq 0$ in a neighborhood of the origin. For example, we have the following result.

Corollary 1**.**

If $\left<D{\bf F}_{0}{\bf x},{\bf x}\right>\,<0$ for all ${\bf x}\in{\mathbb{R}}^{n}\backslash\{{\bf 0}\}$ in addition to the conditions in Main Theorem 1, the potential function $V$ constructed in Main Theorem 1 is a Lyapunov function.

Proof.

As the potential function $V$ constructed in Main Theorem 1 assumes a minimum at the origin, Proposition 1 can be applied. We need only show that $\,{}^{t}D{\bf F}_{\bf 0}\left({\rm Hess}\,V\right)+\left({\rm Hess}\,V\right)D{\bf F}_{\bf 0}$ is negative definite.

As the Hessian of the potential function $V$ constructed in Main Theorem 1 is $\frac{\rho}{n}I,$ we have

[TABLE]

if ${\bf x}\neq 0$. ∎

Remark 2**.**

The condition that $\left<D{\bf F}_{0}{\bf x},{\bf x}\right>\,<0$ for all ${\bf x}\in{\mathbb{R}}^{n}\backslash\{{\bf 0}\}$ is equivalent to the symmetric matrix $D{\bf F}_{0}+\,^{t}D{\bf F}_{0}$ being negative definite. Therefore, the stability of the origin is easily obtained without resort to Lyapunov functions.

Even though a potential function that assumes a minimum at the origin can be chosen under a relatively mild condition, it is not obvious that $\,{}^{t}D{\bf F}_{\bf 0}\left({\rm Hess}\,V\right)+\left({\rm Hess}\,V\right)D{\bf F}_{\bf 0}$ can be chosen so as to be negative definite. Thus, this problem is reformulated as an optimization problem, as shown by the next theorem.

Theorem 2.2** (Main Theorem 2).**

Let ${\bf F}:\bar{\Omega}\to{\mathbb{R}}^{n}$ be a vector field with an equilibrium point at the origin and ${\bf F}=-\nabla V+{\bf u}$ be an HHD, assume $V$ has a minimum at the origin, and let the largest singular value of $D{\bf u}_{\bf 0}$ be $\lambda_{\bf u}$, the smallest singular value of $D{\bf F}_{\bf 0}$ be $\mu_{\bf F}$, and the smallest eigenvalue of ${\rm Hess}\,V$ be $\mu_{V}$. If $\lambda_{\bf u}^{2}-\mu_{V}^{2}<\mu_{\bf F}^{2}$, then $V$ is a (local) Lyapunov function of the origin.

Proof.

The proof is by direct calculation. Using the representation obtained in Proposition 1, we have

[TABLE]

Therefore, $({\rm Hess}\,\dot{V})$ is negative definite. Combining this with Proposition 1, we conclude that ${\dot{V}}$ assumes a maximum value [math] at the origin and $V$ is a Lyapunov function. ∎

Using Main Theorem 2, a Lyapunov function may be constructed by seeking a decomposition that satisfies $\lambda_{\bf u}^{2}-\mu_{V}^{2}<\mu_{\bf F}^{2}$, even though a solution does not necessarily exist. In particular, this method is not applicable to vector fields with degenerate Jacobian matrix $D{\bf F}_{\bf 0}$.

For optimization, a harmonic function $h$ may be constructed so that the condition $\lambda_{{\bf u}+\nabla h}^{2}-\mu_{V+h}^{2}<\mu_{\bf F}^{2}$ is satisfied. This could be carried out by, for example, considering harmonic polynomials or harmonic functions constructed using the Poisson kernel representation. As only the quadratic terms of harmonic functions are relevant to the calculation of $\lambda_{{\bf u}+\nabla h}$ and $\mu_{V+h}$, it will suffice to consider the addition of quadratic harmonic polynomials. However, the implementation of such a procedure requires further investigation.

Example 1**.**

A simple example of the use of Main Theorem 2 is now provided. Let the vector field ${\bf F}$ on ${\mathbb{R}}^{2}$ be given by

[TABLE]

with $a>b>c>0$. Then, there is an HHD of ${\bf F}=-\nabla V+{\bf u}$ given by

[TABLE]

This decomposition is obtained by selecting the “diagonal terms” of the vector field as $-\nabla V$ and the remaining terms as ${\bf u}$. It is now shown that $V$ is a Lyapunov function of the origin. This is easily verified by direct calculation as well. As $\lambda_{\bf u}=\mu_{V}=b$, it suffices to show that $\mu_{\bf F}^{2}>0$. This is equivalent to the matrix $\,{}^{t}\left(D{\bf F}_{\bf 0}\right)\left(D{\bf F}_{\bf 0}\right)$ being regular, which is obviously true because the coefficient matrix is regular. By Main Theorem 2, we conclude that $V$ is a Lyapunov function of the origin.

This analysis can be extended to vector fields involving higher-order terms with the same linear part if they have a decomposition with the lower-order terms of the potential function and the divergence-free vector field in the form given here. For example, let a vector field ${\bf F}$ be given by

[TABLE]

where $p,q,r$ are homogeneous polynomials with ${\rm deg}(p)>2$, ${\rm deg}(q),{\rm deg}(r)>1$, and $\frac{\partial q}{\partial x}+\frac{\partial r}{\partial y}=0$. Then, $V=\frac{a}{2}x^{2}+\frac{b}{2}y^{2}+p(x,y)$ is a Lyapunov function of the origin.

Example 2**.**

Another example of the use of Main Theorem 2 is now presented to illustrate the features of the proposed method, compared to other schemes. The following ordinary differential equation is considered:

[TABLE]

This equation was analyzed as an example of numerical methods for constructing Lyapunov functions in [10]. Samples of solution curves are given in the right panel of Figure 1. As the origin is an asymptotically stable equilibrium point of the vector field, a local Lyapunov function of the origin may be constructed. For example, linear approximation of the vector field yields a local Lyapunov function $V_{0}=\frac{1}{2}\left(x^{2}+y^{2}\right)$, which is strictly decreasing on the disk $\{(x,y)\in{\mathbb{R}}^{2}\mid V_{0}(x,y)\leq 0.09\}$. For the details of the construction via linear approximation, the reader is referred to [10].

Local Lyapunov functions are now constructed by applying Main Theorem 2. A potential function $V$ of an HHD for the vector field in the equation (3) satisfies the following Poisson equation:

[TABLE]

which is obtained by taking the divergence in the definition of HHD. As $-\nabla V_{0}$ gives the linear part of the vector field, it is natural to seek a potential function $V$ with the terms of order less than two given by $V_{0}$, so that the linear approximation of $-\nabla V$ coincides with $-\nabla V_{0}$. If such $V$ is chosen, it takes minimum value at the origin. Furthermore, we have $D{\bf u}_{\bf 0}=0$; therefore, $\lambda_{\bf u}=0$ and the hypotheses of Main Theorem 2 hold trivially.

As no cross term of $x$ and $y$ appears in Equation (4), an obvious solution is given by integrating each non-constant term twice, namely,

[TABLE]

Contours of $V_{1}$ and the sign of $\dot{V_{1}}$ are shown in the left panel of Figure 1. Compared with $V_{0}$, it succeeds in expanding the domain where solutions are determined to be divergent, whereas it fails to improve the estimate for the local attracting basin at the origin.

By retaking HHD, it may be possible to improve it. Another Lyapunov function is now constructed by adding harmonic polynomials to $V_{1}$. As the original vector field is roughly symmetric about the $y$ axis, the addition of, for example, the quartic harmonic functions of the form $a\left(x^{4}-6x^{2}y^{2}+y^{4}\right)$ is considered, where $a$ is a constant. If $a=0.5$, we obtain a Lyapunov function $V_{2}(x,y)=V_{1}(x,y)+0.5\left(x^{4}-6x^{2}y^{2}+y^{4}\right)$. Contours of $V_{1}$ and the sign of $\dot{V_{1}}$ are shown in Figure 2. It should be noted that the estimate for the local attracting basin at the origin is improved.

Only linear approximation is essentially required in the construction above. However, the resulting Lyapunov functions respect the nonlinear nature of the original equation, despite the ease of calculation. This is due to the fact that HHD reflects the nonlinearity of the original vector field. Furthermore, as numerous different Lyapunov functions may be obtained at a relatively low computational cost, it appears possible that a better picture of the local attracting basin could be obtained by comparing them. This points out another method for the computational analysis of the attracting basin.

3. Strictly orthogonal HHD

As was mentioned in Introduction, orthogonality in the $L^{2}$ sense is frequently imposed on HHD. Indeed, such a decomposition is obtained by imposing as a boundary condition that the divergence-free vector field should be orthogonal to the normal vector of the boundary, and HHD with this boundary condition is unique [5]. However, this is not sufficient if the aim is to completely analyze the behavior of vector fields. In particular, equilibrium points are not necessarily respected under this condition, and potential functions may have critical points irrelevant to the flow of the original vector field. These properties result in artifacts, as seen, for example, in [16].

In certain cases, a strictly orthogonal HHD may be obtained, which is sufficient for completely describing the behavior of the vector field.

Definition 3.1**.**

The HHD ${\bf F}=-\nabla V+{\bf u}$ is said to be strictly orthogonal on $D\subset{\mathbb{R}}^{n}$ if ${\bf u}({\bf x})\cdot\nabla V({\bf x})=0$ for all ${\bf x}\in D$.

If an HHD is strictly orthogonal, then it is orthogonal in the $L^{2}$ sense. For example, some vector fields given by relatively simple polynomials obviously have such a decomposition.

Example 3** (Strictly orthogonal HHD).**

For $\mu>0$ and $\omega<0$, the following differential equation is considered:

[TABLE]

The vector field on the right-hand side has a strictly orthogonal HHD, namely,

[TABLE]

[TABLE]

Solution curves of the vector field are illustrated in Figure 3. Solution curves of $-\nabla V$ and ${\bf u}$ given above are shown in Figure 4.

For vector fields with strictly orthogonal HHD, there is a complete analogue to gradient systems, which are the extreme case of vector fields with Lyapunov functions. It is well known that the behavior of gradient vector fields is given in terms of potential functions [14]. The next theorem states that the same holds true if there is a strictly orthogonal HHD; therefore, vector fields with such a decomposition are a generalization of gradient vector fields.

Theorem 3.2** (Main Theorem 3).**

Let $D\subset{\mathbb{R}}^{n}$ be a bounded domain. If a vector field ${\bf F}:{\mathbb{R}}^{n}\to{\mathbb{R}}^{n}$ has a strictly orthogonal HHD on $\bar{D}$ with potential function $V$, then the following hold:

1. (1)

If $D$ is forward invariant, then for all ${\bf x}\in D$,

[TABLE] 2. (2)

If $D$ is backward invariant, then for all ${\bf x}\in D$,

[TABLE]

The proof of this theorem is similar to that of LaSalle’s invariance principle [11]. The following lemma is used in the proof of Main Theorem 3.

Lemma 3.3**.**

Let $D\subset{\mathbb{R}}^{n}$ be a bounded domain. It is assumed that a vector field ${\bf F}:{\mathbb{R}}^{n}\to{\mathbb{R}}^{n}$ has a strictly orthogonal HHD on $\bar{D}$ with potential function $V$ and $D$ is forward invariant. For each point ${\bf x}\in D$, let $V^{*}\left({\bf x},\cdot\right):[0,\infty)\rightarrow{\mathbb{R}}$ be defined by

[TABLE]

Then, the following hold:

1. (1)

$V^{*}\left({\bf x},\cdot\right)$* is well-defined and given explicitly by*

[TABLE]

In particular, $V^{*}\left({\bf x},\cdot\right)$ is monotonically decreasing. 2. (2)

For all ${\bf y}\in\omega({\bf x})$, we have

[TABLE]

Proof.

(1) As $\dot{V}\leq 0$ for all ${\bf y}\in\bigcup_{T\leq t}\left\{\phi^{t}({\bf x})\right\}$, we have

[TABLE]

For all $\tilde{\bf y}\in\overline{\bigcup_{T\leq t}\left\{\phi^{t}({\bf x})\right\}}$, there is a sequence ${\bf y}_{n}$ such that

[TABLE]

The continuity of $V$ implies that

[TABLE]

Therefore,

[TABLE]

Moreover,

[TABLE]

thus, $V^{*}\left({\bf x},T\right)\geq V(\phi^{T}({\bf x}))$. Therefore, $V^{*}\left({\bf x},T\right)=V\left(\phi^{T}({\bf x})\right)$. Monotonicity is obvious.

(2) The compactness of $\bar{D}$ implies that $V$ is bounded on $D$. Therefore, the monotonic function $V^{*}\left({\bf x},\cdot\right)$ is also bounded. Thus, the quantity

[TABLE]

is well-defined. If ${\bf y}\in\omega({\bf x})$, there exists a sequence $t_{n}$ such that

[TABLE]

Using the continuity of $V$ and (1), we conclude that

[TABLE]

Therefore, $V$ is constant on $\omega({\bf x})$ and equal to $V^{*}\left({\bf x},\infty\right)$. ∎

Main Theorem 3 is proved below.

Proof.

The results for $\alpha$-limit sets follow from those for $\omega$-limit sets. Thus, only the latter case is proved here.

Let ${\bf y}\in\omega({\bf x})$ be fixed. As $\omega({\bf x})$ is invariant for all $t\in{\mathbb{R}}$, we have

[TABLE]

By Lemma 3.3, we have for all $t\in{\mathbb{R}}$,

[TABLE]

As

[TABLE]

we conclude that $\nabla V({\bf y})={\bf 0}$. Therefore, $\omega({\bf x})\subset\{{\bf y}\in{\bar{D}}\mid\nabla V({\bf y})={\bf 0}\}$. ∎

Remark 3**.**

The method of Demongeot, Glade, and Forest is essentially an approximation by polynomial vector fields with a strictly orthogonal HHD[4].

The main limitation of Main Theorem 3 is the difficulty in finding a strictly orthogonal HHD. In general, its existence cannot be ensured. Furthermore, this theorem implies that the limit sets for a vector field with a strictly orthogonal HHD are restricted in smooth manifolds, which is not plausible for higher dimensional vector fields.

Nevertheless, if a strictly orthogonal HHD can be obtained, the behavior of the vector field can be rigorously determined, as the next example illustrates.

Example 4**.**

As an application of Main Theorem 3, it is shown that the vector field in Example 3 has a limit cycle. As is easily verified, the vector field has a single equilibrium point at the origin. The set $S=\{{\bf y}\in{\mathbb{R}}^{2}\mid\nabla V({\bf y})={\bf 0}\}$ is the union of the origin and the circle $C=\{{\bf y}\in{\mathbb{R}}^{2}\mid{\bf y}\cdot{\bf y}=\mu\}$. By Main Theorem 3, $\omega({\bf x})\subset S$ for all ${\bf x}\in{\mathbb{R}}^{2}$. As $\omega({\bf x})$ is connected, we either have $\omega({\bf x})=\{{\bf 0}\}$ or $\omega({\bf x})\subset C$. If ${\bf x}\neq{\bf 0}$, the former is impossible because the origin is unstable. Thus, the only possibility is that $\omega({\bf x})\subset C$. By the Poincaré-Bendixon theorem, $\omega({\bf x})$ is a periodic orbit; therefore, $\omega({\bf x})=C$.

The existence of a strictly orthogonal HHD is compatible with a boundary condition commonly used if the domain is properly chosen, even though boundary conditions alone are not sufficient.

Lemma 3.4**.**

Let ${\bf F}$ be a vector field defined on ${\mathbb{R}}^{2}$. Let $D\subset{\mathbb{R}}^{2}$ be a domain and $\partial D$ be a piecewise $C^{1}$ Jordan curve. It is assumed that all ${\bf x}\in\partial D$ are equilibrium points. If there is a strictly orthogonal HHD of ${\bf F}$ on $\bar{D}$, it satisfies the following on $\partial D$:

[TABLE]

where ${\bf n}$ is normal vector of $\partial D$ pointing outward.

Proof.

We have ${\bf u}({\bf y})={\bf F}({\bf y})={\bf 0}$ on $\partial D$. Therefore, ${\bf u}\cdot{\bf n}=0$. ∎

Lemma 3.5**.**

Let ${\bf F}$ be a vector field defined on ${\mathbb{R}}^{2}$. Let $D\subset{\mathbb{R}}^{2}$ be a domain and $\partial D$ be a piecewise $C^{1}$ Jordan curve. It is assumed that there is a point ${\bf x}\in D$ such that $\omega({\bf x})=\partial D$. If there is a strictly orthogonal HHD of ${\bf F}$ on $\bar{D}$, it satisfies the following on $\partial D$:

[TABLE]

where ${\bf n}$ is normal vector of $\partial D$ pointing outward.

Proof.

Main Theorem 3 implies that $\nabla V={\bf 0}$ on $\partial D$. Therefore, ${\bf u}({\bf y})={\bf F}({\bf y})$ on $\partial D$ and ${\bf u}$ is tangent to $\partial D$. Thus, ${\bf u}\cdot{\bf n}=0$. ∎

Under the assumptions in Theorem 3.5, a strictly orthogonal HHD is unique.

Corollary 2**.**

Let ${\bf F}$ be a vector field defined on ${\mathbb{R}}^{2}$. Let $D\subset{\mathbb{R}}^{2}$ be a domain and $\partial D$ be a partially $C^{1}$ Jordan closed curve. It is assumed that there is a point ${\bf x}\in D$ such that $\omega({\bf x})=\partial D$. Then, a strictly orthogonal HHD is unique if it exists.

Proof.

The HHD that satisfies ${\bf u}\cdot{\bf n}=0$ on $\partial D$ is unique (a proof may be found in [5]). ∎

Uniqueness is useful for constructing strictly orthogonal HHD. In Corollary 2, the domain $D$ under consideration is determined in terms of the dynamics of the original vector field and there is no need for optimization. If a strictly orthogonal HHD exists, it will be constructed by applying the boundary condition ${\bf u}\cdot{\bf n}=0$ on $\partial D$, owing to uniqueness.

4. Analysis of planar vector fields

In this section, to study the limitations of the analysis based on the potential functions of HHDs, planar vector fields are considered. In particular, the effect of boundary conditions on the properties of potential functions are studied. The results obtained here indicate that if the derivative of the vector field is degenerate at the equilibrium point, a Lyapunov function cannot be easily constructed by choosing a potential function.

For planar vector fields with isolated equilibrium points, the Poisson equation can be explicitly solved in neighborhoods of such points, and a formula for representing the derivatives of the potential function of the HHD can be obtained in terms of the Fourier coefficients of the boundary condition. Here, the following property is used: the addition of a harmonic function to the potential function of an HHD yields another HHD. For simplicity, in what follows, it is assumed that the vector field has an equilibrium point at the origin. The Poisson kernel and its orbital derivative are first defined.

Definition 4.1**.**

The Poisson kernel on the unit disk is defined by

[TABLE]

where $\theta\in[0,2\pi)$ is a parameter. For a vector field ${\bf F}:=(f_{1},f_{2})$ on the unit disk, let

[TABLE]

In terms of the vector field, $L[{\bf F}]$ is given explicitly as follows.

Lemma 4.2**.**

For a vector field ${\bf F}$ on the unit disk, $L[{\bf F}](x,y,\theta)$ is explicitly given by

[TABLE]

where ${\bf x}=(x,y)$ and ${\bf e}_{\theta}=\left(\cos(\theta),\sin(\theta)\right)$.

Representing a harmonic function on the unit disk using the Poisson kernel, we obtain the following result by interchanging differentiation and integration.

Lemma 4.3**.**

If a harmonic function on the unit disk is given by

[TABLE]

where $\alpha$ is a bounded function on the unit circle $S^{1}$, then

[TABLE]

The addition of harmonic functions changes the choice of the decomposition. If an HHD is given by ${\bf F}=-\nabla V+{\bf u}$, the addition of a harmonic function $h$ yields another decomposition of the form ${\bf F}=-\nabla\left(V+h\right)+\left({\bf u}+\nabla h\right)$. Therefore, it suffices to determine a harmonic function $h$ so that $V+h$ is a Lyapunov function of the origin, where $V$ is the potential function of the initial decomposition. The aim here is to calculate the value of the partial differentials of $V+h$ and ${\dot{V}}+{\dot{h}}$ at the origin using Lemma 4.3 and thereby establish that $V+h$ attains a minimum, whereas ${\dot{V}}+{\dot{h}}$ attains a maximum at the origin. The following result can be shown by direct calculation.

Lemma 4.4**.**

The values of the partial differentials of $K$ at the origin are

[TABLE]

[TABLE]

[TABLE]

[TABLE]

[TABLE]

By Lemma 4.4, we obtain formulas for the values of the derivatives of a harmonic function at the origin by interchanging differentiation and integration.

Proposition 2**.**

Let a harmonic function on the unit disk be given by

[TABLE]

where $\alpha$ has a Fourier series expansion

[TABLE]

Then, the values of the partial differentials of $h$ at the origin are

[TABLE]

[TABLE]

[TABLE]

[TABLE]

[TABLE]

Combining the results above, we obtain a sufficient condition for the Fourier coefficients whereby $V+h$ attains a minimum at the origin.

Proposition 3**.**

Let $V$ be a potential function of an HHD and $h$ be a harmonic function that satisfies the condition in Proposition 2. If $V$ assumes a minimum at the origin and $a_{1}=b_{1}=0$, then the origin is a critical point of $V+h$. Furthermore, if $V$ is the potential function constructed in Main Theorem 1, then the condition

[TABLE]

ensures that $V+h$ also assumes a minimum at ${\bf 0}$.

Proof.

By Proposition 2, we immediately obtain that $\nabla(V+h)({\bf 0})={\bf 0}$ if $a_{1}=b_{1}=0$. Let V be the potential function constructed in Main Theorem 1. The Hessian of $V+h$ at the origin is

[TABLE]

where $\rho={\rm tr}D{\bf F}_{0}={\rm tr}\,{\rm Hess}\,V$. As the trace of this matrix is positive, it is positive definite if its determinant is positive. This is equivalent to $\frac{\rho^{2}}{4}-16\left(a_{2}^{2}+b_{2}^{2}\right)>0$. ∎

${\dot{V}}+{\dot{h}}$ is now considered. By calculations similar to those above, the following results are obtained.

Lemma 4.5**.**

The values of the partial differentials of $L[{\bf F}]$ at the origin are

[TABLE]

[TABLE]

[TABLE]

[TABLE]

[TABLE]

[TABLE]

[TABLE]

Proposition 4**.**

Let a harmonic function on the unit disk be given by

[TABLE]

where $\alpha$ has a Fourier series expansion

[TABLE]

Then, the values of partial differentials of ${\dot{h}}$ at the origin are

[TABLE]

[TABLE]

[TABLE]

[TABLE]

[TABLE]

[TABLE]

Thus, the problem of obtaining a decomposition that yields a Lyapunov function is reformulated as an optimization problem for Fourier coefficients. By Proposition 4, the following result is immediately obtained.

Proposition 5**.**

Let ${\bf F}:{\mathbb{R}}^{2}\to{\mathbb{R}}^{2}$ be a vector field with an equilibrium point at the origin, and let ${\bf F}=-\nabla V+{\bf u}$ be an HHD and $h$ be a harmonic function that satisfies the condition in Proposition 2. If $V$ assumes a minimum at [math] and $a_{1}=b_{1}=0$, then [math] is a critical point of $\dot{V}+\dot{h}$.

Proof.

From Proposition 1, we have $\nabla\dot{V}({\bf 0})=0$. Therefore, it suffices to show that $\nabla\dot{h}({\bf 0})=0$, but this is obvious from Proposition 4. ∎

To obtain a Lyapunov function of the origin, it suffices that $\dot{V}+\dot{h}\leq 0$. Thus, the Fourier coefficients $a_{2},b_{2}$ should satisfy the condition (6) so that the Hessian of $\dot{V}+\dot{h}$ is negative definite. However, this is not always possible because in the calculation above, only the first-order derivatives were considered. For example, the negative definiteness of $\dot{V}+\dot{h}$ cannot be ensured in the case of $D{\bf F}_{\bf 0}=0$. This is a limitation of the construction of Lyapunov functions based on the potential functions of HHDs. In general, higher-order derivatives and, therefore, Fourier coefficients for larger $n$ should be considered for such equilibrium points.

5. Concluding remarks

Even though the results obtained in this study cannot be applied to “real” problems, the construction of Lyapunov functions using HHD was shown to be possible. There are two questions to be answered regarding the general construction of Lyapunov functions of vector fields using HHD:

1. (1)

To what extent is the Lyapunov function globally valid if it is constructed by the proposed method? 2. (2)

Under what conditions is it possible to construct Lyapunov functions using Main Theorem 2? In that case, how can it be efficiently constructed? 3. (3)

Is it possible to construct Lyapunov functions using HHD in the case of vector fields with degenerate Jacobian matrix?

The first question is important in the applications because a global Lyapunov function would yield information about the global behavior of solutions. However, the construction in this study does not ensure the validity outside a neighborhood of the equilibrium point. The second question should be answered if the method proposed here is to be applied to real problems. A certain stability condition is clearly necessary; however, it is not obvious. Furthermore, an optimization problem should be solved to apply Main Theorem 2, and therefore an efficient scheme is required. The third question concerns the superiority of the analysis based on HHD compared to the existing methods (e.g., linearization). An affirmative answer would imply that the stability of an equilibrium point could be determined for a wider class of vector fields.

Moreover, the properties of vector fields with strictly orthogonal HHD should be further studied. In particular, conditions ensuring its existence are particularly important, as such vector fields can be rigorously analyzed in detail. Another problem is to refine Main Theorem 3 so that it provides a characterization of $\alpha$- or $\omega$-limit sets.

For the problem of choosing an HHD that respects the dynamics of vector fields, a decomposition that yields a Lyapunov function is a potential solution. However, this requires further investigation.

Acknowledgment

I would like to express my gratitude to Professor Masashi Kisaka for his patient guidance and useful critique of this research. This study was supported by Grant-in-Aid for JSPS Fellows (17J03931). I would like to thank Editage (www.editage.jp) for English language editing. Further, I would like to thank the anonymous reviewers for their valuable comments and suggestions.