remarkRemark hypothesisHypothesis claimClaim exampleExample Nonholonomic
Herglotz IntegratorE. Maciel, I. Ortiz, and C. E. Schaerer ex_supplement

# A Herglotz-based integrator for nonholonomic mechanical systems

Elias Maciel  Inocencio Ortiz  Christian E. Schaerer<EMAIL_ADDRESS><EMAIL_ADDRESS><EMAIL_ADDRESS>
###### Abstract

We propose a numerical scheme for the time-integration of nonholonomic
mechanical systems, both conservative and nonconservative. The scheme is
obtained by simultaneously discretizing the constraint equations and the
Herglotz variational principle. We validate the method using numerical
simulations and contrast them against the results of standard methods from the
literature.

Polytechnic School - National University of Asuncion

San Lorenzo SL 2160 - Paraguay

††2020 Mathematics Subject Classification. Primary: 37M15, 65D30; Secondary:
70G45††Key words and phrases. Nonholonomic systems, Nonconservative systems,
Herglotz principle, Contact integrator.

##  1 Introduction

Numerical integration of differential equations is a very active research
area, both from a theoretical perspective as well as a more practical or
computational standpoint. This last perspective has had a boost with the
advent of digital computers, and we now have a plethora of methods, some of
general purposes and some of more ad-hoc nature.

In the last category, we find the structure-preserving integrators. These are
integrators designed to take into account specific properties of the dynamical
systems they are meant to solve. Relevant examples of structure-preserving
integrators are the symplectic Euler and the symplectic Runge-Kutta, both of
which are compatible with the symplectic structure naturally associated with
Hamiltonian systems. When the relevant structure is geometric in nature, the
structure-preserving integrators are usually known as geometric integrators.
For more details of these and other developments, see [25].

In the last two decades, the differential geometric description of dynamical
systems, known as geometric mechanics, and the variational formulation of
mechanics, have provided a successful approach for developing geometric
integrators. Arguably, the seminal result in this context was the fact that,
for Hamiltonian systems, discretizing the Hamiltonian principle instead of the
Euler-Lagrange equation of motion automatically yields an integrator that
preserves the symplectic structure naturally associated with the system, i.e.,
it yields a symplectic integrator. This result can be traced back at least to
[63], and ever since, integrators coming from the discretization of some
variational principle are known as variational integrators.

The results in [63] was expanded in [53] to consider nonconservative systems
by discretizing the Lagrange-d’Alembert principle. There, the authors also
showed that any symplectic integrator is a variational integrator. Since then,
in the following two decades, there has been several generalizations in an
attempt to incorporate systems of more general nature, as for instance:
coupled multi-body systems, field theories, electrical systems [22, 45, 60,
56, 55, 50]; systems with more specific features (as symmetries) [36, 59, 21,
52, 35]; interactions with control theory [11, 37, 8, 30, 16]; higher-order
and asynchronous techniques [41, 26, 39]. For more recent developments, see
[46, 49, 23, 48, 61, 38, 10, 12]. For surveys and reviews, see [42, 40, 43]
and the references therein.

In this work, we construct a numerical integrator for nonconservative and
nonholonomically constrained systems. To this end, we combine a criterion for
the discretization of nonholonomic constraints (introduced in [14], in the
context of Lagrange-d’Alembert principle) and a recently developed integrator
for dissipative systems (introduced in [62]) coming from a discrete version of
the Herglotz’s variational principle. Then we validate our integrator by
running simulations of two archetypal mechanical systems in order to contrast
it with the outcome of more standard procedures available in the literature.

The organization of this work is as follows. In Section 2 we give a brief
background on geometric mechanics and variational integrators, mainly to fix
notation and terminology. In Section 3 we provide the technical details of our
methodology to get our integrator. In Section 4 we show the outcome of several
numerical simulations, contrasted with a couple of more standard procedures.
Finally, in Section 5 we offer some concluding remarks. We add an appendix
containing the explicit equations for the simulations.

##  2 Geometric mechanics and Variational integrators

Geometric mechanics is, roughly speaking, the description of mechanics in the
language of differential geometry [28, 29, 27, 1, 51]. In many interesting
cases, the configuration space of a mechanical system with $n\in\mathbb{N}$
degrees of freedom can be modeled as a $n$-dimensional smooth manifold $Q$,
and its state space can be modeled by the tangent bundle $TQ$, of the
configuration space, or a subset of it. Kinetic energy is given by a
Riemannian metric on $Q$, potential energy is given by a function $V\colon
Q\to\mathbb{R}$, forces are given by $1$-forms $F\colon Q\to T^{*}Q$, and the
dynamic is encoded by a Lagrangian function $L\colon TQ\to\mathbb{R}$, which
typically involves the kinetic and potential energy [3, 51].

Constraints of motion are relations among the configuration variables and its
derivatives. If those relations involve only the configuration variables, the
constraint is called holonomic, otherwise it is called nonholonomic. From a
geometric point of view, this dichotomy can be seen as follows: holonomically
constrained systems are those for which the configuration space is a smooth
manifold $Q$, and the state space is the entire tangent bundle $TQ$, which
means that, at each configuration $q\in Q$, the entire tangent space $T_{q}Q$
is allowed for the velocity of the system. On the other hand, nonholonomically
constrained systems are those whose configuration space is a smooth manifold
$Q$, but the state space is a subset ${\mathcal{D}}$ of the tangent bundle
$TQ$, which means that at each configuration $q\in Q$, only a subset
${\mathcal{D}}_{q}\subset T_{q}Q$ is allowed for the velocity of the system
[3, 13, 44]. For many interesting cases of nonholonomic systems, the
constraint subset ${\mathcal{D}}$ is an affine, or even linear, subbundle of
$TQ$ (i.e.: if $Q$ is the configuration manifold, for each point $q\in Q$ the
velocity of the system is constrained to a subspace ${\mathcal{D}}_{q}$
(linear or affine) of the tangent space $T_{q}Q$) [15, 9]. In this work we
will focus on this kind of constraints. Readers interested in more general
cases, including nonlinearity and time dependency may consult the following
references [31, 32, 33, 20, 58, 19].

######  Example 2.1 (Rolling without slipping).

Figure 1 illustrates a solid sphere rolling without slipping on a fixed
horizontal surface. The configuration of this sphere is fully determined by
the point $(x,y)\in\mathbb{R}^{2}$ of contact between the sphere and the
plane, and an element of the special orthogonal group $SO(3)$ (rotations on
$\mathbb{R}^{3}$). Hence, its configuration space is $Q=\mathbb{R}^{2}\times
SO(3)$, but the nonslipping condition in the contact point imposes a
restriction on the velocities; thus, the state space is not the entire tangent
bundle $TQ$. This is an archetypal example of a nonholonomically constrained
mechanical system.

Figure 1: A solid sphere rolling without slipping on a horizontal surface.

Variational formalism: Given a mechanical system with configuration space $Q$
and constraint subset ${\mathcal{D}}\subset TQ$, the dynamic (or equation of
motion) of the system is modeled by a differential equation on $TQ$ subject to
the restriction ${\mathcal{D}}$. In geometric mechanics, it is common to
describe the equation of motion as the extreme of a suitable variational
principle. In what follows, we are going to use some standard constructions
from the calculus of variations, so for the reader’s convenience, let us
introduce them for future reference.

Given two points $q_{0}$ and $q_{f}$ in a smooth manifold $Q$, and
$t_{0}<t_{f}\in\mathbb{R}$, let
$\mathcal{C}(t_{0},t_{f}):=\\{q\colon[t_{0},t_{f}]\to Q\\}$ be the space of
smooth paths satisfying $q(t_{0})=q_{0}$ and $q(t_{f})=q_{f}$. A smooth
variation with fixed endpoints of an element
$\overline{q}\in\mathcal{C}(t_{0},t_{f})$ is a smooth map
$q\colon(-\epsilon,\epsilon)\times[t_{0},t_{f}]\to Q$ such that, for each
$s\in(-\epsilon,\epsilon)\subset\mathbb{R}$, $q_{s}(t):=q(s,t)$ for all
$t\in[t_{0},t_{f}]$ is an element of $\mathcal{C}(t_{0},t_{f})$, and
$q_{0}=\overline{q}$. The variational vector field associated to $q$ is the
vector field along $\overline{q}$ given by

$\delta q(\overline{q}(t)):=\frac{\partial}{\partial
s}\Big{|}_{s=0}q(s,t),\forall t\in[t_{0},t_{f}].$

The simplest scenario for a variational formulation of mechanics is that of a
conservative system with holonomic constraints. For this kind of system, the
equation of motion is given by Hamilton’s Variational Principle, which can be
stated as follows.

######  Definition 2.2 (Hamilton’s Principle).

Let $L\colon TQ\to\mathbb{R}$ be a conservative Lagrangian system. Given two
points $q_{0}$ and $q_{f}$ in the configuration space $Q$, and
$t_{0}<t_{f}\in\mathbb{R}$, the path
$\overline{q}(t)\in\mathcal{C}(t_{0},t_{f})$ followed by the system to go from
$q_{0}$ to $q_{f}$ is the one that extremizes the action

$\int_{t_{0}}^{t_{f}}L(q,\dot{q})dt.$

The optimization mentioned in the Hamilton’s Principle is in the sense that,
for all smooth variation $q(s,t)$ of $\overline{q}$ in
$\mathcal{C}(t_{0},t_{f})$, with $s$ in a sufficiently small interval
$(-\epsilon,\epsilon)\subset\mathbb{R}$, the functional variation of the
action vanishes, i.e.:

$\delta\int_{t_{0}}^{t_{f}}L(q,\dot{q})dt:=\frac{\partial}{\partial
s}\bigg{|}_{s=0}\left(\int_{t_{0}}^{t_{f}}L(q(s,t),\dot{q}(s,t))\right)dt=0.$

It is well known that Hamilton’s principle is not appropriate for systems that
are either nonconservative or with nonholonomic constraints [3, 13]. The most
widely used variational-like principle to obtain the equation of motion for
nonconservative and nonholonomic systems is the Lagrange-d’Alembert principle,
together with the so-called nonholonomic principle. This situation can be
stated as follows.

######  Definition 2.3 (L-A Principle).

Let $L\colon TQ\to\mathbb{R}$ be the Lagrangian of a system, subject to
external forces $F^{e}\colon Q\to T^{*}Q$, and nonholonomic constraints
${\mathcal{D}}\subset TQ$. Given two points $q_{0}$ and $q_{f}$ in $Q$, and
$t_{0}<t_{f}\in\mathbb{R}$, the path $\overline{q}(t)$ followed by the system
to go from $q_{0}$ to $q_{f}$ is the one that solves

(1) $\delta\int_{t_{0}}^{t_{f}}L(q,\dot{q})dt+\int_{t_{0}}^{t_{f}}F^{e}\delta
q=0,$

subject to

(2) $\begin{split}\delta
q(\overline{q}(t))&\in{\mathcal{D}}_{\overline{q}(t)},\forall
t\in[t_{0},t_{f}],\\\
\dot{\overline{q}}(t)&\in{\mathcal{D}}_{\overline{q}(t)},\forall
t\in[t_{0},t_{f}].\end{split}$

######  Remark 2.4.

Equation (2) is known as the nonholonomic principle. It basically says that
both the solution as well as the variational vector field must satisfy the
restriction. However, notice that the curves $q_{s}(t)$, given by the
variation $q(s,t)$, are not asked to satisfy the restrictions. Imposing
restrictions on those curves poses a different problem, namely, one of optimal
control, and in general, it does not provide the equation of motion we are
interested in. Readers interested in this dichotomy may consult [3] and the
references therein.

There is an alternative way to handle an interesting spectrum of
nonconservative systems by using the Herglotz variational principle [24, 64].
As it turns out, Herglotz principle has a close relationship with contact
geometry, which in turn has been found to be a natural arena for several kinds
of dissipative systems [17, 4, 24, 5, 7, 6]. In the context of Herglotz
principle, the nonconservative nature of mechanical systems with configuration
space $Q$ are usually incorporated via a function $L\colon
TQ\times\mathbb{R}\to\mathbb{R}$. Due to the aforementioned relationship with
contact geometry, we will refer to this function as a contact-type Lagrangian,
to distinguish it from the usual Lagrangian.

######  Remark 2.5.

The term action-dependent Lagrangian is also used in the literature for
Lagrangians in the context of Herglotz principle (see [34, 18, 47] and the
references therein).

######  Definition 2.6 (Herglotz Principle).

Let $Q$ be the configuration manifold of a mechanical system whose dynamic is
encoded by a contact-type Lagrangian $L\colon TQ\times\mathbb{R}\to\mathbb{R}$
as follows: given $q_{0},q_{f}\in Q$ and $t_{0}<t_{f}\in\mathbb{R}$, let
${\mathcal{Z}}$ be the map that assigns to each
$\xi\in{\mathcal{C}}(t_{0},t_{f})$, the solution of the initial value problem

