A Coordinate-Free Variational Approach to Fourth-Order Dynamical Systems on Manifolds: A System and Control Theoretic Viewpoint

Abstract

The present paper describes, in a theoretical fashion, a variational approach to formulate fourth-order dynamical systems on differentiable manifolds on the basis of the Hamilton–d’Alembert principle of analytic mechanics. The discussed approach relies on the introduction of a Lagrangian function that depends on the kinetic energy and the covariant acceleration energy, as well as a potential energy function that accounts for conservative forces. In addition, the present paper introduces the notion of Rayleigh differential form to account for non-conservative forces. The corresponding fourth-order equation of motion is derived, and an interpretation of the obtained terms is provided from a system and control theoretic viewpoint. A specific form of the Rayleigh differential form is introduced, which yields non-conservative forcing terms assimilable to linear friction and jerk-type friction. The general theoretical discussion is complemented by a brief excursus about the numerical simulation of the introduced differential model.

1. Introduction

Newton’s laws of motion for a rigid body give rise to dynamical systems that are perhaps the best known at any level. Such mechanical dynamical systems are described by second-order differential equations in variables that carry information about the positional and attitudinal (i.e., orientational) state of a body [1]. The first-order derivative of position and attitude are linear velocity and rotational velocity, respectively, while the second-order derivative of position and attitude are linear acceleration and angular acceleration, respectively.

As a matter of fact, a good deal of mathematical models encountered in applied sciences and engineering are either first-order [2,3] or second-order [4,5]. Often, such mathematical models describe non-linear oscillatory phenomena associated with a chaotic behavior. The category of higher-order differential models is likewise rich. An instance of higher-order differential model is the mathematical model of vibration of a cantilever beam, which is fourth-order [6].

In linear and non-linear electrical circuit theory, there exists no theoretical limitation to the differential order of mathematical models since building a higher-order electrical circuit is a matter of adding memory components (such as inductors, capacitors and memristors) and even higher-order electrical components [7]. Notable is the quality of research work conducted toward the characterization and application of third-order circuits (often termed ‘jerk’ circuits) [8,9,10]. Jerk-type systems are not confined to the field of electrical circuit theory, as they are encountered in mechanics [11,12] and biomechanics [13], to name a few.

Classical research in differential dynamical systems concerns unconstrained systems, such as the ones mentioned above, the dynamics of which occur on flat spaces (typically ${\mathbb{R}}^n$, where the integer n denotes the number of descriptive variables, most often $n=2$ or $n=3$). Unconstrained systems may be studied and characterized using classical analysis and geometry tools. In contrast, whenever the state of a dynamical system is subjected to holonomic constraints, manifold calculus plays a central role in its formulation and analysis [14,15], as well as its design and control [16,17,18]. In particular, the formulation of general dynamical systems on manifolds requires a general principle (and a general formalism) as an anchor point.

In the present paper, in continuity with a well-established practice in mechanics and physics, the Lagrangian formulation, based on the Hamiltonian stationary-action principle, extended by the d’Alembertian virtual work principle, is chosen as a general framework (for an introduction, see [19], Chapter 13). The classical Hamiltonian stationary-action principle is based on the definition of a Lagrangian function as a combination of two fundamental forms of energy, namely, kinetic energy and potential energy. Classical Lagrangians of this form describe so-termed ‘simple mechanical systems’ [20]. The kinetic energy is a quadratic function of the system’s state velocity and describes the internal structure of the system, while the potential energy subsumes all sorts of conservative forces acting on the system. The d’Alembertian extension, based on the principle of virtual work, allows taking into account all sorts of non-conservative forces acting on a dynamical system, such as frictions and thrusts produced by actuators.

In classical mechanics, a Lagrangian is a function of the system’s state and of the state’s velocity only. Generalizations to the classical setting have been developed through the decades. The authors of [21] studied Lagrangians linear in velocity, the authors of [22] explain details about Lagrangians depending on the acceleration, the authors of [23] describe non-conventional discontinuous Lagrangian functions for viscous flows, and the authors of [24,25,26] describe a formal procedure to make the main properties of the classical Lagrangian function extendable to functions of any order of the time-derivatives of the state variable.

In the present paper, a Lagrangian function is introduced, which depends on kinetic energy (quadratic in the tangent velocity, a metric-related term) and acceleration energy (quadratic in the covariant acceleration, a Palais-metric term [27]), as well as a potential energy function. In addition, the present paper introduces the notion of Rayleigh differential form to account for non-conservative forces that might act on a system (including active/passive damping). The corresponding Euler–Lagrange–d’Alembert equation of motion is derived and is shown to be of fourth order. A specific form of the Rayleigh differential form is introduced, which yields linear and third-order friction terms. The paper ends with a short discussion about a simple numerical scheme to numerically integrate the differential equation describing a fourth-order dynamical system on manifold and with a few concluding remarks.

A distinguishing feature of the present formulation is that it is derived and expressed in a coordinate-free fashion, which makes it particularly well suited to perform numerical simulations whenever the base manifold is a large-dimension matrix manifold or matrix Lie group. Such an approach is particularly viable when regarded from an applied system and control viewpoint.

2. Notation and Fundamental Properties

The aim of this section is to recall relevant notations and properties from manifold calculus. For an introduction specifically dedicated to system theory and non-linear control, interested readers might want to consult the open-access-published tutorial papers [28,29] appearing in the journal Symmetry.

