# Energy function of 2D and 3D coarse systems

Jean-Marc Ginoux, Riccardo Meucci, Jaume Llibre, Julien Clinton Sprott

► To cite this version:

Jean-Marc Ginoux, Riccardo Meucci, Jaume Llibre, Julien Clinton Sprott. Energy function of 2D and 3D coarse systems. Chaos, Solitons & Fractals, 2025, 199, pp.116643. 10.1016/j.chaos.2025.116643. hal-05117743

HAL Id: hal-05117743

https://hal.science/hal-05117743v1

Submitted on 18 Jun 2025

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
<figure><img src="image_2.png" /></figure>

<figure><img src="image_3.png" /></figure>

# Energy function of 2D and 3D coarse systems

Jean-Marc Ginoux<sup>a ID,*</sup>, Riccardo Meucci<sup>a</sup>, Jaume Llibre<sup>b ID</sup>, Julien Clinton Sprott<sup>c ID</sup>

## ARTICLE INFO

### Keywords:

Energy function

Coarse systems

Zero-Hopf bifurcation

Van der Pol system

Li-Sprott model

## ABSTRACT

In this work, while using the Flow Curvature Method developed by one of us (JMG), we prove that the energy function of two and three-dimensional coarse systems involving a small parameter $\mu$ can be directly deduced from the curvature of their trajectory curves when $\mu$ tends to zero. Such a result thus confirms the relationship between curvature and energy function for a certain class of differential systems already established in one of our previous contributions. Then, we state that the rate of change of the energy function of such coarse systems is equal to the scalar product of the velocity vector field and its first time derivative, i.e. the acceleration vector field. The comparison of these results with the so-called Frénet frame enables to prove that energy function is proportional to the normal component of the acceleration when $\mu$ tends to zero while the rate of change of the energy function is proportional to the tangential component of the acceleration at first order in $\mu$. Two and three-dimensional examples are then used to emphasize these two main results.

## 1. Introduction

In the 1930's the famous soviet mathematicians Aleksandr' Andronov and his colleague Lev Pontryagin introduced the concept of “coarse systems” also known as “rough systems” or “systèmes grossiers” (in French) in dimension two, see [1]. Here, the terms “coarse systems” and “roughness” are used as synonymous in a broader sense, e.g. to mean merely the preservation of some property of the system under a small perturbation. Thus, in such a case it is preferable to speak about the structural stability of the property in question. According to Andronov et al. [2]:

> “In setting out the differential equations we cannot take account of all the factors that influence in some manner or other the behaviour of the physical system. On the other hand, none of the factors taken into account can remain absolutely constant during a motion of the system, so that when we attribute to the parameters perfectly determined values; this has only a meaning on condition that small variations of the parameters do not substantially vary the character of the motion. A certain number of parameters corresponding to physical parameters of the problem occur in the functions of our system equations, so these functions are never known exactly. Small variations of these parameters must leave unchanged the qualitative structure of the phase portrait. It is natural, therefore, to separate the class of dynamic systems whose topological structure of the

phase paths does not vary for small variations of the differential equations. We call such systems ‘coarse’ or structurally stable ...”

In this work, we will follow such a definition of coarse systems and we will make therefore the assumption that the qualitative nature of the phase portrait of the systems of differential equations (1)–(2) defined below remains unchanged under a small variation of the coefficient $\mu$ in factor of the polynomials functions $f, g$ and $h$.

### 1.1. Two and three-dimensional coarse systems

In dimension two the “coarse systems” here studied are the following:

$$ 
\begin{aligned}
\frac{dx}{dt} &= y + \mu f(x, y, \mu), \\
\frac{dy}{dt} &= -x + \mu g(x, y, \mu),
\end{aligned}
\qquad (1) $$

where $0 < \mu \ll 1$ and $f(x, y, \mu)$ and $g(x, y, \mu)$ are functions of class $C^2$. When $f(x, y, \mu) = x - x^3/3$ and $g(x, y, \mu) = 0$, Eqs. (1) represent the so-called Van der Pol system. Here, in dimension three, we consider

<footnote><sup>a</sup> Université de Toulon, Aix Marseille Univ, CNRS, LIS, Toulon, France</footnote>

<footnote><sup>b</sup> Departament de Matemàtiques, Universitat Autònoma de Barcelona, 08193 Bellaterra, Barcelona, Spain</footnote>

<footnote><sup>c</sup> University of Wisconsin 1150 University Avenue Madison, WI 53706-1390, USA</footnote>

<footnote>* Corresponding author.</footnote>

<footnote>E-mail address: jmginoux@orange.fr (J.-M. Ginoux).</footnote>
<figure><img src="image_4.png" /><figcaption>Fig. 1. Frénet trihedron.</figcaption></figure>

the “coarse systems” as follows:

$$ 
\begin{aligned}
\frac{dx}{dt} &= y + \mu f(x, y, z, \mu), \\
\frac{dy}{dt} &= -x + \mu g(x, y, z, \mu), \\
\frac{dz}{dt} &= \mu h(x, y, z, \mu),
\end{aligned}
\qquad (2) $$

where $0 < \mu \ll 1$ and $f(x, y, z, \mu)$, $g(x, y, z, \mu)$ and $h(x, y, z, \mu)$ are functions of class $C^3$. System (2) has been studied by Llibre and Zhang [3,4].

Very recently, we have proposed a method for obtaining the energy function of differential systems starting from their corresponding equation of motion, i.e. their acceleration equation in dimension two and jerk equation in dimension three [5]. Two years ago, we had established a link in dimension two between the energy function of the generalized Liénard systems which can be considered as “coarse systems” in the particular case where $g(x, y, \mu) = 0$ and the curvature of their trajectory curve [6]. In this work, we extend this result to dimension three and thus we prove that the energy function of two and three-dimensional “coarse systems” can be directly deduced from the curvature of the trajectory curve of such systems. The concept of curvature is one of the foundations of the Flow Curvature Method briefly recalled below.

## 1.2. Flow curvature method

Let us first consider the following set of nonlinear ordinary differential equations:

$$ \frac{d\vec{X}}{dt} = \vec{\mathfrak{F}}(\vec{X}), \qquad (3) $$

with $\vec{X} = [x_1, x_2, \dots, x_n]^t \in E \subset \mathbb{R}^n$ and $\vec{\mathfrak{F}}(\vec{X}) = [f_1(\vec{X}), f_2(\vec{X}), \dots, f_n(\vec{X})]^t \in E \subset \mathbb{R}^n$. The vector $\vec{\mathfrak{F}}(\vec{X})$ defines a velocity vector field $\dot{\vec{X}} = [\dot{x}_1, \dot{x}_2, \dots, \dot{x}_n]^t$ in $E$ whose components $f_i$ which are supposed to be infinitely differentiable with respect to all $x_i$, i.e. are $C^\infty$ functions in $E$ and with values included in $\mathbb{R}$, satisfy the assumptions of the Cauchy-Lipschitz theorem. For more details, see for example Coddington & Levinson [7]. A solution of this system is a trajectory curve $\vec{X}(t)$ tangent (except at the fixed points) to $\vec{\mathfrak{F}}$ whose values define the states of the dynamical system described by the Eq. (3).