(3) $\dot{z}=L(\xi,\dot{\xi},z(t));\quad z(t_{0})=z_{0}.$

Then, the Herglotz principle establishes that the path $\overline{q}$ followed
by the system to go from $q_{0}$ to $q_{f}$ is the one that extremizes the
functional ${\mathcal{Z}}(\xi)(t_{f})$.

Table 1 summarizes some kinds of mechanical systems and the variational
principles to obtain their corresponding equations of motion.

| Holonomic | Nonholonomic
---|---|---
Conservative |  | •Hamilton’s principle
---
$\delta\int_{a}^{b}L(q,\dot{q})\mathop{\mathrm{\mathstrut d}}\nolimits\\!t=0$.
| •Lagrange-d’Alembert principle
---
$\delta\int_{a}^{b}L(q,\dot{q})\mathop{\mathrm{\mathstrut d}}\nolimits\\!t=0$,
$\sum_{k=1}^{n}a_{k}^{j}\delta q^{k}=0$.
Nonconservative |  | •Lagrange-d’Alembert principle
---
$\delta\int_{a}^{b}L(q,\dot{q})\mathop{\mathrm{\mathstrut
d}}\nolimits\\!t+\int_{a}^{b}F\delta q\mathop{\mathrm{\mathstrut
d}}\nolimits\\!t=0$.
•Herglotz principle
$\delta z(b)=0$.
| •Lagrange-d’Alembert principle
---
$\delta\int_{a}^{b}L(q,\dot{q})\mathop{\mathrm{\mathstrut
d}}\nolimits\\!t+\int_{a}^{b}F^{e}\delta q\mathop{\mathrm{\mathstrut
d}}\nolimits\\!t=0$,
$\sum_{k=1}^{n}a_{k}^{j}\delta q^{k}=0$.
•Herglotz principle with constraints
$\delta z(b)=0$,
$\sum_{k=1}^{n}a_{k}^{j}\delta q^{k}=0$.
Table 1: Types of mechanical systems and the corresponding variational
principles from which to obtain their equations of motion.

The variable $z$ in the contact-type Lagrangian for the Herglotz principle
allows to incorporate dissipation directly. This makes the Herglotz principle
very attractive to obtain integrators that are well suited for dissipative
systems, as argued in [62]. Next we review the main ideas behind variational
integrators.

Variational integrators: The main paradigm to get integrators from a
variational formulation was introduced in [63, 53]. There, the authors
essentially propose the following procedure: Given a holonomic system $L\colon
Q\to\mathbb{R}$, subject to external, nonconservative forces $F^{e}\colon
TQ\to T^{*}Q$, discretize the force as $F_{d}^{+},F_{d}^{-}\colon Q\times Q\to
T^{*}Q$, as well as the action functional
$\int_{t_{0}}^{t_{f}}L(q,\dot{q})dt$, namely, let $\\{t_{i}\\}_{i=0}^{N}$ be a
partition of the interval $[t_{0},t_{f}]$, and consider

$L_{d}(q_{j},q_{j+1})\approx\int_{t_{j}}^{t_{j+1}}L(q,\dot{q})dt,$

where $L_{d}(q_{j},q_{j+1})$ is the approximation given by some quadrature of
the integral, and $\\{q_{j}\\}_{j=0}^{N}$ is a sequence of configuration
points which approximate the true trajectory $q(t)$ at the times $t_{j}$, for
$j=0,1,2,\ldots,N$. The sequence $\\{q_{j}\\}_{j=0}^{N}$ is chosen as the one
that satisfies the expression

$\delta\left(\sum_{j=0}^{N-1}L_{d}(q_{j},q_{j+1})\right)+\sum_{j=0}^{N-1}[F_{d}^{-}(q_{j},q_{j+1})\cdot\delta
q_{j}+F_{d}^{+}(q_{j},q_{j+1})\cdot\delta q_{j+1}]=0$

for all variations $\\{\delta q_{j}\\}_{j=0}^{N}$ vanishing at the endpoints.

For holonomic systems, the variations can be taken in arbitrary directions in
each tangent space $T_{q_{j}}Q$, and thus we get the so-called forced discrete
Euler-Lagrange equations:

$D_{1}L_{d}(q_{j},q_{j+1})+D_{2}L_{d}(q_{j-1},q_{j})+F_{d}^{+}(q_{j-1},q_{j})+F_{d}^{-}(q_{j},q_{j+1})=0,$

where $D_{i}$, for $i=1,2$ denotes the derivative with respect to the $i$-th
argument. Under some mild regularity conditions, these equations yield an
integrator

$(q_{j-1},q_{j})\mapsto(q_{j},q_{j+1}).$

One of the main results in [53] is that, for conservative, holonomic systems,
this integrator preserves the symplectic structure naturally associated with
the given Lagrangian system, i.e., it yields a symplectic integrator. The
symplecticity of the integrator results in very good energy conservation
behavior for long-time simulations.

In this article, we are interested in the approach given in [14] to handle
nonholonomic systems. There, the authors construct an integrator for
nonholonomic systems from the Lagrange-d’Alembert principle, where a main
result is a criterion for discretizing the nonholonomic constraint in a
compatible way with the discretization of the Lagrangian. To follow up, let us
observe that the discrete Lagrangian $L_{d}$ can be interpreted as a map
$L_{d}\colon Q\times Q\to\mathbb{R}$, which is obtained by means of a
discretization map $\Psi\colon Q\times Q\to TQ$ as $L_{d}=L\circ\Psi$. Then,
if ${\mathcal{D}}\subset TQ$ is the restriction distribution, they propose to
take $D:=\Psi({\mathcal{D}})\subset Q\times Q$ as the discretized restriction.
Thus, they get the following integrator:

(4)
$\begin{split}D_{1}L_{d}(q_{j},q_{j+1})+D_{2}L_{d}(q_{j-1},q_{j})+F_{d}^{+}(q_{j-1},q_{j})+F_{d}^{-}(q_{j},q_{j+1})&=\lambda_{a}\Phi^{a},\\\
\Phi_{d}^{a}(q_{j},q_{j+1})&=0,\end{split}$

where $\lambda_{a},a\in\\{1,2,...,m\\}$ are Lagrange multipliers and
$\Phi^{a}$ are linearly independent 1-forms defining the annihilator
${\mathcal{D}}^{0}\subset T^{*}Q$ of ${\mathcal{D}}$.

In [62], the authors derive an integrator for dissipative systems with
holonomic constraints by discretizing the Herglotz principle and show that it
is compatible with the contact structure naturally associated to dissipative
systems. This compatibility is completely parallel to the one found between
symplectic integrators and the symplectic structure of Hamiltonian systems,
and thus they named it a contact integrator.

##  3 Construction of the integrator

In this section, we construct an integrator for nonholonomically constrained
dissipative system following ideas from [17], [62], and [14].

###  3.1 Herglotz principle with nonholonomic constraints

Consider a dissipative mechanical system of contact-type given by $L\colon
TQ\times\mathbb{R}\to\mathbb{R}$ (in the setting of Herglotz principle as in
Definition 2.6), subject to a nonholonomic constraint ${\mathcal{D}}\subset
TQ$. Hence, the Herglotz principle is stated as follows: the path
$\overline{q}\in{\mathcal{C}}(t_{0},t_{f})$ followed by the system is the one
that extremizes the functional ${\mathcal{Z}}(\xi)(t_{f})$, among variations
tangent to ${\mathcal{D}}$, and satisfies itself the constraints. This means
that, at each time $t\in[t_{0},t_{f}]$, the velocity $\dot{\overline{q}}(t)$
of the path $\overline{q}$ as well as the value $\delta q(t)$ of the
variational vector field $\delta q$ must belong to
${\mathcal{D}}_{\overline{q}(t)}$.

To express this in a more concrete way, let us consider that the restriction
${\mathcal{D}}$ is (locally) defined by the vanishing of $m$ functions
$\Phi^{j}\colon TQ\to\mathbb{R}$, linear in the velocities. This means that on
a local chart $(q,\dot{q})$ the restrictions are given as the kernel of a
matrix $A(q)$; namely, locally we have

${\mathcal{D}}=\\{(q,\dot{q});A(q)\dot{q}=0\\}.$

Then, it can be shown that $\overline{q}$ solves the nonholonomic Herglotz
principle if, and only if, it satisfies the following equations (see [17] for
details):

(5) $\begin{split}\frac{d}{dt}\frac{\partial
L}{\partial\dot{q}_{i}}-\frac{\partial L}{\partial q_{i}}-\frac{\partial
L}{\partial\dot{q}_{i}}\frac{\partial L}{\partial
z}&=\sum_{j=1}^{m}\lambda_{j}A_{i}^{j}(q),\\\ A(q)\dot{q}&=0,\end{split}$

for some Lagrange multipliers $\lambda_{j}$.

###  3.2 Discrete Herglotz Variational Principle

Consider the discretized Lagrangian $L_{d}\colon
Q^{2}\times\mathbb{R}^{2}\to\mathbb{R}$, and the discrete Herglotz principle,
stated as follows: given a discrete curve $q=(q_{0},q_{1},\ldots,q_{N})\in
Q^{N+1}$, define $z=(z_{0},\ldots z_{N})\in\mathbb{R}^{n}$ by $z_{0}=0$ and

(6) $z_{j+1}-z_{j}=hL_{d}(q_{j},q_{j+1},z_{j},z_{j+1}).$

Then, a discrete curve
$\overline{q}=(\overline{q}_{0},\ldots,\overline{q}_{N})$ is called a solution
of the discrete Herglotz principle if

$\frac{\partial z_{j+1}}{\partial\overline{q_{j}}}=0,$

for all $j\in\\{1,2,\ldots,N\\}$.

It can be shown (see Theorem 1 in [62]), that a discrete curve
$\overline{q}=(\overline{q}_{0},\ldots,\overline{q}_{N})$ solves the discrete
Herglotz principle if, and only if, it satisfies the discrete generalized
Euler-Lagrange equations

(7)
$\begin{split}0&=D_{2}L_{d}(\overline{q}_{j-1},\overline{q}_{j},z_{j-1},z_{j})\frac{1+hD_{3}L_{d}(\overline{q}_{j},\overline{q}_{j+1},z_{j},z_{j+1})}{1-hD_{4}L_{d}(\overline{q}_{j-1},\overline{q}_{j},z_{j-1},z_{j})}\\\
&+D_{1}L_{d}(\overline{q}_{j},\overline{q}_{j+1},z_{j},z_{j+1}).\end{split}$

Assuming the nondegeneracy condition
$D_{1}D_{2}L_{d}(q_{j},q_{j+1},z_{j},z_{j+1})\neq 0$, equation (7) yields an
integrator $(q_{j-1},q_{j})\mapsto(q_{j},q_{j+1})$, for the unconstrained
system $L\colon TQ\times\mathbb{R}\to\mathbb{R}$. In [62], it is proved that
this integrator preserves the natural contact structure associated to the
discrete contact-type Lagrangian $L_{d}$, hence the authors called it a
contact integrator.

######  Remark 3.1.

A simpler (and perhaps more natural) discretization of $TQ\times\mathbb{R}$
would be $Q^{2}\times\mathbb{R}$, as was used in [2]. This would lead to a
discretized action-type Lagrangian $L_{d}$ without the $z_{j+1}$ argument and
without the $D_{4}L_{d}$ term in (7). As pointed out in [62], for the
theoretical development of contact integrators, there is no need for this
extra argument in $L_{d}$. However, the discretization without the $z_{j+1}$
argument is of first order, whereas the one with $z_{j+1}$ is of second order,
and this is relevant for numerical behavior. In our work, for the numerical
experiments, we will use both alternatives.

###  3.3 Discretization of the Constraint

Now, let us bring into play the constraint ${\mathcal{D}}\subset TQ$. Let us
first recall the method introduced in [14] to incorporate a nonholonomic
constraint into the integrator derived from the Lagrange-d’Alembert principle,
given a mechanical system $L\colon TQ\to\mathbb{R}$. Since the discretization
of the Lagrangian yields a map $L_{d}\colon Q\times Q\to\mathbb{R}$, the
discrete constraint must be a submanifold $\Delta_{d}\subset Q\times Q$, which
satisfies $(q,q)\in\Delta_{d}$ for all $q\in Q$, (i.e., $\Delta_{d}$ has to
contain the diagonal of $Q\times Q$). This submanifold imposes a restriction
to the motion in the sense that the discrete curve
$q=(q_{0},q_{1},\ldots,q_{N})$ must satisfy that
$(q_{j},q_{j+1})\in\Delta_{d}$, and the variation $\delta q$ must satisfy the
condition $\delta q(q_{j})\in{\mathcal{D}}_{q_{j}}$, i.e., the allowed
variations are those who respect the continuous constraint
${\mathcal{D}}\subset TQ$. In order to construct $\Delta_{d}$ in a compatible
way with the discrete Lagrangian, it is assumed that there is a discretizing
map $\Psi\colon Q\times Q\to TQ$ such that $L_{d}=L\circ\Psi$, and then,
$\Delta_{d}:=\\{(q,q^{\prime})\in Q\times
Q;\Psi(q,q^{\prime})\in{\mathcal{D}}\\}$. This means that, if ${\mathcal{D}}$
is (locally) defined by the vanishing of the functions $\Phi^{c}\colon
TQ\to\mathbb{R}$, then $\Delta_{d}$ is (locally) defined by the vanishing of
the functions $\Phi^{c}_{d}:=\Phi^{c}\circ\Psi$.