A smooth state space exhibiting the structure of a differentiable manifold is denoted as $\mathbb{M}$. Given a point $x\in \mathbb{M}$, the associated tangent space is denoted as $T_x\mathbb{M}$. The tangent bundle associated with the manifold $\mathbb{M}$ is denoted as $T\mathbb{M}:=\{(x,v)\in \mathbb{M}\times T_x\mathbb{M}\}$. A smooth vector field $v:\mathbb{M}\to T\mathbb{M}$ associates with each point $x\in \mathbb{M}$ a tangent vector $v_x\in T_x\mathbb{M}$. The set of smooth vector fields of the manifold $\mathbb{M}$ is denoted as $\Gamma (T\mathbb{M})$.

Manifolds of interest in applied system theory and control are, for instance, the hypersphere ${\mathbb{S}}^{p-1}$; the general linear group of matrices $\mathbb{GL}(p)$; the special linear group of matrices $\mathrm{Sl}(p)$; the manifold of special orthogonal matrices $\mathbb{SO}(p)$; the compact Stiefel manifold $\mathrm{St}(n,p)$, with $p⩽n$, the real symplectic group $\mathrm{Sp}(2n)$; and the space of symmetric, positive-definite matrices ${\mathbb{S}}^{+}(p)$, where $p,n$ denote positive integers.

Given a smooth function $f:\mathbb{M}\to \mathbb{R}$, its Riemannian gradient at a point $x\in \mathbb{M}$ is denoted as ${\mathrm{grad}}_{x}f\in T_x\mathbb{M}$. Given a smooth path $\rho :[0,1]\to \mathbb{M}$, it holds that

$$\frac{\mathrm{d}}{\mathrm{d}t}(f\circ \rho )={\langle {\mathrm{grad}}_{\rho }f,\dot{\rho }\rangle }_{\rho },$$

(1)

which may be recognized as a chain rule of derivation, where $\dot{\rho }:=\frac{\mathrm{d}\rho }{\mathrm{d}t}$ denotes the tangent velocity field associated with the path at point $\rho (t)\in \mathbb{M}$.

The manifold $\mathbb{M}$ is taken as a Riemannian manifold endowed with a smooth family of inner products ${\langle \star ,\star \rangle }_{\star }:{(T\mathbb{M})}^2\to \mathbb{R}$ and the associated Levi-Civita connection ∇. It pays to recall that a Levi-Civita connection is metric, in the sense that, given a smooth path $\rho :[0,1]\to \mathbb{M}$ and two vector fields $v,w\in \Gamma (T\mathbb{M})$ defined at least on the path $\rho$, it holds that

$$\frac{\mathrm{d}}{\mathrm{d}t}{\langle v,w\rangle }_{\rho }={\langle {\nabla }_{\dot{\rho }}v,w\rangle }_{\rho }+{\langle v,{\nabla }_{\dot{\rho }}w\rangle }_{\rho }.$$

(2)

Furthermore, it is worth recalling that the Levi-Civita connection is torsion-free; hence, the covariant derivative satisfies

$${\nabla }_{v}w={\nabla }_{w}v+[v,w],$$

(3)

where $v,w\in \Gamma (T\mathbb{M})$ and $[\star ,\star ]:{(\Gamma (T\mathbb{M}))}^2\to \Gamma (T\mathbb{M})$ denotes the commutator of two vector fields.

A useful relation that appears as a homolog of integration by parts is the path-integration by parts formula.

Lemma 1.

Given a path $\rho :[0,1]\to \mathbb{M}$ and two smooth vector fields $v,w\in \Gamma (T\mathbb{M})$ defined at least on the path ρ, it holds that

$${\int }_0^1{\langle w,{\nabla }_{\dot{\rho }}v\rangle }_{\rho }\mathrm{d}t={\langle v,w\rangle }_{\rho (1)}-{\langle v,w\rangle }_{\rho (0)}-{\int }_0^1{\langle {\nabla }_{\dot{\rho }}w,v\rangle }_{\rho }\mathrm{d}t.$$

(4)

Proof. 

The proof follows directly from the property (2). □

It is also instrumental to recall the fundamental relationship between the double covariant derivative of a vector field and the curvature endomorphism. Given three smooth vector fields $u,v,w\in \Gamma (T\mathbb{M})$, it holds that

$${\nabla }_u{\nabla }_{v}w={\nabla }_v{\nabla }_{u}w+{\nabla }_{[u,v]}w+{\mathbf{R}}^{u,v}(w),$$

(5)

where $\mathbf{R}:{(\Gamma (T\mathbb{M}))}^3\to \Gamma (T\mathbb{M})$ denotes the curvature endomorphism associated with the connection ∇. The curvature endomorphism exhibits several symmetry properties, two of which are

$$\begin{array}{c}{\mathbf{R}}^{u,v}=-{\mathbf{R}}^{v,u}, (\mathrm{skew} \mathrm{symmetry})\hfill \end{array}$$

(6)

$$\begin{array}{c}\langle {\mathbf{R}}^{u,v}(w),z\rangle =\langle {\mathbf{R}}^{w,z}(u),v\rangle , (\mathrm{interchange} \mathrm{symmetry})\hfill \end{array}$$

(7)

for every quadruple of vector fields $u,v,w,z\in \Gamma (T\mathbb{M})$.

