# Action principles for dissipative, non-holonomic Newtonian mechanics

Amit Acharya Ambar N. Sengupta Department of Civil & Environmental
Engineering, and Center for Nonlinear Analysis, Carnegie Mellon University,
Pittsburgh, PA 15213, email<EMAIL_ADDRESS>Corresponding
author.Department of Mathematics, University of Connecticut, Storrs, CT 06269,
email<EMAIL_ADDRESS>
###### Abstract

A methodology for deriving dual variational principles for the classical
Newtonian mechanics of mass points in the presence of applied forces,
interaction forces, and constraints, all with a general dependence on particle
velocities and positions, is presented. Methods for incorporating constraints
are critically assessed. General theory, as well as explicitly worked out
variational principles for a dissipative system (due to Lorenz) and a system
with anholonomic constraints (due to Pars) are demonstrated. Conditions under
which a (family of) dual Hamiltonian flow(s), as well as a constant(s) of
motion, may be associated with a conservative or dissipative, and possibly
constrained, primal system naturally emerge in this work.

##  1 Introduction

In this paper, we develop a variational/action principle for Newton’s
equations of motion for a system of particles possibly subjected to kinematic
constraints which are nonlinear in the particle velocities, and forces that
may not arise from a potential function of the particle positions. Along the
way, we also critically consider the various prevalent formulations in the
literature of constraint forces. Although a classical subject, the subtleties
involved, even in the setting of anholonomic, linear-in-velocities
constraints, and the nature of unresolved issues can be appreciated from the
works [CR09, LM95]. The primary applications that we envisage of such
variational principles is to facilitate path integral based statistical
formulations of dissipative particle systems and for computing periodic orbits
of general ODE systems displaying bounded long-time behavior.

A scheme for generating action principles for dissipative initial value
problems (ordinary differential equations (ODE)) of classical physics has been
considered by [Gal13]. A special family of Hamiltonians for thermostatted
molecular dynamic systems with conservative applied forces is proposed and
reviewed in [DM96], [MD98, Sec. IV]. These works are different in concept and
scope from ours, which is also applicable to partial differential equations
(PDE) [Ach22, Ach23] and related to the line of thought advanced in [Bre18,
Bre20]. We also mention at the outset that the proposed methodology is
different from a ‘Least-Squares’ approach, and in no way hinges on the
proposed functional attaining a particular value for its validity; the
essential idea of our method can be understood in the finite dimensional
setting, as explained in [Ach23, Sec. 2]. For a review of the classical
mechanics-related aspects of our work, see [Del20].

An outline of the paper is as follows: in Sec. 2 we critically discuss the
assumptions behind the formulation of constrained particle mechanics. Sec. 3
details the main achievements of the work before presenting the fundamental
formalism for producing dual variational principles for conservative and/or
dissipative classical mechanics with general nonlinear constraints on
positions and velocities. This includes an explicitly worked out example on
the dissipative Lorenz system of ODE as well as remarks on potential use of
our methods for computing periodic orbits of general ODE systems. Sec. 4
continues with further demonstration of our scheme in the setting of a key
system with anholonomic constraints due to Pars [Par54], and its
generalization. Sec. 5 contains a proposal, free of ad-hoc assumptions, for
resolving the non-uniqueness in the form of the constraint forces in particle
systems with nonlinear constraints. Finally, Sec. 6 contains some concluding
remarks.

Holonomic systems were identified and studied by Hertz [Her10]. Non-holonomic
systems, though realistic and common in practice, have been discussed in the
literature over the years ([Fla05, CR09, BMZ05] to mention a few), but there
seems to be no standardized treatment that goes beyond constraints linear in
velocities. In the latter case, there is the coincidental advantage that the
velocities of the solution path satisfy the fixed-time variation constraints.
Griffiths [Gri83] has developed a framework that realizes the variational
equations for a given Lagrangian, subject to constraints, as describing
integral curves of exterior differential systems. Our task is different, in
that we start with the equations of motion and develop a (corresponding family
of) Lagrangian(s) (22) purely on the dual space. The Euler-Lagrange equations
of the corresponding action functionals are the given equations of motion,
interpreted through a mapping or a change of variables, of the dual variables
and their first order derivatives, defining the primal variables appearing in
those equations of motion. An additional feature of our method is the role
played by a function $H$ that effectively parametrizes classes of solutions to
the primal problem.

##  2 On the definition of forces maintaining constraints

It will suffice for us to consider all particle positions, velocities, and
forces to be expressed in components w.r.t a fixed rectangular Cartesian
coordinate system. A superposed dot will represent a time derivative and two
dots, the second time derivative. A twice repeated index in a monomial will
indicate summation over the range of the index, if unaccompanied by an
explicit summation symbol.

Consider a mechanical system of $N$ particles without any kinematic
constraints, each subjected to a net force; with
$x_{A}=(x_{1A},\ldots,x_{dA})\in{\mathbb{R}}^{d}$ denoting the configuration
coordinates of the $A$-th particle, and $f_{A}=(f_{1A},\ldots,f_{dA})$
denoting the force divided by the mass of the $A$-th particle, the equation of
motion for the system is:

$\ddot{x}=f;\qquad(\ddot{x}_{iA}=f_{iA}),\quad i=1,\dots,d;A=1,\dots,N,$

where $x=(x_{A})_{A=1,\ldots,N}$ and $f=(f_{A})_{A=1,\ldots,N}$, and, as
noted, we have absorbed the particle masses into the ‘forces’ (which now have
physical dimensions of Force/Mass = acceleration) $f_{A}$. For simplicity, it
suffices for this preliminary discussion to even consider the forces $f$ as a
given functions of time. Suppose now the system is subjected to a single
holonomic kinematic constraint given by

$g(x)=0,$

where $g:{\mathbb{R}}^{dN}\to{\mathbb{R}}$ is a smooth function. Then a
physically reasonable proposition, and one that is universally accepted, is
that the constrained dynamics is defined by the system

$\displaystyle\ddot{x}$
$\displaystyle=f-f^{(c)};\qquad\left(\ddot{x}_{iA}=f_{iA}-f^{(c)}_{iA}\right),$
(1a) $\displaystyle 0$ $\displaystyle=g(x),$ (1b)

where $f^{(c)}$ is the constraint force that arises to maintain the kinematic
constraint(s) (1b). As such, the task of the theory of Newtonian particle
mechanics is to determine a trajectory

$t\mapsto x(t)\in{\mathbb{R}}^{dN},$

with $t$ running over an interval in ${\mathbb{R}}$ containing the initial
time $0$, that satisfies (1), subject to prescribed initial conditions

$x(0)=x^{0};\qquad\dot{x}(0)=v^{0}.$

Uniqueness of solutions $t\mapsto x(t)$ satisfying (1), at least for short
times, is desirable. In the process, it is also desirable to have an idea of
the constraint forcing

$t\mapsto f^{(c)}(t)$

required in the problem.

It is clear that the above system is formally underdetermined as it has, in
general, more variables than equations. Differentiating (1b) w.r.t time along
a trajectory yields

$\nabla g(x(t))\cdot\dot{x}(t)=0,$ (2)

and if one now invokes the additional requirement that the constraint forces
expend no power along actual trajectories, i.e., for functions $x(\cdot)$ that
satisfy the equations of motion and constraints

$f^{(c)}\cdot\dot{x}=0,$ (3)

then it is clear that for $N\geq 1,d>2$, while

$f^{(c)}(t):=\lambda(t)\,\nabla g(x(t))$ (4)

for $t\mapsto\lambda(t)\in\mathbb{R}$ is a sufficient condition for the
satisfaction of (3), it is certainly not necessary.

If it is demanded that the requirement that ‘constraints do no work’ (3) hold
for _all_ possible trajectories consistent with only the constraint and not
necessarily (1a) (this is also a consequence of a constrained Hamiltonian
description of the dynamics when the applied forces, $f$, arise from a
potential), then (4) becomes a necessary condition, but it is to be realized
that this is, in essence, an extraneous physical assumption as it seems
physically natural to demand/know that a constraint force arises to maintain
the constraint only along an actual trajectory that also satisfies Newton’s
second law (1a).

Of course, the above assumption (4), while physically dubious in the above
sense, has the major advantage that it immediately makes the dynamical
description at least formally determined, that is having the right number of
equations and unknowns.

Viewed in the above manner, the argument remains essentially unchanged in the
case of anholonomic, homogeneous constraints linear in the velocities
characterized by

$a(x)\cdot\dot{x}=0,$ (5)