In our situation, as we have a contact-type Lagrangian $L\colon
TQ\times\mathbb{R}\to\mathbb{R}$, we might consider a discretizing map of the
form $\Psi=\Psi_{1}\otimes\Psi_{2}\colon Q^{2}\times\mathbb{R}^{2}\to
TQ\times\mathbb{R}$, such that
$\Psi(q,q^{\prime},z,z^{\prime})=(\Psi_{1}(q,q^{\prime}),\Psi_{2}(z,z^{\prime}))$.
Thus, we want our discrete restriction to be compatible with the discretizing
map $\Psi_{1}\colon Q^{2}\to TQ$, in the same way as discussed in the previous
paragraph.

More concretely, suppose as before that on a chart $(q,\dot{q})$ of $TQ$ our
restriction ${\mathcal{D}}$ is given in the form

${\mathcal{D}}=\\{(q,\dot{q})\in TQ;A(q)\dot{q}=0\\},$

for some matrix $A(q)$. Let $(q_{d},\dot{q}_{d})\in TQ$ be the image of a
point $(q,q^{\prime})\in Q^{2}$ under the map $\Psi_{1}$. Then we define the
discrete restriction map as $A_{d}(q,q^{\prime})=A(q_{d})\dot{q}_{d}$, and the
discrete restriction as

$\Delta_{d}:=\\{(q,q^{\prime})\in Q^{2};A_{d}(q,q^{\prime})=0\\}.$

Therefore, the discrete generalized Euler-Lagrange equations with nonholonomic
constraints linear in the velocities are

(8)
$\begin{split}D_{1}L_{d}(q_{j},q_{j+1},z_{j},z_{j+1})+D_{2}L_{d}(q_{j-1},q_{j},z_{j-1},z_{j})\frac{1+hD_{3}L_{d}(q_{j},q_{j+1},z_{j},z_{j+1})}{1-hD_{4}L_{d}(q_{j-1},q_{j},z_{j-1},z_{j})}&=\lambda
A(q_{j}),\\\ A_{d}(q_{j},q_{j+1})&=0,\end{split}$

which correspond to the equations of a contact integrator with nonholonomic
constraints linear in the velocities.

##  4 Numerical simulations and comparison

In this section, we will apply our integrator to two nonholonomic mechanical
systems, namely the Foucault pendulum and the falling rolling disk.

###  4.1 Foucault pendulum

We analyze the Foucault pendulum with Rayleigh dissipation. We derive its
Lagrangian, both for the Lagrange-d’Alembert as well as for the Herglotz
formulation, in order to compare our integrator with the one coming from the
Lagrange-d’Alembert principle.

General setting and simplifications: The Foucault pendulum consists of a
pendulum of length $l$ and mass $m$ located at latitude $\beta$ on Earth’s
surface. We may consider here a Rayleigh dissipation with parameter $\alpha$.
The plane of oscillation does not rotate in a reference $(X,Y,Z)$ fixed in
space, hence as Earth rotates, this plane rotates with respect to a reference
$(x,y,z)$ attached to Earth (see Figure 2).

Figure 2: Foucault pendulum and its inertial and noninertial frames.

To model this problem, we consider the inertial frame $(X,Y,Z)$ with its
origin at the center of the Earth and $Z$ passing through the north pole.
Hence, the angular velocity of Earth, $\Omega$, points along $Z$. On the other
hand, the noninertial frame $(x,y,z)$ is such that $x$ points along a meridian
in the south direction, $y$ points to the east along the parallel $\beta$ and
$z$ coincides with the vertical at the pendulum location. Hence, the vector
position $r=(x,y,z)$ of a particle in the noninertial frame satisfies the
relation $\omega=\frac{r\times\dot{r}}{\|r\|^{2}}$, where
$\omega=(\Omega\cos\beta,0,-\Omega\sin\beta)$ is the angular velocity of the
particle in the noninertial frame.

As the pendulum forms a small angle $\phi$ in its oscillatory motion, the
coordinates $x,y$ are of order $l\phi$, whereas $z$ is of order $l\phi^{2}$
and so it is negligible. Thus, we may consider the movement of the mass
pendulum in the plane ${z=0}$ and take $q=(x,y)$ as generalized coordinates.
With these considerations, the kinetic energy $K(\dot{q})$ and potential
energy $V(q)$, in terms of the mass $m$, the length $l$ and the gravitational
acceleration $g$, are given by:

(9) $K(\dot{q})=\frac{1}{2}m(\dot{x}^{2}+\dot{y}^{2})\qquad\text{and}\qquad
V(q)=\frac{1}{2}m\frac{g}{l}(x^{2}+y^{2}),$

while the relation $\omega=\frac{r\times\dot{r}}{\|r\|^{2}}$ reads

(10) $-y\dot{x}+x\dot{y}+\Omega\sin\beta(x^{2}+y^{2})=0,$

which is a nonholonomic constraint for the system.

Lagrange-d’Alembert description: In this context, the Lagrangian is
$L(q,\dot{q})=K(\dot{q})-V(q)$, which according to our previous computations
reads

(11)
$L(q,\dot{q})=\frac{1}{2}m(\dot{x}^{2}+\dot{y}^{2})-\frac{1}{2}m\frac{g}{l}(x^{2}+y^{2}),$

subject to the constraint (10). The damping is modeled as an external force
$F(q,\dot{q})=-\alpha m\dot{q}$.

To derive the integrator for this case, we use equation (4) with a linear-
order quadrature for both the Lagrangian and the external force, while the
discrete constraint is obtained using the discretizing map $\Psi$ on equation
(10) (see [14]). The resulting equations are

(12)
$\begin{split}\frac{-x_{j+1}+2x_{j}-x_{j-1}}{h}-h\frac{g}{l}x_{j}-\alpha(x_{j+1}-x_{j})+\lambda_{1}\frac{y_{j}}{m}&=0\\\
\frac{-y_{j+1}+2y_{j}-y_{j-1}}{h}-h\frac{g}{l}y_{j}-\alpha(y_{j+1}-y_{j})-\lambda_{1}\frac{x_{j}}{m}&=0\\\
-y_{j}\frac{x_{j+1}-x_{j}}{h}+x_{j}\frac{y_{j+1}-y_{j}}{h}+\Omega\sin\beta(x_{j}^{2}+y_{j}^{2})&=0,\end{split}$

where $\lambda_{1}$ is a Lagrange multiplier.

Herglotz description: Here, the Lagrangian is taken as
$L(q,\dot{q},z)=K(\dot{q})-V(q)-\alpha z$, which, according to our previous
analysis becomes

(13)
$L(q,\dot{q},z)=\frac{m}{2}(\dot{x}^{2}+\dot{y}^{2})-\frac{mg}{2l}\left(x^{2}+y^{2}\right)-\alpha
z,$

also subject to the nonholonomic constraint (10).

To derive the Contact integrator for this system, we use a linear-order
approximation in equations (8), i.e.

(14) $z_{j+1}-z_{j}=hL(x_{j},x_{j+1},z_{j},z_{j+1}),$

and the constraints are discretized in the same way as in [14]. The resulting
equations are

(15)
$\begin{split}\frac{-x_{j+1}+2x_{j}-x_{j-1}}{h^{2}}-\frac{g}{l}x_{j}-\alpha\left(\frac{x_{j}-x_{j-1}}{h}-\frac{h}{2}\frac{g}{l}x_{j}\right)+\lambda_{1}\frac{y_{j}}{m}&=0\\\
\frac{-y_{j+1}+2y_{j}-y_{j-1}}{h^{2}}-\frac{g}{l}y_{j}-\alpha\left(\frac{y_{j}-y_{j-1}}{h}-\frac{h}{2}\frac{g}{l}y_{j}\right)-\lambda_{1}\frac{x_{j}}{m}&=0\\\
-y_{j}\frac{x_{j+1}-x_{j}}{h}+x_{j}\frac{y_{j+1}-y_{j}}{h}+\Omega\sin\beta(x_{j}^{2}+y_{j}^{2})&=0,\end{split}$

where $\lambda_{1}$ is a Lagrange multiplier.

As can be noted, the discretization of the constraints in both integrators LA
and Contact are the same.

###  4.2 Falling disk

The rolling falling disk is one of the archetypal examples of nonholonomically
constrained mechanical systems, frequently used as a benchmark for numerical
methods [3, 54]. It consists of an idealized (without thickness) homogeneous
disk rolling without slipping over a fixed horizontal plane, subject to
gravity (see Figure 3). Depending on the research interest, we may also
consider external forces acting upon the disk, either dissipative forces or
control forces. To derive its Lagrangian, we will follow the formulation given
in [57].

Figure 3: Falling disk and its state variables.

The configuration space of this system can be identified as
$Q=\mathbb{R}^{2}\times SO(3)$. As generalized coordinates, we can take a
point $(X,Y)\in\mathbb{R}^{2}$ describing the orthogonal projection of the
center of the disk over the plane $\\{Z=0\\}$ (see Figure 3), and an element
of $SO(3)$ in terms of Eulerian angles $(\theta,\phi,\psi)$, defined as
follows: we consider a right-handed noninertial orthogonal reference frame
$(x,y,z)$ whose origin is attached to the center of the disk, the $x$ axis
being orthogonal to the plane of the disk and the $y$ axis remaining in the
plane of the disk and parallel to the plane $\\{Z=0\\}$. Then, $\theta$ is the
angle between the axes $z$ and $Z$, while $\phi$ is the angle between the axes
$y$ and $X$, and $\psi$ is the angle of rotation of the disk with respect to
the $x$ axis.

In the Lagrange-d’Alembert formulation, the Lagrangian for this system is
$L=T-V$, where $T$ is the kinetic energy and $V$ is the potential energy. For
the kinetic energy, we have the expression

(16)
$T=\frac{1}{2}m\hat{v}^{2}+\frac{1}{2}[I_{A}\omega_{x}^{2}+I_{T}(\omega_{y}^{2}+\omega_{z}^{2})],$

where $\hat{v}$ is the linear velocity of the center of the disk, $I_{A}$ and
$I_{T}$ are the moment of inertia with respect to $x$ axis, and the $y$ and
$z$ axes, respectively, while $\omega_{i}$, for $i\in\\{x,y,z\\}$ are the
angular velocity with respect to the $i$-th axis. In terms of the generalized
velocities associated to the Euler angles, the angular velocities are given by

(17) $\begin{split}\omega_{x}&=-\dot{\psi}+\dot{\phi}\sin\theta,\\\
\omega_{y}&=-\dot{\theta},\\\ \omega_{z}&=\dot{\phi}\cos\theta.\end{split}$

On the other hand, the potential energy is given by

$V=mgR\cos\theta,$

where $m$ is the mass of the disk, $g$ is the gravitational acceleration, and
$R$ is the radius of the disk. Hence, the complete Lagrangian is

(18)
$\begin{split}L=T(\dot{q})-V(q)&=\frac{1}{2}m[\dot{X}^{2}+\dot{Y}^{2}+R^{2}\sin^{2}\theta\dot{\theta}^{2}]\\\
&+\frac{1}{2}\left[I_{A}\left(\dot{\psi}-\dot{\phi}\sin\theta\right)^{2}+I_{T}\left(\dot{\theta}^{2}+\dot{\phi}^{2}\cos^{2}\theta\right)\right]-mgR\cos\theta.\end{split}$

Finally, the non-slipping condition imposes restrictions on the velocities
given by

(19)
$\begin{split}\dot{X}&=-R\cos\theta\sin\phi\dot{\theta}-R\sin\theta\cos\phi\dot{\phi}+R\cos\phi\dot{\psi},\\\
\dot{Y}&=R\cos\theta\cos\phi\dot{\theta}-R\sin\theta\sin\phi\dot{\phi}+R\sin\phi\dot{\psi}.\end{split}$

Let us consider that the motion of the disk undergoes a Rayleigh dissipation
and also that there is some external nonconservative force applied upon the
disk. The Herglotz description allows us to incorporate these features into
the Lagrangian of the system by modifying the expression given in (18) in the
following way

(20) $L_{H}(t,q,\dot{q},z)=L-\alpha z+F(t)q,$

where $L$ is the Lagrangian given by (18), $\alpha$ is the dissipation
parameter, and $F(t)$ is the nonconservative external force. The restriction
of the system, given by (19) are unchanged in the Herglotz formulation.

###  4.3 Numerical Experiments: Foucault Pendulum