A curve $\rho :[0,1]\to \mathbb{M}$ that satisfies the zero covariant acceleration property ${\nabla }_{\dot{\rho }}\dot{\rho }=0$ is a geodesic. A solution to the geodesic equation corresponding to the initial conditions $\rho (0)=x\in \mathbb{M}$ and $\dot{\rho }(0)=v\in T_x\mathbb{M}$ is denoted as ${\rho }_{x,v}(t)$.

A canonical map associated with a geodesic curve is the exponential map $exp:T\mathbb{M}\to \mathbb{M}$, defined by ${exp}_x(v):={\rho }_{x,v}(1)$. A further canonical map between tangent spaces is parallel transport $\mathbb{P}:T\mathbb{M}\times \mathbb{M}\to T\mathbb{M}$. Parallel transport moves a tangent vector $v\in T_x\mathbb{M}$ from tangent space $T_x\mathbb{M}$ to tangent space $T_y\mathbb{M}$ along the geodesic arc connecting point x to point y. The resulting vector is denoted as ${\mathbb{P}}^{x\to y}(v)$.

The covariant derivative along a curve $\rho$ is related to parallel transport by the formula

$${({\nabla }_{\dot{\rho }}w)}_{\rho (t)}=\underset{h\to 0}{lim}\frac{{\mathbb{P}}^{\rho (t+h)\to \rho (t)}{w}_{\rho (t+h)}-w_{\rho (t)}}{h},$$

(8)

where the vector field $w\in \Gamma (T\mathbb{M})$ is defined and smooth at least over the path $\rho$. Such a relation proves relevant when dealing with computer-based simulations of dynamical systems on a manifold.

3. Lagrangian Formulation of Fourth-Order Dynamical Systems on Manifold

The formulation of a fourth-order dynamical system on a smooth manifold is based on the Hamiltonian principle of stationary action, modified through the d’Alembertian principle of virtual work to accommodate for non-conservative generalized forces. According to such principles, the state of a system whose dynamics is governed by a Lagrangian function and is subjected to an exogenous forcing term follows the unique trajectory on the state manifold corresponding to a stationary action/virtual work.

The motion of the dynamical system is a continuous set of pairs of the form $(t,x)\in [0,1]\times \mathbb{M}$, where the variable t denotes absolute time and ranges conventionally in the interval $[0,1]$ (since reparameterization is always possible), and the variable x denotes the state of the system belonging to the Riemannian state manifold $\mathbb{M}$. The function $t\mapsto x(t)$ is termed the trajectory of the system.

In order to formulate the Hamilton–d’Alembert variational principle, it pays to define a smooth family of trajectories $c:[-ϵ,ϵ]\times [0,1]\to \mathbb{M}$, where $ϵ>0$, among which the fiducial path is sought. A member $c(s,\star )$ of the family of curves is individuated by the index s, while a specific point on a specific curve is denoted as $c(s,t)$ through the affine parameter t. The family member corresponding to the fiducial trajectory has index $s=0$; namely, it is denoted as $c(0,t)\equiv x(t)$. In systems and control, the endpoints of a sought trajectory are fixed; therefore, it is assumed that $c(\star ,0)=x(0)$ is known and $c(\star ,1)=x(1)$ is known, meaning that all curves originate from the same point and terminate at the same point.

To ease notation, the velocity field associated with the fiducial trajectory is denoted as $v:={\left.\frac{\partial c}{\partial t}\right|}_{s=0}$, while the variation from the fiducial trajectory is denoted through the ‘transverse’ vector field $\xi :={\left.\frac{\partial c}{\partial s}\right|}_{s=0}$. The vector field $v\in T_x\mathbb{M}$ and the transverse variation $\xi \in T_x\mathbb{M}$ are tangent vector fields. Using the smoothness of the family of curves c, it follows that $\frac{{\partial }^{2}c}{\partial s\partial t}=\frac{{\partial }^{2}c}{\partial t\partial s}$; hence, the vector fields v and $\xi$ commute, namely, $[v,\xi ]=0$.

The dynamics of the system being constructed are governed by an (extended) Lagrangian function

$${\mathcal{L}}_x(v):=\frac{\sigma }{2}{\langle {\nabla }_{v}v,{\nabla }_{v}v\rangle }_x+\frac{\phi }{2}{\langle v,v\rangle }_x-{\mathcal{V}}_x,$$

(9)

where $\sigma ,\phi \in \mathbb{R}$ are arbitrary constants, ${({\nabla }_{v}v)}_x$ denotes the covariant acceleration of the trajectory at point x and $\mathcal{V}:\mathbb{M}\to \mathbb{R}$ denotes a potential energy function (only depending on the point over the trajectory) that subsumes all conservative forces acting on the system. The term $\langle {\nabla }_{v}v,{\nabla }_{v}v\rangle$ in the Lagrangian may be seen as a generalization of the Appell–Gibbs ‘acceleration energy’ [30].

The Hamilton–d’Alembert principle states that a dynamical system will follow a trajectory that makes a certain functional stationary, with such a functional being a sum of the total action of the system and the total virtual work of the exogenous forces. In the present context, such a principle will be expressed as

$$0={\left\{\left.\frac{\mathrm{d}}{\mathrm{d}s}{\int }_0^1{\mathcal{L}}_c\left(\frac{\partial c}{\partial t}\right)\mathrm{d}t+{\int }_0^1{\mathcal{R}}_c\left(\frac{\partial c}{\partial s}\right)\mathrm{d}t\right\}\right|}_{s=0},$$

(10)