where the real-valued function $a$, defined on the configuration space
${\mathbb{R}}^{dN}$, takes the place of $\nabla g$, and it is d’Alembert’s
principle, an assumption, that the constraint force is required to do no work
on all trajectories consistent with the constraint (5) (not necessarily
satisfying (1a)). The above arguments also hold, in essence, in the presence
of more than one constraint, involving then the same number of Lagrange
multiplier fields $\lambda$ (the treatments are standard and can be found in
any text book on Analytical Dynamics (cf. [Ros91, Gol57]).

Within this context, it is interesting to assess the controversy regarding the
treatment of linear-in-velocity, homogeneous anholonomic constraints by
d’Alembert’s principle and by the standard Lagrange multiplier rule of the
Calculus of Variations, as described in [CR09, Fla05]. It suffices to discuss
the case of a single anholonomic constraint given by

$a(x)\cdot\dot{x}=0.$ (6)

The Calculus of Variations treatment gives the constraint force (see, e.g.,
[CR09]), written in terms of generalized coordinates as (and, which, in this
specific context can be written in terms of
$(A,i)\mapsto\alpha(i,A)=(A-1)d+i$, $x_{ia}=q_{\alpha(i,A)}$, e.g.)

$f^{(c)}_{\alpha}=-\dot{\mu}\,a_{\alpha}(q)-\mu(\partial_{\beta}\,a_{\alpha}-\partial_{\alpha}\,a_{\beta})\dot{q}_{\beta},$
(7)

whereas the constraint force from the d’Alembert treatment could be
interpreted simply as

$f^{(c)}_{\alpha}=-\dot{\mu}\,a_{\alpha}(q).$ (8)

Cronström and Raita [CR09] prove that the solutions of the dynamical equations
using these two descriptions are not identical for identical, prescribed
initial conditions. While it is certainly possible to argue the (de)merits of
any particular description, we simply note that the constraint force
description (7) does satisfy all requirements of the d’Alembert specification,
namely, that all trajectories consistent with the constraint
$a_{\alpha}\dot{q}_{\alpha}=0$ expend no power when acted on by the constraint
force.

This further underlines the indeterminacy inherent in the constrained initial
value problem (ivp) of analytical dynamics, even under the best of situations
when the corresponding unconstrained problem may have unique solutions to the
ivp.

In the same spirit, of obtaining equations of motion for constrained systems
that have equal number of unknowns and equations by making special
assumptions, is the principle of Least Constraint of Gauss [Gau29], as clearly
explained in [EHF+83, Sec. II] (also see, [Del20]). Consider a constraint with
a completely general, nonlinear in the velocities, dependence of the form

$n(x,\dot{x},t)\cdot\ddot{x}+w(x,\dot{x},t)=0,$ (9)

where $n$ and $w$ are smooth real-valued functions on
${\mathbb{R}}^{dN}\times{\mathbb{R}}^{dN}\times{\mathbb{R}}$; this
encompasses, holonomic and anholonomic constraints considered earlier. The
constraint applies to any trajectory satisfying the equation of motion (1a);
hence,

$n\cdot\left(f-f^{(c)}\right)+w=0\Longrightarrow n\cdot f^{(c)}=(n\cdot f+w).$

_Assuming_ the constraint force to be of the form

$f^{(c)}:=\lambda\,n$

for $\lambda$ a scalar valued function, one obtains the representation

$f^{(c)}=\frac{(n\cdot f+w)}{n\cdot n}\,a.$ (10)

It can be checked that a trajectory satisfying (1a) with the prescription (10)
also satisfies the constraint (9). Thus, Gauss’s principle of Least Constraint
is yet another choice of producing a constrained dynamics involving equal
number of equations and unknowns in the constrained problem. Of note here is
also the work of [FKNU05] which proposes a generalization of the Gauss
principle of Least Constraint. It is perhaps important to observe that Gauss’
Least constraint principle, given full knowledge of the arrays $(f,n,w)$,

$\displaystyle\min_{b}(b_{i}-f_{i})(b_{i}-f_{i})$ $\displaystyle\mbox{subject
to }n^{\alpha}_{i}b_{i}+w^{\alpha}=0,$

is an additional, eminently physical, assumption which does not imply a
Lagrange multiplier formulation of the problem. When a Lagrange multiplier
formulation is employed to solve the optimization problem, the critical point
is characterized by solutions of

$b_{i}=f_{i}-\lambda^{\alpha}n^{\alpha}_{i},$

identifying, in the process, the nature of the constraint force specific to
this protocol. Moreover, when utilizing the Lagrange multiplier formulation of
the problem, the uniqueness of the (instantaneous) values of the Lagrange
multipliers rests on the invertibility of the matrix
$n^{\alpha}_{i}n^{\beta}_{i}$ in

$n^{\alpha}_{i}n^{\beta}_{i}\lambda^{\beta}=w^{\alpha}+n^{\alpha}_{i}f_{i},$

and this makes it clear that even with the Lagrange multiplier formulation of
Gauss’ principle, the constraint forces may not be unique and further
conditions may then be required to induce even local-in-time uniqueness of the
solutions $t\mapsto x(t)$ of the system (1) or its anholonomic counterpart.

Ad-hoc as the above assumptions on the form of the constraint force may seem,
it is also extremely important to note that some assumption(s) beyond simply
invoking the presence of a constraint force in the equations of motion and
requiring that it expend no power when the system experiences internal
kinematic constraints is necessary on physical grounds. As a simple example,
consider the free (without any external applied force) motion of a single
particle constrained to move on the plane $x_{3}=0$. If the only requirement
of the constraint forces is that they maintain the constraint on the motion of
the particle and expend no power, circular, in-plane motions would clearly be
allowed, with the (in-plane component of the) constraint force pointing from
the particle to the center of the circle - this is physical nonsense, of
course.

In what follows, we demonstrate a family of variational principles for the
equations of constrained analytical dynamics - dissipative or conservative -
without ad-hoc assumptions, leaving enough flexibility for the use of specific
models as described above, as well as for the general mathematical study of
solutions to the equations, including that of the physically relevant
conditions needed for inducing uniqueness of solutions for the initial value
problem (ivp) of Newtonian classical mechanics of mass points.

##  3 Dual variational principle for constrained particle dynamics

None of the arguments and assumptions in Sec. 2 shed light on when the final
set of (un)constrained equations of motion can be expected to arise as the
Euler-Lagrange equations of an action principle defined on trajectories/paths
of an appropriately defined set of variables. Indeed, even in the
unconstrained case, if the applied forces, $f$, were not assumed to arise as
the gradient of a potential function in $x$, then it is generally believed
that an action principle does not exist for the dynamics.

With the above considerations in mind, we would now like to propose a scheme
that generates a family of variational principles for each of the forms of
(un)constrained dynamics considered above. In fact, we will also allow for
constraints that can do non-negative work along trajectories. As a natural
result, the Euler-Lagrange equations of our formulation always involve the
same number of dual variables as the number of primal equations ((1a) and
constraints), even when the primal system is formally underdetermined, as
discussed earlier in this Section. Furthermore, unlike Hamilton’s Principle
that has the slightly unpleasant feature of requiring information on final-
time boundary conditions on the primal position and velocities (also see
discussion in [Gal13]), a specification inconsistent with solving the well-set
ivp arising from Newton’s Laws, our dual variational principle exactly
recovers as its Euler-Lagrange equations and natural boundary conditions the
Newtonian ivp. Interestingly, our methodology can associate a family of dual
action functionals with the primal Euler-Lagrange equations arising from a
Hamiltonian structure. Our work is a natural application of ideas presented in
[Ach22, Ach23, KA23], the last of which shows a computation of (approximate)
solutions to Euler’s equations for the the dynamics of a rigid body
constrained at a single point under no torque as well as a dissipative torque
(among other PDE examples like the linear heat and transport equations). Among
the nontrivial PDE examples solved by our duality-based method are nonlinear
elasticity and the inviscid Burgers equation [SGA24, KA24].

In this section we develop a general framework for systems that are subject to
possibly non-holonomic time-dependent constraints and dissipative forces. We
demonstrate the power of the method by applying it to the case of the
dissipative and nonlinear Lorenz system (29), obtaining the dual Lorenz action
explicitly (31). The reader may also consult Sec. 4 for another example.

As before, we are considering a system with $N$ constituent components, such
as particles, with the $A$-th component having coordinates
$(x_{1A},\ldots,x_{dA})\in{\mathbb{R}}^{d}$ and velocity denoted
$(v_{1A},\ldots,v_{dA})$. The following constrained differential-algebraic
system is considered, with $f,g,W$ being suitably smooth given functions of
their arguments:

$\displaystyle\dot{x}_{iA}=v_{iA},\qquad
i=1,\ldots,d;A=1,\ldots,N\qquad\mbox{(position-velocity relationship)}$ (11)
$\displaystyle\dot{v}_{iA}=f_{iA}(x,v,Q,R,t),\qquad\mbox{(Newton's Second
Law)}$ $\displaystyle
0=g^{\alpha}(x,v,t)\qquad\alpha=1,\ldots,m<d.N,\qquad\mbox{(constraints)}$
$\displaystyle\dot{Q}_{\Gamma}=R_{\Gamma};\qquad\Gamma=1,\ldots,\bar{\Gamma},\quad
m\leq\bar{\Gamma}\leq d.N\qquad\mbox{(evolution of constraint parameters)}$
$\displaystyle 0=W(x,v,Q,R,s,t),\qquad\mbox{(power expended by constraint
forces)}$ $\displaystyle t\mapsto s(t)\in\mathbb{R}\ \mbox{is a slack variable
to convert an inequality constraint into an equality}$ $\displaystyle
x^{0}_{iA}=x_{iA}(0)\qquad\mbox{(initial condition on position)}$
$\displaystyle v^{0}_{iA}=v_{iA}(0)\qquad\mbox{(initial condition on
velocity)},$

and it is required that the prescribed initial data $(x^{0},v^{0})$ be
consistent with the constraints at the initial time, i.e.,
$0=g(x^{0},v^{0},0)$. The domain of the ${\mathbb{R}}^{dN}$-valued function
$f=(f_{iA})$ is an open subset of
${\mathbb{R}}^{dN}\times{\mathbb{R}}^{dN}\times{\mathbb{R}}^{{\bar{\Gamma}}}\times{\mathbb{R}}^{{{\bar{\Gamma}}}}\times[0,\infty)$,
and $g$ is a smooth ${\mathbb{R}}^{m}$-valued function on
${\mathbb{R}}^{2dN+1}$. Further, $W$ is a smooth real-valued function on
${\mathbb{R}}^{2dN+2{\bar{\Gamma}}+2}$. Solving the system (11) means finding
a smooth mapping

$t\mapsto\bigl{(}x(t),v(t),Q(t),R(t),s(t)\bigr{)}\in{\mathbb{R}}^{dN}\times{\mathbb{R}}^{dN}\times{\mathbb{R}}^{\bar{\Gamma}}\times{\mathbb{R}}^{\bar{\Gamma}}\times{\mathbb{R}},$

with $t\in[0,T]$, for some positive $T$, such that (11) holds when both sides
are evaluated at each time instant $t$.

We assume the following structure for the forces

$f_{iA}(x,v,Q,R,t)=f^{(a)}_{iA}(x,v,t)-f^{(c)}_{iA}(x,v,Q,R,t),$ (12)

where $f^{(c)}$ represents the constraint force.

Examples of the constraint forces for the various models discussed are

$\displaystyle f^{(c)}_{iA}$ $\displaystyle=Q^{(c)}_{iA}\qquad\mbox{no
assumption made on the form of the constraint force}$ (13a) $\displaystyle
f^{(c)}_{iA}$
$\displaystyle=-Q_{\alpha}a^{\alpha}_{iA}(x,t)\qquad\mbox{d'Alembert; }a\mbox{
is a known function of its arguments}$ (13b) $\displaystyle f^{(c)}_{iA}$
$\displaystyle=-Q_{\alpha}n^{\alpha}_{iA}(x,v,t)\qquad\mbox{Gauss; }n\mbox{ is
a known function of its arguments}$ (13c) $\displaystyle f^{(c)}_{iA}$
$\displaystyle=-\dot{Q}_{\alpha}a^{\alpha}_{iA}(x)-Q_{\alpha}\left(\frac{\partial
a^{\alpha}_{iA}}{\partial x_{jB}}-\frac{\partial a^{\alpha}_{jB}}{\partial
x_{iA}}\right)v_{jB}$ (13d) $\displaystyle\qquad\mbox{Hamiltonian formalism
for a homogeneous, linear-in-velocity constraint}.$

An example of constraint forces expending non-negative power, which also shows
the use of a slack variable $s$, is

$W(x,v,Q,R,s,t):=\sum_{i,a}m_{A}f^{(c)}_{iA}(x,v,Q,R,t)v_{iA}-s^{2}=0.$ (14)

where $m_{A}$ is the mass of particle $A$; note that, by our notational
convention, $f_{iA}$ is the $i$-th component of the force on the $A$-th
particle divided by the mass $m_{A}$.

When $f^{(c)}(x,v,Q,R,t)=Q$ then the constraint parameter $Q$ has the meaning
of the usual Lagrange multiplier. The form of the constraint force in the
Hamiltonian formalism with constraints dictates that the constraint forces can
depend on $\dot{Q}$; motivated by this, we allow for a general nonlinear
dependence on $\dot{Q}$ of the constraint force, which requires, within our
formalism, to introduce a new variable which we denote by $R:=\dot{Q}$ and
include it within our governing set of equations - this is no different, in
principle, from the introduction of a new variable $v:=\dot{x}$ in classical
mechanics.

The totality of the _primal_ physical functions whose evolution is of interest
are

$U:=(x,v,Q,R,s);$ (15)

the components of $U$ can be thought of as coordinate functions on an open
subset of ${\mathbb{R}}^{2dN+2{\bar{\Gamma}}+1}$, and we are interested in
motions $t\mapsto U(t)$ that satisfy the conditions (11).

We introduce an array of _dual_ functions of time

$\displaystyle D:=(\rho,\lambda,\mu,\Lambda,\kappa)$ (16)
$\displaystyle\rho_{iA},\qquad i=1,\ldots,d;A=1,\ldots,N$
$\displaystyle\lambda_{iA},\qquad i=1,\ldots,d;A=1,\ldots,N$
$\displaystyle\mu_{\alpha},\qquad\alpha=1,\ldots,m<d.N$
$\displaystyle\kappa_{\Gamma},\qquad\Gamma=1,\ldots,\bar{\Gamma},\quad
m\leq\bar{\Gamma}\leq d.N$

to be used subsequently, as well as the array

$\mathcal{D}:=(\dot{\rho},\rho,\dot{\lambda},\lambda,\mu,\Lambda,\dot{\kappa},\kappa).$
(17)

Technically, the components of $D$ can be thought of as coordinate functions
on ${\mathbb{R}}^{2dN+m+1+{\bar{\Gamma}}}$. Then define the pre-dual
functional

$\displaystyle\widehat{S}_{H}[x,v,Q,\rho,\lambda,\mu,\Lambda,\kappa]=-\rho(0)\cdot
x^{0}-\lambda(0)\cdot
v^{0}+\int_{0}^{T}\mathcal{L}_{H}(U(t),\mathcal{D}(t),t)\,dt,$ (18)

where on the left the argument $(x,v,Q,\rho,\lambda,\mu,\Lambda,\kappa)$ is
any smooth function
$[0,T]\to{\mathbb{R}}^{2dN+{\bar{\Gamma}}+2dN+m+{\bar{\Gamma}}}$, and the
integrand on the right is

$\mathcal{L}_{H}(U,\mathcal{D},t):=-\dot{\rho}\cdot x-\rho\cdot
v-\dot{\lambda}\cdot v-\lambda\cdot f+\mu\cdot g+\Lambda W-\dot{\kappa}\cdot
Q-\kappa\cdot R+H(U,\overline{U}(t)),$ (19)

and

$\overline{U}(t):=(\bar{x}(t),\bar{v}(t),\bar{Q}(t),\bar{R}(t),\bar{s}(t))$
(20)

is a collection of the arbitrarily specified functions of time displayed, _and
$H$ is any smooth function of its arguments that enables the construction of a
function_

$U^{(H)}(\mathcal{D},t)=\left(x^{(H)}(\mathcal{D},t),v^{(H)}(\mathcal{D},t),Q^{(H)}(\mathcal{D},t),R^{(H)}(\mathcal{D},t),s^{(H)}(\mathcal{D},t)\right)$

_such that_

$\frac{\partial\mathcal{L}_{H}}{\partial
U}\left(U^{(H)}(\mathcal{D},t),\mathcal{D},t\right)=0\qquad\forall(\mathcal{D},t),$
(21)

_i.e., $\frac{\partial\mathcal{L}_{H}}{\partial U}(U,\mathcal{D},t)=0$ should
be solvable for $U$ in terms of $(\mathcal{D},t)$._ To be more precise, here
$H$ is a real-valued smooth function on
$({\mathbb{R}}^{2dN+2{\bar{\Gamma}}+1})^{2}$. (Note that $H$ does not refer to
a Hamiltonian.)

_We refer to the function $U^{(H)}$ as the dual-to-primal (DtP) mapping._

The functional $\widehat{S}_{H}$ is obtained by taking appropriate scalar
products of (11) with the dual fields (Lagrange multipliers) $D$, integrating
by parts over the time interval $[0,T]$, and applying the primal boundary
conditions available and simply ignoring the remaining boundary terms at this
point.

Then define the _dual_ functional

$\displaystyle S_{H}[\rho,\lambda,\mu,\Lambda,\kappa]=-\rho(0)\cdot
x^{0}-\lambda(0)\cdot
v^{0}+\int_{0}^{T}\mathcal{L}_{H}\left(U^{(H)}(\mathcal{D}(t),t),\mathcal{D}(t),t\right)\,dt$
(22) $\displaystyle\mbox{with Dirichlet boundary
conditions}\quad\lambda(T),\rho(T),\kappa(T),\kappa(0)\quad\mbox{specified
arbitrarily.}$

Note that, on the left, $[\rho,\lambda,\mu,\Lambda,\kappa]$ is any smooth
function $[0,T]\to{\mathbb{R}}^{2dN+m+1+{\bar{\Gamma}}}$.

Utilizing the requirement (21), and the boundary conditions

$\delta\lambda(T)=0,\delta\rho(T)=0,\delta\kappa(T)=0,\delta\kappa(0)=0,$

the first variation of $S_{H}$ is given by (using the short-hand and slight
abuse of notation $U^{(H)}(t)=U^{(H)}(\mathcal{D}(t),t)$ and
$f^{(H)}(t)=f(U^{(H)}(t))$, $g^{(H)}(t)=g(U^{(H)}(t))$,
$W^{(H)}(t)=W(U^{(H)}(t))$ )

$\displaystyle\delta S_{H}$ $\displaystyle=-\delta\rho(0)\cdot
x^{0}-\delta\lambda(0)\cdot
v^{0}+\int_{0}^{T}\frac{\partial\mathcal{L}_{H}}{\partial\mathcal{D}}\left(U^{(H)}(\mathcal{D}(t),t),\mathcal{D}(t),t\right)\cdot\delta\mathcal{D}\,dt$
(23)
$\displaystyle=\left(x^{(H)}(0)-x^{0}\right)\cdot\delta\rho(0)+\int_{0}^{T}\delta\rho(t)\cdot\left(\dot{\overline{x^{(H)}}}(t)-v^{(H)}(t)\right)\,dt$
$\displaystyle\quad+\left(v^{(H)}(0)-v^{0}\right)\cdot\delta\rho(0)+\int_{0}^{T}\delta\lambda(t)\cdot\left(\dot{\overline{v^{(H)}}}(t)-f^{(H)}(t)\right)\,dt$
$\displaystyle\quad+\int^{T}_{0}\delta\mu(t)\cdot
g^{(H)}(t)\,dt+\int^{T}_{0}\delta\Lambda(t)W^{(H)}(t)\,dt+\int^{T}_{0}\delta\kappa(t)\cdot\left(\dot{\overline{Q^{(H)}}}(t)-R^{(H)}(t)\right)\,dt.$

Apart from a direct computation, this may also be understood by observing that
$\mathcal{L}_{H}(U,\mathcal{D},t)$ is affine in $\mathcal{D}$ by construction.
We note that the Euler-Lagrange equations of the dual functional $S_{H}$,
_parametrized by the function $H$_, is the set of primal equations (11) with
the replacement

$(U,f,g,W)\rightarrow\left(U^{(H)},f^{(H)},g^{(H)},W^{(H)}\right).$

To understand some salient properties of our variational principle, let us now
consider a special case where no assumptions are made on the form of the
constraint forces (13a) so that

$f(x,v,Q,t)=f^{(a)}(x,v,t)-Q.$

For simplicity in conveying ideas, let us also assume that the power of
constraint forces statement is presented as an equality so that the slack
variable $s$ is not required, and that the functions $W$ is independent of
$R$. Furthermore, we assume that the constraint function $g$ does not depend
on $Q$ (and $s$).

Then it suffices to choose the function $H$ to be a ‘shifted quadratic’ form
given by

$H(x,v,Q,t)=\frac{1}{2}\left(c_{x}|x-\bar{x}(t)|^{2}+c_{v}|v-\bar{v}(t)|^{2}+c_{Q}|Q-\bar{Q}(t)|^{2}\right),$
(24)

where $(x,v,Q)$ runs over ${\mathbb{R}}^{3dN}$ here.

The key requirement (21) now becomes the condition that the following
algebraic system of equations

$\displaystyle\frac{\partial\mathcal{L}_{H}}{\partial x_{iA}}$
$\displaystyle:\qquad
c_{x}(x-\bar{x})_{iA}-\dot{\rho}_{iA}-\lambda_{jB}\frac{\partial
f^{(A)}_{jB}}{\partial x_{iA}}(x,v,t)+\mu_{\alpha}\frac{\partial
g^{\alpha}}{\partial x_{iA}}(x,v,t)+\Lambda\frac{\partial W}{\partial
x_{iA}}(x,v,t)=0$ (25) $\displaystyle\frac{\partial\mathcal{L}_{H}}{\partial
v_{iA}}$ $\displaystyle:\qquad
c_{v}(v-\bar{v})_{iA}-\rho_{iA}-\dot{\lambda}_{iA}-\lambda_{jB}\frac{\partial
f^{(A)}_{jB}}{\partial v_{iA}}(x,v,t)+\mu_{\alpha}\frac{\partial
g^{\alpha}}{\partial v_{iA}}(x,v,t)+\Lambda\frac{\partial W}{\partial
v_{iA}}(x,v,t)=0$ $\displaystyle\frac{\partial\mathcal{L}_{H}}{\partial
Q_{iA}}$ $\displaystyle:\qquad
c_{Q}(Q-\bar{Q})_{iA}+\lambda_{iA}+\Lambda\frac{\partial W}{\partial
Q_{iA}}(x,v,Q,t)=0$

be solvable for $(x,v,Q)$ in terms of

$\mathcal{D}=(\dot{\rho},\rho,\dot{\lambda},\lambda,\mu,\Lambda)\qquad\mbox{and}\qquad(\bar{x}(t),\bar{v}(t),\bar{Q}(t))$
(26)

to define the functions

$\displaystyle(\mathcal{D},t)$ $\displaystyle\mapsto x^{(H)}(\mathcal{D},t)$
(27) $\displaystyle(\mathcal{D},t)$ $\displaystyle\mapsto
v^{(H)}(\mathcal{D},t)$ $\displaystyle(\mathcal{D},t)$ $\displaystyle\mapsto
Q^{(H)}(\mathcal{D},t).$

Suppose, now, this can be done (and definitely in a perturbative sense, by
choosing $c_{x},c_{v},c_{q}\to+\infty$). Then the following remarks are in
order:

  1. 1.

Let ($t\mapsto\bar{x}(t)$, $t\mapsto\bar{v}(t)$, $t\mapsto\bar{Q}(t)$) be a
solution to the primal problem (11) in this special case for given initial
data. Then, for the quadratic $H$ defined in (24) by these specified
$\bar{(\cdot)}$ functions, _the dual variational principle has a critical
point, explicitly given by_

$\left(t\mapsto\rho(t)=0\in\mathbb{R}^{d.N},\quad
t\mapsto\lambda(t)=0\in\mathbb{R}^{d.N},\quad
t\mapsto\mu(t)=0\in\mathbb{R}^{m},\quad
t\mapsto\Lambda(t)=0\in\mathbb{R}\right).$

This is an existence results for our proposed scheme. More importantly, it
makes it reasonable to expect solutions to the dual problem to exist when the
functions $(\bar{x},\bar{v},\bar{Q})$ chosen to define $H$ are close to actual
solutions. Of course, even in other situations vastly different from these
choices, solutions very likely exist, as motivated in [KA23] through the
computation of controlled approximations. Rigorous existence results for very
weak and weak solutions in more general PDE contexts exist, e.g., the Euler
equations and inviscid Burgers equation, see [Bre18, SGA24].

  2. 2.

We note that if the primal problem has a unique solution for given initial
data and the variational dual problem has solutions for a collection of $H$’s
and Dirichlet b.c.s, then the primal solutions defined by the differing class
of $U^{(H)}$ are one and the same, regardless of the choice of $H$ (and b.c.s)
defining them. This may be considered as a special type of gauge invariance of
our procedure.

  3. 3.

For $c_{x},c_{v},c_{Q}$ large (and definitely when the primal system is
linear),

$x^{(H)}\sim\bar{x}+\frac{1}{c_{x}}\dot{\rho},\qquad
v^{(H)}\sim\bar{v}+\frac{1}{c_{v}}\dot{\lambda}$ (28)

and the dual Euler-Lagrange equation for $x^{(H)},v^{(H)}$ become second-order
ODE and the scheme proposes to solve it as a boundary value problem by
arbitrarily specifying final-time conditions on ($\rho(T),\lambda(T)$). It is
then reasonable to ask whether such a specification can act as an obstruction
to describing the correct solution of the primal _initial value problem_ at
time $T$, for solutions defined through the DtP mapping. The answer to the
question is that such an obstruction does not arise as the mapping (28) shows
that $(x^{(H)}(T),v^{(H)}(T))$ are functions of
$(\dot{\rho}(T),\dot{\lambda}(T))$, respectively, and specifying the values of
$(\rho(T),\lambda(T))$ leaves the derivatives adjustable to the demands of
solving the primal ivp.

  4. 4.

Defining a function

$\mathbb{L}(D,\dot{D},t):=\mathcal{L}_{H}\left(U^{(H)}(\mathcal{D},t),\mathcal{D},t\right)$

(since $\mathcal{D}$ is a function of $(D,\dot{D},t)$), if

$\frac{\partial\mathbb{L}}{\partial\dot{D}}(D,\dot{D},t)=P$

is _uniquely_ solvable for

$\dot{D}=\mathcal{R}(D,P,t),$

then it is an elementary result of Hamiltonian mechanics that a Hamiltonian

$\mathbb{H}(D,P,t):=P\cdot\mathcal{R}(D,P,t)-\mathbb{L}(D,\mathcal{R}(D,P,t),t)$

can be constructed by a Legendre transform as shown, and the corresponding
Hamiltonian dynamics is essentially equivalent to our dual Lagrangian
dynamics, i.e.,

$\frac{dD}{dt}=\frac{\partial\mathbb{H}}{\partial
P}\,;\quad\frac{dP}{dt}=-\frac{\partial\mathbb{H}}{\partial
D}\qquad\mbox{`}\Longleftrightarrow\mbox{'}\qquad\frac{d}{dt}\left(\frac{\partial\mathbb{L}}{\partial\dot{D}}\right)-\frac{\partial\mathbb{L}}{\partial
D}=0.$

Moreover, when $\mathbb{L}$ is independent of $t$ (certainly when the $H$
chosen to define $\mathcal{L}_{H}$ is independent of $t$, e.g., by not
involving any of the overhead $\ \bar{}\ $ functions), then

$\mathbb{H}(D,P)$

is a constant of the dual evolution and, hence, of the primal, possibly
dissipative, evolution - since every dual evolution satisfying the dual E-L
equation defines a primal motion by our basic scheme.

Since our scheme allows the association of, potentially, entire families of
dual evolutions with a single primal system depending on the choice of the
potential $H$, this suggests that there is the possibility of associating many
constants of motion, in this sense, with a primal evolution, conservative or
dissipative. This is an intriguing prospect whose implication needs to be
explored in future work.

  5. 5.

The dual action generated above for a dissipative system by our scheme and, in
general, for a system of ODE of the type

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

including the particle systems that are the main subject of this paper, can be
useful for computing (approximate) periodic orbits by looking for critical
points of the functionals in the class of periodic functions (cf. [Rab78]),
complementary to the ‘Least-Squares’ approaches in [BFL+11, LC04]. This is a
potential application of our methods that is of considerable practical
importance.

For $T$ (a possible period of a periodic orbit to be determined), introducing
the change of the time variable

$s=\frac{2\pi t}{T};\qquad\mathcal{P}:=\frac{T}{2\pi}$

one seeks to find $2\pi$-periodic solutions $(x,\mathcal{P})$ of the system

$\frac{dx}{ds}=\mathcal{P}F(x);\qquad\frac{d\mathcal{P}}{ds}=0.$

The problem can now be approached as looking for critical points of dual
functionals developed by our scheme on the fixed interval $[0,2\pi]$ within
the class of functions satisfying periodic boundary conditions on the dual
fields. The flexibility provided by our scheme through the choice of the
function $H$ which, in turn, allows the further use of the ‘base states’
$\bar{x}$ may be expected to benefit the search for identifying a substantial
class of periodic orbits.

The proposed methodology is general and applies to a large collection of ODE
systems (converted to first-order form). To explicitly demonstrate the
‘algorithmic’ nature of the scheme, we apply our method to the dissipative and
nonlinear Lorenz system [Lor63] written as (with $A,R,B>0$ as given
parameters)

$\displaystyle\dot{x}-A(y-x)=0;\qquad\dot{y}-x(R-z)+y=0;\qquad\dot{z}-xy+Bz=0$
(29) $\displaystyle\mbox{with initial conditions }\quad x(0)=x^{0};\quad
y(0)=y^{0};\quad z(0)=z^{0}.$

It is known that the Lorenz dynamics is dissipative (and chaotic for a certain
parameter range) and long-time solutions are restricted to a small
neighborhood of a bounded attracting set of phase space after, possibly, an
initial transient; we will assume that initial conditions are chosen within a
bounded region of phase space, which contains the aforementioned attracting
set.

Next, with the choice

$H(x,y,z,t)=\frac{1}{2}c\left((x-\bar{x}(t))^{2}+(y-\bar{y}(t))^{2}+(z-\bar{z}(t))^{2}\right),$

where the functions $(\bar{x},\bar{y},\bar{z})$ of time are free to choose as
is the constant $c>0$, form the pre-dual functional

$\displaystyle\widehat{S}_{H}[U,D]$
$\displaystyle=\int_{0}^{T}\mathcal{L}_{H}\left(U,\mathcal{D}\right)\,dt-\lambda(0)x^{0}-\mu(0)y^{0}-\gamma(0)z^{0}\mbox{
where }$ $\displaystyle\mathcal{L}_{H}\left(U,\mathcal{D}\right)$
$\displaystyle=-x\dot{\lambda}-\lambda A(y-x)-y\dot{\mu}-\mu x(R-z)+\mu
y-z\dot{\gamma}-\gamma xy+\gamma Bz+H(U,t),$

where

$U=(x,y,z);\qquad
D=(\lambda,\mu,\gamma);\qquad\mathcal{D}=(\lambda,\dot{\lambda},\mu,\dot{\mu},\gamma,\dot{\gamma}).$

The DtP mapping $U^{(H)}(\mathcal{D})$ is generated from solving the equations

$\left.\begin{aligned} \begin{aligned} \frac{\partial\mathcal{L}_{H}}{\partial
x}=0&:\qquad-\dot{\lambda}+A\lambda-\mu R-\gamma y+\mu z+c(x-\bar{x})=0\\\
\frac{\partial\mathcal{L}_{H}}{\partial y}=0&:\qquad-\lambda
A-\dot{\mu}+\mu-\gamma x+c(y-\bar{y})=0\\\
\frac{\partial\mathcal{L}_{H}}{\partial z}=0&:\qquad\mu x-\dot{\gamma}+\gamma
B+c(z-\bar{z})=0\end{aligned}\end{aligned}\right\\}\qquad\mathbb{A}\big{|}_{D}\left(U^{(H)}-\bar{U}\right)=p\big{|}_{(\mathcal{D},\bar{U})}$

with the matrix defined as

$\mathbb{A}\big{|}_{D}=c\begin{bmatrix}1&-\frac{\gamma}{c}&\frac{\mu}{c}\\\
-\frac{\gamma}{c}&1&0\\\ \frac{\mu}{c}&0&1\end{bmatrix}\quad\mbox{ with
}\quad\left(\mathbb{A}\big{|}_{D}\right)^{-1}=\frac{1}{c\left(1-\frac{\gamma^{2}}{c^{2}}-\frac{\mu^{2}}{c^{2}}\right)}\begin{bmatrix}1&\frac{\gamma}{c}&-\frac{\mu}{c}\\\
\frac{\gamma}{c}&\left(1-\frac{\mu^{2}}{c^{2}}\right)&-\frac{\gamma\mu}{c^{2}}\\\
-\frac{\mu}{c}&-\frac{\gamma\mu}{c^{2}}&\left(1-\frac{\gamma^{2}}{c^{2}}\right)\end{bmatrix}=:\mathbb{B}\big{|}_{D},$

and the vectors as

$U^{(H)}-\bar{U}=\begin{bmatrix}x^{(H)}-\bar{x}\\\ y^{(H)}-\bar{y}\\\
z^{(H)}-\bar{z}\end{bmatrix};\qquad
p\big{|}_{(\mathcal{D},\bar{U})}=\begin{bmatrix}\gamma\bar{y}-\mu\bar{z}+\dot{\lambda}-A\lambda+R\mu\\\
\gamma\bar{x}+\dot{\mu}-\mu+A\lambda\\\
-\mu\bar{x}+\dot{\gamma}-B\gamma\end{bmatrix},$

with solution given by

$U^{(H)}-\bar{U}=\mathbb{B}\big{|}_{D}\ p\big{|}_{(\mathcal{D},\bar{U})}.$
(30)

The inverse matrix exists by making a choice of $c$ large enough (using the
boundedness of the flow). Computation, with suitable algebraic grouping, shows
that

$\displaystyle\widehat{S}_{H}[U^{(H)},D]$
$\displaystyle=\int_{0}^{T}-(x^{(H)}-\bar{x})\left(\dot{\lambda}-\lambda
A+R\mu+\gamma\bar{y}-\mu\bar{z}\right)\,dt-\lambda(0)x^{0}$
$\displaystyle\qquad+\int_{0}^{T}-(y^{(H)}-\bar{y})\left(\lambda
A+\dot{\mu}-\mu+\gamma\bar{x}\right)\,dt-\mu(0)y^{0}$
$\displaystyle\qquad+\int_{0}^{T}-(z^{(H)}-\bar{z})\left(\dot{\gamma}+B\gamma-\mu\bar{x}\right)\,dt-\gamma(0)z^{0}$
$\displaystyle\qquad-\int_{0}^{T}\bar{U}\cdot
p\big{|}_{(\mathcal{D},\bar{U})}\,dt+\int_{0}^{T}\mu\left(x^{(H)}-\bar{x}\right)\left(z^{(H)}-\bar{z}\right)-\gamma\left(x^{(H)}-\bar{x}\right)\left(y^{(H)}-\bar{y}\right)\,dt$
$\displaystyle\qquad+\int_{0}^{T}-\mu\bar{x}\bar{z}+\gamma\bar{x}\bar{y}+H\left(U^{(H)},t\right)\,dt.$

Noting, further, that

$\displaystyle\frac{1}{2}(U-\bar{U})\cdot\mathbb{A}(U-\bar{U})$
$\displaystyle=\frac{1}{2}c\left(\left(x-\bar{x}\right)^{2}+\left(y-\bar{y}\right)^{2}+\left(z-\bar{z}\right)^{2}\right)$
$\displaystyle\quad+\mu\left(x-\bar{x}\right)\left(z-\bar{z}\right)-\gamma\left(x-\bar{x}\right)\left(y-\bar{y}\right),$

we obtain the corresponding _dual Lorenz action_ as:

$\displaystyle S_{H}[D]$
$\displaystyle=\int_{0}^{T}\left(-\frac{1}{2}p\big{|}_{(\mathcal{D},\bar{U})}\cdot\mathbb{B}\big{|}_{D}\
p\big{|}_{(\mathcal{D},\bar{U})}-\bar{U}\cdot
p\big{|}_{(\mathcal{D},\bar{U})}-\mu\bar{x}\bar{z}+\gamma\bar{x}\bar{y}\right)\,dt$
(31) $\displaystyle\quad-\lambda(0)x^{0}-\mu(0)y^{0}-\gamma(0)z^{0}$
$\displaystyle\mbox{ with }\lambda(T),\mu(T),\gamma(T)\mbox{ specified
arbitrarily}.$

Our formalism and results guarantee that the corresponding E-L equation and
natural boundary conditions lead to

$\displaystyle\frac{d{x}^{(H)}}{dt}-A\left(y^{(H)}-x^{(H)}\right)=0;\quad\frac{d{y}^{(H)}}{dt}-x^{(H)}\left(R-z^{(H)}\right)+y^{(H)}=0;\quad\frac{d{z}^{(H)}}{dt}-x^{(H)}y^{(H)}+Bz^{(H)}=0$
$\displaystyle\quad x^{(H)}\big{|}_{t=0}=x^{0};\quad
y^{(H)}\big{|}_{t=0}=y^{0};\quad z^{(H)}\big{|}_{t=0}=z^{0}.$

We note from (30) that if
$(t\mapsto\bar{x}(t),t\mapsto\bar{y}(t),t\mapsto\bar{z}(t))$ is a solution of
the Lorenz system then the corresponding dual Lorenz action has a critical
point given by $(t\mapsto\lambda(t)=0,t\mapsto\mu(t)=0,t\mapsto\gamma(t)=0)$.
In passing, we also note that for $(\bar{x}=0,\bar{y}=0,\bar{z}=0)$ the dual
Lagrangian is a negative semi-definite form for $c\gg 1$ (assuming boundedness
of the dual flow).

Based on what has been said in remark 4. above, it can be checked that when
there is a solution to the dual Lorenz problem, defining (with
$D=(\lambda,\mu,\gamma)$)

$\displaystyle\mathbb{L}(D,\dot{D},t)$
$\displaystyle=-\frac{1}{2}p\big{|}_{(D,\dot{D},\bar{U}(t))}\cdot\mathbb{B}\big{|}_{D}\
p\big{|}_{(D,\dot{D},\bar{U}(t))}-\bar{U}\big{|}_{t}\cdot
p\big{|}_{(D,\dot{D},\bar{U}(t))}-\mu(\bar{x}\bar{z})\big{|}_{t}+\gamma(\bar{x}\bar{y})\big{|}_{t}$
$\displaystyle\tilde{p}\big{|}_{(D,\bar{U})}$
$\displaystyle:=p\big{|}_{(D,\dot{D},\bar{U})}-\dot{D}$
$\displaystyle\mathcal{R}(D,P,t)$
$\displaystyle=-\mathbb{A}\big{|}_{D}\left(P+\bar{U}(t)\right)-\tilde{p}\big{|}_{(D,\bar{U}(t))},$

a Hamiltonian

$\mathbb{H}(D,P):=\mathcal{R}(D,P,t)\cdot
P-\mathbb{L}(D,\mathcal{R}(D,P,t),t)$

can be associated with the dissipative Lorenz evolution, and it is a constant
of the Lorenz motion when $\bar{U}$ is constant (in the sense of the remark
4).

##  4 Generalization of the example of Pars [Par54]

We consider the system (with all particle masses assumed equal to $1$ in
appropriate units)

$\displaystyle\dot{x}_{i}-v_{i}$ $\displaystyle=0;\qquad
i=1,\dots,M\in\mathbb{Z}^{+}$ (32a) $\displaystyle\dot{v}_{i}+Q_{i}$
$\displaystyle=0$ (32b) $\displaystyle(b_{i}+L_{ij}x_{j})v_{i}$
$\displaystyle=0$ (32c) $\displaystyle Q_{i}v_{i}-\frac{1}{2}s^{2}$
$\displaystyle=0$ (32d) with initial conditions $\displaystyle x_{i}(0)$
$\displaystyle=x^{0}_{i}$ (32e) $\displaystyle v_{i}(0)$
$\displaystyle=v^{0}_{i}$ (32f)
$\displaystyle\left(b_{i}+L_{ij}x^{0}_{j}\right)v^{0}_{i}$ $\displaystyle=0.$
(32g)

Here $L,b$ are a given constant matrix and a vector, respectively, and the
initial data $(x^{0},v^{0})$ is assumed to be consistent with (32g). The
system above can represent holonomic, anholonomic, and non-integrable
constraints, including the example of Pars [Par54, CR07].

It should be noted that (32d, 32b) imply that along any solution of (32) the
rate of change of kinetic energy is non-positive, i.e.,

$\dot{K}\leq 0,\qquad K=\frac{1}{2}v_{i}v_{i}.$ (33)

In case external forces and forces of interaction arising from a potential
$E(x)$ were to be involved to give the equation of motion

$\dot{v}+Q=f^{ext}-\partial_{x}E,$

the corresponding power of constraints statement would be an expression of
non-negative mechanical dissipation given by

$Q\cdot v=f^{ext}\cdot v-\frac{d}{dt}(E+K)\geq 0.$

Denoting

$U=(x,v,Q,s);\qquad
D=(\rho,\lambda,\mu,\Lambda);\qquad\mathcal{D}=(\dot{\rho},\rho,\dot{\lambda},\lambda,\mu,\Lambda),$

we have

$\widehat{S}_{H}[x,v,Q,s,\rho,\lambda,\mu,\Lambda]=\int^{T}_{0}\mathcal{L}_{H}(U(t),\mathcal{D}(t),t)\,dt-\rho(0)\cdot
x^{0}-\lambda(0)\cdot v^{0},$

where

$\displaystyle\mathcal{L}_{H}(U,\mathcal{D},t)$
$\displaystyle=-x\cdot\dot{\rho}-v\cdot\rho-v\cdot\dot{\lambda}+\lambda\cdot
Q+\mu v\cdot(Lx)+\mu b\cdot v+\Lambda Q\cdot v-\frac{1}{2}\Lambda s^{2}$ (34)
$\displaystyle\quad+\frac{1}{2}\left(c_{x}|x-\bar{x}(t)|^{2}+c_{v}|v-\bar{v}(t)|^{2}+c_{Q}|Q-\bar{Q}(t)|^{2}+c_{s}(s-\bar{s}(t))^{2}\right).$

The DtP mapping is obtained from solving the system

$\displaystyle\frac{\partial\mathcal{L}_{H}}{\partial x_{i}}$
$\displaystyle:\qquad-\dot{\rho}_{i}+\mu
L_{ki}v_{k}+c_{x}(x_{i}-\bar{x}_{i})=0$ (35)
$\displaystyle\frac{\partial\mathcal{L}_{H}}{\partial v_{i}}$
$\displaystyle:\qquad-\rho_{i}-\dot{\lambda}_{i}+\mu b_{i}+\mu
L_{ij}x_{j}+\Lambda Q_{i}+c_{v}(v_{i}-\bar{v}_{i})=0$
$\displaystyle\frac{\partial\mathcal{L}_{H}}{\partial Q_{i}}$
$\displaystyle:\qquad\lambda_{i}+\Lambda v_{i}+c_{Q}(Q_{i}-\bar{Q}_{i})=0$
$\displaystyle\frac{\partial\mathcal{L}_{H}}{\partial s}$
$\displaystyle:\qquad-\Lambda s+c_{s}(s-\bar{s})=0.$

We now _assume_ that the functions $\mu^{2}$ and $\Lambda^{2}$ are bounded
above in $[0,T]$ so that $c_{x},c_{v},c_{Q}$ can be chosen large enough to
make the matrix

$\mathbb{M}_{ik}(\mu^{2},\Lambda^{2}):=\delta_{ik}-\frac{1}{c_{x}c_{v}}\mu^{2}L_{ij}L_{kj}-\frac{1}{c_{Q}c_{v}}\Lambda^{2}\delta_{ik}$

invertible.

Then the solution of (35) is given by

$\displaystyle v_{k}^{(H)}(\mathcal{D},t)$
$\displaystyle=(\mathbb{M}^{-1})_{ki}\bar{v}_{i}(t)+\frac{1}{c_{v}}(\mathbb{M}^{-1})_{ki}\left(-\mu
L_{ij}\bar{x}_{j}(t)-\Lambda\bar{Q}_{i}(t)+\rho_{i}+\dot{\lambda}_{i}-\mu
b_{i}-\frac{1}{c_{x}}\mu
L_{ij}\dot{\rho}_{j}+\frac{1}{c_{Q}}\Lambda\lambda_{i}\right)$ (36)
$\displaystyle x^{(H)}_{k}(\mathcal{D},t)$
$\displaystyle=\bar{x}_{i}(t)+\frac{1}{c_{x}}\dot{\rho}_{i}-\frac{1}{c_{x}}\mu
L_{ki}v^{(H)}_{k}(\mathcal{D},t)$ $\displaystyle Q^{(H)}_{i}(\mathcal{D},t)$
$\displaystyle=\bar{Q}_{i}(t)-\frac{1}{c_{Q}}\lambda_{i}-\frac{1}{c_{Q}}\Lambda
v^{(H)}_{i}(\mathcal{D},t)$ $\displaystyle s^{(H)}(\mathcal{D},t)$
$\displaystyle=\frac{c_{s}\bar{s}(t)}{c_{s}-\Lambda}.$

The dual functional is then given by

$S_{H}[D]=\int^{T}_{0}\mathcal{L}_{H}\left(U^{(H)}(\mathcal{D}(t),t),\mathcal{D}(t),t\right)\,dt-\rho(0)\cdot
x^{0}-\lambda(0)\cdot v^{0},$ (37)

with

$\rho(T),\lambda(T)\quad\mbox{specified (arbitrarily)}.$

Since

$\frac{\partial\mathcal{L}_{H}}{\partial
U}\left(U^{(H)}(\mathcal{D},t),\mathcal{D},t\right)=0\qquad\mbox{for all
}(\mathcal{D},t)\ \mbox{by design},$

we have

$\displaystyle\delta
S_{H}=\int^{T}_{0}\frac{\partial\mathcal{L}_{H}}{\partial\mathcal{D}}\left(U^{(H)}(\mathcal{D},t),\mathcal{D},t\right)\cdot\delta\mathcal{D}\,dt-\delta\rho(0)\cdot
x^{0}-\delta\lambda(0)\cdot v^{0}$

and repeating the general arguments of Sec. 3 it can be checked that the
primal system (32) is recovered with the replacement

$(x(t),v(t),Q(t),s(t)):=\left(x^{(H)}(\mathcal{D}(t),t),v^{(H)}(\mathcal{D}(t),t),Q^{(H)}(\mathcal{D}(t),t),s^{(H)}(\mathcal{D}(t),t)\right).$

###  4.1 The example of Pars: variations on the theme

A dynamical system is holonomic if its evolution is specified by applied
forces and constraints, possibly time-dependent, on the coordinates specifying
the constituents of the system. If $\\{x_{i}\\}$ are coordinates specifying a
state of the system, holonomic dynamics is given by an equation for the
accelerations of the form ${\ddot{x}}_{i}=F_{i}^{(a)}(x,t)$ and constraints
$g^{\alpha}(x,t)=0$, for $\alpha$ running over a finite set of indices, and
$F_{i}^{(a)}$ are the applied forces.

Our interest is mainly in non-holonomic, possibly dissipative, systems,
specifically those where the constraints are of the form $g^{\alpha}(x,v,t)=0$
and forces of the form $F_{i}^{(a)}(x,v,t)$, where $v$ runs over the
velocities of the constituents of the system. For a brief history and a modern
perspective, as well as other references, on non-holonomic dynamics we refer
to [BMZ05].

For holonomic systems, the traditional method is the d’Alembert principle
which, in its basic form, states that the work done for virtual displacements
by the constraint forces is zero (the d’Alembert principle also applies to
systems with constraints that are linear, but not affine, in the velocities);
technically, this means that, with $F^{(a)}_{i}$ denoting the $i$-th
constituent component of the applied non-constraint force,

$\sum_{i}(m_{i}{\ddot{x}}_{i}-F^{(a)}_{i})\delta x_{i}=0$ (38)

for all variations $\\{\delta x_{i}\\}$ for which the constraints are
maintained at each fixed time $t$;

$\sum_{i}\frac{\partial g^{\alpha}(x,t)}{\partial x_{i}}\delta x_{i}=0.$ (39)

In terms of Lagrange multipliers $\lambda_{\alpha}(t,x)$, one for each
constraint $g_{\alpha}$, we have then

$m_{i}{\ddot{x}}_{i}-F^{(a)}_{i}-\lambda_{\alpha}\frac{\partial
g^{\alpha}(x,t)}{\partial x_{i}}=0,$ (40)

with summation over the repeated index $\alpha$. Thus the solution $t\mapsto
x(t)$ solves these equations and satisfies the constraints. In the case where
the constraints are _time-independent_ , the power expended by the constraint
forces is then

$(m_{i}{\ddot{x}}_{i}-F^{(a)}_{i}){\dot{x}}_{i}=\lambda_{\alpha}\frac{dg^{\alpha}\bigl{(}x(t)\bigr{)}}{dt},$
(41)

which is $0$ if the solution $t\mapsto x(t)$ satisfies the constraints.
Assuming that the applied force comes from a potential, the equation (41) can
be obtained as the Euler-Lagrange equations from an action of the form $\int
L\bigl{(}t,x(t),{\dot{x}}(t),\lambda(t)\bigr{)}\,dt$. For non-holonomic
dynamics, with velocity-dependent constraints (and/or velocity dependent
applied forces) this classical method generally fails. As already discussed,
our framework provides an action functional, explicitly involving the initial
configuration of the system, whose E-L equations coincide with the equations
of motion. In this section we consider a non-holonomic system, with equations
of motion given in (42), and work through our framework and obtain the dual
action functional (49) for this system.

Pars [Par54] studied two examples of non-holonomic systems in relation to
variational principles. One such system consists of the motion of a single
particle in 3-d space without any applied forces but subject to the
requirement that the path $t\mapsto x(t)$ must respect the constraint:

$x_{3}d{x}_{1}-d{x}_{2}=0.$

In particular, the velocity $\dot{x}$ and position coordinates $x$ must
satisfy:

$x_{3}\dot{x}_{1}-\dot{x}_{2}=0.$

The differential form $\eta=x_{3}dx_{1}-dx_{2}$ gives a non-holonomic
constraint, in the sense that it cannot be expressed as $\mu\,dg$, for any
constraint function $g$ and ‘intergrating factor’ $\mu$, as any such form
satisfies

$(\mu\,dg)\wedge d(\mu\,dg)=0$

whereas $\eta\wedge d\eta=-dx_{2}\wedge dx_{3}\wedge dx_{1}\neq 0$.

In the following, we will assume $C^{2}$ differentiability of all functions on
the time interval $[0,T]$. Thus, the governing system of equations is

$\displaystyle\dot{x}_{i}-v_{i}=0$ (42a) $\displaystyle\dot{v}_{i}+Q_{i}=0$
(42b) $\displaystyle x_{3}v_{1}-v_{2}=0$ (42c) $\displaystyle
x_{i}(0)=x_{i}^{0};\qquad v_{i}(0)=v^{0}_{i};\qquad
x_{3}^{0}v^{0}_{1}-v^{0}_{2}=0;\quad\mbox{(initial data consistent with
constraint)}.$ (42d)

The constraint can be written as $a_{i}v_{i}=0$, where $a=(x_{3},-1,0)$, as
well as

$n_{i}\ddot{x}_{i}+w=0$ (43)

with

$n=a=(x_{3},-1,0);\qquad w=\dot{x}_{3}\dot{x}_{1}.$

In terms of the matrix $L$ and vector $b$, $L_{13}=1,b_{2}=-1$, respectively,
with all other components of the arrays being $0$. Gauss’ principle of Least
Constraint can be used to reduce this system to a system of three equations in
three unknowns whose solutions satisfy the constraint (42c), by _assuming_

$Q_{i}:=\lambda^{*}n_{i}=\lambda^{*}a_{i}$

with $\lambda^{*}$ evaluated from substituting $Q_{i}=\lambda^{*}n_{i}$ and
(42b) into the constraint written as

$0=n_{i}\dot{v}_{i}+w=-n_{i}Q_{i}+w=-\lambda^{*}n\cdot
n+w=-\lambda^{*}x_{3}^{2}-\lambda^{*}+\dot{x}_{3}\dot{x}_{1}\Longrightarrow\lambda^{*}=\frac{v_{3}v_{1}}{(1+x_{3}^{2})}.$

Then the reduced equations of motion satisfying the constraint are:

$\ddot{x}_{1}=\frac{-\dot{x}_{1}x_{3}\dot{x}_{3}}{(1+x_{3}^{2})};\qquad\ddot{x}_{2}=\frac{\dot{x}_{1}\dot{x}_{3}}{(1+x_{3}^{2})};\qquad\ddot{x}_{3}=0.$
(44)

Given the homogeneous, linear-in-velocity constraint $a\cdot v=0$ and the
choice of the form of the constraint force, it is clear that the constraint
force of the formalism expends no power along solutions of (42), i.e.,

$Q\cdot v=0$

as well as all trajectories satisfying only the constraint (42c).

Suppose we now define a solution to (44) as
$t\mapsto\bar{x}(t),t\mapsto\bar{v}(t)$ and the corresponding constraint force
$t\mapsto\lambda^{*}(t)a(t)=\bar{Q}(t)$. For $M=3$ in (32), defining
$\mathcal{L}_{H}$ in (34) and the dual action in (37) using these barred
functions with $t\mapsto\bar{s}(t)=0$ and the functions
$x^{(H)},v^{(H)},Q^{(H)},s^{(H)}$ defined in (36) replacing $x,v,Q,s$, we have
the result that the Euler-Lagrange equations of the action (37) recover (42)
(and $Q\cdot v=0$) with the change of variables $x\to x^{(H)},v\to
v^{(H)},Q\to,Q^{(H)}$ and that a solution to this system of equations, viewed
as equations for determining the dual state, is given by
$t\mapsto\rho(t)=0,t\mapsto\lambda(t)=0,t\mapsto\mu(t)=0,t\mapsto\Lambda(t)=0$.
In turn, $x^{(H)},v^{(H)},Q^{(H)},s^{(H)}$ in (36) defined in terms of this
dual set of functions solve the Pars system (42).

Let us now consider a constraint force of the form

$Q=\lambda^{*}\,a(x)+a(x)\times v+\nu\,v;\qquad\nu\geq 0,\mbox{ a scalar
constant}.$ (45)

Then, it is reasonable to expect solutions $(t\mapsto(x(t),v(t),Q(t))$ to the
system (42) with this particular form of $Q$ by considering the system

$\displaystyle\dot{x}_{i}-v_{i}=0$ (46a)
$\displaystyle\dot{v}_{i}=-\lambda^{*}n_{i}(x)-\varepsilon_{ijk}n_{j}(x)v_{k}-\nu
v_{i}$ (46b) $\displaystyle x_{3}v_{1}-v_{2}=0$ (46c)
$\displaystyle(\lambda^{*}n_{i}(x)+\varepsilon_{ijk}n_{j}(x)v_{k}+\nu
v_{i})v_{i}-\frac{1}{2}s^{2}=0$ (46d) $\displaystyle x_{i}(0)=p_{i}^{0};\qquad
v_{i}(0)=v^{0}_{i};\qquad x_{3}(0)v^{0}_{1}-v^{0}_{2}=0;\quad\mbox{(initial
data consistent with constraint)},$ (46e)

and obtaining a solution $(t\mapsto(x(t),v(t),\lambda^{*}(t))$ to the
constrained ivp given by (46a, 46b, 46c, 46e), as well as a solution to (46)
by using that solution to assign the value of $s^{2}$ from (46d).

Indeed, a solution can be generated as follows: with $n\times v$ given by

$n\times v=(-\dot{x}_{3},-x_{3}\dot{x}_{3},x_{3}v_{2}+v_{1}),$

substitute (46b) (corresponding to the constraint force (45)) into (43) to
solve for $\lambda^{*}$ as

$\lambda^{*}(x,v;\nu)=\frac{v_{3}v_{1}+\nu(v_{2}-x_{3}v_{1})}{x_{3}^{2}+1}.$

With this expression of $\lambda^{*}(x,v;\nu)$, solve the system of ODE (46a,
46b, 46e). It can be checked that, by construction, the obtained solution
satisfies the constraint (46c), since the initial data
$\left(x^{0},v^{0}\right)$ is required to satisfy (46e).

Now, (46c) is simply the expression $n_{i}(x(t))v_{i}(t)=0$, so that (46d)
becomes

$2\nu\,v_{i}v_{i}=s^{2}$

and assigning $t\mapsto s(t)$ to be any of $\pm\sqrt{2\nu\,v(t)\cdot v(t)}$
satisfies (46d).

For $\nu=0$, the constraint force expends no power but does not satisfy the
d’Alembert requirement of being in the direction of $a(x)$; moreover, the
equations of motion are given by

$\ddot{x}_{1}=\frac{-\dot{x}_{1}x_{3}\dot{x}_{3}}{(1+x_{3}^{2})}+\dot{x}_{3};\qquad\ddot{x}_{2}=\frac{\dot{x}_{1}\dot{x}_{3}}{(1+x_{3}^{2})}+x_{3}\dot{x}_{3};\qquad\ddot{x}_{3}=-x_{3}\dot{x}_{2}-\dot{x}_{1}.$
(47)

These equations cannot yield the same motion as (44) from identical, but
generic, initial conditions, even for short times (it suffices to note the
evolution of $x_{3}$ for the two cases).

Of course, for all choices of the constraint force, our scheme presented in
Sec. 3-4 generates an action principle for the corresponding system of
governing constrained ODE. For example, consider (44) and denote

$\displaystyle
x_{3}(t)=v^{0}_{3}t+x_{3}^{0}=:z(t);\qquad\dot{v}_{1}-\frac{v_{1}zv_{3}^{0}}{1+z^{2}}=0;\qquad\dot{v}_{2}-\frac{v_{1}v_{3}^{0}}{1+z^{2}}=0,$
(48) $\displaystyle v_{1}(0)=v_{1}^{0};\qquad v_{2}(0)=v_{2}^{0}.$

Next, with the choice
$H(v_{1},v_{2})=\frac{1}{2}\left(v_{1}^{2}+v_{2}^{2}\right)$, form the pre-
dual functional

$\displaystyle\widehat{S}_{H}[v,\lambda]$
$\displaystyle=\int_{0}^{T}\mathcal{L}_{H}\left(v,\lambda,\dot{\lambda}\right)\,dt-\lambda_{1}(0)v_{1}^{0}-\lambda_{2}(0)v_{2}^{0}$
$\displaystyle\mbox{where
}\mathcal{L}_{H}\left(v,\lambda,\dot{\lambda}\right)$
$\displaystyle=-v_{1}\dot{\lambda}_{1}-\frac{\lambda_{1}v_{1}zv_{3}^{0}}{1+z^{2}}-v_{2}\dot{\lambda}_{2}-\frac{\lambda_{2}v_{1}v_{3}^{0}}{1+z^{2}}+\frac{1}{2}\left(v_{1}^{2}+v_{2}^{2}\right).$

The DtP mapping $v^{(H)}\left(\lambda,\dot{\lambda}\right)$ is generated from
solving the equations

$\displaystyle\frac{\partial\mathcal{L}_{H}}{\partial v_{1}}=0$
$\displaystyle:\qquad
v_{1}^{(H)}\left(\lambda,\dot{\lambda}\right)=\dot{\lambda}_{1}+\frac{\lambda_{1}zv_{3}^{0}}{1+z^{2}}+\frac{\lambda_{2}v_{3}^{0}}{1+z^{2}}$
$\displaystyle\frac{\partial\mathcal{L}_{H}}{\partial v_{2}}=0$
$\displaystyle:\qquad
v_{2}^{(H)}\left(\lambda,\dot{\lambda}\right)=\dot{\lambda}_{2}.$

Then the _dual action_ is given by
$S_{H}[\lambda]=S_{H}[v^{(H)}[\lambda],\lambda]$:

$\displaystyle S_{H}[\lambda]$
$\displaystyle=-\frac{1}{2}\int_{0}^{T}\left(\left(v^{(H)}_{1}\right)^{2}+\left(v^{(H)}_{1}\right)^{2}\right)\,dt-\lambda_{1}(0)v_{1}^{0}-\lambda_{2}(0)v_{2}^{0}$
(49)
$\displaystyle=-\frac{1}{2}\int_{0}^{T}\left(\left(\dot{\lambda}_{1}+\frac{\lambda_{1}zv_{3}^{0}}{1+z^{2}}+\frac{\lambda_{2}v_{3}^{0}}{1+z^{2}}\right)^{2}+\left(\dot{\lambda}_{2}\right)^{2}\right)\,dt-\lambda_{1}(0)v_{1}^{0}-\lambda_{2}(0)v_{2}^{0},$
$\displaystyle\qquad\mbox{with }\lambda_{1}(T),\lambda_{2}(T)\mbox{ specified
arbitrarily,}$

along with the guarantee that its E-L equation and natural boundary conditions
are

$\qquad\dot{\overline{{v}^{(H)}_{1}}}-\frac{v^{(H)}_{1}zv_{3}^{0}}{1+z^{2}}=0;\qquad\dot{\overline{{v}^{(H)}_{2}}}-\frac{v^{(H)}_{1}v_{3}^{0}}{1+z^{2}}=0;\qquad
v^{(H)}_{1}(0)=v_{1}^{0};\qquad v^{(H)}_{2}(0)=v_{2}^{0},$

with the function $z$ defined in (48).

The closed set of dual E-L systems can be written down for (47) as well as the
case when $Q$ is left free, but this is not instructive. The main point of the
above two examples is to show through explicit computation of a non-trivial
example that, in a formulation in which an additive constraint force is added
to the equations of unconstrained motion to account for the presence of some
kinematic constraints subject only to the power of constraint forces being
non-negative, uniqueness of solutions, even for short times, is not to be
expected.

##  5 A resolution

The apparent non-uniqueness in the nature of the motion (as well as the
constraint force), as evidenced from the examples considered in this paper, is
certainly disconcerting (but not entirely unanticipated, cf. Sec. 2).
Formulations with equal number of variables and equations have been
demonstrated, but uniqueness of solutions can still be an issue, even with
restrictions imposed on working of constraint forces that ensure they are
dissipative or conservative. Within this general setting, some formulations
may deliver uniqueness but then the question of which formulation is ‘best,’
and why, needs to be explored. In this context, the variational principles we
suggest allows the use of different choices of $H$ \- both through its
functional form and its dependence on the ‘base states,’ the collection of
specified functions of time specified by overhead bars - as a selection
criterion for solutions to the primal problem of constrained particle
dynamics.

In closing, we consider the following conceptually minimal generalities free
of ad-hoc physical choices, but not necessarily implementable in practical
terms, especially for analytical work in terms of explicit formulae. For the
constrained system of $2d.N$ degrees of freedom (including both position and
velocities) with $M$ constraints written in the form

$\displaystyle\dot{x}_{i}$ $\displaystyle=v_{i},\qquad i=1,\ldots,d.N$ (50)
$\displaystyle\dot{v}_{i}$ $\displaystyle=F_{i}(x,v,t)+F_{i}^{(c)},\qquad
i=1,\ldots,d.N$ $\displaystyle 0$
$\displaystyle=g^{\alpha}(x,v,t),\qquad\alpha=1,\ldots,M$ $\displaystyle
v_{i}^{0}$ $\displaystyle=v_{i}(0);x_{i}^{0}=x_{i}(0)\qquad\mbox{specified
satisfying }0=g^{\alpha}\left(x^{0},v^{0},0\right),$

assume, as an application of the implicit function theorem, that, in a small
neighborhood of the initial condition at least, the array $v$ can be split
into two parts and written, by renumbering if necessary, as

$v=\left(v^{(r)},v^{(s)}\right)$

such that there exists a function

$\hat{v}^{(s)}\left(x,v^{(r)},t\right)$

satisfying

$g^{\alpha}\left(x,v^{(r)},\hat{v}^{(s)}\left(x,v^{(r)},t\right),t\right)=0.$

The number of degrees of freedom in the splitting of $v$ depends on the
maximal rank of the matrix $\frac{\partial g^{\alpha}}{\partial
v_{i}}(x^{0},v^{0},0)$; let this maximal rank be $K\leq M$. We include
holonomic constraints in the treatment by formally differentiating such
constraints w.r.t time and considering them in ‘rate’ form. The array
$v^{(s)}$ may be indexed as $v^{(s)}_{i_{s}},i_{s}=1,\ldots,K$. Assuming we
have reordered the list $((x_{i},v_{i}),i=1,\ldots,d.M)$ as

$\left(x^{(s)},x^{(r)},v^{(s)},v^{(r)}\right):=\left(\left(x_{i_{s}},x_{i_{r}},v_{i_{s}},v_{i_{r}}\right)|\quad
i_{s}=1,\ldots K,\quad i_{r}=K+1,\ldots,d.M\right),$

and using the notation

$x=\left(x^{(s)},x^{(r)}\right),\quad F=\left(F^{(s)},F^{(r)}\right),\quad
F^{(c)}=\left(F^{(c)(s)},F^{(c)(r)}\right),$

one solves the reduced system

$\displaystyle\dot{x}^{(r)}(t)$ $\displaystyle=v^{(r)}(t)$ (51)
$\displaystyle\dot{x}^{(s)}(t)$
$\displaystyle=\hat{v}^{(s)}\left(x(t),v^{(r)}(t),t\right)$
$\displaystyle\dot{v}^{(r)}(t)$
$\displaystyle=F^{(r)}\left(x(t),\hat{v}^{(s)}\left(x(t),v^{(r)}(t),t\right),v^{(r)}(t),t\right)$
$\displaystyle\left(x,v^{(s)},v^{(r)}\right)(0)$
$\displaystyle=\left(x^{0},v^{(s)0},v^{(r)0}\right)\quad\mbox{specified
satisfying }$ $\displaystyle g^{\alpha}\left(x^{0},v^{(s)0},v^{(r)0},0\right)$
$\displaystyle=0,\qquad\alpha=1,\ldots,M$

and obtains the solution to the full system by the following evaluation of the
constraint forces:

$\displaystyle F^{(c)(r)}(t)$ $\displaystyle=0$ $\displaystyle F^{(c)(s)}(t)$
$\displaystyle:=\left.\frac{d}{dt}\left(\hat{v}^{(s)}\left(x(\cdot),v^{(r)}(\cdot),\cdot\right)\right)\right|_{t}-F^{(s)}\left(x(t),\hat{v}^{(s)}\left(x(t),v^{(r)}(t),t\right),v^{(r)}(t),t\right)$

(in the neighborhood where $\hat{v}^{(s)}$ is defined). In order to compute
orbits involving states not contained in the domain, say $\Omega_{0}$, of the
function $\hat{v}^{s}$ corresponding to the specified initial condition in
(50), the ‘initial condition’ can be reset to the attained state when orbits
reach states at or near the boundary of $\Omega_{0}$, and the scheme above
repeated.

The above considerations show the minimal constraint reactions that can be in
play in a constrained particle system. No extraneous, ad-hoc, assumptions have
been made here about the nature of the constraint forces and, under $C^{1}$
smoothness of the functions $F$ and $\hat{v}^{(s)}$, local in time uniqueness
of solutions is also expected.

We make two observations:

  * •

Generating the function $\hat{v}^{(s)}$ is a non-trivial matter, but may not
be impenetrable in a computational approximation scheme, although definitely
computationally expensive. Within a computational setting, the local nature of
the definition of $\hat{v}^{(s)}$ is also not a fundamental barrier to
computation of approximate global trajectories of the constrained system.

  * •

The methodology for generating dual action principles we have proposed applies
seamlessly to the system (51).

##  6 Concluding remarks

We have developed a method that takes as input a broad class of equations of
motion, including those with kinematical constraints that are nonlinear in the
component velocities as well as those involving non-conservative forces, and
transforms them into a dual system that is governed by a Lagrangian depending
on a functional parameter. The dual system can, in many cases, be reformulated
in a Hamiltonian framework, giving rise to constants of motion. Our method is
applicable to PDE systems and, we expect, should make it possible to extend
the method of path integrals to a wider class of physical systems. More
speculatively, we hope to connect with astrophysical applications as in
[GR06].

## Acknowledgments

This work was supported by the Simons Pivot Fellowship grant # 983171. We
thank G. J. Wang for bringing the paper [DM96] to our attention.

## References

  * [Ach22] Amit Acharya.  Variational principles for nonlinear PDE systems via duality.  Quarterly of Applied Mathematics, LXXXI:127–140, 2023, Article electronically published on September 26, 2022.
  * [Ach23] Amit Acharya.  A dual variational principle for nonlinear dislocation dynamics.  Journal of Elasticity, https://doi.org/10.1007/s10659-023-09998-5, 2023.
  * [BFL+11] B. M. Boghosian, L. M. Fazendeiro, J. Laett, H. Tang, and P. V. Coveney.  New variational principles for locating periodic orbits of differential equations.  Philosophical Transactions of the Royal Society A: Mathematical, Physical and Engineering Sciences, 369(1944):2211–2218, 2011.
  * [BMZ05] Anthony M. Bloch, J. E. Marsden, and D. V. Zenkov.  Nonholonomic dynamics.  Notices Amer. Math. Soc., 52(3):324–333, 2005.
  * [Bre18] Y. Brenier.  The initial value problem for the Euler equations of incompressible fluids viewed as a concave maximization problem.  Communications in Mathematical Physics, 364(2):579–605, 2018.
  * [Bre20] Y. Brenier.  Examples of hidden convexity in nonlinear PDEs.  https://hal.science/hal-02928398/document, 2020.
  * [CR07] C. Cronström and T. Raita.  On the existence of Hamiltonians for non-holonomic systems.  arXiv preprint arXiv:0711.3997, 2007.
  * [CR09] C. Cronström and T. Raita.  On nonholonomic systems and variational principles.  Journal of Mathematical Physics, 50(4):042901, 2009.
  * [Del20] D. H. Delphenich.  Selected papers on Gauss’s Principle of Least Constraint.  Neo-Classical Press, Spring Valley, Ohio, 2020.
  * [DM96] C. P. Dettmann and G. P. Morriss.  Hamiltonian formulation of the Gaussian isokinetic thermostat.  Physical Review E, 54(3):2495, 1996.
  * [EHF+83] D. J. Evans, W. G. Hoover, B. H. Failor, W. Moran, and A. J. C. Ladd.  Nonequilibrium molecular dynamics via Gauss’s principle of least constraint.  Physical Review A, 28(2):1016, 1983.
  * [FKNU05] Y. Y. Fan, R. E. Kalaba, H. H. Natsuyama, and F.E. Udwadia.  Reflections on the Gauss principle of least constraint.  Journal of Optimization Theory and Applications, 127:475–484, 2005\.
  * [Fla05] M. R. Flannery.  The enigma of nonholonomic constraints.  American Journal of Physics, 73(3):265–272, 2005.
  * [Gal13] C. R. Galley.  Classical mechanics of nonconservative systems.  Physical Review Letters, 110(17):174301, 2013.
  * [Gau29] Carl Friedrich Gauß.  Über ein neues allgemeines grundgesetz der mechanik.  Journal für die reine und angewandte Mathematik, 4:232–235, 1829\.
  * [Gol57] H. Goldstein.  Classical Mechanics, $5^{th}$ ed. Addison-Wesley Publsihing Company, Inc., 1957.
  * [GR06] W. D. Goldberger and Ira Z. Rothstein.  Dissipative effects in the worldline approach to black hole dynamics.  Physical Review D, 73(10):104030, 2006.
  * [Gri83] Phillip A. Griffiths.  Exterior differential systems and the calculus of variations. Birkhäuser, Boston, Mass., 1983.
  * [Her10] H. Hertz.  Die Prinzipien der Mechanik in neuem Zusammenhange dargestellt.  J. A. Barth, 1910, Germany, 1910.
  * [KA23] Uditnarayan Kouskiya and Amit Acharya.  Hidden convexity in the heat, linear transport, and Euler’s rigid body equations: A computational approach.  to appear in Quarterly of Applied Mathematics, https://arxiv.org/abs/2304.09418, 2023.
  * [KA24] Uditnarayan Kouskiya and Amit Acharya.  Inviscid burgers as a degenerate elliptic problem.  to appear in Quarterly of Applied Mathematics, https://arxiv.org/abs/2401.08814,2024.
  * [LC04] Y. Lan and P. Cvitanović.  Variational method for finding periodic orbits in a general flow.  Physical Review E, 69(1):016217, 2004.
  * [LM95] A. D. Lewis and R. M. Murray.  Variational principles for constrained systems: theory and experiment.  International Journal of Non-Linear Mechanics, 30(6):793–815, 1995\.
  * [Lor63] E. N. Lorenz.  Deterministic nonperiodic flow.  Journal of Atmospheric Sciences, 20(2):130–141, 1963.
  * [MD98] G. P. Morriss and C. P. Dettmann.  Thermostats: analysis and application.  Chaos: An Interdisciplinary Journal of Nonlinear Science, 8(2):321–336, 1998.
  * [Par54] L. A. Pars.  Variational principles in dynamics.  Journal of Mechanics and Applied Mathematics, 7(3):338–351, 1954\.
  * [Rab78] P. H. Rabinowitz.  A variational method for finding periodic solutions of differential equations.  In Nonlinear evolution equations (Proc. Sympos., Univ. Wisconsin, Madison, Wis., 1977), volume 40 of Publ. Math. Res. Center Univ. Wisconsin, pages 225–251. Academic Press, New York-London, 1978\.
  * [Ros91] R. Rosenberg.  Analytical dynamics of discrete systems.  Plenum Press, New York and London, 1991.
  * [SGA24] Siddharth Singh, Janusz Ginster, and Amit Acharya.  A hidden convexity of nonlinear elasticity.  https://arxiv.org/abs/2401.08538, 2024.