For the Foucault pendulum, we perform two simulations, using as reference the
$4^{th}$ order Runge-Kutta-Fehlberg method, to compare the behavior of both
integrators. In both simulations, we use the configuration of the original
experiment held in the Observatory of Paris, in 1851, namely, we consider a
pendulum of mass $m=28$ kg, length $l=67$ m, and latitude $\beta=49^{\circ}$.
We consider two different values of the dissipation parameter $\alpha$.

Experiment 1: Here we consider $\alpha=1\times 10^{-3}$ for both the Contact
and the LA integrators. For this setting, both integrators display
indistinguishable behavior in terms of the trajectory (see Figure 4 and Figure
5) as well as in energy dissipation, as can be seen in part b) of Figure 6.
However, we observe that the LA integrator do have a better outcome in terms
of the error $\|q_{tj}-q_{ref}(t_{j})\|_{2}$, where $q_{t_{j}}$ is the
approximation given by the integrator and $q_{ref}(t_{j})$ is the reference
solution, both at time $t=t_{j}$. This can be seen in the part a) of Figure 6.

Figure 4: Simulation 1 for the Foucault Pendulum with initial conditions
$q(0)=(0,l/100)$, $\dot{q}(0)=(0,0)$. Dissipation parameter $\alpha=1\times
10^{-3}$. The graphics show the values of the $x$ and $y$ coordinates as
functions of time for the last $100$ seconds of simulation.

Figure 5: Simulation 1 for the Foucault Pendulum with initial conditions
$q(0)=(0,l/100)$, $\dot{q}(0)=(0,0)$. Dissipation parameter $\alpha=1\times
10^{-3}$. The graphics show the projection over the horizontal plane of the
oscillation plane in the interval $0\leq t\leq 3600$ seconds.

(a) Trajectory error

(b) Energy

Figure 6: Simulation 1 for the Foucault Pendulum with initial conditions
$q(0)=(0,l/100)$, $\dot{q}(0)=(0,0)$. Dissipation parameter $\alpha=1\times
10^{-3}$. The graphics show the errors of the trajectory and the energy for
the last $100$ seconds of simulation.

Experiment 2: In order to simulate a more realistic situation, where
dissipation is due to friction with the air, here we consider $\alpha=1\times
10^{-4}$. In this setting, the contact integrator outcome is closer to the
reference, both in terms of the error $\|q_{t_{j}}-q_{ref}(t_{j})\|$ as well
as in term of energy dissipation, as can be seen in Figure 7. More
interestingly, in this setting, the LA integrator displays anomalous behavior.
Concretely the plane of oscillation varies discontinuously at a given time
(see part b) of Figure 8. On the other hand, the contact integrator does not
show this anomaly, as can be seen in part a) of Figure 8.

Figure 7: Simulation 2 for the Foucault Pendulum with initial conditions
$q(0)=(0,l/100)$, $\dot{q}(0)=(0,0)$. Dissipation parameter $\alpha=1\times
10^{-4}$. The graphics show the errors of the trajectory and the energy for
the last $100$ seconds of simulation.

(a) Contact (1st order)

(b) LA (1st order)

(c) Reference (RKF45)

Figure 8: Foucault pendulum with initial conditions $q(0)=(0,l/100)$,
$\dot{q}(0)=(0,0)$. Dissipation parameter $\alpha=1\times 10^{-4}$. The
graphics show the trajectories followed by the pendulum on the horizontal
plane, for an oscillation time of 3600 s.

###  4.4 Numerical Experiments: Falling rolling disk

For the falling rolling disk, we perform $4$ numerical experiments, in each
case, using the MATLAB solver ode15i as the reference, with a step-size 10
times smaller than the one used for the proposed integrator. In all the four
experiments, we use the following configuration for the system: $m=5$ kg,
$R=0.5$ m, $I_{A}=\frac{1}{2}mR^{2}$ kgm2, $I_{T}=\frac{1}{4}mR^{2}$ kgm2,
$g=9.81$ m/s2. The numerical parameters were set as: step-size $h=0.1$, and
tolerance $\epsilon=1\times 10^{-6}$, both for the ode15i and for the modified
multivariate Newton-Raphson method. The four sets of numerical experiments
thus were used to explore different initial conditions and different values of
the dissipation parameter $\alpha$. In each case, we use the MATLAB function
“decic” to set consistent initial conditions for the solver ode15i.

Experiment 1: Here we consider a disk starting in a vertical position and an
external force $F^{\psi}=\frac{1}{2}$ N which forces it to roll. Under these
conditions, we perform two simulations, with dissipation parameter
$\alpha=0.005$ and $\alpha=0.1$. The outcome of the proposed integrator and
the reference are indistinguishable in these settings.

Experiment 2: For the second set of simulations we consider a disk starting in
an inclined position, precisely, with $\theta_{0}=\frac{\pi}{36}$ rad, and
initial rolling velocity $\dot{\psi}_{0}=2\pi$. Under these conditions, we
perform three simulations corresponding to $\alpha=0$, $\alpha=0.005$ and
$\alpha=0.1$. For the first two cases, the reference and the proposed
integrator have indistinguishable behaviors. For the third case, we observe
something interesting. In this case, the solver ode15i was not able to
continue the simulation after $t\approx 11.1$ s, concretely displaying the
following message.

Warning: Failure at t=1.112335e+01. Unable to meet integration tolerances
without reducing the step size below the smallest value allowed (3.951806e-14)
at time t.

On the other hand, the proposed integrator was able to continue the simulation
during the entire preset time interval. This phenomena can be observed in
Figure 9, Figure 10, and Figure 11, displaying the evolution of the five
generalized coordinates and the total energy of the system. It is worth
noticing that the warning message displayed by the solver ode15i is not
because the disk has reached a configuration corresponding to a completely
fallen disk since, as can be seen in Figure 12, for the value of $t\approx
11.1$ s, the value of $\theta_{ref}$ is far from $\pi/2$, which is the value
corresponding to a completely fallen disk.

Figure 9: Experiment 2.3 of the falling disk with initial conditions
${q(0)=(0,0,\pi/36,0,0)}$, ${\dot{q}(0)=(\pi,0,0,0,2\pi)}$. Dissipation
${\alpha=0.1}$ and ${F(t)=(0,0,0,0,0)}$. The graphics display the coordinate
functions $(X,Y,\theta,\phi,\psi)$ of the system, obtained by the contact
integrator and by ode15i, which is the reference method. For each coordinate,
the $ref$ subindex indicates the solution given by the reference method.

Figure 10: Experiment 2.3 of the falling disk with initial conditions
${q(0)=(0,0,\pi/36,0,0)}$, ${\dot{q}(0)=(\pi,0,0,0,2\pi)}$. Dissipation
${\alpha=0.1}$ and ${F(t)=(0,0,0,0,0)}$. The graphics display the velocity
functions $(\dot{X},\dot{Y},\dot{\theta},\dot{\phi},\dot{\psi})$ of the
system, obtained by the contact integrator and by ode15i which is the
reference method. For each velocity, the $ref$ subindex indicates the solution
given by the reference method. Figure 11: Experiment 2.3 of the falling disk
with initial conditions ${q(0)=(0,0,\pi/36,0,0)}$,
${\dot{q}(0)=(\pi,0,0,0,2\pi)}$. Dissipation ${\alpha=0.1}$ and
${F(t)=(0,0,0,0,0)}$. The graphic displays the energy functions obtained by
the contact integrator and by ode15i which is the reference method. Figure 12:
Experiment 2.3 of the falling disk with initial conditions
${q(0)=(0,0,\pi/36,0,0)}$, ${\dot{q}(0)=(\pi,0,0,0,2\pi)}$. Dissipation
${\alpha=0.1}$ and ${F(t)=(0,0,0,0,0)}$. The graphic displays the function of
the configuration variable $\theta$ obtained by the contact integrator and
ode15i, which is the reference method. It can be observed that the disk does
not fall completely for either of the two integrators in $t=11.1$ s, as the
disk falls completely when $\theta=\pi/2$.

Experiment 3: For this set of simulations, we consider again a disk starting
in a vertical position, with initial rolling velocity $\dot{\psi}_{0}=\pi$.
Additionally, here we consider two applied forces $F^{\psi};\frac{t}{16}$ N
and $F^{\phi}=\frac{t}{16}$ N, changing the rolling speed and the direction,
respectively. We perform simulations for three values of the dissipation
parameter, namely: $\alpha=0$, $\alpha=0.005$ and $\alpha=0.1$. The
interesting case now corresponds to $\alpha=0$, where once again, the
reference solver ode15i fails to complete the simulation, displaying the same
warning message as before for $t\approx 7.7$ s. As before, this failure is not
due to a configuration corresponding to a completely fallen disk, as can be
seen in Figure 16. On the other hand, the proposed integrator does not suffer
from this issue, as can be observed in Figure 13, Figure 14, and Figure 15.
For the other values of $\alpha$, the outcomes of both integrators are
qualitatively well-behaved.

Figure 13: Experiment 3.1 of the falling rolling disk with initial condition
$q(0)~{}=~{}(0,0,0,0,0)$, $\dot{q}(0)=(\pi/2,0,0,0,\pi)$. Dissipation
$\alpha=0$ and $F(t)=(0,0,0,t/16,t/16)$. The graphics show the time evolution
of the generalized coordinates of the system for both integrators.

Figure 14: Experiment 3.1 of the falling rolling disk with initial condition
$q(0)~{}=~{}(0,0,0,0,0)$, $\dot{q}(0)=(\pi/2,0,0,0,\pi)$. Dissipation
$\alpha=0$ and $F(t)=(0,0,0,t/16,t/16)$. The graphics show the time evolution
of the generalized velocities of the system for both integrators. Figure 15:
Experiment 3.1 of the falling rolling disk with initial condition
$q(0)~{}=~{}(0,0,0,0,0)$, $\dot{q}(0)=(\pi/2,0,0,0,\pi)$. Dissipation
$\alpha=0$ and $F(t)=(0,0,0,t/16,t/16)$. The graphic shows the energy function
for both integrators. Figure 16: Experiment 3.1 of the falling rolling disk
with initial condition $q(0)~{}=~{}(0,0,0,0,0)$,
$\dot{q}(0)=(\pi/2,0,0,0,\pi)$. Dissipation $\alpha=0$ and
$F(t)=(0,0,0,t/16,t/16)$. The graphic displays the function of the
configuration variable $\theta$ obtained by the contact integrator and ode15i,
which is the reference method. It can be observed that the disk does not fall
completely for either of the two integrators in $t\approx 7.7$ s, as the disk
falls completely when $\theta=\pi/2$.

Experiment 4: This set of simulations considers a condition under which the
disk follows a stable circular path. This evolution is possible if the
inclination angle $\theta=\theta_{0}$, the precession
$\dot{\phi}=\dot{\phi}_{0}$, and the rolling velocity
$\dot{\psi}=\dot{\psi}_{0}$ satisfy the following restriction

$(I_{T}-I_{A}-mR^{2})^{2}\sin(\theta_{0})\dot{\theta}_{0}^{2}-(I_{A}+mR^{2})\tan(\theta_{0})\dot{\theta}_{0}\dot{\psi}_{0}-mgR=0.$

We perform again three simulations, corresponding to $\alpha=0$,
$\alpha=0.005$ and $\alpha=0.1$. In these simulations, both integrators
display, qualitatively, the same outcome, but once again, for $\alpha=0.1$,
the solver ode15i was unable to complete the simulation, displaying the same
warning message as in the previous situations, while the proposed integrator
completed the simulation satisfactorily. We summarize this observations with
Figures 17 to 20.

Figure 17: Experiment 4.3 of the falling rolling disk with initial conditions
$q(0)~{}=~{}(0,0,20\pi/180,0,0)$,
$\dot{q}(0)=(\pi/2,0,0,-3\pi/10,\dot{\psi}_{0})$, with
$\dot{\psi}_{0}=((I_{T}-I_{A}-mR^{2})\sin(\theta_{0})\dot{\phi}_{0}^{2}-mgR)/((I_{A}+mR^{2})\tan(\theta_{0})\dot{\phi}_{0})$.
Dissipation $\alpha=0.1$ and $F(t)=(0,0,0,0,0)$. The graphics show the time
evolution of the generalized coordinates for both integrators.