where $\mathcal{R}\in \Gamma (T^{\star }\mathbb{M})$ denotes a differential form (or covector field) that subsumes all non-conservative phenomena and may depend on the trajectory (here, the symbol $T^{\star }\mathbb{M}$ denotes the dual tangent bundle associated with the manifold $\mathbb{M}$). As a covector, it acts like $\mathcal{R}:T\mathbb{M}\to \mathbb{R}$. In systems and control, such quantity accounts for frictions and typically forces/torques exerted by actuators. The form $\mathcal{R}$ in the variational formulation (10) may be regarded as a generalization of the Rayleigh energy dissipation function [30]. The integral $\mathcal{A}:=\int \mathcal{L}\mathrm{d}t$ denotes the total action of the system, while the integral $\mathcal{W}:=\int \mathcal{R}(\xi )\mathrm{d}t$ denotes the virtual work of the exogenous forces along a variation.

Canonically to the variational formulation, the variation $\xi$ from the fiducial trajectory is assumed to be arbitrary, except that, at the endpoints of the fiducial trajectory, it must hold that ${\xi }_{x(0)}={\xi }_{x(1)}=0$ and ${({\nabla }_v\xi )}_{x(0)}={({\nabla }_v\xi )}_{x(1)}=0$.

Theorem 1.

The general fourth-order differential equation of dynamics arising from the variational principle (10) reads

$$\sigma {({\nabla }_v{\nabla }_v{\nabla }_{v}v)}_x-\phi {({\nabla }_{v}v)}_x+\sigma {\mathbf{R}}_x^{{({\nabla }_{v}v)}_x,v_x}(v_x)={\mathrm{grad}}_x\mathcal{V}-f_x.$$

(11)

Proof. 

In order to derive the differential equation governing the time evolution of the state of a dynamical system with Lagrangian $\mathcal{L}$, it is necessary to evaluate the variation in the three addenda in (9), which make up the action $\mathcal{A}$. The term $\mathcal{W}$ accounting for non-conservative phenomena will be detailed in Section 4.

The term including the potential function may be written as

$${\left\{\left.\frac{\mathrm{d}}{\mathrm{d}s}{\int }_0^1\mathcal{V}\mathrm{d}t\right\}\right|}_{s=0}={\left\{\left.{\int }_0^1\frac{\partial {\mathcal{V}}_c}{\partial s}\mathrm{d}t\right\}\right|}_{s=0}={\int }_0^1{\langle {\mathrm{grad}}_x\mathcal{V},\xi \rangle }_x\mathrm{d}t.$$

(12)

The term including the velocity field along the family of curves may be written explicitly as

$${\left\{\left.\frac{\mathrm{d}}{\mathrm{d}s}{\int }_0^1\frac{1}{2}{〈\frac{\partial c}{\partial t},\frac{\partial c}{\partial t}〉}_c\mathrm{d}t\right\}\right|}_{s=0}={\int }_0^1{\langle {\nabla }_{\xi }v,v\rangle }_x\mathrm{d}t.$$

(13)

From Property (3) and the recalled property $[v,\xi ]=0$, it follows that

$${\left\{\left.\frac{\mathrm{d}}{\mathrm{d}s}{\int }_0^1\frac{1}{2}{〈\frac{\partial c}{\partial t},\frac{\partial c}{\partial t}〉}_c\mathrm{d}t\right\}\right|}_{s=0}={\int }_0^1{\langle {\nabla }_v\xi ,v\rangle }_x\mathrm{d}t.$$

(14)

Thanks to the path integration via parts formula (4), the above expression may be recast as

$${\left\{\left.\frac{\mathrm{d}}{\mathrm{d}s}{\int }_0^1\frac{1}{2}{〈\frac{\partial c}{\partial t},\frac{\partial c}{\partial t}〉}_c\mathrm{d}t\right\}\right|}_{s=0}=-{\int }_0^1{\langle {\nabla }_{v}v,\xi \rangle }_x\mathrm{d}t$$

(15)

since, by assumption, ${\langle v_{x(0)},{\xi }_{x(0)}\rangle }_{x(0)}={\langle v_{x(1)},{\xi }_{x(1)}\rangle }_{x(1)}=0$.

The last term to be taken into account is the one including covariant acceleration that, upon taking variation as in the previous two terms, may be written explicitly as

$${\int }_0^1{\langle {\nabla }_{\xi }{\nabla }_{v}v,{\nabla }_{v}v\rangle }_x\mathrm{d}t.$$

(16)

By recalling Expression (5) and invoking again the fact that $[v,\xi ]=0$, it can be immediately concluded that

$${\int }_0^1{\langle {\nabla }_{\xi }{\nabla }_{v}v,{\nabla }_{v}v\rangle }_x\mathrm{d}t={\int }_0^1{\langle {\nabla }_v{\nabla }_{\xi }v+{\mathbf{R}}^{\xi ,v}(v),{\nabla }_{v}v\rangle }_x\mathrm{d}t.$$

(17)

From the right-hand side of the above relation, there emerge two addenda that are worth analyzing separately. The rightmost addendum may be rewritten as

$$\begin{array}{c}\hfill {\int }_0^1{\langle {\mathbf{R}}^{\xi ,v}(v),{\nabla }_{v}v\rangle }_x\mathrm{d}t=-{\int }_0^1{\langle {\mathbf{R}}^{v,\xi }(v),{\nabla }_{v}v\rangle }_x\mathrm{d}t=\\ \hfill -{\int }_0^1{\langle {\mathbf{R}}^{v,{\nabla }_{v}v}(v),\xi \rangle }_x\mathrm{d}t={\int }_0^1{\langle {\mathbf{R}}^{{\nabla }_{v}v,v}(v),\xi \rangle }_x\mathrm{d}t,\end{array}$$

