Action principles for dissipative, non-holonomic Newtonian mechanics

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.


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.

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 , . . ., x dA ) ∈ R d denoting the configuration coordinates of the A-th particle, and f A = (f 1A , . . ., f dA ) denoting the force divided by the mass of the A-th particle, the equation of motion for the system is: where x = (x A ) A=1,...,N and f = (f A ) A=1,...,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 : R dN → 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 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 with t running over an interval in R containing the initial time 0, that satisfies (1), subject to prescribed initial conditions Uniqueness of solutions t → 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 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 and if one now invokes the additional requirement that the constraint forces expend no power along actual trajectories, i.e., for functions x(•) that satisfy the equations of motion and constraints then it is clear that for N ≥ 1, d > 2, while for t → λ(t) ∈ 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 where the R dN -valued function a, defined on the configuration space R dN , takes the place of ∇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 λ (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 linearin-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) • ẋ = 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) → α(i, A) = (A − 1)d + i, x ia = q α(i,A) , e.g.) whereas the constraint force from the d'Alembert treatment could be interpreted simply as 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 α qα = 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 where n and w are smooth real-valued functions on R dN × R dN × R; this encompasses, holonomic and anholonomic constraints considered earlier.The constraint applies to any trajectory satisfying the equation of motion (1a); hence, Assuming the constraint force to be of the form for λ a scalar valued function, one obtains the representation 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), is an additional, eminently physical, assumption.When a Lagrange multiplier formulation is employed to solve the optimization problem, the critical point is characterized by solutions of 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 α i n β i in 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 → 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.

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 nonholonomic 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 , . . ., x dA ) ∈ R d and velocity denoted (v 1A , . . ., v dA ).The following constrained differential-algebraic system is considered, with f, g, W being suitably smooth given functions of their arguments: (power expended by constraint forces) t → s(t) ∈ R is a slack variable to convert an inequality constraint into an equality 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 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 where f (c) represents the constraint force.
Examples of the constraint forces for the various models discussed are no assumption made on the form of the constraint force (13a) 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 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 Q; motivated by this, we allow for a general nonlinear dependence on Q of the constraint force, which requires, within our formalism, to introduce a new variable which we denote by R := Q and include it within our governing set of equations -this is no different, in principle, from the introduction of a new variable v := ẋ in classical mechanics.
The totality of the primal physical functions whose evolution is of interest are the components of U can be thought of as coordinate functions on an open subset of R 2dN +2 Γ +1 , and we are interested in motions t → U (t) that satisfy the conditions (11).
We introduce an array of dual functions of time to be used subsequently, as well as the array Technically, the components of D can be thought of as coordinate functions on R 2dN +m+1+ Γ .Then define the pre-dual functional where on the left the argument (x, v, Q, ρ, λ, µ, Λ, κ) is any smooth function [0, T ] → R 2dN + Γ +2dN +m+ Γ , and the integrand on the right is and 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 i.e., ∂L H ∂U (U, D, t) = 0 should be solvable for U in terms of (D, t).To be more precise, here H is a real-valued smooth function on (R 2dN +2 Γ +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 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 with Dirichlet boundary conditions λ(T ), ρ(T ), κ(T ), κ(0) specified arbitrarily. ( Note that, on the left, [ρ, λ, µ, Λ, κ] is any smooth function [0, T ] → R 2dN +m+1+ Γ .Utilizing the requirement (21), and the boundary conditions δλ(T ) = 0, δρ(T ) = 0, δκ(T ) = 0, δκ(0) = 0, the first variation of S H is given by (using the short-hand and slight abuse of notation Apart from a direct computation, this may also be understood by observing that L H (U, D, t) is affine in 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 H) , W (H) .
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 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 where (x, v, Q) runs over R 3dN here.
The key requirement (21) now becomes the condition that the following algebraic system of equations to define the functions (D, t) → x (H) (D, t) Suppose, now, this can be done (and definitely in a perturbative sense, by choosing c x , c v , c q → +∞).
Then the following remarks are in order: 1. Let (t → x(t), t → v(t), t → 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 (•) functions, the dual variational principle has a critical point, explicitly given by 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 (x, v, 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. 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. For c x , c v , c Q large (and definitely when the primal system is linear), 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 (ρ(T ), λ(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 ( ρ(T ), λ(T )), respectively, and specifying the values of (ρ(T ), λ(T )) leaves the derivatives adjustable to the demands of solving the primal ivp.

Defining a function
then it is an elementary result of Hamiltonian mechanics that a Hamiltonian H(D, P, t) := P • R(D, P, t) − L(D, 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., Moreover, when L is independent of t (certainly when the H chosen to define L H is independent of t, e.g., by not involving any of the overhead ¯functions), then 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. The dual action generated above for a dissipative system by our scheme and, in general, for a system of ODE of the type ẋ = 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 = 2πt T ; P := T 2π one seeks to find 2π-periodic solutions (x, P) of the system The problem can now be approached as looking for critical points of dual functionals developed by our scheme on the fixed interval [0, 2π] 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' 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) 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.
The DtP mapping U (H) (D) is generated from solving the equations with the matrix defined as and the vectors as with solution given by 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 Noting, further, that we obtain the corresponding dual Lorenz action as: with λ(T ), µ(T ), γ(T ) specified arbitrarily. (31) Our formalism and results guarantee that the corresponding E-L equation and natural boundary conditions lead to x (H) t=0 = x 0 ; y (H) t=0 = y 0 ; z (H) t=0 = z 0 .We note from (30) that if (t → x(t), t → ȳ(t), t → z(t)) is a solution of the Lorenz system then the corresponding dual Lorenz action has a critical point given by (t → λ(t) = 0, t → µ(t) = 0, t → γ(t) = 0).In passing, we also note that for (x = 0, ȳ = 0, z = 0) the dual Lagrangian is a negative semi-definite form for c ≫ 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 = (λ, µ, γ)) a Hamiltonian H(D, P ) := R(D, P, t) can be associated with the dissipative Lorenz evolution, and it is a constant of the Lorenz motion when Ū is constant (in the sense of the remark 4).

Generalization of the example of Pars [Par54]
We consider the system (with all particle masses assumed equal to 1 in appropriate units) with initial conditions 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., In case external forces and forces of interaction arising from a potential E(x) were to be involved to give the equation of motion the corresponding power of constraints statement would be an expression of non-negative mechanical dissipation given by we have where The DtP mapping is obtained from solving the system We now assume that the functions µ 2 and Λ 2 are bounded above in [0, T ] so that c x , c v , c Q can be chosen large enough to make the matrix Then the solution of (35) is given by (36) The dual functional is then given by with ρ(T ), λ(T ) specified (arbitrarily).

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, defined 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 ẍi = F (a) i (x, t) and constraints g α (x, t) = 0, for α running over a finite set of indices, and F (a) i are the applied forces.Our interest is mainly in non-holonomic, possibly dissipative, systems, specifically those where the constraints are of the form g α (x, v, t) = 0 and forces of the form F (a) i (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, i for all variations {δx i } for which the constraints are maintained at each fixed time t; In terms of Lagrange multipliers λ α (t, x), one for each constraint g α , we have then with summation over the repeated index α.Thus the solution t → 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 which is 0 if the solution t → 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 L t, x(t), ẋ(t), λ(t) 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 → x(t) must respect the constraint: In particular, the velocity ẋ and position coordinates x must satisfy: The differential form η = x 3 dx 1 − dx 2 gives a non-holonomic constraint, in the sense that it cannot be expressed as µ dg, for any constraint function g and 'intergrating factor' µ, as any such form satisfies (µ dg) ∧ d(µ dg) = 0 whereas η ∧ dη = −dx 2 ∧ dx 3 ∧ dx 1 = 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 ẋi The constraint can be written as a i v i = 0, where a = (x 3 , −1, 0), as well as with n = a = (x 3 , −1, 0); w = ẋ3 ẋ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 with λ * evaluated from substituting Q i = λ * n i and (42b) into the constraint written as .
Then the reduced equations of motion satisfying the constraint are: Given the homogeneous, linear-in-velocity constraint a • 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 • v = 0 as well as all trajectories satisfying only the constraint (42c).
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 Next, with the choice The DtP mapping v (H) λ, λ is generated from solving the equations Then the dual action is given by S along with the guarantee that its E-L equation and natural boundary conditions are v(H) 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.

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 ẋi = v i , i = 1, . . ., d.N vi = F i (x, v, t) + F (c) i , i = 1, . . ., d.N 0 = g α (x, v, t), α = 1, . . ., M v 0 i = v i (0); x 0 i = x i (0) specified satisfying 0 = g α x 0 , v 0 , 0 , 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 = v (r) , v (s)   such that there exists a function v(s) x, v (r) , t satisfying g α x, v (r) , v(s) x, v (r) , t , t = 0.
The number of degrees of freedom in the splitting of v depends on the maximal rank of the matrix ∂g α ∂v i (x 0 , v 0 , 0); let this maximal rank be K ≤ 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) is , i s = 1, . . ., K. Assuming we have reordered the list ((x i , v i ), i = 1, . . ., d.M ) as x (s) , x (r) , v (s) , v (r) := ((x is , x ir , v is , v ir ) | i s = 1, . . .K, i r = K + 1, . . ., d.M ) , and using the notation x = x (s) , x (r) , F = F (s) , F (r) , F (c) = F (c)(s) , F (c)(r) , one solves the reduced system ẋ(r) (t) = v (r) (t) ẋ(s) (t) = v(s) x(t), v (r) (t), t v(r) (t) = F (r) x(t), v(s) x(t), v (r) (t), t , v (r) (t), t x, v (s) , v (r) (0) = x 0 , v (s)0 , v (r)0 specified satisfying g α x 0 , v (s)0 , v (r)0 , 0 = 0, α = 1, . . ., M and obtains the solution to the full system by the following evaluation of the constraint forces: s) x(t), v(s) x(t), v (r) (t), t , v (r) (t), t (in the neighborhood where v(s) is defined).In order to compute orbits involving states not contained in the domain, say Ω 0 , of the function vs 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 Ω 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 v(s) , local in time uniqueness of solutions is also expected.
We make two observations: • Generating the function 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 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).

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].