Figure 18: Experiment 4.3 of the falling rolling disk with initial conditions
$q(0)~{}=~{}(0,0,20\pi/180,0,0)$,
$\dot{q}(0)=(\pi/2,0,0,-3\pi/10,\dot{\psi}_{0})$, with
$\dot{\psi}_{0}=((I_{T}-I_{A}-mR^{2})\sin(\theta_{0})\dot{\phi}_{0}^{2}-mgR)/((I_{A}+mR^{2})\tan(\theta_{0})\dot{\phi}_{0})$.
Dissipation $\alpha=0.1$ and $F(t)=(0,0,0,0,0)$. The graphics show the time
evolution of the generalized velocities for both integrators. Figure 19:
Experiment 4.3 of the falling rolling disk with initial conditions
$q(0)~{}=~{}(0,0,20\pi/180,0,0)$,
$\dot{q}(0)=(\pi/2,0,0,-3\pi/10,\dot{\psi}_{0})$, with
$\dot{\psi}_{0}=((I_{T}-I_{A}-mR^{2})\sin(\theta_{0})\dot{\phi}_{0}^{2}-mgR)/((I_{A}+mR^{2})\tan(\theta_{0})\dot{\phi}_{0})$.
Dissipation $\alpha=0.1$ and $F(t)=(0,0,0,0,0)$. The graphic shows the time
evolution of the energy of the system for both integrators. Figure 20:
Experiment 4.3 of the falling rolling disk with initial conditions
$q(0)~{}=~{}(0,0,20\pi/180,0,0)$,
$\dot{q}(0)=(\pi/2,0,0,-3\pi/10,\dot{\psi}_{0})$, with
$\dot{\psi}_{0}=((I_{T}-I_{A}-mR^{2})\sin(\theta_{0})\dot{\phi}_{0}^{2}-mgR)/((I_{A}+mR^{2})\tan(\theta_{0})\dot{\phi}_{0})$.
Dissipation $\alpha=0.1$ and $F(t)=(0,0,0,0,0)$. The graphic displays the
function of the configuration variable $\theta$ obtained by the contact
integrator and ode15i, which is the reference method. It can be observed that
the disk does not fall completely for either of the two integrators in
$t\approx 22.8$ s, as the disk falls completely when $\theta=\pi/2$.

##  5 Concluding remarks

In this work, we propose an integrator for nonholonomic, nonconservative
mechanical systems. To do so, we combine previous works dealing with a
derivation of a geometric integrator based on Herglotz variational principle,
as well as one incorporating nonholonomic constraints for integrators derived
from the Lagrange-d’Alembert principle. For many interesting mechanical
systems, the use of Herglotz principle allows to incorporate dissipation
phenomena directly into the contact-type Lagrangian, providing a more
intrinsic treatment of these kinds of mechanical systems.

We validate our integrator by numerical experiments performed over two
mechanical systems, namely the Foucault pendulum with dissipation and the
rolling falling disk with dissipation and a forcing term. The validation was
done by comparing the numerical outcome of our integrator with standard
methods such as the fourth-order Runge-Kutta-Fehlberg (for the Foucault
pendulum) and the ode15i solver available in Matlab (for the rolling falling
disk). Also, for the Foucault pendulum, we contrast our integrator against the
one coming from the Lagrange-d’Alembert principle, which is somehow the
“standard” variational-like approach to nonholonomic systems.

The numerical experiments show that our integrator has good qualitative
behavior, matching the outcome of the references and, in some situations,
outperforming them, at least in the qualitative sense. This is worth
mentioning since one of the main paradigms of geometric integration is to
develop numerical schemes with good qualitative behavior in long-time
simulations. One potential drawback of the proposed integrator is the
execution time which is significantly larger than the ode15i. However, we
point out that this is mainly because the ode15i is an optimized package
designed for Matlab, while the implementation of the proposed integrator makes
use of a particular subroutine implementing a modified multivariate Newton-
Raphson method, which has not been optimized for the Matlab platform.

It is also worth noticing that the anomaly displayed by the
Lagrange-d’Alembert integrator (regarding Figure 8) disappears when we refine
the time-step by a factor of 20, which certainly implies a computational cost
that might be undesirable. This raises the question of whether this anomaly is
geometric in nature or is it a purely numerical issue.

## Acknowledgments

C.E.S. acknowledges the support given by PRONII. I.O. and E.M. acknowledge
partial support given by FEEI-CONACYT-PROCIENCIA.

## References

  * [1] R. Abraham and J. Marsden, Foundations of Mechanics, American Mathematical Society, 2008.
  * [2] A. Anahory Simoes, D. Martín de Diego, M. Lainz Valcázar, and M. de León, On the geometry of discrete contact mechanics, Journal of Nonlinear Science, 31 (2021), pp. 1–30.
  * [3] A. M. Bloch, J. E. Marsden, P. E. Crouch, and J. Baillieul, Nonholonomic Mechanics and Control, Springer, 2015.
  * [4] A. Bravetti, Contact hamiltonian dynamics: the concept and its use, Entropy, 19 (2017), p. 535.
  * [5] A. Bravetti, Contact geometry and thermodynamics, International Journal of Geometric Methods in Modern Physics, 16 (2019).
  * [6] A. Bravetti, C. Hans, and D. Tapias, Contact hamiltonian mechanics, Annals of Physics, 376 (2017), pp. 17–89.
  * [7] A. Bravetti, C. Lopez-Monsalvo, and F. Nettel, Contact symmetries and hamiltonian thermodynamics, Annals of Physics, 361 (2015), pp. 377–400.
  * [8] C. M. Campos, S. Ober-Blöbaum, and E. Trélat, High order variational integrators in the optimal control of mechanical systems, arXiv preprint arXiv:1502.00325, (2015).
  * [9] J. F. Carinena and M. F. Ranada, Lagrangian systems with constraints: a geometric approach to the method of lagrange multipliers, Journal of Physics A: Mathematical and General, 26 (1993), p. 1335.
  * [10] L. Colombo, M. G. Fernández, and D. M. de Diego, Variational integrators for non-autonomous systems with applications to stabilization of multi-agent formations, arXiv preprint arXiv:2202.01471, (2022).
  * [11] L. Colombo, F. Jiménez, and D. M. de Diego, Variational integrators for underactuated mechanical control systems with symmetries, arXiv preprint arXiv:1209.6315, (2012).
  * [12] L. J. Colombo and H. G. de Marina, A variational integrator for the distance-based formation control of multi-agent systems, IFAC-PapersOnLine, 51 (2018), pp. 76–81.
  * [13] J. Cortés, Geometric, Control and Numerical Aspects of Nonholonomic Systems, Springer, 2002.
  * [14] J. Cortés and S. Martínez, Non-holonomic integrators, Nonlinearity, 14 (2001), pp. 1365–1392.
  * [15] M. de León and D. M. de Diego, On the geometry of non-holonomic lagrangian systems, Journal of Mathematical Physics, 37 (1996), pp. 3389–3414.
  * [16] M. de León, D. M. de Diego, and A. Santamaría-Merino, Discrete variational integrators and optimal control theory, Advances in Computational Mathematics, 26 (2007), pp. 251–268.
  * [17] M. de León, V. M. Jiménez, and M. L. Valcázar, Contact hamiltonian systems with nonholonomic constraints, arXiv preprint arXiv:1911.05409, (2019).
  * [18] M. de León and M. Lainz, A review on contact hamiltonian and lagrangian systems, arXiv preprint arXiv:2011.05579, (2020).
  * [19] M. de León, J. Marrero, and D. Martin De Diego, New Developments in Differential Geometry, Budapest 1996, Springer, Dordrecht, 1999, ch. Time-Dependent Mechanical Systems With Non-Linear Constraints.
  * [20] M. de León, J. C. Marrero, and D. M. de Diego, Mechanical systems with nonlinear constraints, International Journal of Theoretical Physics, 36 (1997), pp. 979–995.
  * [21] F. Demoures, F. Gay-Balmaz, S. Leyendecker, S. Ober-Blöbaum, T. S. Ratiu, and Y. Weinand, Discrete variational lie group formulation of geometrically exact beam dynamics, Numerische Mathematik, 130 (2015), pp. 73–123.
  * [22] J. Ding and Z. Pan, Higher order variational integrators for multibody system dynamics with constraints, Advances in Mechanical Engineering, 6 (2014), p. 383680.
  * [23] V. Duruisseaux, J. Schmitt, and M. Leok, Adaptive hamiltonian variational integrators and applications to symplectic accelerated optimization, SIAM Journal on Scientific Computing, 43 (2021), pp. A2949–A2980, https://doi.org/10.1137/20M1383835.
  * [24] B. Georgieva, The Variational Principle of Herglotz and Reletad Results, Geometry, Integrability and Quantization, 12 (2011), pp. 214–225.
  * [25] E. Hairer, C. Lubich, and G. Wanner, Geometric numerical integration: structure-preserving algorithms for ordinary differential equations, vol. 31, Springer Science & Business Media, 2006.
  * [26] J. Hall and M. Leok, Spectral variational integrators, Numerische Mathematik, 130 (2015), pp. 681–740.
  * [27] D. Holm, Geometric Mechanics and Symmetry: from finite o infinite dimensions, Oxford University Press, 2009.
  * [28] D. Holm, Geometric Mechanics - Part I: Dynamics and Symmetry, World Scientific Publishing Company, 2011.
  * [29] D. Holm, Geometric Mechanics - Part II: Rotating, Translating and Rolling, World Scientific Publishing Company, 2011.
  * [30] I. I. Hussein, M. Leok, A. K. Sanyal, and A. M. Bloch, A discrete variational integrator for optimal control problems on so (3), in Proceedings of the 45th IEEE Conference on Decision and Control, IEEE, 2006, pp. 6636–6641.
  * [31] O. Krupková, Mechanical systems with nonholonomic constraints, Journal of Mathematical Physics, 38 (1997), pp. 5098–5126.
  * [32] O. Krupková, The nonholonomic variational principle, Journal of Physics A: Mathematical and Theoretical, 42 (2009), p. 185201.
  * [33] O. Krupková, Geometric mechanics on nonholonomic submanifolds, Communications in Mathematics, 18 (2010), pp. 51–77.
  * [34] M. J. Lazo, J. Paiva, J. T. Amaral, and G. S. Frederico, An action principle for action-dependent lagrangians: Toward an action principle to non-conservative systems, Journal of Mathematical Physics, 59 (2018), p. 032902.
  * [35] T. Lee, M. Leok, and N. H. McClamroch, Lie group variaional integrator for the full body problem in orbital mechanics, Celestial Mechanics and Dynamical Astronomy, 98 (2007), pp. 121–144.
  * [36] T. Lee, M. Leok, and N. H. McClamroch, Lie group variational integrators for the full body problem in orbital mechanics, Celestial Mechanics and Dynamical Astronomy, 98 (2007), pp. 121–144.
  * [37] M. Leok, An overview of lie group variational integrators and their applications to optimal control, in International conference on scientific computation and differential equations, The French National Institute for Research in Computer Science and Control, 2007, p. 1.
  * [38] M. Leok, Variational discretizations of gauge field theories using group-equivariant interpolation, Foundations of Computational Mathematics, 19 (2019), pp. 965–989.
  * [39] M. Leok and T. Shingel, General techniques for constructing variational integrators, Frontiers of Mathematics in China, 7 (2012), pp. 273–303.
  * [40] A. Lew, J. E. Marsden, and M. Ortiz, Variational time integrators, International Journal for Numerical Methods in Engineering, 60 (2004), pp. 153–212.
  * [41] A. Lew, J. E. Marsden, M. Ortiz, and M. West, Asynchronous variational integrators, Archive for Rational Mechanics and Analysis, 167 (2003), pp. 85–146.
  * [42] A. Lew, J. E. Marsden, M. Ortiz, and M. West, An overview of variational integrators, International Center for Numerical Methods in Engineering (CIMNE), 2004.
  * [43] A. J. Lew and P. Mata A, A Brief Introduction to Variational Integrators, Springer International Publishing, Cham, 2016, pp. 201–291.
  * [44] A. Lewis, Lagrangian mechanics, dynamics, and control, Preprint available online at http://penelope. mast. queensu. ca/~ andrew, 237 (2003).
  * [45] S. Leyendecker, J. E. Marsden, and M. Ortiz, Variational integrators for constrained dynamical systems, ZAMM-Journal of Applied Mathematics and Mechanics/Zeitschrift für Angewandte Mathematik und Mechanik: Applied Mathematics and Mechanics, 88 (2008), pp. 677–708.
  * [46] D. J. Limebeer, S. Ober-Blöbaum, and F. H. Farshi, Variational integrators for dissipative systems, IEEE Transactions on Automatic Control, 65 (2020), pp. 1381–1396, https://doi.org/10.1109/TAC.2020.2965059.
  * [47] A. López-Gordón, L. Colombo, and M. de León, Nonsmooth herglotz variational principle, arXiv preprint arXiv:2208.02033, (2022).
  * [48] S. Man, Q. Gao, and W. Zhong, Variational integrators in holonomic mechanics, Mathematics, 8 (2020), https://www.mdpi.com/2227-7390/8/8/1358.
  * [49] Z. Manchester, N. Doshi, R. J. Wood, and S. Kuindersma, Contact-implicit trajectory optimization using variational integrators, The International Journal of Robotics Research, 38 (2019), pp. 1463–1476, https://doi.org/10.1177/0278364919849235.
  * [50] J. E. Marsden, S. Pekarsky, S. Shkoller, and M. West, Variational methods, multisymplectic geometry and continuum mechanics, Journal of Geometry and Physics, 38 (2001), pp. 253–284.
  * [51] J. E. Marsden and T. S. Ratiu, Introduction to Mechanics and Symmetry: a basic exposition of classical mechanical systems, vol. 17, Springer-Verlag New York, 2 ed., 1999.
  * [52] J. E. Marsden and J. M. Wendlandt, Mechanical systems with symmetry, variational principles, and integration algorithms, in Current and future directions in applied mathematics, Springer, 1997, pp. 219–261.
  * [53] J. E. Marsden and M. West, Discrete mechanics and variational integrators, Acta Numerica, 10 (2001), pp. 357–514.
  * [54] K. Modin and O. Verdier, What makes nonholonomic integrators work?, Numerische Mathematik, 145 (2020), pp. 405–435.
  * [55] S. Ober-Blöbaum, Galerkin variational integrators and modified symplectic Runge-Kutta methods, IMA Journal of Numerical Analysis, 37 (2016), pp. 375–406.
  * [56] S. Ober-Blöbaum, M. Tao, M. Cheng, H. Owhadi, and J. E. Marsden, Variational integrators for electric circuits, Journal of Computational Physics, 242 (2013), pp. 498–530.
  * [57] P. C. Paris and L. Zhang, A disk rolling on a horizontal surface without slip, Mathematical and computer modelling, 36 (2002), pp. 855–860.
  * [58] M. F. Rañada, Time-dependent lagrangian systems: A geometric approach to the theory of systems with constraints, Journal of Mathematical Physics, 35 (1994).
  * [59] X. Shen and M. Leok, Lie group variational integrators for rigid body problems using quaternions, arXiv preprint arXiv:1705.04404, (2017).
  * [60] A. Stern, Y. Tong, M. Desbrun, and J. E. Marsden, Geometric computational electrodynamics with variational integrators and discrete differential forms, in Geometry, mechanics, and dynamics, Springer, 2015, pp. 437–475.
  * [61] B. Tran and M. Leok, Multisymplectic hamiltonian variational integrators, International Journal of Computer Mathematics, 99 (2022), pp. 113–157.
  * [62] M. Vermeeren, A. Bravetti, and M. Seri, Contact variational integrators, Journal of Physics A: Mathematical and Theoretical, 52 (2019), p. 445206.
  * [63] J. M. Wendlandt and J. E. Marsden, Mechanical integrators derived from a discrete variational principle, Physica D: Nonlinear Phenomena, 106 (1997), pp. 223–246.
  * [64] Y. Zhang and X. Tian, Conservation laws of nonconservative nonholonomic system based on Herglotz variational problem, Physics Letters A, 383 (2019), pp. 691–696.