(18)

by applying the skew-symmetry and interchange-symmetry properties of the curvature endomorphism. Likewise, the leftmost addendum in the relation (18) may be rewritten as

$$\begin{array}{c}\hfill {\int }_0^1{\langle {\nabla }_{v}v,{\nabla }_v{\nabla }_{\xi }v\rangle }_x\mathrm{d}t={\int }_0^1{\langle {\nabla }_{v}v,{\nabla }_v{\nabla }_v\xi \rangle }_x\mathrm{d}t=\\ \hfill -{\int }_0^1{\langle {\nabla }_v{\nabla }_{v}v,{\nabla }_v\xi \rangle }_x\mathrm{d}t={\int }_0^1{\langle {\nabla }_v{\nabla }_v{\nabla }_{v}v,\xi \rangle }_x\mathrm{d}t,\end{array}$$

(19)

by applying twice the path integration by parts formula and invoking the non-commutativity of the vector–field pair $(v,\xi )$ and the boundary conditions on the variation $\xi$.

By definition, the differential form $\mathcal{R}$ may be written as $\mathcal{R}(\star )=\langle f,\star \rangle$, where $f_x\in T_x\mathbb{M}$ denotes a forcing term, in such a way that ${\mathcal{R}}_x(\xi )={\langle f_x,\xi \rangle }_x$. Collecting the above partial expressions leads to the relation:

$${\int }_0^1{\langle \sigma \{{\nabla }_v{\nabla }_v{\nabla }_{v}v+{\mathbf{R}}^{{\nabla }_{v}v,v}(v)\}-\phi {\nabla }_{v}v+f-\mathrm{grad}\mathcal{V},\xi \rangle }_x\mathrm{d}t=0.$$

(20)