Then, let us recall the definition of the so-called Frénet frame or trihedron. The Frénet frame [8] built from the space curve $\vec{X}(t)$ and directed towards the motion of mobile M consists in a unit tangent vector $\vec{\tau}$ to the space curve in M, a unit normal vector $\vec{n}$, i.e. the principal normal in M directed towards the interior of the concavity of the space curve and a unit binormal vector $\vec{\beta}$ to the space curve in M so that the trihedron $(\vec{\tau}, \vec{n}, \vec{\beta})$ is direct since the velocity vector field $\vec{V} = \dot{\vec{X}}$ is tangent to any point M to the space curve $\vec{X}(t)$ (see Fig. 1 below).

By taking the time derivative of the velocity vector field (3), one can deduce the acceleration vector field $\ddot{\vec{X}}$ of the trajectory curve $\vec{X}(t)$. In the Frénet frame, this acceleration has two components according to the unit tangent vector $\vec{\tau}$ and to the unit normal vector $\vec{n}$ to the trajectory curve: $\ddot{\vec{X}} = \ddot{\vec{X}}_r + \ddot{\vec{X}}_n$, where

* the tangential acceleration $\ddot{\vec{X}}_r$ reads:

$$ \ddot{\vec{X}}_r = \frac{\dot{\vec{X}} \cdot \ddot{\vec{X}}}{\|\dot{\vec{X}}\|} \vec{\tau}, \qquad (4) $$

where here and in what follows $\|\cdot\|$ represents the Euclidean norm. Let us recall that the time derivative of the norm of any vector field $\vec{U}$ reads:

$$ \frac{d}{dt} \|\vec{U}\| = \frac{\vec{U} \cdot \dot{\vec{U}}}{\|\vec{U}\|}. $$

It follows that the tangential acceleration and, as we will see below, the rate of change of the energy function is proportional to the time derivative of the norm of the velocity vector field.

* the normal acceleration $\ddot{\vec{X}}_n$ reads:

$$ \ddot{\vec{X}}_n = \frac{\|\dot{\vec{X}} \wedge \ddot{\vec{X}}\|}{\|\dot{\vec{X}}\|} \vec{n}. \qquad (5) $$

According to the Flow Curvature Method developed by Ginoux et al. [9,10] and Ginoux [11], the flow curvature manifold is defined as the location of the points where the curvature of the flow, i.e., the curvature of the trajectory curves of any n-dimensional differential system (3) vanishes. Then, it reads:

$$ \kappa(\vec{X}) = \| \dot{\vec{X}} \wedge \ddot{\vec{X}} \| = 0. \qquad (6) $$

The curvature is a measure of the rate of change of direction of a mobile point M evolving on the trajectory curve. By comparing Eqs. (5) and (6), it follows that curvature expresses the magnitude of the normal acceleration of such a point. Let us notice that in dimension two, $\kappa(\vec{X}) = |\det(\dot{\vec{X}}, \ddot{\vec{X}})| = 0$.

## 1.3. Energy function and its rate of change

Usually, to obtain the energy function and its rate of change one computes the corresponding equation of motion of two and three-dimensional “coarse systems” (1)–(2), i.e. respectively their acceleration and jerk equation. In dimension two, it is far well-known that energy function of differential systems can be simply obtained by multiplying its corresponding second-order ordinary differential equation, i.e., its equation of motion by the first time derivative of its state variable. In dimension three, we have proved in a previous work [5] that the jerk equation of such differential systems must be multiplied by the second time derivative of the state variable and not by the first like in dimension two to obtain the energy function. In this work, we will see below that in the case of three-dimensional “coarse systems” (2), its jerk equation reduces to its acceleration equation like for two-dimensional “coarse systems” (1).

However, more than twenty years ago, Sarasola et al. [12], Wang et al. [13], Ma et al. [14] and Yu et al. [15] to name but a few proposed to use a method based on the Helmholtz’s theorem to deduce the energy function of any n-dimensional differential systems (3). This method is briefly recalled below. According to Helmholtz’s theorem, such differential system (3) can be decomposed into a gradient and a rotational field, i.e., as the sum of a conservative and a dissipative vector field $\vec{\mathfrak{F}}(\vec{X}) = \vec{\mathfrak{F}}_c(\vec{X}) + \vec{\mathfrak{F}}_d(\vec{X})$. Thus, it can be expressed in a generalized Hamiltonian form:

$$ \frac{d\vec{X}}{dt} = [J(\vec{X}) + R(\vec{X})] \nabla H(\vec{X}), \qquad (7) $$
where $\nabla H$ is the gradient vector of the smooth energy function $H(\vec{X})$, $J(\vec{X})$ is a skew-symmetric matrix and $R(\vec{X})$ is a symmetric matrix. The Hamiltonian energy function can thus be calculated by:

$$ \frac{dH}{dt} = \nabla H^T R(\vec{X}) \nabla H, \qquad (8) \\ \nabla H^T J(\vec{X}) \nabla H = 0. $$

Since the vector field (3) can be decomposed into a gradient and a rotational field, it follows that the energy function and its rate of change can be deduced by:

$$ \frac{dH}{dt} = \nabla H^T R(\vec{X}) \nabla H = \nabla H^T \vec{\nabla}_d(\vec{X}), \qquad (9) \\ \nabla H^T J(\vec{X}) \nabla H = 0 = \nabla H^T \vec{\nabla}_c(\vec{X}). $$

In this work, we use the Flow Curvature Method, i.e. Eq. (6) and we compute the curvature of the trajectory curves of each of these two “coarse systems”. Then, by comparing with the energy function given by the two other methods (classical method and Helmholtz’s theorem), we prove that when $\mu \to 0$ the curvature of “2D & 3D coarse systems” tend to be equal to their corresponding energy function. Finally, we prove that the rate of change of the energy function of “2D & 3D coarse systems” is equal to the scalar product of its vector field and its first time derivative. Interpretations of these two main results are given in the discussion.

## 2. Acceleration equation and energy function of 2D “coarse systems”

By computing the time derivative of $\dot{x}$ and by replacing $\dot{y}$ by its expression (1), we obtain the equation of motion, i.e., the acceleration equation of the two-dimensional “coarse system”:

$$ \dot{x} + x = \mu (g(x, y, \mu) + \dot{f}(x, y, \mu)). \qquad (10) $$

### 2.1. Classical method

A classical method for obtaining the energy function $H$ of a two-dimensional differential “coarse system” consists in multiplying its acceleration equation (10) by $\dot{x}$. This directly provides its rate of change which reads:

$$ \frac{dH}{dt} = \frac{d}{dt} \left( \frac{\dot{x}^2}{2} + \frac{\dot{x}^2}{2} \right) = \mu \left( g(x, y, \mu) + \dot{f}(x, y, \mu) \right) \dot{x}, \qquad (11) $$

where the energy function is

$$ H = \frac{\dot{x}^2}{2} + \frac{\dot{x}^2}{2}. \qquad (12) $$

The energy function $H$ is analogous to Hamiltonian function. Following this idea, $\frac{\dot{x}^2}{2}$ is considered as a kinetic energy while $\frac{\dot{x}^2}{2}$ is the potential energy. This definition holds for dimension three.

### 2.2. Helmholtz’s theorem

Then, by using Helmholtz’s theorem recalled above the 2D coarse system (1) can be decomposed into a gradient and a rotational field, i.e., as the sum of conservative and dissipative vector fields as follows:

$$ \begin{pmatrix} \dot{x} \\ \dot{y} \end{pmatrix} = \begin{pmatrix} 0 & 1 \\ -1 & 0 \end{pmatrix} \begin{pmatrix} x \\ y \end{pmatrix} + \mu \begin{pmatrix} f(x, y, \mu) \\ g(x, y, \mu) \end{pmatrix} \qquad (13) $$

But in order to facilitate comparison between the energy function (12) given by the acceleration equation (10) and that given by Helmholtz’s theorem, it is necessary to use the first time derivative of system (13) and not the original system. As we will see below, it will also enable the comparison with the energy function given by the curvature. So, we have:

$$ \begin{pmatrix} \dot{x} \\ \dot{y} \end{pmatrix} = \begin{pmatrix} 0 & 1 \\ -1 & 0 \end{pmatrix} \begin{pmatrix} \dot{x} \\ \dot{y} \end{pmatrix} + \mu \begin{pmatrix} \dot{f}(x, y, \mu) \\ \dot{g}(x, y, \mu) \end{pmatrix} \qquad (14) $$

Then, according to Helmholtz’s theorem, the energy function $H$ associated with the 2D coarse system (1) must satisfy the partial differential equation:

$$ \left( \frac{\partial H}{\partial \dot{x}}, \frac{\partial H}{\partial \dot{y}} \right) \left( \dot{y} \mid \dot{x} \right) = 0 \quad \Leftrightarrow \quad \left\{ \begin{array}{l} \frac{\partial H}{\partial \dot{x}} = \dot{x}, \\ \frac{\partial H}{\partial \dot{y}} = \dot{y}, \end{array} \right. \qquad (15) $$

This leads to the energy function:

$$ H = \frac{\dot{x}^2}{2} + \frac{\dot{y}^2}{2}. \qquad (16) $$

By replacing in Eq. (16) $\dot{y}$ by its expression (1), it is easy to demonstrate that both Eqs. (12) and (16) are identical when $\mu \to 0$. From (9) and (15), the rate of change of the energy function is:

$$ \frac{dH}{dt} = \left( \frac{\partial H}{\partial \dot{x}}, \frac{\partial H}{\partial \dot{y}} \right) \mu \left( \dot{f}(x, y, \mu) \mid \dot{g}(x, y, \mu) \right) = \mu \left( \dot{x} \dot{f}(x, y, \mu) + \dot{y} \dot{g}(x, y, \mu) \right) \qquad (17) $$

Let us notice that we would have found the same result as above (17) by taking the time derivative of the energy function $H$ given either by Eq. (12) or by (16).

### 2.3. Flow curvature method

By using the Flow Curvature Method, i.e. Eq. (6), we compute the curvature of the trajectory curves of the 2D coarse system (1). We have:

$$ \kappa(x, y, \mu) = |\det(\dot{\vec{X}}, \ddot{\vec{X}})| = \dot{x}^2 + \dot{y}^2 + \mu (\dot{y}\dot{f}(x, y, \mu) - \dot{x}\dot{g}(x, y, \mu)) = 0. \qquad (18) $$

Obviously, when $\mu \to 0$, we find that:

$$ \kappa(x, y, 0) = \dot{x}^2 + \dot{y}^2 = 2H. \qquad (19) $$

Now, we compute the scalar product of the vector field $\dot{\vec{X}}$ of (1) and its time derivative $\ddot{\vec{X}}$, we obtain:

$$ \dot{\vec{X}} \cdot \ddot{\vec{X}} = \mu (\dot{x} \dot{f}(x, y, \mu) + \dot{y} \dot{g}(x, y, \mu)). \qquad (20) $$

Then, it follows that the right hand side of Eq. (20) is identical with that of Eq. (17).

A comparison of the energy function and its rate of change given by the three methods is presented in the Table 1 below.

Although the energy function $H$ and its rate of change $dH/dt$ given by the classical method (11)–(12) and those given by both Helmholtz’s theorem (16)–(17) and Flow Curvature Method (19)–(20) seem to be different, they are absolutely identical (see the proof in Appendix A). Concerning the additional term in the expression of $H$ given by the Flow Curvature Method (18), its expression and importance are given and discussed in Appendix B.

In summary, we proved the next result:

**Theorem 1.** For the 2D coarse systems (1), the energy function is given by:

$$ H = \lim_{\mu \to 0} \frac{1}{2} \kappa(x, y, \mu) = \frac{\dot{x}^2}{2} + \frac{\dot{x}^2}{2} = \frac{\dot{x}^2}{2} + \frac{\dot{y}^2}{2}, $$

and its rate of change is given by:

$$ \frac{dH}{dt} = \dot{\vec{X}} \cdot \ddot{\vec{X}}. $$

## 3. Jerk equation and energy function of three-dimensional “coarse systems”

By computing the second time derivative of $\dot{x}$ and by replacing $\dot{y}$ by its expression deduced from the time derivative of (2), we obtain the jerk equation of the three-dimensional “coarse system”:

$$ \ddot{x} + \dot{x} = \mu (\dot{g}(x, y, z, \mu) + \dot{f}(x, y, z, \mu)). \qquad (21) $$
<table><caption>**Table 1**
Comparison of the energy function and its rate of change given by the three methods.</caption><thead><tr><th></th><th>Classical method</th><th>Helmholtz's theorem</th><th>Flow curvature method</th></tr></thead><tbody><tr><td>H</td><td>$\frac{\dot{x}^2}{2} + \frac{x^2}{2}$</td><td>$\frac{\dot{x}^2}{2} + \frac{\dot{y}^2}{2}$</td><td>$x^2 + \dot{y}^2 + \mu (\dot{y}\dot{f}(x, y, \mu) - \dot{x}\dot{g}(x, y, \mu))$</td></tr><tr><td>$\frac{dH}{dt}$</td><td>$\mu (g(x, y, \mu) + \dot{f}(x, y, \mu)) \dot{x}$</td><td>$\mu (\dot{x}\dot{f}(x, y, \mu) + \dot{y}\dot{g}(x, y, \mu))$</td><td>$\mu (\dot{x}\dot{f}(x, y, \mu) + \dot{y}\dot{g}(x, y, \mu))$</td></tr></tbody></table>

But, in the case of three-dimensional “coarse systems” (2), its *jerk equation* reduces to its *acceleration equation* by a simple and direct integration. Thus, we find the following equation which is the same as for two-dimensional “coarse systems” (1):

$$ \dot{x} + x = \mu (g(x, y, z, \mu) + \dot{f}(x, y, z, \mu)). \quad (22) $$

### 3.1. Classical method

Then, by using the classical method, i.e. by multiplying this *acceleration equation* (22) by $\dot{x}$, we obtain the rate of change of the *energy function*:

$$ \frac{dH}{dt} = \frac{d}{dt} \left( \frac{\dot{x}^2}{2} + \frac{x^2}{2} \right) = \mu (g(x, y, z, \mu) + \dot{f}(x, y, z\mu)) \dot{x}, \quad (23) $$

where the *energy function* is

$$ H = \frac{\dot{x}^2}{2} + \frac{x^2}{2}. \quad (24) $$

### 3.2. Helmholtz's theorem

As previously recalled, another method consists in using Helmholtz's theorem. Thus, the 3D coarse system (2) is decomposed into a gradient and a rotational field, i.e., as the sum of conservative and dissipative vector fields as follows:

$$ \begin{pmatrix} \dot{x} \\ \dot{y} \\ \dot{z} \end{pmatrix} = \begin{pmatrix} 0 & 1 & 0 \\ -1 & 0 & 0 \\ 0 & 0 & 0 \end{pmatrix} \begin{pmatrix} x \\ y \\ z \end{pmatrix} + \mu \begin{pmatrix} f(x, y, z, \mu) \\ g(x, y, z, \mu) \\ h(x, y, z, \mu) \end{pmatrix}. \quad (25) $$

As above, in order to facilitate comparison between the *energy function* (24) given by the *acceleration equation* (22) and that given by Helmholtz's theorem, it is necessary to use the first time derivative of system (25) and not the original system. Then, it will also enable the comparison with the *energy function* given by the *curvature*. We take the first time derivative of this Eq. (25), and we have:

$$ \begin{pmatrix} \dot{x} \\ \dot{y} \\ \dot{z} \end{pmatrix} = \begin{pmatrix} 0 & 1 & 0 \\ -1 & 0 & 0 \\ 0 & 0 & 0 \end{pmatrix} \begin{pmatrix} \dot{x} \\ \dot{y} \\ \dot{z} \end{pmatrix} + \mu \begin{pmatrix} \dot{f}(x, y, z, \mu) \\ \dot{g}(x, y, z, \mu) \\ \dot{h}(x, y, z, \mu) \end{pmatrix} \quad (26) $$

Then, according to the Helmholtz's theorem, the *energy function* $H$ associated with the 3D coarse system (2) must satisfy the partial differential equation:

$$ \left( \frac{\partial H}{\partial \dot{x}}, \frac{\partial H}{\partial \dot{y}}, \frac{\partial H}{\partial \dot{z}} \right) \begin{pmatrix} \dot{y} \\ -\dot{x} \\ 0 \end{pmatrix} = 0 \quad \Leftrightarrow \quad \begin{cases} \frac{\partial H}{\partial \dot{x}} = \dot{x}, \\ \frac{\partial H}{\partial \dot{y}} = \dot{y}, \\ \frac{\partial H}{\partial \dot{z}} = 0. \end{cases} \quad (27) $$

This leads to the *energy function*:

$$ H = \frac{\dot{x}^2}{2} + \frac{\dot{y}^2}{2}. \quad (28) $$

Then, by replacing in Eq. (28) $\dot{y}$ by its expression (2), it is easy to demonstrate that both Eqs. (24) and (28) are identical when $\mu \to 0$. From (9) and (27), the rate of change of the *energy function* is given by:

$$ \begin{aligned} \frac{dH}{dt} &= \left( \frac{\partial H}{\partial \dot{x}}, \frac{\partial H}{\partial \dot{y}}, \frac{\partial H}{\partial \dot{z}} \right) \mu \begin{pmatrix} \dot{f}(x, y, z, \mu) \\ \dot{g}(x, y, z, \mu) \\ \dot{h}(x, y, z, \mu) \end{pmatrix} \\ &= \mu (\dot{x}\dot{f}(x, y, z, \mu) + \dot{y}\dot{g}(x, y, z, \mu)). \end{aligned} \quad (29) $$

Let us notice that we would have found the same result as above (29) by taking the time derivative of the *energy function* $H$ given either by Eq. (24) or by (28).

### 3.3. Flow curvature method

Now, by using the *Flow Curvature Method*, i.e. Eq. (6), we compute the *curvature* of the trajectory curves of the 3D coarse system (2), and we have:

$$ \kappa(x, y, z, \mu) = \sqrt{[\dot{x}^2 + \dot{y}^2 + \mu (\dot{y}\dot{f}(x, y, z, \mu) - \dot{x}\dot{g}(x, y, z, \mu))]^2 + O(\mu^2)} = 0. \quad (30) $$

Thus, obviously, when $\mu \to 0$, we find that:

$$ \kappa(x, y, z, 0) = \dot{x}^2 + \dot{y}^2 = 2H. \quad (31) $$

Now, we compute the scalar product of the vector field $\dot{X}$ of (2) and its time derivative $\ddot{X}$, we obtain:

$$ \dot{X} \cdot \ddot{X} = \mu (\dot{x}\dot{f}(x, y, z, \mu) + \dot{y}\dot{g}(x, y, z, \mu)) + \mu^2 h(x, y, z, \mu)h(x, y, z, \mu). \quad (32) $$

Then, it follows that the right hand side of Eq. (32) is identical with that of Eq. (29) provided the second order terms in $\mu$ are neglected. A comparison of the *energy function* and its rate of change given by the three methods would have lead to the same results as those presented in Table 1 (except for the arguments of functions). In fact, here again, the *energy function* $H$ and its rate of change $dH/dt$ given by the classical method (23)–(24) and those given by both Helmholtz's theorem (28)–(29) and *Flow Curvature Method* (31)–(32) seem to be different. Nevertheless, as highlighted in Appendix A, they are absolutely identical. Concerning the additional term in the expression of $H$ given by the *Flow Curvature Method* (30), its expression and importance are given and discussed in Appendix C. In summary, we proved the next result:

**Theorem 2.** For the 3D coarse systems (2) the *energy function* is:

$$ H = \lim_{\mu \to 0} \frac{1}{2} \kappa(x, y, z, \mu) = \frac{\dot{x}^2}{2} + \frac{\dot{y}^2}{2}, $$

and its rate of change is:

$$ \frac{dH}{dt} = \dot{X} \cdot \ddot{X}, $$

at first order in $\mu$.

### 4. Applications

In this section we propose two applications of these main results to two and three-dimensional “coarse systems”.

#### 4.1. Van der Pol system

The Van der Pol system [16] can be written as follows:

$$ \begin{aligned} \frac{dx}{dt} &= y + \mu f(x), \\ \frac{dy}{dt} &= -x, \end{aligned} \quad (33) $$

where $f(x, y, \mu) = f(x) = x - x^3/3$ is a cubic function, $g(x, y, \mu) = 0$ and $\mu$ is a small real positive parameter. While using either the classical
<figure><img src="image_5.png" /><figcaption>Fig. 2. Phase portraits, curvature and energy of Van der Pol model (33) for various values $\mu$.</figcaption></figure>

method (11)–(12) or Helmholtz’s theorem (16)–(17), the energy function and its rate of change read:

$$H = \frac{\dot{x}^2}{2} + \frac{\dot{y}^2}{2} \quad ; \quad \frac{dH}{dt} = \mu (1 - x^2) \dot{x}^2. \qquad (34)$$

According to the Flow Curvature Method, the curvature of the trajectory curves (18) of the Van der Pol “coarse system” (33) provides the following result.

$$\kappa(x, y, \mu) = |\det(\dot{\bar{X}}, \ddot{\bar{X}})| = \dot{x}^2 + \dot{y}^2 + \mu (x^2 - 1) \left(y + \mu \left(x - \frac{x^3}{3}\right)\right) x. \quad (35)$$

where we denote the additional term in the curvature (see Appendix B) as:

$$\Delta(x, y, \mu) = \mu (x^2 - 1) \left(y + \mu \left(x - \frac{x^3}{3}\right)\right) x \qquad (36)$$

Thus, when $\mu \to 0$, $\kappa(x, y, 0) = 2H$. Then, the scalar product of the vector field $\dot{\bar{X}}$ of (33) by its time derivative $\ddot{\bar{X}}$ reads:

$$\dot{\bar{X}} \cdot \ddot{\bar{X}} = \mu (1 - x^2) \dot{x}^2. \qquad (37)$$

Identity between Eqs. (34)–(37) is obvious.

In Figs. 2, have been plotted in red the Poincaré’s limit cycle solution of the Van der Pol system (33), in blue the curvature (35) and in black the energy function (34) for various values of $\mu$. We observe in Fig. 2a that for $\mu = 0.01$, these three curves are superimposed. In Fig. 2b, for $\mu = 0.1$, we notice that curvature slightly deviates from the energy function. In Fig. 2c, for $\mu = 0.5$, these two curves are moving very far apart. When $\mu = 5$, we observe in Fig. 2d that both energy function and curvature are further and further apart. We also notice that for such a large value of $\mu$ Van der Pol’s limit cycle tends asymptotically to a portion of the left and right part of the curvature which then corresponds to the so-called slow invariant manifold. Moreover, according to Diliberto [17], “a periodic solution is stable, if at each point it is in a region where the curvature is negative”. In the case of the Van der Pol system (33), the meshed region (in left and right parts of Fig. 2d) delimited by the inside of limit cycle and the vertical asymptotes $x = \pm 1$ corresponds to the location of the points where the curvature is negative, i.e. according to Diliberto where the periodic solution is stable. We recall that it has been established by Ginoux et al. [9,10], by Ginoux [11,18] and then by Ginoux et al. [5,6] that the curvature of the flow provides the slow invariant manifold equation of Van der Pol’s slow-fast system. Then, it follows that outside the meshed region (plotted
<figure><img src="image_6.png" /><figcaption>Fig. 3. Bifurcation diagram of the difference (36) between curvature and energy function of the Van der Pol system (33).</figcaption></figure>

in Fig. 2d) the curvature is therefore positive. It corresponds to the fast evolution of the trajectory curve which is therefore unstable. Moreover, according to (37), the rate of change of energy is positive inside the meshed region in Fig. 2d. This confirms a result already claimed by Minorsky [19] and established by Ginoux et al. [6], i.e., “each time the system absorbs energy the curvature of its trajectory decreases and vice versa”. Thus, we can summarize this two main results as follows: the location of the points where the curvature is negative and where its rate of change is positive corresponds to the location of the points where the periodic solution is stable and where the system absorbs energy.

Then, in order to evaluate the robustness of the curvature-energy relationship between the energy function (34) and the curvature (35) of Van der Pol system (33), we built the bifurcation diagram (presented below in Fig. 3). In this diagram we have plotted the difference Δ (36) between the curvature (35) and the energy function (34) for the maximal values of x(t) and y(t) as a function of μ (see also Eqs. (B.1)–(B.4) in Appendix B).

We observe in Figs. 3 that the discrepancies between the curvature (35) and energy function (34) of the Van der Pol system (33) is increasing according to μ. We notice that when μ becomes greater than μ ≈ 0.1 the difference (36) becomes more and more important.

## 4.2. Modified Li-Sprott model

As a first example of three-dimensional “coarse system”, we have slightly modified a differential system originally designed by Li & Sprott [20] and we obtained:

$$ 
\begin{aligned}
\frac{dx}{dt} &= y + \mu yz, \\
\frac{dy}{dt} &= -x + \mu (yz - xz), \\
\frac{dz}{dt} &= \mu (-z - xy + b),
\end{aligned} 
\quad (38) $$

where 0 < μ ≪ 1 and b a positive real parameter. We notice that for μ = 1/10 and b = 3 and for two different set of initial conditions (x₀, y₀, z₀) = (-0.5, 0.1, 0.7) and (x₀, y₀, z₀) = (0.5, 0.1, 1.7), this model (38) enables to highlight the coexistence of several attractors (see Fig. 4).

While using either the classical method (23)–(24) or Helmholtz’s theorem (28)–(29), the energy function and its rate of change read:

$$ H = \frac{\dot{x}^2}{2} + \frac{\dot{y}^2}{2}, \quad \frac{dH}{dt} = \mu x^2 z + O(\mu^2). \quad (39) $$

<figure><img src="image_7.png" /><figcaption>Fig. 4. Chaotic attractor of the modified Li-Sprott model (38).</figcaption></figure>

According to the Flow Curvature Method, the curvature of the trajectory curves (30) of the Li-Sprott “coarse system” (38) provides the following result.

$$ \kappa(x, y, z, \mu) = \sqrt{[\dot{x}^2 + \dot{y}^2 + \mu(\dot{y}\dot{f}(x, y, z, \mu) - \dot{x}\dot{g}(x, y, z, \mu))]^2 + O(\mu^2)} = 0. \quad (40) $$

where μ (ẏf(x, y, z, μ) - ẋg(x, y, z, μ)) = Δ(x, y, z, μ) and

$$ \Delta(x, y, z, \mu) = \mu \left[ (x + y + \mu xz)(xz + \mu(x - y)z^2 + \mu y(z + xy - b)) + y(1 + \mu z)(yz + \mu yz^2 - \mu xz - \mu x^2 y + \mu bx) \right] \quad (41) $$

Thus, when μ → 0, κ(x, y, z, 0) = 2H. Then, the scalar product of the vector field of (38), i.e. Ẋ and its first time derivative, i.e. Ẋ reads:

$$ \dot{X} \cdot \ddot{X} = \mu x^2 z + O(\mu^2). \quad (42) $$
<figure><img src="image_8.png" /><figcaption>Fig. 5. Bifurcation diagram of the difference (41) between curvature and energy function of the Li-Sprott model (38).</figcaption></figure>

Identity between Eqs. (39)–(42) is again obvious. Then, in order to evaluate the robustness of the curvature-energy relationship between the energy function (39) and the curvature (40) of Li-Sprott model (38), we built the bifurcation diagram (presented below in Fig. 5). In this diagram we have plotted the difference Δ (41) between the curvature (40) and the energy function (39) for the maximal values of x(t), y(t), z(t) as a function of μ (see also Eqs. (C.1)–(C.3) in Appendix C)

We observe in Figs. 5 that the discrepancies between the curvature (40) and energy function (39) of the Li-Sprott model (38) is increasing exponentially according to μ. We notice that for small μ regime, i.e. for μ ≪ 1 the difference (41) stay at a low level such that both curvature and energy can be considered as identical.

Finally, we recall that a zero-Hopf equilibrium of the differential system (38) is an equilibrium point such the Jacobian matrix of the differential system evaluated at it has eigenvalues 0 and ±ωi with ω > 0. The unique equilibrium point (0, 0, h) of the differential system (38) is zero-Hopf when μ = 0, because its eigenvalues are 0 and ±i.

The next result shows that there is a zero-Hopf bifurcation at the zero-Hopf equilibrium (0, 0, h), where an unstable periodic orbit bifurcates from this equilibrium.

**Theorem 3.** For $|\mu|$ sufficiently small and $h\mu < 0$ a periodic orbit $(x(t, \mu), y(t, \mu), z(t, \mu))$ bifurcates from the zero-Hopf equilibrium $(0, 0, h)$ of the differential system (38) satisfying that $(x(0, \mu), y(0, \mu), z(0, \mu)) = (O(\mu), O(\mu^2), h + O(\mu))$.

**Proof.** We write the differential system (38) in cylindrical coordinates $(r, \theta, z)$ where $x = r \cos \theta$, $y = r \sin \theta$, and we get

$$ 
\begin{aligned}
\dot{r} &= \mu r z \sin^2 \theta, \\
\dot{\theta} &= -1 + \mu z(\sin \theta \cos \theta - 1), \\
\dot{z} &= \mu(h - z - r^2 \sin \theta \cos \theta).
\end{aligned} 
\quad (43) $$

Taking as the new independent variable $\theta$ the differential system (43) becomes

$$ 
\begin{aligned}
\frac{dr}{d\theta} &= r' = -\mu r z \sin^2 \theta + O(\mu^2) = \mu F_1(\theta, r, z) + O(\mu^2), \\
\frac{dz}{d\theta} &= z' = -\mu(h - z - r^2 \sin \theta \cos \theta) + O(\mu^2) = \mu F_2(\theta, r, z) + O(\mu^2).
\end{aligned} 
\quad (44) $$

Now we apply to the differential system (44) the averaging theory of first order described in the Theorems 11.5 and 11.6 of [21]. Then, the first averaged function $f(r, \theta) = (f_1(r, \theta), f_2(r, \theta))$ is

$$ 
\begin{aligned}
f_1(r, z) &= \frac{1}{2\pi} \int_0^{2\pi} F_1(\theta, r, z) d\theta = -\frac{1}{2}rz, \\
f_2(r, z) &= \frac{1}{2\pi} \int_0^{2\pi} F_2(\theta, r, z) d\theta = z - h,
\end{aligned} 
$$

<figure><img src="image_9.png" /><figcaption>Fig. 6. A periodic orbit of Theorem 3.</figcaption></figure>

The unique solution of the system $f_1(r, z) = 0, f_2(r, z) = 0$, is $(r, z) = (0, h)$. The Jacobian matrix of the function $f(r, z)$ evaluated at the solution $(0, h)$ is

$$ A = \begin{pmatrix} -1/2 & 0 \\ 0 & 1 \end{pmatrix}. $$

Since the determinant of A is non-zero, by Theorem 11.5 of [21], the solution $(0, h)$ provides for $|\mu|$ sufficiently small a periodic orbit $(r(\theta, \mu), z(\theta, \mu))$ of the differential system (44) such that $(r(0, \mu), z(0, \mu)) = (O(\mu), h + O(\mu))$.

The periodic orbit $(r(\theta, \mu), z(\theta, \mu))$ of the differential system (44) in the differential system (43) becomes the periodic orbit $(x(t, \mu), y(t, \mu), z(t, \mu))$ such that $(x(0, \mu), y(0, \mu), z(0, \mu)) = (O(\mu), O(\mu^2), h + O(\mu))$. By Theorem 11.6 of [21], since one of the eigenvalues of the matrix A is negative, the periodic orbit $(x(t, \mu), y(t, \mu), z(t, \mu))$ is unstable. This completes the proof of the theorem. □

In Fig. 6 we plot the periodic orbit $(x(t, \mu), y(t, \mu), z(t, \mu))$ of Theorem 3 for the values of $\mu = 1/500$ and $h = 1$ computed numerically from the initial conditions $(x(0, \mu), y(0, \mu), z(0, \mu)) = (O(\mu), O(\mu^2), h + O(\mu))$.

In Fig. 7 we plot the periodic orbit of model (38) in the case of Theorem 3 for the values of $\mu = 1/20, b = 0$ and $h = 1$ computed numerically from the initial conditions $(x(0), y(0), z(0)) = (0.01, 0.01, h)$.

## 5. Discussion

In this work we proved that the energy function of two and three-dimensional coarse systems can be directly deduced from the so-called
<figure><img src="image_10.png" /><figcaption>Fig. 7. A periodic orbit of model (33).</figcaption></figure>

flow curvature manifold, i.e. the curvature when $\mu \to 0$. Then, we prove that the rate of change of the energy function of such coarse systems is equal to the scalar product of the vector field (velocity) and its first time derivative (acceleration) at first order in $\mu$.

The comparison between both tangential and normal components of the acceleration has enabled to provide an interpretation of the two main relationships established in this work. We thus proved for 2D & 3D coarse systems that when $\mu \to 0$, curvature, i.e., energy function expresses the magnitude of normal acceleration. We also proved that the rate of change of the energy function is proportional to the tangential acceleration and expresses the rate of change of the norm of the velocity vector field.

Thus, these two main results could be extended and generalized to “coarse systems” of higher dimensions than three. However, such a generalization would imply to define as an example a 4D coarse system by adding a fourth equation to the third original ones (2). But on one hand, the definition of this last equation will not be necessary the unique one, and on the other hand, the number of remaining terms in the curvature (6) will increase according to the dimension as already highlighted through the two examples of 2D & 3D coarse systems. Thus, the convergence will be more difficult to establish. However, by considering that the last equation can be written as $\dot{u} = \mu k(x, y, z, \mu)$ where $k$ is a function of class $C^4$, the generalization will then pose no difficulty.

In this work, the considered 2D & 3D coarse systems are not a generalized version of multi scale problems where the diffusion term or the high error derivatives are also having parameters which are arbitrarily small. Hence the present class of systems is having limited applications as they are not of singularly perturbed nature where the solution varies rapidly in a small region and smoothly outside of it. This induces a limitation for proving the energy function for higher dimensional systems. Singularly perturbed Differential equations are extremely difficult for numerical approximations as well as for the stability analysis for nonlocal cases including well posedness properties (see the works of Kumar et al. [22–26] or those of Sarkar et al. [27] and Saw et al. [28]).

### CRediT authorship contribution statement

**Jean-Marc Ginoux:** Writing – review & editing, Writing – original draft. **Riccardo Meucci:** Writing – review & editing, Writing – original draft. **Jaume Llibre:** Writing – review & editing, Writing – original draft. **Julien Clinton Sprott:** Writing – review & editing, Writing – original draft.

### Declaration of competing interest

The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.

### Acknowledgments

Jaume Llibre is partially supported by the Agencia Estatal de Investigación of Spain grant PID2022-136613NB-100, AGAUR (Generalitat de Catalunya) grant 2021SGR00113, and by the Reial Acadèmia de Ciències i Arts de Barcelona.

## Appendix A. Identity between energy functions and their rate of change

In this first appendix, we prove the identity between the energy function $H$ and its rate of change $dH/dt$ given by the classical method (11)–(12) and those given by both Helmholtz’s theorem (16)–(17) and Flow Curvature Method (19)–(20). For sake of simplicity and since the following proof holds for both 2D and 3D coarse systems (1) and (2), we omit to express the arguments of functions $f$ and $g$ and their time derivatives in what is presented below. Let us start from the rate of change of the energy function given either by the Helmholtz’s theorem (17) or by the Flow Curvature Method (20):

$$ \frac{dH}{dt} = \frac{d}{dt} \left( \frac{\dot{x}^2}{2} + \frac{\dot{y}^2}{2} \right) = \mu (\dot{x}f + \dot{y}g) \quad (A.1) $$

Then, according to (1), $\dot{y} = -x + \mu g$. Thus, we have:

$$ \frac{\dot{y}^2}{2} = \frac{x^2}{2} - \mu xg + \frac{\mu^2}{2} g^2 \quad (A.2) $$

By replacing this expression (A.2) in Eq. (A.1) leads to:

$$ \frac{d}{dt} \left( \frac{\dot{x}^2}{2} + \frac{x^2}{2} - \mu xg + \frac{\mu^2}{2} g^2 \right) = \mu (\dot{x}f + \dot{y}g) \quad (A.3) $$

This gives:

$$ \frac{d}{dt} \left( \frac{\dot{x}^2}{2} + \frac{x^2}{2} \right) - \mu \dot{x}g - \mu x \dot{g} + \mu^2 g \dot{g} = \mu (\dot{x}f + \dot{y}g) \quad (A.4) $$

After factorization, we have:

$$ \frac{d}{dt} \left( \frac{\dot{x}^2}{2} + \frac{x^2}{2} \right) - \mu \dot{x}g + \mu (-x + \mu g) \dot{g} = \mu (\dot{x}f + \dot{y}g) \quad (A.5) $$

Finally, we obtain:

$$ \frac{d}{dt} \left( \frac{\dot{x}^2}{2} + \frac{x^2}{2} \right) = \mu (f + g) \dot{x} \quad (A.6) $$

Then, the proof is stated.

## Appendix B. Additional term in the energy function $H$ of 2D coarse systems

According to Eq. (18), the curvature reads:

$$ \kappa = \dot{x}^2 + \dot{y}^2 + \mu (\dot{y}f - \dot{x}g) = 0. \quad (B.1) $$

Let us compute the expression of the third term in the right hand side. To this aim, let us replace $\dot{f}$ and $\dot{g}$ by:

$$ \begin{aligned} \dot{f} &= \frac{\partial f}{\partial x} \dot{x} + \frac{\partial f}{\partial y} \dot{y}, \\ \dot{g} &= \frac{\partial g}{\partial x} \dot{x} + \frac{\partial g}{\partial y} \dot{y}. \end{aligned} \quad (B.2) $$

Thus, we have:

$$ \mu (\dot{y}f - \dot{x}g) = \mu \left( \frac{\partial f}{\partial x} \dot{x} \dot{y} + \frac{\partial f}{\partial y} \dot{y}^2 - \frac{\partial g}{\partial x} \dot{x}^2 - \frac{\partial g}{\partial y} \dot{x} \dot{y} \right) \quad (B.3) $$
By replacing $\dot{x}$ and $\dot{y}$ by their expressions given by (1), we find:

$$ \mu (\dot{y}f - \dot{x}\dot{g}) = \mu \left( \frac{\partial f}{\partial x} x^2 - \left( \frac{\partial f}{\partial x} - \frac{\partial g}{\partial y} \right) xy - \frac{\partial g}{\partial y} y^2 \right) \\ - \mu^2 \left( 2xg \frac{\partial f}{\partial y} - \left( \frac{\partial f}{\partial x} - \frac{\partial g}{\partial y} \right) (yg - xf) + 2yf \frac{\partial g}{\partial x} \right) \\ + \mu^3 \left( g^2 \frac{\partial f}{\partial y} + \left( \frac{\partial f}{\partial x} - \frac{\partial g}{\partial y} \right) fg - f^2 \frac{\partial g}{\partial x} \right) \quad (B.4) $$

From Eq. (B.4), we observe that for small $\mu$ regime the whole expression tends to zero provided that $x, y$, the functions $f, g$ and their partial derivatives are bounded. This condition might be satisfied either for a periodic motion or for a chaotic attractor. This is the case for the Van der Pol system (33) the solution of which is a limit cycle and for the Li-Sprott model (38) which exhibits a bounded chaotic attractor. When $\mu$ is not considered “small” the identity between energy function and curvature is no more valid.

## Appendix C. Additional term in the energy function H of 3D coarse systems

According to Eq. (30), the curvature reads:

$$ \kappa(x, y, z, \mu) = \sqrt{[\dot{x}^2 + \dot{y}^2 + \mu(\dot{y}f(x, y, z, \mu) - \dot{x}\dot{g}(x, y, z, \mu))]^2 + O(\mu^2)} = 0. \quad (C.1) $$

Let us compute the expression of the third term in the right hand side. To this aim, let us compute the coordinates of the vector $\ddot{\vec{x}} \wedge \ddot{\vec{x}}$ appearing in the definition of the curvature (6). For sake of simplicity, we omit to express the arguments of functions $f$ and $g$ and their time derivatives in what is presented below. We find:

$$ \dot{\vec{x}} \wedge \ddot{\vec{x}} = \begin{pmatrix} -\dot{x}\dot{z} + \mu(\dot{z}\dot{g} - \dot{y}h) \\ -\dot{y}\dot{z} + \mu(\dot{x}h - \dot{z}f) \\ \dot{x}^2 + \dot{y}^2 + \mu(\dot{y}f - \dot{x}\dot{g}) \end{pmatrix} \quad (C.2) $$

Since the curvature (6) is defined by the Euclidean norm of the above vector (C.2), the additional term in the energy function (30) reads:

$$ (-\dot{x}\dot{z} + \mu(\dot{z}\dot{g} - \dot{y}h))^2 + (-\dot{y}\dot{z} + \mu(\dot{x}h - \dot{z}f))^2 \quad (C.3) $$

Here again, from Eq. (C.3), we observe that for small $\mu$ regime the whole expression of the curvature (6) tends to $\dot{x}^2 + \dot{y}^2$ provided that $x, y, z$, the functions $f, g, h$ and their partial derivatives are bounded. This condition might be satisfied either for a periodic motion or for a chaotic attractor. This is the case for the Li-Sprott model (38) which exhibits a limit cycle or a bounded chaotic attractor. When $\mu$ is not considered “small” the identity between energy function and curvature is no more valid.

## Data availability

No data was used for the research described in the article.

## References

[1] Andronov AA, Pontryagin LS. Systèmes grossiers (rough systems). Dokl Akad Nauk SSSR 1937;14(5):247–50, (in French).

[2] Andronov AA, Khaikin SE, Vitt AA. Theory of oscillators I. Moscow, Princeton, N. J: Princeton Univ. Press; 1937, Engl. transl. 1949.

[3] Llibre J, Zhang X. Hopf bifurcation in higher dimensional systems via the averaging method. Pacific J Math 2009;240(2):321–41.

[4] Llibre J, Makhlouf A, Badi S. 3-dimensional Hopf bifurcation via averaging theory of second order. Discret Contin Dyn Syst 2009;25:1287–95.

[5] Ginoux JM, Meucci R, Llibre J, Sprott JC. Energy function of 2D and 3D dynamical systems. Chaos Solitons Fractals 2025;190:115768, 6.

[6] Ginoux J-M, Lebiedz D, Meucci R, Llibre J, Sprott JC. Flow curvature manifold and energy of generalized Liénard systems. Chaos Solitons Fractals 2022;161:112354, 7.

[7] Coddington EA, Levinson N. Theory of ordinary differential equations. New York: Mac Graw Hill; 1955.

[8] Frénet F. Sur les courbes à double courbure, thèse toulouse. Résumé Dans J Math 1847;17, 1847.

[9] Ginoux JM, Rossetto B. Differential geometry and mechanics applications to chaotic dynamical systems. Int J Bifurc Chaos 2006;4(16):887–910.

[10] Ginoux JM, Rossetto B, Chua LO. Slow invariant manifolds as curvature of the flow of dynamical systems. Int J Bifurc Chaos 2008;11(18):3409–30.

[11] Ginoux JM. Differential geometry applied to dynamical systems. World scientific series on nonlinear science, series A, vol. 66, Singapore: World Scientific; 2009.

[12] Sarasola C, Torrealdea FJ, d'Anjou A, Moujahid A, Grana M. Energy balance in feedback synchronization of chaotic systems. Phys Rev E 2004;69:011606, 12.

[13] Wang CN, Wang Y, Ma Y. Calculation of hamilton energy function of dynamical system by using Helmholtz theorem. Acta Phys Sin 2016;65(24):240501, 6.

[14] Ma J, Wu F, Jin W, Zhou P, Hayat T. Calculation of hamilton energy and control of dynamical systems with different types of attractors. Chaos 2017;27:053108.

[15] Yu J, Njitacke ZT, Jiang D, Wu J, Awrejcewicz J. Energy function and complex dynamics from a jerk system. Phys Scr 2024;99:015245, 11.

[16] Van der Pol B. On 'relaxation-oscillations'. Phil Mag 1926;7(2):978–92.

[17] Diliberto SP. On systems of ordinary differential equations. In: Contributions to the theory of nonlinear oscillations. Ann. math. stud., vol. 20, 1950, p. 1–38.

[18] Ginoux J-M. Slow invariant manifolds of slow-fast dynamical systems. Int J Bifurc Chaos 2021;31(7):2150112, 17.

[19] Minorsky N. Théorie des oscillations. Méml. Sci Math. Fasc 1967;163.

[20] Li C, Sprott JC. Multistability in a butterfly flow. Int J Bifurc Chaos 2013;23(12):1350199, 10.

[21] Verhulst F. Nonlinear differential equations and dynamical systems. Universitex, Berlin: Springer-Verlag; 1990, 1990.

[22] Kumar S, Ishwariya R, Das P. Impact of mixed boundary conditions and non-smooth data on layer originated non-premixed combustion problems: Higher order convergence analysis. Stud Appl Math 2024;153(4):e12763.

[23] Kumar S, Kumar S, Das P. Second order a priori and a posteriori error estimations for integral boundary value problems of nonlinear singularly perturbed parameterized form. Numer Algorithms 2024. https://link.springer.com/article/10.1007/s11075-024-01918-5.

[24] Kumar S, Das P, Kumar K. Adaptive mesh based efficient approximations for Darcy scale precipitation-dissolution models in porous media. Int J Numer Methods Fluids (1D & 2D Problems) 2024;96(8):1415–44.

[25] Kumar S, Das P. A uniformly convergent analysis for multiple scale parabolic singularly perturbed convection-diffusion coupled systems: Optimal accuracy with less computational time. Appl Numer Math 2025;207:534–57.

[26] Kumar S, Das P. Analysis of an efficient parameter uniform domain decomposition approach for singularly perturbed gierer-meinhardt type nonlinear coupled systems of parabolic problems. Int J Numer Anal Model 2025;22(4):459–82.

[27] Sarkar D, Kumar S, Das P, Ramos H. Higher-order convergence analysis for interior and boundary layers originated semi-linear reaction diffusion system networked by a k-star graph with non-smooth source terms. Netw Heterog Media 2024;19(3):1085–115.

[28] Saw V, Das P, Srivastava HM. Investigations of a class of Liouville-Caputo fractional order Pennes bioheat flow partial differential equations through orthogonal polynomials on collocation points. Bull Sci Math. 2025;201:103637, 2025.