##  Appendix A Discrete equations of the contact integrator with constraints
for the falling disk

The following equations enumerated from (21) to (25) correspond to the
resulting equations of the integrator based on the discretization of the
Herglotz variational principle with constraints that we propose in this work,
applied to the problem of the falling disk described in Section 4.2.

(21)
$\begin{split}&R\cos(\phi_{j})\left[\frac{F^{X}_{j}}{2}-\frac{1}{\frac{\alpha
h}{2}+1}\left(\frac{F^{X}_{j}}{2}-\frac{m(X_{j-1}-X_{j})}{h^{2}}\right)\left(\frac{\alpha
h}{2}-1\right)+\frac{m(X_{j}-X_{j+1})}{h^{2}}\right]\\\
&+R\sin(\phi_{j})\left[\frac{F^{Y}_{j}}{2}-\frac{1}{\frac{\alpha
h}{2}+1}\left(\frac{F^{Y}_{j}}{2}-\frac{m(Y_{j-1}-Y_{j})}{h^{2}}\right)\left(\frac{\alpha
h}{2}-1\right)+\frac{m(Y_{j}-Y_{j+1})}{h^{2}}\right]\\\
&-\frac{1}{\frac{\alpha h}{2}+1}\left(\frac{\alpha
h}{2}-1\right)\left[\frac{F^{\psi}_{j}}{2}-\frac{1}{h}I_{A}\left(\frac{\psi_{j-1}-\psi_{j}}{h}-\sin\left(\frac{\theta_{j-1}+\theta_{j}}{2}\right)\frac{\phi_{j-1}-\phi_{j}}{h}\right)\right]\\\
&+\frac{1}{h}I_{A}\left(\frac{\psi_{j}-\psi_{j+1}}{h}-\sin\left(\frac{\theta_{j}+\theta_{j+1}}{2}\right)\frac{\phi_{j}-\phi_{j+1}}{h}\right)+\frac{F^{\psi}_{j}}{2}=0.\end{split}$

(22)
$\begin{split}&R\cos\left(\frac{\phi_{j}+\phi_{j+1}}{2}\right)\frac{\psi_{j}-\psi_{j+1}}{h}-\frac{X_{j}-X_{j+1}}{h}\\\
&-R\cos\left(\frac{\phi_{j}+\phi_{j+1}}{2}\right)\sin\left(\frac{\theta_{j}+\theta_{j+1}}{2}\right)\frac{\phi_{j}-\phi_{j+1}}{h}\\\
&-R\cos\left(\frac{\theta_{j}+\theta_{j+1}}{2}\right)\sin\left(\frac{\phi_{j}+\phi_{j+1}}{2}\right)\frac{\theta_{j}-\theta_{j+1}}{h}=0.\end{split}$

(23)
$\begin{split}&R\sin\left(\frac{\phi_{j}+\phi_{j+1}}{2}\right)\frac{\psi_{j}-\psi_{j+1}}{h}-\frac{Y_{j}-Y_{j+1}}{h}\\\
&-R\sin\left(\frac{\phi_{j}+\phi_{j+1}}{2}\right)\sin\left(\frac{\theta_{j}+\theta_{j+1}}{2}\right)\frac{\phi_{j}-\phi_{j+1}}{h}\\\
&+R\cos\left(\frac{\phi_{j}+\phi_{j+1}}{2}\right)\cos\left(\frac{\theta_{j}+\theta_{j+1}}{2}\right)\frac{\theta_{j}-\theta_{j+1}}{h}=0.\end{split}$

(24)
$\begin{split}&\frac{F^{\theta}_{j}}{m}+\frac{I_{T}}{m}\left(\frac{2(\theta_{j}-\theta_{j+1})}{h^{2}}-\frac{\sin(\theta_{j}+\theta_{j+1})}{2}\frac{(\phi_{j}-\phi_{j+1})^{2}}{h^{2}}\right)+Rg\sin(\theta_{j})\\\
&+R^{2}\left(\sin^{2}\left(\frac{\theta_{j}+\theta_{j+1}}{2}\right)\frac{2(\theta_{j}-\theta_{j+1})}{h^{2}}+\frac{\sin(\theta_{j}+\theta_{j+1})}{2}\frac{(\theta_{j}-\theta_{j+1})^{2}}{h^{2}}\right)\\\
&+\frac{\alpha h-2}{\alpha
h+2}\Bigg{\\{}\frac{I_{T}}{m}\left[\frac{2(\theta_{j-1}-\theta_{j})}{h^{2}}+\frac{\sin(\theta_{j-1}+\theta_{j})}{2}\frac{(\phi_{j-1}-\phi_{j})^{2}}{h^{2}}\right]-\frac{F^{\theta}_{j}}{m}-Rg\sin(\theta_{j})\\\
&+R^{2}\left[\sin^{2}\left(\frac{\theta_{j-1}+\theta_{j}}{2}\right)\frac{2(\theta_{j-1}-\theta_{j})}{h^{2}}-\frac{\sin(\theta_{j-1}+\theta_{j})}{2}\frac{(\theta_{j-1}-\theta_{j})^{2}}{h^{2}}\right]\\\
&+\frac{I_{A}(\phi_{j-1}-\phi_{j})}{m}\cos\left(\frac{\theta_{j-1}+\theta_{j}}{2}\right)\left[\frac{\psi_{j-1}-\psi_{j}}{h^{2}}-\sin\left(\frac{\theta_{j-1}+\theta_{j}}{2}\right)\frac{\phi_{j-1}-\phi_{j}}{h^{2}}\right]\Bigg{\\}}\\\
&+2R\cos(\phi_{j})\cos(\theta_{j})\left[\frac{F^{Y}_{j}}{2m}-\frac{\alpha
h-2}{\alpha
h+2}\left(\frac{F^{Y}_{j}}{2m}-\frac{Y_{j-1}-Y_{j}}{h^{2}}\right)+\frac{Y_{j}-Y_{j+1}}{h^{2}}\right]\\\
&-2R\cos(\theta_{j})\sin(\phi_{j})\left[\frac{F^{X}_{j}}{2m}-\frac{\alpha
h-2}{\alpha
h+2}\left(\frac{F^{X}_{j}}{2m}-\frac{X_{j-1}-X_{j}}{h^{2}}\right)+\frac{X_{j}-X_{j+1}}{h^{2}}\right]\\\
&-\frac{I_{A}(\phi_{j}-\phi_{j+1})}{m}\cos\left(\frac{\theta_{j}+\theta_{j+1}}{2}\right)\left[\frac{\psi_{j}-\psi_{j+1}}{h^{2}}-\sin\left(\frac{\theta_{j}+\theta_{j+1}}{2}\right)\frac{\phi_{j}-\phi_{j+1}}{h^{2}}\right]=0.\end{split}$

(25)
$\begin{split}&\frac{F^{\phi}_{j}}{2}+I_{T}\cos^{2}\left(\frac{\theta_{j}+\theta_{j+1}}{2}\right)\frac{\phi_{j}-\phi_{j+1}}{h^{2}}\\\
&-\frac{\alpha h-2}{\alpha
h+2}\Bigg{\\{}\frac{F^{\phi}_{j}}{2}-I_{T}\cos^{2}\left(\frac{\theta_{j-1}+\theta_{j}}{2}\right)\frac{\phi_{j-1}-\phi_{j}}{h^{2}}\\\
&+I_{A}\sin\left(\frac{\theta_{j-1}+\theta_{j}}{2}\right)\left[\frac{\psi_{j-1}-\psi_{j}}{h^{2}}-\sin\left(\frac{\theta_{j-1}+\theta_{j}}{2}\right)\frac{\phi_{j-1}-\phi_{j}}{h^{2}}\right]\Bigg{\\}}\\\
&-R\cos(\phi_{j})\sin(\theta_{j})\left[\frac{F^{X}_{j}}{2}-\frac{\alpha
h-2}{\alpha
h+2}\left(\frac{F^{X}_{j}}{2}-\frac{m(X_{j-1}-X_{j})}{h^{2}}\right)+\frac{m(X_{j}-X_{j+1})}{h^{2}}\right]\\\
&-R\sin(\phi_{j})\sin(\theta_{j})\left[\frac{F^{Y}_{j}}{2}-\frac{\alpha
h-2}{\alpha
h+2}\left(\frac{F^{Y}_{j}}{2}-\frac{m(Y_{j-1}-Y_{j})}{h^{2}}\right)+\frac{m(Y_{j}-Y_{j+1})}{h^{2}}\right]\\\
&-I_{A}\sin\left(\frac{\theta_{j}+\theta_{j+1}}{2}\right)\left[\frac{\psi_{j}-\psi_{j+1}}{h^{2}}-\sin\left(\frac{\theta_{j}+\theta_{j+1}}{2}\right)\frac{\phi_{j}-\phi_{j+1}}{h^{2}}\right]=0.\end{split}$

E-mail address:
E-mail address:
E-mail address:


# Deep scattering network for speech emotion recognition

Premjeet Singh1, Goutam Saha1, Md Sahidullah2 1Dept of Electronics and ECE,
Indian Institute of Technology Kharagpur, Kharagpur, India
2Université de Lorraine, CNRS, Inria, LORIA, F-54000, Nancy, France
<EMAIL_ADDRESS><EMAIL_ADDRESS><EMAIL_ADDRESS>
###### Abstract

This paper introduces scattering transform for speech emotion recognition
(SER). Scattering transform generates feature representations which remain
stable to deformations and shifting in time and frequency without much loss of
information. In speech, the emotion cues are spread across time and localised
in frequency. The time and frequency invariance characteristic of scattering
coefficients provides a representation robust against emotion irrelevant
variations e.g., different speakers, language, gender etc. while preserving
the variations caused by emotion cues. Hence, such a representation captures
the emotion information more efficiently from speech. We perform experiments
to compare scattering coefficients with standard mel-frequency cepstral
coefficients (MFCCs) over different databases. It is observed that frequency
scattering performs better than time-domain scattering and MFCCs. We also
investigate layer-wise scattering coefficients to analyse the importance of
time shift and deformation stable scalogram and modulation spectrum
coefficients for SER. We observe that layer-wise coefficients taken
independently also perform better than MFCCs.

###### Index Terms:

Deep convolutional networks, Deep scattering transform, EmoDB, IEMOCAP,
RAVDESS, Shift invariance, Speech emotion recognition

##  I Introduction