Since the variation $\xi$ is arbitrary in the open interval $]0,1[$, the integrand must vanish to zero, which leads to the differential system (11), and needs four boundary conditions to be solved. A common set of boundary conditions is $x(0)=x_0$, $x(1)=x_1$, $\dot{x}(0)=v_0$ and $\dot{x}(1)=v_1$. □

In the differential Equation (11) on a tangent bundle, the third term on the left-hand side is perceived as a virtual force owing to the curved nature of the state manifold, while the second term on the right-hand side denotes the resultant of all non-conservative forces acting on the system.

Fourth-order (or jounce) dynamical systems on manifolds and their control were studied in detail in a series of contributions, among which we can cite [31,32,33]. We also mention that the special case, where $\sigma =1$, $\phi =0$, $\mathrm{grad}\mathcal{V}=0$ and $f=0$, is related to the theory of ‘cubic polynomials’ (for a review and connections with non-linear control, readers might consult [34]). The corresponding differential equation reads

$${\nabla }_v{\nabla }_v{\nabla }_{v}v=-{\mathbf{R}}^{{\nabla }_{v}v,v}(v),$$

(21)

the solutions of which are related to generalizations of geodesic arcs termed ‘cubic polynomials on manifolds’, defined as the solutions of third-order differential equations on manifolds of the form ${\nabla }_v{\nabla }_v{\nabla }_{v}v=0$. (The analysis of the differences between curvature-compliant cubic polynomials arising from Equation (21) and the standard cubic polynomials would be, per se, an endeavor of interest.)

4. Non-Conservative Forcing Terms and Reduced-Order Systems

The non-conservative forcing term f may arise as the resultant of several contributions that range from friction-type damping terms to active thrusts owing to propellers in actuated mechanical systems. Friction is often associated with energy dissipation, while the active forcing terms coincide essentially with driving and stabilizing forces/torques as computed by a control algorithm and exerted by an actuator of sorts. Nonetheless, friction may be active, hence injecting energy into a system. It is interesting to discuss, in detail, some kinds of friction-type terms.

The general principle (10) fits well with the amply documented notion of Rayleigh potential [35,36], which, in the present setting, is replaced by the Rayleigh differential form $\mathcal{R}$. Such differential form gives rise to the force

$$f:={\mathcal{R}}^{♯},$$

(22)

where superscript${ }^{♯}:T^{\star }\mathbb{M}\to T\mathbb{M}$ denotes the sharp isomorphism. (It is perhaps worth recalling that, often, forcing terms are identified with covectors rather than being identified with vectors, which, in fact, may be thought of as ‘coordinates’ of covectors. This is tantamount to saying that ${\mathcal{R}}_x(\xi )=\langle \langle {\mathcal{R},\xi \rangle \rangle }_x$, where the symbol $\langle \langle \star ,\star \rangle \rangle$ denotes a so-called ‘pairing’. In order to comply with the above notation, it suffices to observe that $\langle \langle {\mathcal{R},\xi \rangle \rangle }_x\equiv {\langle {\mathcal{R}}^{♯},\xi \rangle }_x$.)

An instance of Rayleigh form is $\mathcal{R}(\star ):=-\mu {\langle v,v\rangle }^{\frac{\epsilon }{2}}\langle v,\star \rangle$, giving rise to the friction term $-\mu \parallel v_x{\parallel }_x^{\epsilon }{v}_x$, which represents standard viscous damping, provided the constant $\mu$ takes a positive value [37]. The constant exponential $\epsilon \in \mathbb{R}$ denotes the degree of damping introduced in the system. The value $\epsilon =0$ corresponds to linear damping, namely,

$$F_x:=-\mu v_x.$$

(23)

A further instance of damping term is jerk-type friction. In a flat space, such a term has been mentioned, e.g., in [38]. In the present setting, we associate with such kind of friction the differential form $-\tau \langle {\nabla }_v{\nabla }_{v}v,\star \rangle$, where $\tau \in \mathbb{R}$ is constant, to which corresponds jerk-type friction

$$J_x:=-\tau {({\nabla }_v{\nabla }_{v}v)}_x.$$

(24)

Such a friction term may quantify the attitude of a system to avoid ‘shaken’ trajectories or any sort of active/passive force proportional to jerk. An instance of such kinds of forces is the Abraham–Lorentz friction term that represents the recoil action on an accelerating charged particle caused by the particle emitting electromagnetic radiation [39]. (Such force is proportional to the square of a moving particle’s charge times its jerk. In a cyclotron, where the jerk points opposite to the velocity, Abraham–Lorentz force is directed opposite to the velocity of the particle, hence exerting a braking action.)

We henceforth defined the Rayleigh differential form as $\mathcal{R}:=\langle -\mu v-\tau {\nabla }_v{\nabla }_{v}v+\tilde{f},\star \rangle$, where the vector field $\tilde{f}\in \Gamma (T\mathbb{M})$ represents the resultant of any additional forcing term not explicitly accounted for by the other terms. Having established such an expression for the non-conservative resultant force, the general equation of motion (11) becomes

$$\sigma {({\nabla }_v{\nabla }_v{\nabla }_{v}v)}_x-\tau {({\nabla }_v{\nabla }_{v}v)}_x-\phi {({\nabla }_{v}v)}_x+\sigma {\mathbf{R}}_x^{{({\nabla }_{v}v)}_x,v_x}(v_x)-\mu v_x={\mathrm{grad}}_x\mathcal{V}-{\tilde{f}}_x,$$

(25)

The derivative ${\nabla }_v{\nabla }_{v}v$ is termed covariant jerk, while the derivative ${\nabla }_v{\nabla }_v{\nabla }_{v}v$ is termed covariant jounce [29].

On the basis of the definition of the curvature-compliant cubic polynomials recalled in Section 3, it can immediately be seen that all terms in Equation (25), besides those with coefficient $\sigma$, tend to make the solution path deviate from a pure cubic polynomial. In addition, it is interesting to note that if the coefficient $\sigma =0$ while the coefficient $\tau \ne 0$, the above mathematical model collapses to a third-order jerk-type system, namely,

$${\nabla }_v{\nabla }_{v}v=-\frac{\phi }{\tau }{\nabla }_{v}v-\frac{\mu }{\tau }v-\frac{1}{\tau }(\mathrm{grad}\mathcal{V}-\tilde{f}).$$

(26)

Here, it is perhaps the right place to highlight the fact that the introduction of a third-order friction term is the only way to obtain a third-order system since a purely conservative Lagrangian approach would yield exclusively second- and fourth-order terms.

If one takes both $\sigma =0$ and $\tau =0$, while the coefficient $\phi \ne 0$, the system (25) reduces to the second-order dynamical model

$${\nabla }_{v}v=-\frac{\mu }{\phi }v-\frac{1}{\phi }(\mathrm{grad}\mathcal{V}-\tilde{f}),$$

(27)

studied in detail in [29], while if $\phi =0$, but $\mu \ne 0$, the system (25) further reduces to the first-order system

$$\dot{x}=\frac{1}{\mu }({\tilde{f}}_x-{\mathrm{grad}}_x\mathcal{V}),$$

(28)

which, in fact, is non-dynamical since it prescribes the velocity of a motion rather than its acceleration, jerk or jounce.

The first-order system (28) is classical in gradient-based optimization on Euclidean spaces, where the potential function $\mathcal{V}$ represents a criterion function, and its extremal points (either local or global) are sought after. The second-order system (27) represents an improvement to the first-order optimization method (28) and traces back to Polyak’s heavy-ball method [40], which, however, came as a discrete-time iterative algorithm. The heavy-ball method (later known as Nesterov accelerated gradient optimization method) may indeed be studied through an associated continuous-time differential equation, as shown in, e.g., [41,42,43]. The dynamical system (27) would constitute a natural generalization of the second-order differential equation associated with Nesterov’s iteration, under the proviso that one allows the friction coefficient $\mu$ to decrease in time as fast as $\frac{1}{t}$. A main finding in [41] is that the shape of the friction coefficient must indeed be exactly $\mu =\frac{3}{t}$ to ensure the fastest convergence, at least in the special Euclidean case $\mathbb{M}={\mathbb{R}}^n$.

A variety of alternative formulations may be found in the specialized literature. We mention that the authors of [44] studied forcing terms arising from non-holonomic constraints on the trajectory of a fourth-order dynamical system on manifold. It is also worth mentioning the alternative method to account for some sorts of non-conservative forces that may be found in the work of Herglotz (see, for example, [45]) based on action-dependent Lagrangian functions.

5. An Invariant of Motion for Fourth-Order Conservative Systems

The conservative version of the general fourth-order dynamical system (11) admits an invariant, namely, a function $\mathcal{H}$ such that $\dot{\mathcal{H}}=0$ whenever the resultant of the non-conservative forcing terms is zero. Such invariant takes the form of an energy

$$\mathcal{H}:=\frac{\sigma }{2}\left\{\langle {\nabla }_{v}v,{\nabla }_{v}v\rangle -2\langle {\nabla }_v{\nabla }_{v}v,v\rangle \right\}+\frac{\phi }{2}\langle v,v\rangle +\mathcal{V}.$$

(29)

Such an invariant depends on kinetic energy, acceleration energy and potential, as well as jerk-type force. In fact, it holds that $\mathcal{H}=\mathcal{L}-2\mathcal{V}-\sigma \langle {\nabla }_v{\nabla }_{v}v,v\rangle$.

Theorem 2.

Expression (29) is an invariant of the system (11) in the absence of non-conservative forces.

Proof. 

In order to show that the energy function $\mathcal{H}$ is indeed an invariant of System (11) in the absence of non-conservative forces, it pays to apply the rule of derivation (2), which leads to

$$\dot{\mathcal{H}}=\langle \phi {\nabla }_{v}v+\mathrm{grad}\mathcal{V}-\sigma {\nabla }_v{\nabla }_v{\nabla }_{v}v,v\rangle .$$

(30)

From the relation (11), it follows that $\sigma {\nabla }_v{\nabla }_v{\nabla }_{v}v-\phi {\nabla }_{v}v-\mathrm{grad}\mathcal{V}=-f-\sigma {\mathbf{R}}^{{\nabla }_{v}v,v}(v)$; therefore,

$$\dot{\mathcal{H}}=\langle \sigma {\mathbf{R}}^{{\nabla }_{v}v,v}(v)+f,v\rangle .$$

(31)

From the properties of the Riemannian curvature endomorphism (interchange symmetry and skew symmetry), it follows that ${\mathbf{R}}^{{\nabla }_{v}v,v}(v)\perp v$, irrespective of the chosen metric, whereby

$$\dot{\mathcal{H}}=\langle f,v\rangle .$$

(32)

Such a simple expression is essentially due to the fact that the forcing term arising from curvature ${\mathbf{R}}^{{\nabla }_{v}v,v}(v)$ is immaterial. The right-hand side of the expression (32) represents, in fact, the power (or energy rate) absorbed (if positive) or emitted (if negative) by the system.

The net consequence intrinsic to the relationship (32) is that, for a conservative system in which $f\equiv 0$, it holds that the function $\mathcal{H}$ stays constant to its initial value along any trajectory of the dynamical system (11). □

It is interesting to notice that, from the above discussion, it emerges that the term due to curvature does not produce work. In real-world settings, however, it was shown recently (both theoretically and experimentally) that curvature may indeed be responsible for motion. In particular, the authors of [46] showed that a precision apparatus confined to a spherical surface without a solid substrate is able to self-propel without environmental momentum exchange thanks to shape changes. This phenomenon is due to the lack of commutativity of state transitions on curved spaces.

For a non-conservative system of the type (25), where we take the unspecified term to be $\tilde{f}\equiv 0$, it rather holds that

$$\dot{\mathcal{H}}=-\mu \langle v,v\rangle -\tau \langle {\nabla }_v{\nabla }_{v}v,v\rangle .$$

(33)

The first term on the right-hand side is certainly non-positive and makes the energy $\mathcal{H}$ decrease over time, while the second term in the sum may take any value, depending on the orientation of the covariant jerk with respect to the velocity (For example, in an Abraham–Lorenz setting, if $\tau >0$, then $\langle J,v\rangle ⩽0$, and the Abraham–Lorenz jerk-type force J describes a dissipative phenomenon, as mentioned in Section 4).

We mention that in the special case in which $\sigma =1$, $\phi =0$, $\mathrm{grad}\mathcal{V}=0$ and $f=0$, the function (29) coincides with the invariant described in Theorem 3.4 of paper [34].

6. Numerical Recipe to Simulate a Fourth-Order Dynamical System on Manifold

Numerical simulations of dynamical systems are currently an integral part of system modeling and control as they afford predictions on the quality of a mathematical model as well as of related control schemes. In the present section, we shall briefly tackle the problem of numerically simulating the behavior of a fourth-order system of the kind (11) for illustrative purposes. We shall underline that, even though the mathematical model (11) does not include a jerk-dependent term, third-order covariant derivatives nonetheless enter the numerical recipe; hence, their numerical approximations are discussed as well.

It is convenient to preliminarily recall two formulas to numerically approximate the covariant derivative of a smooth vector field $w\in \Gamma (T\mathbb{M})$ at a point $x(t)\in \mathbb{M}$ of a smooth trajectory along the tangent direction $\dot{x}(t)\in T_{x(t)}\mathbb{M}$, namely,

$$\begin{array}{c}\hfill \mathrm{Forward} \mathrm{Euler} \mathrm{method} (\mathrm{fEul}):{({\nabla }_{\dot{x}}w)}_{x(t)}\approx \frac{{\mathbb{P}}^{x(t+h)\to x(t)}{w}_{x(t+h)}-w_{x(t)}}{h},\end{array}$$

(34)

$$\begin{array}{c}\hfill \mathrm{Backward} \mathrm{Euler} \mathrm{method} (\mathrm{bEul}):{({\nabla }_{\dot{x}}w)}_{x(t)}\approx \frac{w_{x(t)}-{\mathbb{P}}^{x(t-h)\to x(t)}{w}_{x(t-h)}}{h},\end{array}$$

(35)

which play the role of path-dependent approximate covariant derivative formulas. The constant $h>0$ denotes a numerical stepsize, and even though the above formulas represent numerical approximations, they are nonetheless manifold-calculus-compliant. (By this, we mean that, albeit being numerical approximations, they produce results that belong to the correct tangent spaces.) By employing the above two approximation rules, it is possible to write a numerical recipe to numerically simulate the behavior of the system (11).

Let us note down the following approximation formulas for the covariant jerk and the covariant jounce that hold true by virtue of the bEul approximation method:

$$\begin{array}{c}{({\nabla }_v{\nabla }_{v}v)}_{x(t)}\approx \frac{{({\nabla }_{v}v)}_{x(t)}-{\mathbb{P}}^{x(t-h)\to x(t)}{({\nabla }_{v}v)}_{x(t-h)}}{h},\hfill \end{array}$$

(36)

$$\begin{array}{c}{({\nabla }_v{\nabla }_v{\nabla }_{v}v)}_{x(t)}\approx \frac{{({\nabla }_v{\nabla }_{v}v)}_{x(t)}-{\mathbb{P}}^{x(t-h)\to x(t)}{({\nabla }_v{\nabla }_{v}v)}_{x(t-h)}}{h}.\hfill \end{array}$$

(37)

From relationship (37), it follows that

$${({\nabla }_v{\nabla }_{v}v)}_{x(t)}={\mathbb{P}}^{x(t-h)\to x(t)}{({\nabla }_v{\nabla }_{v}v)}_{x(t-h)}+h{({\nabla }_v{\nabla }_v{\nabla }_{v}v)}_{x(t)}+\mathcal{O}(h^2),$$

(38)

while from relationship (36), it follows that

$$\begin{array}{cc}\hfill {({\nabla }_{v}v)}_{x(t)}& ={\mathbb{P}}^{x(t-h)\to x(t)}{({\nabla }_{v}v)}_{x(t-h)}+h{({\nabla }_v{\nabla }_{v}v)}_{x(t)}+\mathcal{O}(h^2)\hfill \\ \hfill & ={\mathbb{P}}^{x(t-h)\to x(t)}{({\nabla }_{v}v)}_{x(t-h)}+h{\mathbb{P}}^{x(t-h)\to x(t)}{({\nabla }_v{\nabla }_{v}v)}_{x(t-h)}\hfill \\ \hfill & +h^2{({\nabla }_v{\nabla }_v{\nabla }_{v}v)}_{x(t)}+\mathcal{O}(h^3),\hfill \end{array}$$

(39)

where the Landau symbol $\mathcal{O}$ is characterized by the property that ${lim}_{h\to 0}{h}^{1-n}\mathcal{O}(h^n)=0$ for $n\in \{2,3\}$.

Let us now assume that a numerical approximation of the solution to the differential equation $\dot{x}=v_x$ on a tangent bundle $T\mathbb{M}$ is sough. Here, $v:\mathbb{M}\to T\mathbb{M}$ denotes a known vector field, and x denotes the path on the base manifold $\mathbb{M}$ that satisfies such a first-order differential equation. A fEul-type stepping method on manifold to advance the solution to such a differential equation starting from the initial condition $x(0)=x_0\in \mathbb{M}$ reads

$$x(t+h)\approx {exp}_{x(t)}(h{v}_{x(t)}),t\ge 0.$$

(40)

Similarly, a fEul-type approximation method on $T\mathbb{M}$ to recover a vector field from its covariant derivative and from the initial condition $v_{x(0)}=v_0\in T_{x(0)}\mathbb{M}$ reads

$$v_{x(t+h)}={\mathbb{P}}^{x(t)\to x(t+h)}(v_{x(t)})+h{\mathbb{P}}^{x(t)\to x(t+h)}({({\nabla }_{v}v)}_{x(t)})+\mathcal{O}(h^2).$$

(41)

The above approximations afford estimating the trajectory of a dynamical system step by step according to the chosen stepsize value. We are not diving into further details here, while other resources from the present author are certainly more exhaustive.

7. Conclusions

The body of the present paper described an approach to formulate fourth-order dynamical systems on manifold on the basis of the Hamilton–d’Alembert stationary action/virtual work principle.

In particular, the present paper introduced a Lagrangian function that depends on the kinetic energy and the covariant acceleration energy associated with an abstract dynamical system whose state evolves on a possibly curved Riemannian manifold, as well as a potential energy function. The corresponding fourth-order equation of motion is derived through a variational approach. Depending on the values of the coefficients introduced in the equations, the obtained system may also be of the third order (i.e., jerk-type), second order (i.e., Newton-type) or first order (i.e., kinetic).

Moreover, the present paper introduced the notion of Rayleigh differential form that affords taking into account non-conservative exogenous forces, such as braking effects and active control, into the devised general equation of dynamics. The paper ended with a short discussion of a numerical recipe to approximately simulate the introduced differential-type system.

As a distinguishing feature of the present paper, following the general line of research of the present author, the devised fourth-order system, as well as the presented recipe, to numerically simulate its behavior, were written in coordinate-free terms. Such a mindset makes the devised notation prone to seamless tailoring to arbitrary Riemannian state manifolds of any dimension.

We trust that the presented mathematical framework may be seamlessly extended to deal with applications that require higher orders of covariant differentiation, as in minimum jerk trajectory generation [47] or minimum jounce trajectory generation [48].