The introduction of deep learning has caused rapid development in different
speech signal processing tasks. One such domain is _speech emotion
recognition_ (SER) which finds applications in health-care systems, sentiment
analysis, and many other human-computer interaction applications [1, 2, 3].
Even after two decades of research, SER is still considered as a challenging
task [2]. One of the major challenges is the interpersonal and intrapersonal
variability in emotion expression, e.g., different speaking styles, cultural
background, language, context, speaker’s mood etc [1, 3]. This causes SER
systems to face difficulty in generalization across different emotion speech
samples, especially when the train and test data samples belong to different
domains [4]. This downside motivates the search for a robust feature that can
extract emotion-specific information irrespective of different variabilities.

The deep learning based SER approaches generally addresses the issue of
variabilities by using deep convolutional neural networks (CNN), which are
considered efficient feature extractors in terms of invariances learned
against signal shifts [5, 6]. However, the complexity and _lack of
explainability_ of deep CNNs come as an added disadvantage to automatic
feature extraction ability [7].

To obtain an end-to-end framework, some SER studies use 1D CNNs to learn the
relevant features directly from raw speech [8]. Such studies observe that deep
networks learn representations similar to handcrafted features [9]. However,
due to inherent instability and requirement of large data samples for
training, suitability of deep neural networks becomes uncertain [10].

Following similar lines, we propose the use of _scattering transform_ for SER.
Scattering transform was introduced in [11, 12] as a deep convolutional
network involving convolution with _predefined kernel_ instead of automatic
kernel learning. It also introduces stability to both temporal shift and
deformations in the feature representation of signals. As emotion cues spread
temporally in speech, such characteristics of scattering transform should
provide a representation robust to both time-shifting of cues and emotion
irrelevant temporal variations. In [13], authors compute scattering
coefficients across the log-frequency domain for frequency transposition
invariance to obtain speaker-independent representation for speech
recognition. Such characteristics may also help in learning speaker-
independent emotion cues.

Scattering transform is used in both 1-D and 2-D data processing. Regarding
1-D signals, authors in [14] introduce joint time-frequency scattering that
incorporates multiscale frequency energy distribution over time-invariant
representation for various audio classification frameworks. Authors in [15]
use two-layer scattering coefficients with CNN layers to obtain a stable
descriptor of speaker information from raw speech. In [16], authors compute
different moments of scattering coefficients layers and found that such
moments contain enough information for decent voice synthesis quality. Authors
in [17] used scattering transform for urban sound classification and report a
marginal performance improvement but with reduced training data. They
attribute their finding to better background information learning ability of
temporal modulations. Along similar lines, [18] also applies to scatter
coefficients for environmental sound classification. In [19] a joint time-
frequency based scattering representation is used to analyze the timbral
similarity between different acoustic instruments. Our work is the first to
apply scattering transform for SER. We conduct experiments on three different
SER datasets using time and frequency scattering coefficients and achieve
noticeable improvement over standard mel-frequency cepstral coefficients
(MFCCs). The main contributions of this work are summarized below:

  1. 1.

Optimization of scattering transform parameters and analysing its performance
for SER,

  2. 2.

Analysis of time-domain and frequency-domain scattering coefficients over
different databases, and

  3. 3.

Layer-wise analysis of scattering coefficients for SER.

##  II Scattering Transform

For 1D signals, the idea behind scattering transform is to obtain a robust
feature that remains invariant to the location of information cues
(translation invariance) and time warping (diffeomorphism) of information
across different instances of the signal.

###  II-A Notion of Lipschitz continuity and information loss in MFCC

Let the signal $x(t)$ be translated by $c$. Then we have, $x_{c}(t)=x(t-c)$.
Taking its Fourier transform, we obtain, $X_{c}(\omega)=e^{-j\omega
c}X(\omega)$. Taking modulus on both sides, we obtain,
$|X_{c}(\omega)|=|X(\omega)|$. In short time Fourier transform (STFT),
translation $c$ is localized to the time frame duration $T$. If $|c|\ll T$,
the STFT representation is already invariant to such translation. However, a
robust feature representation also requires stability to time deformations
appearing in the signal. Let $x(t)$ is warped (deformed) in time by a factor
$\tau$, i.e., $x_{\tau}(t)=x(t-\tau(t))$. A transformation $\phi$ of $x(t)$ is
said to be stable to deformation $\tau$ if, $||\phi(x)-\phi(x_{\tau})||\leq
C\underset{t}{sup}|\tau^{{}^{\prime}}(t)|||x||$, i.e., the distance between
the non-deformed ($\phi(x)$) and deformed ($\phi(x_{\tau})$) feature in
transformation space $\phi$ should be less than a factor $C$ of maximum
amplitude/size of deformation ($\underset{t}{sup}|\tau^{{}^{\prime}}(t)|$).
This is also known as Lipschitz continuity condition. Applying this to Fourier
transform, we have, $|||X(\omega)|-|X_{\tau}(\omega)|||\leq C\epsilon||x||$,
where, time deformation $\tau(t)=\epsilon t$. Now the Fourier transform of
$x_{\tau}(t)=x(t-\epsilon t)$ is given as
$X_{\tau}(\omega)=\frac{1}{1-\epsilon}X(\frac{\omega}{1-\epsilon})$. Hence,
the frequency components at $\omega$ are shifted to $\frac{1}{1-\omega}$. This
effect is more prominent at high frequencies. Hence, modulus of Fourier
transform is not Lipschitz continuous. In MFCC, filters with non-linearly
varying bandwidths are applied over the STFT. As output of every mel-filter is
averaged to get a single coefficient, filter application can be viewed as
averaging performed over frequency domain which counters the effect of
instability. The mel-filters applied around high frequency regions have higher
bandwidth which reduces the effect of deformation instability. However, this
averaging in frequency also leads to information loss, especially at higher
frequencies. Now, consider $x_{t}(u)$ to be the signal frame at time $t$. Then
$x_{t}(u)=x(u)\phi(u-t)$ where $\phi$ is a window of length $T$. The
application of mel-filters over Fourier transform of $x_{t}(u)$, can be
written in frequency as,

$Mx(t,\lambda)=\frac{1}{2\pi}\int|x_{t}(\omega)|^{2}|\psi_{\lambda}(\omega)|^{2}d\omega$
(1)

where $x_{t}(\omega)$ is the frequency response of $x_{t}(u)$ and
$\psi_{\lambda}(\omega)$ is the mel-filter with support $\lambda$. Since
multiplication in frequency becomes convolution in time, Eq. 1 becomes,

$\displaystyle\begin{split}Mx(t,\lambda)=\int|x_{t}*\psi_{\lambda}(v)|^{2}dv\\\
=\int\left|\int x(u)\phi(u-t)\psi_{\lambda}(v-u)du\right|^{2}dv\end{split}$
(2)

If $T$ is much greater than the filter support $\lambda$ in time, $\phi(t)$ is
constant over $\lambda$. Hence we can write,
$\phi(u-t)\psi_{\lambda}(v-u)\approx\phi(v-t)\psi_{\lambda}(v-u)$, which
gives,

$\displaystyle\begin{split}M_{x}(t,\lambda)\approx\int\left|\int
x(u)\psi_{\lambda}(v-u)du\right|^{2}|\phi(v-t)|^{2}dv\\\
=|x*\psi_{\lambda}|^{2}*|\phi|^{2}(t)\end{split}$ (3)

Eq. 3 shows that averaging in frequency domain is equivalent to time-domain
averaging of $|x*\psi_{\lambda}|^{2}$ with window duration $T$ and becomes the
basis for scattering transform. Hence, $M_{x}(t,\lambda)$ remains invariant to
time shifts smaller than $T$. In MFCC, $T$ generally equals $20$ms which
limits the capturing of long time scale information. If $T$ is increased, loss
of high frequency information takes place due to averaging [13].

###  II-B Layer-wise scattering coefficients

To prevent the information loss appearing in MFCC while, at the same time,
capturing high frequency information, scattering transform computes second
layer coefficients or the modulation spectrogram of signal. This includes
computing the frequency response of the time-series of scalogram frequency
bins. This is followed by averaging of both scalogram and modulation spectrum
coefficients for stability against deformation. Scattering transform uses
constant quality factor wavelet filters over complete frequency range to
compute scalogram and modulation spectrogram. If $Q$ defines the number of
wavelets per octave, the center frequencies of wavelets becomes
$\lambda=2^{\frac{k}{Q}}$, where $k$ is an integer and the bandwidth is
proportional to $1/Q$. Every wavelet has a support of bandwidth $\lambda/Q$
around center frequency $\lambda$. The same wavelet in time provide a spread
of $\frac{2\pi Q}{\lambda}$. To make sure that this spread is less than $T$
(because $T$ is the time window over which time averaging will be computed for
stable feature representation) wavelets are only defined for
$\lambda\geq\frac{2\pi Q}{T}$. For frequencies less than $\frac{2\pi Q}{T}$
the frequency interval is covered with wavelets of equal bandwidth given by
$\frac{2\pi}{T}$, hence covering complete frequency range. The mel-filter
operation and filtering due to time window ($T$) in MFCC is replaced with
wavelet-based filters in scattering transform.

An issue with using wavelets is that they ”commute with translation” [20].
However, invariance to translation can generally be introduced by averaging.
Following this, $\int x(t)*\psi_{\lambda}(t)dt$ should be translation
invariant. But due to symmetric nature of wavelet, $\int
x(t)*\psi_{\lambda}(t)dt=0$. Hence, a modulus non-linear operator is applied
to prevent this vanishing of integral, i.e. $\int|x(t)*\psi_{\lambda}(t)|dt$.
The reason behind choosing modulus non-linearity is because modulus operation
is contractive and preserves the norm of the vector to which it is applied
[14]. The modulus of Fourier transform coefficients are then low-pass filtered
with filter $\phi$ to obtain feature representation robust against time-
warping induced deformation instability. These features are then invariant to
translations smaller than $2^{J}$, where $J$ is the scale of the low-pass
wavelet filter ($\phi$). Therefore, any $m$th order scattering coefficients
are given as,

$S_{J_{m}}x(u)=U_{m}x*\phi_{2^{J}}(u)=\int U_{m}x(v)\phi_{2^{J}}(u-v)dv$ (4)

where,

$U_{m}x=U[\lambda_{m}]...U[\lambda_{2}]U[\lambda_{1}]x=|||x*\psi_{\lambda
1}|*\psi_{\lambda 2}|...|*\psi_{\lambda m}|$

To obtain invariance to frequency transpositions, wavelet transform is first
computed over the log-frequency axis (e.g., $log({\lambda_{1}}$)) of time-
domain scattering coefficients followed by wavelet averaging. This is similar
to the DCT operation in MFCC [13]. Such averaging introduces invariance to
shifts in frequency defined by the scale of the averaging wavelet. The
frequency scattering coefficients are cascaded with time scattering
coefficients to generate a feature representation of the signal.

Figure 1: Graphical representation of scattering transform. Here, $x$ is the
input signal and $\phi$ is the low-pass filter used for stability to
deformations. $\psi_{1}$ and $\psi_{2}$ are wavelet filter banks corresponding
to $1$st and $2$nd layer of scattering transform. The figure describes only
$2$ layer decomposition of signal $x$.

##  III Experimental Details

###  III-A Dataset description

Three different speech corpora used in the experiments are described in Table
I. The speech samples are downsampled at 16 kHz sampling frequency when
required. EmoDB database is a German-language database whereas, both IEMOCAP
and RAVDESS contain speech samples in the English language. For IEMOCAP, we
use only four emotions (Happy, Angry, Sad, and Neutral) following other works
with this dataset [4, 21].

TABLE I: Summary of the speech corpora used in the experiments. (F=Female,
M=Male)

Databases Speakers Emotions Berlin Emotion Database (EmoDB) [22] 10 (5 F, 5 M)
7 (Anger, Sad, Boredom, Fear, Happy, Disgust and Neutral) Ryerson Audio-Visual
Database of Emotional Speech and Song (RAVDESS) [23] 24 (12 F, 12 M) 8 (Calm,
Happy, Sad, Angry, Neutral, Fearful, Surprise, and Disgust) Interactive
Emotional Dyadic Motion Capture Database (IEMOCAP) [24] 10 (5 F, 5 M) 4
(Happy, Angry, Sad and Neutral)

###  III-B Experimental evaluation & Methodology

We first optimize the scattering transform parameters over the EmoDB database.
The optimized parameter set is then used for SER evaluation over other
databases. We perform only two-level decomposition as higher level
coefficients contain a very small fraction of signal energy [14]. We optimize
the number of _wavelets per octave_ ($Q$) and _maximal wavelet length_ or
_averaging scale_ ($T$) while using an empirically chosen fixed value of input
signal duration ($N$). For frequency scattering, we select the maximal wavelet
length as 32. We choose this to obtain invariance over five octaves, found
optimum for EmoDB database. Also, following [13], for the first layer time-
scattering coefficients, we compute only the frequency-domain wavelet
decomposition. The required optimum frequency averaging is computed and
applied by the classifier itself [13]. We use Morlet wavelet for both time and
frequency scattering coefficients computation [13]. For implementation, we use
the ScatNet toolkit111https://github.com/scatnet/scatnet.

For classification purposes, we use a radial basis function (RBF) kernel-based
_support vector machine (SVM)_ , as the feature representations obtained are
utterance-level representations. We use the leave-one-speaker-out (LOSO)
cross-validation strategy and report average performance obtained over
different train-test pairs. For every experiment, we keep speech samples of
one speaker in validation, one in test and remaining in training set. We
optimize the hyperparameters of the SVM classifier on the validation set.

For comparison, we also compute MFCCs for speech utterances of duration $N$
and evaluate the performance with the same LOSO cross-validation strategy. We
choose the MFCC feature because of its characteristic of mimicking human sound
perception and its versatility in various speech processing domains. For
performance analysis, we use both _accuracy_ and _unweighted average recall_
(UAR) performance metrics. The UAR is defined as,

$\mathrm{UAR}=\frac{1}{K}\sum_{i=1}^{K}\frac{A_{ii}}{\sum_{j=1}^{K}A_{ij}}.$
(5)

Here $A$ is the contingency matrix, $A_{ij}$ is the number of samples in class
$i$ classified into class $j$ and $K$ is the number of classes. UAR is better
suited as a performance metric in class imbalance situations as compared to
standard accuracy [25].

Figure 2: Comparison of various time-domain scattering transform parameters
over EmoDB database. $Q$ is number of wavelets per octaves and $T$ is the
averaging scale in samples. Improved performance over lower values of $Q$
shows the importance of coarse level frequency-domain information for SER.

Figure 3: Confusion Matrices of time-domain scattering coefficients and MFCC
features over EmoDB database. Although scattering coefficients show higher
classification rate for every emotion class as compared to MFCC, confusion
across emotion classes of similar arousal characteristics, such as, Happy &
Angry, Neutral & Boredom and Fear & Happy can be observed.

##  IV Results

###  IV-A Parameter optimization over EmoDB

Figure 2 shows the results obtained by varying $Q$ and $T$ parameters of
scattering transform over EmoDB database. The optimum value of signal duration
($N$) is empirically at $51000$ samples (or $3.18$ seconds at $16$ kHz
sampling rate). We start with the lowest value of $T$ to be $4096$ samples,
which correspond to $256$ ms at $16$KHz, a duration considered optimum to
convey sufficient emotion cues [26]. $T$ is then increased by a factor of $2$
until its value is below $N$. Interestingly, $T=4096$ shows poor performance,
whereas, $T=32768$ shows better performance over different values of $Q$. In
terms of parameter $Q$, we observed a general increase in performance for
lower values of wavelets per octave. This indicates the importance of coarse-
level frequency domain information for SER. This result is also comparable to
the experiments done using constant-Q transform (CQT) for SER [27]. We select
$T=16384$ and $Q=5$ as the optimum value of parameters for further
experiments. Figure 3 shows the confusion matrices for scattering coefficients
and MFCCs over EmoDB database. Scattering coefficients can better classify
speech samples of different emotion classes in EmoDB as compared to MFCC.

###  IV-B Evaluation on other datasets

Table II shows the performance obtained with the optimised scattering
parameters (ScatNet) over different databases and its comparison with standard
MFCC. We compute 13 MFCCs over window length of $20$ ms, $10$ ms hop and $512$
frequency bins over same utterance duration $N$ ($3.18$ seconds). Mean and
standard deviation of MFCCs is computed over all the frames to obtain a vector
representation for every utterance. The scattering coefficients are observed
to outperform MFCCs over every database. We also compare performance with
frequency domain scattering transform (F-ScatNet) with similar optimized
parameters ($Q$, $T$ and $N$). The frequency transposition invariance obtained
from frequency domain scattering introduces speaker invariance hence improving
the performance. However, such improvement is not very prominent in the EmoDB
database. One probable reason for such observation could be increased
redundancy in frequency scattering coefficients because of small database size
and low number of speakers.

TABLE II: Performance comparison between frequency scattering (F-ScatNet), time-domain scattering (ScatNet) and MFCC over different SER databases. Given values are in percentages. Database | F-ScatNet | ScatNet | MFCC
---|---|---|---
| Acc. | UAR | Acc. | UAR | Acc. | UAR
EmoDB | 74.59 | 70.27 | 74.40 | 71.30 | 58.39 | 54.03
RAVDESS | 51.81 | 50.52 | 50.00 | 48.50 | 36.74 | 34.77
IEMOCAP | 61.55 | 51.00 | 60.41 | 50.40 | 55.54 | 47.19

###  IV-C Layer-wise analysis

To analyse emotion relevance of different layers of scattering transform, we
compare SER performance by using the first and second layer coefficients
separately in Table III. The first layer coefficients are averaged scalogram
coefficients similar to MFCC. However, the averaging provided by filter $\phi$
introduces time shift invariance which results in better performance than
MFCC. The second layer coefficients further perform better than first layer
coefficients showing the relevance of modulation spectrum coefficients from
SER perspective.

TABLE III: SER performance comparison with first and second layer time-domain scattering coefficients taken independently. Given values are in percentages. Database | First layer | Second layer
---|---|---
| Accuracy | UAR | Accuracy | UAR
EmoDB | 62.18 | 57.94 | 68.64 | 61.22
RAVDESS | 38.06 | 36.72 | 48.61 | 47.20
IEMOCAP | 56.48 | 47.30 | 58.93 | 48.80

##  V Discussion & Conclusion

We observe that the scattering transform coefficients, with optimized
parameters, provide feature representations that are better suited for SER.
The optimized wavelet averaging scale provides sufficient invariance and
stability to irrelevant temporal variations to capture the emotion cues in the
time domain. The low value of wavelets per octave generates coarse-level
frequency domain information improving SER performance. Invariance to
frequency transposition in frequency scattering further reduces the speaker-
dependent variations enhancing the system performance. The performed layer-
wise analysis shows the importance of time-domain averaging over the typical
scalogram/mel-spectrogram coefficients. Results obtained with second layer
scattering coefficients indicate the relevance of amplitude modulation of time
series of different scalogram frequency bins.

We conclude that deep convolutional scattering coefficients capture more
relevant emotion-related information than the standard mel-filterbank based
feature. However, the performance with optimized parameter set on EmoDB varies
noticeably across databases, indicating the lack of generalization. We also
observe that the system faces difficulty in classifying emotions which have
similar arousal characteristics (e.g., Happy & Angry). This hints towards the
requirement of further analysis on improving the scattering coefficient based
feature representation for SER. Our work is a preliminary study that
introduces scattering transform for SER. In future, we will explore other
back-end classifiers like multi-layer perceptron, gaussian mixture models etc.
with scattering coefficients. To further evaluate the generalisation across
samples from different domains, cross-corpus analysis for SER can also be
evaluated.

## References

  * [1] M. Shah Fahad, A. Ranjan, J. Yadav, and A. Deepak, “A survey of speech emotion recognition in natural environment,” _Digital Signal Processing_ , vol. 110, p. 102951, 2021.
  * [2] M. B. Akçay and K. Oğuz, “Speech emotion recognition: Emotional models, databases, features, preprocessing methods, supporting modalities, and classifiers,” _Speech Communication_ , vol. 116, 2020.
  * [3] M. El Ayadi, M. S. Kamel, and F. Karray, “Survey on speech emotion recognition: Features, classification schemes, and databases,” _Pattern Recognition_ , vol. 44, no. 3, pp. 572–587, 2011.
  * [4] J. Parry, D. Palaz, G. Clarke, P. Lecomte, R. Mead, M. Berger, and G. Hofer, “Analysis of deep learning architectures for cross-corpus speech emotion recognition,” in _Proc. INTERSPEECH 2019_ , 2019, pp. 1656–1660.
  * [5] O. Abdel-Hamid, A. Mohamed, H. Jiang, L. Deng, G. Penn, and D. Yu, “Convolutional neural networks for speech recognition,” _IEEE/ACM Transactions on Audio, Speech, and Language Processing_ , vol. 22, no. 10, pp. 1533–1545, 2014.
  * [6] Q. Mao, M. Dong, Z. Huang, and Y. Zhan, “Learning salient features for speech emotion recognition using convolutional neural networks,” _IEEE Transactions on Multimedia_ , vol. 16, no. 8, pp. 2203–2213, 2014.
  * [7] Z. C. Lipton, “The mythos of model interpretability,” _Queue_ , vol. 16, no. 3, pp. 31–57, 2018.
  * [8] G. Trigeorgis, F. Ringeval, R. Brueckner, E. Marchi, M. A. Nicolaou, B. Schuller, and S. Zafeiriou, “Adieu features? End-to-end speech emotion recognition using a deep convolutional recurrent network,” in _Proc. ICASSP_ , 2016, pp. 5200–5204.
  * [9] H. Muckenhirn, V. Abrol, M. Magimai-Doss, and S. Marcel, “Understanding and visualizing raw waveform-based CNNs,” in _Proc. INTERSPEECH_ , 2019, pp. 2345–2349.
  * [10] D. Rolnick, A. Veit, S. Belongie, and N. Shavit, “Deep learning is robust to massive label noise,” _arXiv preprint arXiv:1705.10694_ , 2017.
  * [11] S. Mallat, “Recursive interferometric representation,” in _Proc. of EUSICO conference, Denmark_ , 2010.
  * [12] ——, “Group invariant scattering,” _Communications on Pure and Applied Mathematics_ , vol. 65, no. 10, pp. 1331–1398, 2012.
  * [13] J. Andén and S. Mallat, “Deep scattering spectrum,” _IEEE Transactions on Signal Processing_ , vol. 62, no. 16, pp. 4114–4128, 2014.
  * [14] J. Andén, V. Lostanlen, and S. Mallat, “Joint time–frequency scattering,” _IEEE Transactions on Signal Processing_ , vol. 67, no. 14, pp. 3704–3718, 2019.
  * [15] W. Ghezaiel, L. Brun, and O. Lézoray, “Hybrid network for end-to-end text-independent speaker identification,” in _International Conference on Pattern Recognition_ , Milan (virtual), Italy, Jan. 2021.
  * [16] J. Bruna and S. Mallat, “Audio texture synthesis with scattering moments,” _arXiv e-prints_ , p. arXiv:1311.0407, Nov. 2013.
  * [17] J. Salamon and J. P. Bello, “Feature learning with deep scattering for urban sound analysis,” in _Proc. EUSIPCO_ , 2015, pp. 724–728.
  * [18] C. Baugé, M. Lagrange, J. Andén, and S. Mallat, “Representing environmental sounds using the separable scattering transform,” in _Proc. ICASSP_ , 2013, pp. 8667–8671.
  * [19] V. Lostanlen, C. El-Hajj, M. Rossignol, G. Lafay, J. Andén, and M. Lagrange, “Time–frequency scattering accurately models auditory similarities between instrumental playing techniques,” _EURASIP Journal on Audio, Speech, and Music Processing_ , vol. 2021, no. 1, pp. 1–21, 2021.
  * [20] J. Bruna and S. Mallat, “Invariant scattering convolution networks,” _IEEE Transactions on Pattern Analysis and Machine Intelligence_ , vol. 35, no. 8, pp. 1872–1886, 2013.
  * [21] D. Luo, Y. Zou, and D. Huang, “Investigation on joint representation learning for robust feature extraction in speech emotion recognition,” in _Proc. INTERSPEECH_ , 2018, pp. 152–156.
  * [22] F. Burkhardt, A. Paeschke, M. Rolfes, W. F. Sendlmeier, and B. Weiss, “A database of German emotional speech,” in _Proc. INTERSPEECH_ , 2005, pp. 1517–1520.
  * [23] S. R. Livingstone and F. A. Russo, “The Ryerson Audio-Visual database of emotional speech and song (RAVDESS): A dynamic, multimodal set of facial and vocal expressions in North American English,” _PLOS One_ , vol. 13, no. 5, 2018.
  * [24] C. Busso _et al._ , “IEMOCAP: Interactive emotional dyadic motion capture database,” _Language Resources and Evaluation_ , vol. 42, no. 4, p. 335, 2008.
  * [25] A. Rosenberg, “Classifying skewed data: Importance weighting to optimize average recall,” in _Proc. INTERSPEECH_ , 2012, pp. 2242–2245.
  * [26] S. Zhang, S. Zhang, T. Huang, and W. Gao, “Speech emotion recognition using deep convolutional neural network and discriminant temporal pyramid matching,” _IEEE Transactions on Multimedia_ , vol. 20, no. 6, pp. 1576–1590, 2018.
  * [27] P. Singh, G. Saha, and M. Sahidullah, “Non-linear frequency warping using constant-q transformation for speech emotion recognition,” in _2021 International Conference on Computer Communication and Informatics (ICCCI-2021)_ , 2021.