format: text

text:
# MECHANICAL SYSTEMS SUBJECTED TO GENERALIZED NONHOLONOMIC CONSTRAINTS

## 1 Introduction

Mechanical systems subjected to nonholonomic constraints have received a lot of attention in recent years in the literature of Geometric Mechanics (see  and references therein). Indeed, the dynamics of nonholonomic mechanics have been described from several approaches: Hamiltonian, Lagrangian and even Poisson methods have been used.

The constraints which are usually considered in the literature (both linear and nonlinear) satisfy a certain regularity condition. That is, they are given by a set of independent nonholonomic constraint functions, or, in a global description, by a distribution on the configuration manifold in the linear case or a submanifold of its tangent bundle in the case of nonlinear constraints.

However, there is an increasing interest in engineering and robotics in the motion of special mechanical systems as, for example, dynamical devices that locomote with the enviroment via impacts, sudden changes of phase space, etc. In many cases, the jump of the system’s velocity is produced by an impulse that enforces new constraints on the system. In some cases, these systems admit a nice mathematical modelling.

In this paper, we are interested in the following situation. In a local description, given a set of constraints $`\{\mathrm{\Phi }_1,\mathrm{},\mathrm{\Phi }_m\}`$, one assumes that they become linearly dependent at some points. In a global picture, the constraints are given by a generalized codistribution with variable rank. One could think of simple examples that exhibit this kind of behaviour. For instance, imagine a rolling ball on a surface which is rough on some parts but smooth on the rest. On the rough parts, it will roll without slipping and, hence, nonholonomic linear constraints will be present. However, when the sphere reaches a smooth part, these constraints will disappear.

The first (and, up to our knowledge, the unique) reference in a geometrical context for such kind of systems is . In that paper, the constraints are provided by a set of global 1-forms on the configuration manifold and, using the Frobenius theorem, the authors gave a classification of them according to the existence of some special sets that can exert a big influence on the trajectories of the system. In particular, the existence of an integral manifold gives a sort of partial holonomicity with strong implications. The authors were mainly motivated by problems in motion planning. However, we are interested, at least in this first approach, in the geometrical and topological aspects of the problem. In other words, we are concerned with obtaining the dynamical laws that govern the motion of the system.

In consequence, we consider nonholonomic constraints given by a generalized codistribution $`D`$, that is, a codistribution which does not necessarily have the same rank at all points in the configuration manifold. This approach leads us to the definition of the concepts of regular and singular points. It should be noticed that our definitions are slightly different from those in . Indeed, the regular points are those where the codistribution has locally constant rank. In this sense, the generalized codistribution is a regular codistribution, as is commonly understood, on the connected components of the set of regular points. The singular points are those where the codistribution changes its rank. From a dynamical perspective, the situation on the regular points is already known: we can derive the equations of motion following d’Alembert principle and treat them making use of the well-developed theory for nonholonomic Lagrangian systems.

However, on the singular points the matter is essentially different. The classical derivation of the equations of motion no longer works and we must solve the problem with other methods. Here, we have adopted a point of view strongly inspired by the theory of impulsive mechanics  and we use Newton’s second law in its integral form . Analyzing the trajectories which cross the singular set, we have found that, in certain cases, the constraint force can transmit an impulse to the motion. It is precisely the sudden appearence of new constraints (that is, the change of rank of the codistribution) which induces this impulsive character. More precisely, given a motion $`q(t)`$ crossing the singular set at time $`t_0`$, we define two vector subspaces of $`T_{q(t_0)}^{}Q`$ as follows: $`D_{q(t_0)}^{}`$ is the limit of all the 1-forms in the codistribution based on $`q(t)`$, $`t<t_0`$, when $`tt_0^{}`$. $`D_{q(t_0)}^+`$ is defined analogously. Our conclusion is that there exists a jump of momentum only if $`D_{q(t_0)}^+`$ is not contained in $`D_{q(t_0)}^{}`$ and the “pre-impact” momentum $`p(t_0)_{}`$ does not satisfy the constraints imposed by $`D_{q(t_0)}^+`$. In such a case, we propose that the jump is determined by $`\mathrm{\Delta }p(t_0)D_{q(t_0)}^+`$ and the condition that the “post-impact” momentum $`p(t_0)_+`$ must satisfy the constraints imposed by $`D_{q(t_0)}^+`$.

To find out about the relation between the theory developed here and the Hamiltonian theory of impact  applied to this problem, which seems to be a promising possibility, is the object of current research.

The paper is organized as follows. In Section 2, we introduce the notion of generalized codistribution, which is just the geometrization of constraints with non-constant rank. In Section 3 we review the theory of impulsive forces and impulsive constraints. The ideas exposed here will be helpful in understanding the developments of Section 4, which constitutes the main contribution of this paper, where we study the equations of motion for mechanical systems subjected to generalized constraints and we derive the equations describing the jump of momenta. Finally, in Section 5, some examples are discussed with detail.

## 2 Generalized codistributions

We introduce here the notion of a generalized codistribution. This notion will be helpful in subsequent sections to model geometrically the dynamical systems under consideration, that is, systems subjected to constraints which can “degenerate” at certain points. All the results in this section are adapted from the ones stated for generalized distributions in .

By a generalized codistribution we mean a family of linear subspaces $`D=\{D_q\}`$ of the cotangent spaces $`T_q^{}Q`$. Such a codistribution is called differentiable if $`q\text{Dom}D`$, there is a finite number of differentiable local 1-forms $`\omega _1,\mathrm{},\omega _l`$ defined on some open neighbourhood $`U`$ of $`q`$ such that $`D_q^{}=\text{span}\{\omega _1(q^{}),\mathrm{},\omega _l(q^{})\}`$ for all $`q^{}U`$.

We define the rank of $`D`$ at $`q`$ as $`\rho (q)=\text{dim}D_q`$. Given $`q_0Q`$, if $`D`$ is differentiable, it is clear that $`\rho (q)\rho (q_0)`$ in a neighbourhood of $`q_0`$. Therefore, $`\rho `$ is a lower semicontinuous function. If $`\rho `$ is a constant function, then $`D`$ is a codistribution in the usual sense.

For a generalized differentiable codistribution $`D`$, a point $`qQ`$ will be called regular if $`q`$ is a local maximum of $`\rho `$, that is, $`\rho `$ is constant on an open neighbourhood of $`q`$. Otherwise, $`q`$ will be called a singular point of $`D`$. The set $`R`$ of the regular points of $`D`$ is obviously open. But, in addition, it is dense, since if $`q_0S=QR`$, and $`U`$ is a neighbourhood of $`q_0`$, $`U`$ necessarily contains regular points of $`D`$ ($`\rho _{|U}`$ must have a maximum because it is integer valued and bounded). Consequently, $`q_0\overline{R}`$.

Note that in general $`R`$ will not be connected, as the following simple example shows:

###### Example 2.1

Let us consider $`Q=\text{}^2`$ and the general differentiable codistribution $`D_{(x,y)}=\text{span}\{\varphi (x)(dxdy)\}`$, where $`\varphi (x)`$ is defined by

$$\varphi (x)=\{\begin{array}{cc}0& x0\\ e^{\frac{1}{x^2}}& x>0\end{array}$$

The singular points are those of the $`y`$-axis, and the connected components of $`R`$ are the half-planes $`x>0`$ (where the rank is $`1`$) and $`x<0`$ (where the rank is $`0`$).

###### Remark 2.2

We note that the notion of singular point defined here is different from the one considered in . In that paper, the authors treat the case of generalized constraints given by a globally defined set of 1-forms, $`\omega _1,\mathrm{},\omega _l`$. Then, they consider the $`l`$-form

$$\mathrm{\Omega }=\omega _1\mathrm{}\omega _l.$$

The singular set consists of the points for which $`\mathrm{\Omega }(q)=0`$, that is, the points, $`q`$, such that $`\{\omega _1(q),\mathrm{},\omega _l(q)\}`$ are linearly dependent. Applying this notion to the former example, the set of singular points would be the half-plane $`\{x0\}`$.

Given a generalized codistribution, $`D`$, we define its annihilator, $`D^o`$, as the generalized distribution given by

$$\begin{array}{cccc}\hfill D^o:& \hfill \text{Dom}DQ& & TQ\hfill \\ & \hfill q& & D_q^o=(D_q)^o.\hfill \end{array}$$

Remark that if $`D`$ is differentiable, $`D^o`$ is not differentiable, even continuous, in general (the corresponding rank function of $`D^o`$ will not be lower semicontinuous). In fact, $`D^o`$ is differentiable if and only if $`D`$ is a regular codistribution.

We will call $`M`$ an integral submanifold of $`D`$ if $`T_mM`$ is annihilated by $`D_m`$ at each point $`mM`$. $`M`$ will be an integral submanifold of maximal dimension if

$$T_mM^o=D_m,mM.$$

In particular, this implies that the rank of $`D`$ is constant along $`M`$. A leaf $`L`$ of $`D`$ is a connected integral submanifold of maximal dimension such that every connected integral manifold of maximal dimension of $`D`$ which intersects $`L`$ is an open submanifold of $`L`$. $`D`$ will be a partially integrable codistribution if for every regular point $`qR`$, there exists one leaf passing through $`q`$. $`D`$ will be a completely integrable codistribution if there exists a leaf passing through $`q`$, for every $`qQ`$. In the latter case, the set of leaves defines a general foliation of $`Q`$. Obviously, any completely integrable codistribution is partially integrable.

$`M`$ being an integral submanifold of $`D`$ is exactly the same as being an integral submanifold of its annihilator $`D^o`$, and so on.

In Example 2.1, the leaves of $`D`$ are the half-plane $`\{x<0\}`$ and the half-lines of slope 1 in the half-plane $`\{x>0\}`$. Given any singular point, there is no leaf passing through it. Consequently, $`D`$ is not a completely integrable codistribution, but it is partially integrable.

## 3 Impulsive forces

In this section, we discuss classical mechanical systems with impulsive forces . This field has traditionally been studied by a rich variety of methods (analytical, numerical and experimental), being a meeting place among physicists, mechanical engineers and mathematicians (for an excellent overview on the subject, see ). Recently, such systems have been brought into the context of Geometric Mechanics . We will give here a brief review of the classical approach. These ideas will be useful in understanding the behaviour of the constraint forces acting on mechanical systems subject to generalized constraints. Both situations are not the same, but have many points in common, as we will see in the following.

Consider a system of $`n`$ particles in $`\text{}^3`$ such that the particle $`r`$ has mass $`m_r`$. Introducing coordinates $`(q^{3r2},q^{3r1},q^{3r})`$ for the particle $`r`$, we denote by $`Q`$ the configuration manifold $`R^{3n}`$ and by $`F_r=(F^{3r2},F^{3r1},F^{3r})`$ the resultant of all forces acting on the $`r^{th}`$ particle.

The motion of the particle $`r`$ in an interval $`[t,t^{}]`$ is determined by the system of integral equations

$$m_r(\dot{q}^k(t^{})\dot{q}^k(t))=_t^t^{}F^k(\tau )𝑑\tau ,$$
(1)
where $`3r2k3r`$ and $`k`$ is an integer. The integrals of the right-hand side are the components of the impulse of the force $`F_r`$. Equation (1) establishes the relation between the impulse and the momentum change, i.e. “impulse is equal to momentum change”. Equation (1) is a generalized writing of Newton’s second law, stated in integral form in order to allow us to consider the case of velocities with finite jump discontinuities. This is precisely the case of impulsive forces, which generate a finite non-zero impulse at some time instants.

If $`F`$ is impulsive there exists an instant $`t_0`$ such that

$$\underset{tt_0}{lim}_{t_0}^tF(\tau )𝑑\tau =P0.$$
(2)
Equation (2) implies that the impulsive force has an infinite magnitude at the point $`t_0`$, but we are assuming that its impulse $`P`$ is well defined and bounded. The expression $`P\delta (t_0)=\underset{tt_0}{lim}F(t)`$ can be mathematically seen as a Dirac delta function concentrated at $`t_0`$.

Now, we will derive the equations for impulsive motion following the discussion in . In the sequel, the velocity vector of the $`r^{th}`$ particle, $`(\dot{q}^{3r2},\dot{q}^{3r1},\dot{q}^{3r})`$, will be denoted by $`\dot{q}^r`$. Then, the system of integral equations (1) can be written as

$$m_r(\dot{q}^r(t_0+ϵ)\dot{q}^r(t_0ϵ))=_{t_0ϵ}^{t_0+ϵ}F_r(\tau )𝑑\tau .$$

If we multiply this expression by the virtual displacements at the point $`q(t_0)`$, we obtain

$$(p_r(t_0+ϵ)p_r(t_0ϵ))\delta q^r=_{t_0ϵ}^{t_0+ϵ}F_r(\tau )𝑑\tau \delta q^r.$$

For the entire system, one has

$$\underset{r=1}{\overset{n}{}}\left\{p_r(t_0+ϵ)p_r(t_0ϵ)_{t_0ϵ}^{t_0+ϵ}F_r^{}(\tau )𝑑\tau \right\}\delta q^r=\underset{r=1}{\overset{n}{}}_{t_0ϵ}^{t_0+ϵ}F_r^{\prime \prime }(\tau )𝑑\tau \delta q^r,$$
(3)
where $`F_r^{}`$ and $`F_r^{\prime \prime }`$ are, respectively, the resultant of the given forces and of the constraint reaction forces acting on the $`r^{th}`$ particle at time $`\tau `$.

Now, take a local chart $`(q^A)`$, $`1A3n`$ on a neighbourhood $`U`$ of $`q(t_0)`$ and consider the identification $`T_qQR^{3n}`$, which maps each $`v_qT_qQ`$ to $`(v^A)`$, such that $`v_q=v^A\left({\displaystyle \frac{}{q^A}}\right)_q`$, for each $`qU`$. Let us suppose that the constraints are given on $`U`$ by the 1-forms $`\omega _i=\mu _{iA}dq^A`$, $`1im`$. Then, we have that $`\mu _{iA}(q(t))=\mu _{iA}(q(t_0))+O(tt_0)`$ along the trajectory $`q(t)`$. As the virtual displacements at the point $`q(t)`$ satisfy by definition

$$\mu _{iA}(q(t))(\delta q(t))^A=0,\mathrm{\hspace{0.33em}1}im,$$

we conclude that $`\left(\mu _{iA}(q(t_0))(\delta q(t))^A+O(tt_0)\right)=0`$. Therefore, we have that

$$\delta q^r(t)=\delta q^r(t_0)+O(ϵ),t[t_0ϵ,t_0+ϵ],$$

that is, the virtual displacements at $`q(t)`$ can be approximated by the virtual displacements at $`q(t_0)`$. As a consequence, in the right-hand side of (3) we can write

$$_{t_0ϵ}^{t_0+ϵ}F_r^{\prime \prime }(\tau )𝑑\tau \delta q^r=_{t_0ϵ}^{t_0+ϵ}F_r^{\prime \prime }(\tau )\delta q^r𝑑\tau =_{t_0ϵ}^{t_0+ϵ}F_r^{\prime \prime }(\tau )\delta q^r(\tau )𝑑\tau +O(ϵ).$$

The first term after the last equality is the virtual work done by the constraint forces along the trajectory, and this work is zero since we are considering ideal constraints. The second one goes to zero as $`ϵ`$ tends to zero.

In the presence of given impulsive forces acting on $`m`$ particles, say, at time $`t_0`$, we have

$$\underset{tt_0}{lim}_{t_0}^tF_r^{}(\tau )𝑑\tau =P_r^{}0,\mathrm{\hspace{0.33em}1}r^{}m.$$

Then, taking the limit $`ϵ0`$ in (3), we obtain the equation for impulsive motion

$$\underset{r=1}{\overset{n}{}}\left\{p_r(t_0)_+p_r(t_0)_{}P_r\right\}\delta q^r=0.$$
(4)
An example in which equation (4) can be applied is when we strike with a cue a billiard ball which is initially at rest. In that case we are exerting an impulsive force that puts the ball into motion. But what happens when the ball collides with the edge of the billiard? What we see is that it bounces, i.e. it suffers again a discontinuous jump in its velocity. The constraint imposed by the wall of the billiard exerts an impulsive force on the ball. When the impulsive force is caused by constraints, such constraints are called impulsive constraints. There is a number of different situations in which they can appear. In the following, we examine them.

In the presence of linear constraints of type $`\mathrm{\Psi }=0`$, where $`\mathrm{\Psi }=b_k(q)\dot{q}^k`$ (a situation which covers the case of unilateral holonomic constraints, such as the impact against a wall, and more general types of constraints such as instantaneous nonholonomic constraints), the constraint force, $`F=F_kdq^k`$, is given by

$$F_k=\mu b_k,$$

where $`\mu `$ is a Lagrange multiplier. Then the constraint is impulsive if and only if

$$\underset{tt_0}{lim}_{t_0}^t\mu b_k𝑑\tau =P_k0,$$

for some $`k`$. The impulsive force may be caused by different circumstances: the function $`b_k`$ is discontinuous at $`t_0`$, the Lagrange multiplier $`\mu `$ is discontinuous at $`t_0`$ or both.

The presence of such constraints does not invalidate equation (4). It merely means that the virtual displacements $`\delta q^r`$ must satisfy certain additional conditions, which are just those imposed by the constraints. So, in the abscence of impulsive external forces and in the presence of impulsive constraints, we would have

$$\underset{r=1}{\overset{n}{}}\mathrm{\Delta }p_r(t_0)\delta q^r=0,$$
(5)
where $`\mathrm{\Delta }p_r(t_0)=p_r(t_0)_+p_r(t_0)_{}`$.

###### Remark 3.1

In general, equation (5) is not enough to determine the jump of the momentum. One usually needs additional physical hypothesis, related with elasticity, plasticity, etc. to obtain the post-impact momentum. In this respect, there are two classical approaches, the Newtonian approach and the Poisson approach . The Newtonian approach relates the normal component of the rebound velocity to the normal component of the incident velocity by means of an experimentally determined coefficient of restitution $`e`$, where $`0e1`$. Poisson approach divides the impact into compression and decompression phases and relates the impulse in the restitution phase to the impulse in the compression phase.

###### Remark 3.2

It could happen that impulsive constraints and impulsive forces to be present at the same time. For example, in the collision between a rigid lamina and an immobile plane surface, we must take into account not only the normal component of the contact force, but also the friction force associated to the contact. It is not innocuous the way the friction is entered into the picture. In fact, the Newton and Poisson approaches have been revealed to be physically inconsistent in certain situations. On the one hand, Newton approach can show energy gains . On the other hand, Poisson’s rule is not satisfactory since non-frictional dissipation does not vanish for perfectly elastic impacts . This surprising consequence of the impact laws is only present when the velocity along the impact surface (slip) stops or reverses during collision, due precisely to the friction. Stronge  proposed a new energetically consistent hypothesis for rigid body collisions with slip and friction. It should be noticed that the three approaches are equivalent if slip does not stop during collision and in the perfectly inelastic case ($`e=0`$).

Recently, a new Newton-style model of partly elastic impacts has been proposed  which, interestingly, always dissipates energy, unlike the classical formulation of the Newtonian approach discussed in .

In the frictionless case, one can prove the following

###### Theorem 3.3 (Carnot’s theorem)

() The energy change due to impulsive constraints is always a loss of energy.

## 4 Mechanical systems subjected to generalized constraints

In this section, we study the equations of motion for mechanical systems subjected to generalized constraints. Let us consider a mechanical system with Lagrangian function $`L:TQ\text{}`$, $`L(v)=\frac{1}{2}g(v,v)(U\tau _Q)(v)`$, where $`g`$ is a Riemannian metric on $`Q`$ and $`U`$ is a function on the configuration space $`Q`$ (the potential). Suppose, in addition, that the system is subjected to a set of constraints given by a generalized differentiable codistribution $`D`$ on $`Q`$, that is, we assume that $`\tau _Q(D)=Q`$. The motions of the system are forced to take place satisfying the constraints imposed by $`D`$.

We know that the codistribution $`D`$ induces a decomposition of $`Q`$ into regular and singular points. We write

$$Q=RS.$$

Let us fix $`R_c`$, a connected component of $`R`$. We can consider the restriction of the codistribution to $`R_c`$, $`D_c=D_{|R_c}:R_cQT^{}Q`$. Obviously, we have that $`D_c`$ is a regular codistribution, that is, it has constant rank.

Then, let us denote by $`D_c^o:R_cTQ`$ the annihilator of $`D_c`$. Now, we can consider the dynamical problem with regular Lagrangian $`L`$, subjected to the regular codistribution $`D_c^o`$ and apply the well-developed theory for nonholonomic Lagrangian systems .

Consequently, our problem is solved on each connected component of $`R`$. The situation changes radically if the motion reaches a singular point. The rank of the constraint codistribution can vary suddenly and the classical derivation of the equations of motion for nonholonomic Lagrangian systems is no longer valid. Let us explore the behaviour of the system when such a thing occurs.

Consider a trajectory of the system, $`q(t)`$, which reaches a singular point at time $`t_0`$, i.e. $`q(t_0)S`$, such that $`q(t_0ϵ,t_0)R`$ and $`q(t_0,t_0+ϵ)R`$ for sufficiently small $`ϵ>0`$. The motion along the trajectory $`q(t)`$ is governed by the following equation, which is, as in the impulsive case, an integral writing of Newton’s second law, to consider possible finite jump discontinuities in the velocities (or the momenta). That is, at each component

$$p_A(t^{})p_A(t)=_t^t^{}F_A(\tau )𝑑\tau ,$$
(6)
on any interval $`tt^{}<\mathrm{}`$, where $`F`$ is the resultant of all the forces action on the trajectory $`q(t)`$. In our case, the unique forces acting are the constraint reaction forces.

The nature of the force can become impulsive because of the change of rank of the codistribution $`D`$. We summarize the situations that can be found in Table 1. On entering the singular set, the rank of the codistribution $`D`$ at the singular point $`q(t_0)`$ can be the same as at the preceding points (Case 1) or can be lower (Cases 2 and 3). In these two latter situations, the constraints have collapsed at $`q(t_0)`$ and this induces a finite jump in the constraint force. As the magnitude of the force is not infinite, there is no abrupt change in the momenta. Consequently, in all cases, we find no momentum jumps on entering the singular set.

On leaving the singular set, the rank of $`D`$ at the posterior points can be the same as at $`q(t_0)`$ (Case 2) or can be higher (Cases 1 and 3). In Case 2 nothing special occurs. In Cases 1 and 3, the trajectory must satisfy, immediately after the point $`q(t_0)`$, additional constraints which were not present before. It is in this sense that we affirm that the constraint force can become impulsive: if the motion which passes through the singular set and tries to enter the regular one again does not satisfy the new constraints, then it experiences a jump of its momentum, due to the presence of the constraint force. In this way, the new values of the momentum satisfy the constraints. But one has to be careful: the impulsive force will act just on leaving $`S`$, on the regular set. Consequently, we must take into account the virtual displacements associated to the posterior regular points. The underlying idea of the mathematical derivation of the momentum jumps in Section 4.1 is the following: to take an infinitesimal posterior point $`q(t)`$ to $`q(t_0)`$, to forget for a moment the presence of the constraints on the path $`q(t_0,t)`$ and to derive the momentum jump at $`q(t)`$ due to the appearence of the additional constraints. Afterwards, to make a limit process $`tt_0`$, cancelling out the interval $`(t_0,t)`$ where we “forgot” the constraints. In any case, we will make the convention that the jump happens at $`q(t_0)`$.

We illustrate the above discussion in the following example.

###### Example 4.1

Consider a particle in the plane subjected to the constraints imposed by the generalized codistribution in Example 2.1. The Lagrangian function is

$$\begin{array}{cccc}\hfill L:& \hfill T\text{}^2& & \text{}\hfill \\ & \hfill (x,y,\dot{x},\dot{y})& & \frac{1}{2}m(\dot{x}^2+\dot{y}^2).\hfill \end{array}$$

On the half-plane $`R_1=\{x<0\}`$ the codistribution is zero and the motion is free. Consequently the trajectories are

$`x`$ $`=`$ $`\dot{x}_0t+x_0,`$
$`y`$ $`=`$ $`\dot{y}_0t+y_0.`$
If the particle starts its motion with initial conditions $`x_0=1`$, $`y_0=1`$, $`\dot{x}_0=1`$, $`\dot{y}_0=0`$, after a time $`1`$, it reaches the singular set $`S=\{x=0\}`$. If the motion crosses the $`y`$-axis, something abrupt occurs on entering the half-plane $`R_2=\{x>0\}`$, where the codistribution is no longer zero and, indeed, imposes the additional constraint $`\dot{x}=\dot{y}`$ (Case 1). We know that the integral manifolds of $`D`$ on $`R_2`$ are half-lines of slope 1, so the particle suffers a finite jump in the velocity on going through the singular part in order to adapt its motion to the prescribed direction.

If, on the contrary, the particle starts on $`R_2`$, say with initial conditions $`x_0=1`$, $`y_0=1`$, $`\dot{x}_0=1`$, $`\dot{y}_0=1`$, after a certain time, it reaches the set $`S`$. On crossing it, nothing special happens, because the particle finds less contraints to fulfill, indeed, there are no constraints (Case 2). Its motion on $`R_1`$ is free, on a straight line of slope 1 and with constant velocity equal to the one at the singular point of crossing.

### 4.1 Momentum jumps

Now, we derive a formula, strongly inspired by the theory of impulsive motion, for the momentum jumps which can occur due to the changes of rank of the codistribution $`D`$ in Cases 1 and 3.

At $`q(t_0)`$ we define the following vector subspaces of $`T_{q(t_0)}^{}Q`$

$`D_{q(t_0)}^{}`$ $`=`$ $`\{\alpha T_{q(t_0)}^{}Q/\stackrel{~}{\alpha }:(t_0ϵ,t_0)T^{}Q,\stackrel{~}{\alpha }(t)D_{q(t)}\text{and}\underset{tt_0^{}}{lim}\stackrel{~}{\alpha }(t)=\alpha \}`$
$`D_{q(t_0)}^+`$ $`=`$ $`\{\alpha T_{q(t_0)}^{}Q/\stackrel{~}{\alpha }:(t_0,t_0+ϵ)T^{}Q,\stackrel{~}{\alpha }(t)D_{q(t)}\text{and}\underset{tt_0^+}{lim}\stackrel{~}{\alpha }(t)=\alpha \}`$
From the definition of $`D_{q(t_0)}^{}`$ and $`D_{q(t_0)}^+`$ we have that

$$(D_{q(t_0)}^{})^{}=\underset{tt_0^{}}{lim}(D_{q(t)})^{}\text{ and }(D_{q(t_0)}^+)^{}=\underset{tt_0^+}{lim}(D_{q(t)})^{}$$

where  denotes the orthogonal complement with respect to the bilinear form induced by the metric $`g`$ on the cotangent space $`T_{q(t_0)}^{}Q`$, and the limits $`(D^{})^{}`$ and $`(D^{})^+`$ are defined as in the case of $`D^{}`$ and $`D^+`$. In the following, we shall not make a notational distinction between the metric $`g`$ and the induced bilinear form on $`T_{q(t_0)}^{}Q`$. In each case, the precise meaning should be clear from the context.

Since $`D`$ is a differentiable codistribution then

$$D_{q(t_0)}D_{q(t_0)}^{}\text{ and }D_{q(t_0)}D_{q(t_0)}^+.$$

Along the interval $`[t_0,t]`$, we have

$$p_A(t)p_A(t_0)=_{t_0}^tF_A(\tau )𝑑\tau .$$

Multiplying by the virtual displacements at the point $`q(t)`$, we obtain

$$(p_A(t)p_A(t_0))\delta q_{|q(t)}^A=_{t_0}^tF_A(\tau )𝑑\tau \delta q_{|q(t)}^A.$$
(7)
Summing in $`A`$, we get

$$\underset{A=1}{\overset{n}{}}(p_A(t)p_A(t_0))\delta q_{|q(t)}^A=\underset{A=1}{\overset{n}{}}_{t_0}^tF_A(\tau )𝑑\tau \delta q_{|q(t)}^A.$$
(8)
Since we are dealing with ideal constraints, the virtual work vanishes, that is

$$\underset{A=1}{\overset{n}{}}_{t_0}^tF_A(\tau )\delta q_{|q(\tau )}^A𝑑\tau =0.$$

If $`t`$ is near $`t_0`$, then $`\tau `$ is close to $`t`$, and $`q(\tau )`$ is near $`q(t)`$, so $`\delta q`$ remains both nearly constant and nearly equal to its value at time $`t`$ throughout the time interval $`(t_0,t]`$, in the same way we exposed in Section 3. Therefore,

$$\underset{A=1}{\overset{n}{}}\left(_{t_0}^tF_A(\tau )𝑑\tau \right)\delta q_{|q(t)}^A=\underset{A=1}{\overset{n}{}}_{t_0}^tF_A(\tau )\delta q_{|q(\tau )}^A𝑑\tau +O(tt_0)=O(tt_0).$$

Consequently, equation (8) becomes

$$\underset{A=1}{\overset{n}{}}(p_A(t)p_A(t_0))\delta q_{|q(t)}^A=O(tt_0).$$
(9)
Taking limits we obtain

$$\underset{tt_0^+}{lim}\left(\underset{A=1}{\overset{n}{}}(p_A(t)p_A(t_0))\delta q_{|q(t)}^A\right)=0$$

which implies

$$\underset{A=1}{\overset{n}{}}(p_A(t_0)_+p_A(t_0))\underset{tt_0^+}{lim}\delta q_{|q(t)}^A=0$$
(10)
or, in other words,

$$(p_A(t_0)_+p_A(t_0))dq^A\underset{tt_0^+}{lim}D_{q(t)}=D_{q(t_0)}^+.$$
(11)
#### Conclusion:

Following the above discussion, we will deduce the existence of jump of momenta depending on the relation between $`D_{q(t_0)}^{}`$ and $`D_{q(t_0)}^+`$. The possible cases are shown in Table 2.

In the second case in Table 2, we have a jump of momenta if the “pre-impact” momentum $`p(t_0)_{}=p(t_0)`$ does not satisfy the constraints imposed by $`D_{q(t_0)}^+`$, that is

$$p_A(t_0)_{}dq^A(D_{q(t_0)}^+)^{}.$$

Our proposal for the equations which determine the jump is then

$$\{\begin{array}{c}(p_A(t_0)_+p_A(t_0)_{})dq^AD_{q(t_0)}^+\hfill \\ p_A(t_0)_+dq^A(D_{q(t_0)}^+)^{}.\hfill \end{array}$$

The first equation has been derived above (cf. (11)) from the generalized writing of Newton’s second law (1). The second equation simply encodes the fact that the “post-impact” momentum must satisfy the new constraints imposed by $`D_{q(t_0)}^+`$.

###### Remark 4.2

In Cases 1 and 3, the virtual displacements at $`q(t_0)`$ are radically different from the ones at the regular posterior points, because of the change of rank. From a dynamics point of view, these are the “main” ones, since it is on the regular set where an additional constraint reaction force acts. As we have seen, the momentum jump happens on just leaving $`S`$, due to the presence of this additional constraint force on the regular set. Note that with the procedure we have just derived, we are taking into account precisely the virtual displacements at the regular posterior points, and not those of $`q(t_0)`$. If we took the virtual displacements at $`q(t_0)`$ and multiply by them in (7), we would obtain non-consistent jump conditions. This is easy to see, for instance, in Example 4.1.

An explicit derivation of the momentum jumps for Cases 1 and 3 would be as follows. Let $`m`$ be the maximum between $`\rho =r`$, the rank at the regular preceding points, and $`\rho =s`$, the rank at the regular posterior points. Then there exists a neighbourhood $`U`$ of $`q(t_0)`$ and 1-forms $`\omega _1,\mathrm{},\omega _m`$ such that

$$D_q=\text{span}\{\omega _1(q),\mathrm{},\omega _m(q)\},qU.$$

Let us suppose that $`\omega _1,\mathrm{},\omega _s`$ are linearly independent at the regular posterior points (if not, we reorder them). Obviously, at $`q(t_0)`$, these $`s`$ 1-forms are linearly dependent. In the following, we will denote by $`\omega _i`$ the 1-form evaluated at $`q(t)`$, ($`t`$ time immediately posterior to $`t_0`$) i.e. $`\omega _i\omega _i(q(t))`$, in order to simplify notation.

Since the Lagrangian is of the form $`L=TU`$, where $`T`$ is the kinetic energy of the Riemannian metric $`g`$, that is, $`L=\frac{1}{2}g_{AB}\dot{q}^A\dot{q}^BU(q)`$, then we have that

$$\omega _{jA}(q(t))\dot{q}^A(t)=\underset{A,B}{}\omega _{jA}g^{AB}p_B(t)=0,j=1,\mathrm{},s.$$
(12)
Using the metric $`g`$ we have the following decomposition

$$T_q^{}Q=D_qD_q^{},qQ.$$

The two complementary projectors associated to this decomposition are

$`𝒫_q:T_q^{}QD_q^{}`$
$`𝒬_q:T_q^{}QD_q`$
The projector $`𝒫_q`$ is given by

$$𝒫_q(\alpha _q)=\alpha _q𝒞^{ij}\alpha _q(Z_i)\omega _j,\alpha _qT^{}Q$$

where

$$Z_i=g^{AB}\omega _{iB}\left(\frac{}{q^A}\right)_q$$

and $`𝒞^{ij}`$ are the entries of the inverse matrix of $`𝒞`$, the symmetric matrix with entries $`𝒞_{ij}=\omega _{iA}g^{AB}\omega _{jB}`$, or $`𝒞=\omega g^1\omega ^T`$ with the obvious notations.

By definition

$$p_A(t_0)_+dq^A|_{q(t_0)}=\underset{tt_0^+}{lim}(p_A(t)dq^A|_{q(t)}).$$

From (12), $`𝒫_{q(t)}(p_A(t)dq^A|_{q(t)}))=p_A(t)dq^A|_{q(t)}`$ and then

$$p_A(t_0)_+dq^A|_{q(t_0)}=\left(\underset{tt_0^+}{lim}𝒫_{q(t)}\right)(p_A(t_0)_+dq^A|_{q(t_0)})(D_{q(t_0)}^+)^{}.$$
(13)
Combining (11) and (13), we obtain

$$p_A(t_0)_+dq^A|_{q(t_0)}=\left(\underset{tt_0^+}{lim}𝒫_{q(t)}\right)\left[p_A(t_0)_{}dq^A|_{q(t_0)}\right].$$

In coordinates, this can be expressed as

$$p_A(t_0)_+=p_A(t_0)_{}\underset{tt_0^+}{lim}\left(\underset{i,j,A,B}{}𝒞^{ij}\omega _{jB}g^{BC}\omega _{iA}\right)|_{q(t)}p_C(t_0)_{},A=1,\mathrm{},n.$$
(14)
Equation (14) can be written in matrix form as follows

$$p(t_0)_+=\left(Id\underset{tt_0^+}{lim}(\omega ^T𝒞^1\omega g^1)_{|q(t)}\right)p(t_0)_{}.$$
(15)
With the derived jump rule, we are able to prove the following version of Carnot’s theorem for generalized constraints.

###### Theorem 4.3

The kinetic energy will only decrease by the application of the jump rule (15).

Proof: We have that

$`g(p(t_0)_+,p(t_0)_+)`$ $`=`$ $`g(\left(\underset{tt_0^+}{lim}𝒫_{q(t)}\right)p(t_0)_{},p(t_0)_{}\left(\underset{tt_0^+}{lim}𝒬_{q(t)}\right)p(t_0)_{})`$
$`=`$ $`g(\left(\underset{tt_0^+}{lim}𝒫_{q(t)}\right)p(t_0)_{},p(t_0)_{})`$
$`=`$ $`g(p(t_0)_{},p(t_0)_{})g(\left(\underset{tt_0^+}{lim}𝒬_{q(t)}\right)p(t_0)_{},p(t_0)_{}).`$
Since $`g(\left(\underset{tt_0^+}{lim}𝒬_{q(t)}\right)p(t_0)_{},p(t_0)_{})=g(\left(\underset{tt_0^+}{lim}𝒬_{q(t)}\right)p(t_0)_{},\left(\underset{tt_0^+}{lim}𝒬_{q(t)}\right)p(t_0)_{})0`$, we can conclude that

$$\frac{1}{2}g(p(t_0)_+,p(t_0)_+)\frac{1}{2}g(p(t_0)_{},p(t_0)_{}).$$

$`QED`$

In fact, the jump rule (15) has the following alternative interpretation. Let $`pD_{q(t_0)}^+`$ and observe that

$`g(pp(t_0)_{},pp(t_0)_{})`$ $`=`$ $`g(p(t_0)_{},p(t_0)_{})+g(p,p2p(t_0)_{})`$
$`=`$ $`g(p(t_0)_{},p(t_0)_{})+g(p,p2\left(\underset{tt_0^+}{lim}𝒫_{q(t)}\right)p(t_0)_{}).`$
Now, note that the covector $`p=\left(lim_{tt_0^+}𝒫_{q(t)}\right)p(t_0)_{}D_{q(t_0)}^+`$ is such that the expression $`g(pp(t_0)_{},pp(t_0)_{})`$ is minimized among all the covectors belonging to $`D_{q(t_0)}^+`$. Therefore, the derived jump rule (15) can be stated as follows: the “post-impact” momenta $`p(t_0)_+`$ is such that the kinetic energy corresponding to the difference of the “pre-impact” and “post-impact” momentum is minimized among all the covectors satisfying the constraints. This is an appropriate version for generalized constraints of the well-known jump rule for perfectly inelastic collisions . This is even more clear in the holonomic case, as is shown in Subsection 4.2.

###### Remark 4.4

So far, we have been dealing with impulsive constraints. More generally, we can consider the presence of external impulsive forces associated to external inputs or controls. Then, equation (11) must be modified as follows

$$(p_A(t_0)_+p_A(t_0)_{}P_A^{}(t_0))dq_{|q(t_0)}^A\underset{tt_0^+}{lim}D_{q(t)}=D_{q(t_0)}^+$$
(16)
where $`P_A^{}(t_0)`$, $`1An`$, are the external impulses at time $`t_0`$. Observe that if $`q(t_0)`$ is a regular point then

$$p_A(t_0)_+=p_A(t_0)_{}+P_A^{}(t_0)\underset{i,j,A,B}{}𝒞^{ij}\omega _{jB}g^{BC}\omega _{iA}P_C^{}(t_0),A=1,\mathrm{},n,$$
(17)
and, if $`q(t_0)`$ is a singular point, we have

$$p_A(t_0)_+=p_A(t_0)_{}+P_A^{}(t_0)\underset{tt_0^+}{lim}\left(\underset{i,j,A,B}{}𝒞^{ij}\omega _{jB}g^{BC}\omega _{iA}\right)|_{q(t)}(p_C(t_0)_{}+P_C^{}(t_0)),A=1,\mathrm{},n.$$
(18)
### 4.2 The holonomic case

We show in this section a meaningful interpretation of the proposed jump rule (15) in case the codistribution $`D`$ is partially integrable.

Let us consider a trajectory $`q(t)Q`$ which reaches a singular point $`q(t_0)S`$ and falls in either Case 1 or Case 3. Since $`Q=\overline{R}`$, we have that $`q(t_0)\overline{L}`$, where $`L`$ is the leaf of $`D`$ which contains the regular posterior points of the trajectory $`q(t)`$. On leaving $`q(t_0)`$, we have seen that the trajectory suffers a finite jump in its momentum in order to satisfy the constraints imposed by $`D`$, which in this case implies that the trajectory after time $`t_0`$ belongs to the leaf $`L`$. Consequently, the jump can be interpreted as a perfectly inelastic collision against the “wall” represented by the leaf $`L`$!

Let us see it revisiting Example 4.1.

###### Example 4.5

Consider again the situation in Example 4.1. If the motion of the particle starts on the left half-plane going towards the right one, then it is easy to see that $`D_{(0,y)}^{}=\{0\}`$ and $`D_{(0,y)}^+=\text{span}\{dxdy\}`$. As $`D_{(0,y)}^+D_{(0,y)}^{}`$, a jump of momenta is possible. In fact, if the “pre-impact” velocity $`(\dot{x}_0,\dot{y}_0)`$ does not satisfy $`\dot{x}_0=\dot{y}_0`$, the jump occurs and is determined by $`\mathrm{\Delta }v(t_0)D_{(0,y)}^+`$ and $`\dot{x}(t_0^+)=\dot{y}(t_0^+)`$. Consequently, we obtain

$`\dot{x}(t_0^+)`$ $`=`$ $`{\displaystyle \frac{\dot{x}_0+\dot{y}_0}{2}},`$
$`\dot{y}(t_0^+)`$ $`=`$ $`{\displaystyle \frac{\dot{x}_0+\dot{y}_0}{2}}.`$
We would have obtained the same result if we had considered that our particle hits, in a perfectly inelastic collision, against the “wall” represented by the half-line of slope 1 contained in $`\{x>0\}`$ passing through the point $`(0,y)`$.

If the particle starts on the right half-plane towards the left one, the roles are reversed and $`D_{(0,y)}^{}=<dxdy>`$, $`D_{(0,y)}^+=\{0\}`$. We have that $`D_{(0,y)}^+D_{(0,y)}^{}`$ and therefore there is no jump.

## 5 Examples

Next, we are going to develop two examples illustrating the above discussion. First, we treat a variation of the classical example of the rolling sphere . Secondly, we take one example from Chen, Wang, Chu and Chou .

### 5.1 The rolling sphere

Consider a homogeneous sphere rolling on a plane. The configuration space is $`Q=\text{}^2\times SO(3)`$: $`(x,y)`$ denotes the position of the center of the sphere and $`(\phi ,\theta ,\psi )`$ denote the Eulerian angles.

Let us suppose that the plane is smooth if $`x<0`$ and absolutely rough if $`x>0`$. On the smooth part, we assume that the motion of the ball is free, that is, the sphere can slip. But if it reaches the rough half-plane, the sphere begins rolling without slipping, because of the presence of the constraints imposed by the roughness. We are interested in knowing the trajectories of the sphere and, in particular, the possible changes in its dynamics because of the crossing from one half-plane to the other.

The kinetic energy of the sphere is

$$T=\frac{1}{2}\left(\dot{x}^2+\dot{y}^2+k^2(\omega _x^2+\omega _y^2+\omega _z^2)\right),$$

where $`\omega _x`$, $`\omega _y`$ and $`\omega _z`$ are the angular velocities with respect to the inertial frame, given by

$`\omega _x`$ $`=`$ $`\dot{\theta }\mathrm{cos}\psi +\dot{\phi }\mathrm{sin}\theta \mathrm{sin}\psi ,`$
$`\omega _y`$ $`=`$ $`\dot{\theta }\mathrm{sin}\psi \dot{\phi }\mathrm{sin}\theta \mathrm{cos}\psi ,`$
$`\omega _z`$ $`=`$ $`\dot{\phi }\mathrm{cos}\theta +\dot{\psi }.`$
The potential energy is not considered here since it is constant.

The condition of rolling without sliding of the sphere when $`x>0`$ implies that the point of contact of the sphere and the plane has zero velocity

$`\varphi ^1`$ $`=`$ $`\dot{x}r\omega _y=0,`$
$`\varphi ^2`$ $`=`$ $`\dot{y}+r\omega _x=0,`$
where $`r`$ is the radius of the sphere.

Following the classical procedure, we introduce quasi-coordinates “$`q_1`$”, “$`q_2`$” and “$`q_3`$” such that “$`\dot{q}_1`$$`=\omega _x`$, “$`\dot{q}_2`$$`=\omega _y`$ and “$`\dot{q}_3`$$`=\omega _z`$. These last expressions only have a symbolic meaning where we interpret “$`dq_i`$” and “$`\frac{}{q_i}`$”, $`1i3`$, as adequate combinations of the differentials and partial derivatives, respectively, of the eulerian angles. Note that $`\{\frac{}{x},\frac{}{y},\frac{}{z}`$, “$`\frac{}{q_1}`$”, “$`\frac{}{q_2}`$”, “$`\frac{}{q_3}`$$`\}`$ and $`\{dx,dy,dz,`$$`dq^1`$”, “$`dq^2`$”, “$`dq^3`$$`\}`$ are dual bases. Moreover, observe that, from the nonintegrability of the constraints, the differential forms “$`dq^i`$” do not represent exact differentials.

The non-holonomic generalized differentiable codistribution $`D`$ is given by

$$D_{(x,y,\varphi ,\theta ,\psi )}=\{\begin{array}{cc}\{0\},& \text{if}x0,\\ \text{span}\{dxrdq^2,dy+rdq^1\},& \text{if}x>0.\end{array}$$

The intersection of the regular set of the generalized codistribution and the $`(x,y)`$-plane has two connected components, the half-planes $`R_1=\{x<0\}`$ and $`R_2=\{x>0\}`$. The line $`\{x=0\}`$ belongs to the singular set of $`D`$.

On $`R_1`$ the codistribution is zero, so the motion equations are

$`m\ddot{x}`$ $`=`$ $`0,`$
$`m\ddot{y}`$ $`=`$ $`0,`$
$`mk^2\dot{\omega _x}`$ $`=`$ $`0,`$ (19)
$`mk^2\dot{\omega _y}`$ $`=`$ $`0,`$
$`mk^2\dot{\omega _z}`$ $`=`$ $`0.`$
On $`R_2`$ we have to take into account the constraints to obtain the following equations of motion

$`m\ddot{x}`$ $`=`$ $`\lambda _1,`$
$`m\ddot{y}`$ $`=`$ $`\lambda _2,`$
$`mk^2\dot{\omega _x}`$ $`=`$ $`r\lambda _2,`$ (20)
$`mk^2\dot{\omega _y}`$ $`=`$ $`r\lambda _1,`$
$`mk^2\dot{\omega _z}`$ $`=`$ $`0,`$
with the constraint equations $`\dot{x}r\omega _y=0`$ and $`\dot{y}+r\omega _x=0`$. One can compute the Lagrange multipliers by an algebraic procedure described in .

Suppose that the sphere starts its motion at a point of $`R_1`$ with the following initial conditions at time $`t=0`$: $`x_0<0`$, $`y_0`$, $`\dot{x}_0>0`$, $`\dot{y}_0`$, $`(\omega _x)_0`$, $`(\omega _y)_0`$ and $`(\omega _z)_0`$. Integrating equations (5.1) we have that if $`x(t)<0`$

$`x(t)`$ $`=`$ $`\dot{x}_0t+x_0,`$
$`y(t)`$ $`=`$ $`\dot{y}_0t+y_0,`$
$`\omega _x(t)`$ $`=`$ $`(\omega _x)_0,`$ (21)
$`\omega _y(t)`$ $`=`$ $`(\omega _y)_0,`$
$`\omega _z(t)`$ $`=`$ $`(\omega _z)_0.`$
At time $`\overline{t}=x_0/\dot{x}_0`$ the sphere finds the rough surface of the plane, where the codistribution is no longer zero and it is suddenly forced to roll without sliding (Case 1). Following the discussion in Section 4, we calculate the instantaneous change of velocity (momentum) at $`x=0`$.

First of all we compute the matrix $`𝒞`$

$`𝒞`$ $`=`$ $`\left(\begin{array}{ccccc}1& 0& 0& r& 0\\ 0& 1& r& 0& 0\end{array}\right)\left(\begin{array}{ccccc}1& 0& 0& 0& 0\\ 0& 1& 0& 0& 0\\ 0& 0& k^2& 0& 0\\ 0& 0& 0& k^2& 0\\ 0& 0& 0& 0& k^2\end{array}\right)\left(\begin{array}{cc}1& 0\\ 0& 1\\ 0& r\\ r& 0\\ 0& 0\end{array}\right)`$
$`=`$ $`(1+r^2k^2)\left(\begin{array}{cc}1& 0\\ 0& 1\end{array}\right).`$
Next, a direct computation shows that the projector $`𝒫`$ does not depend on the base point

$$𝒫=\left(\begin{array}{ccccc}\frac{r^2}{r^2+k^2}& 0& 0& \frac{r}{r^2+k^2}& 0\\ 0& \frac{r^2}{r^2+k^2}& \frac{r}{r^2+k^2}& 0& 0\\ 0& \frac{rk^2}{r^2+k^2}& \frac{k^2}{r^2+k^2}& 0& 0\\ \frac{rk^2}{r^2+k^2}& 0& 0& \frac{k^2}{r^2+k^2}& 0\\ 0& 0& 0& 0& 1\end{array}\right).$$

Therefore, we have

$`(p_x)_+`$ $`=`$ $`{\displaystyle \frac{r^2(p_x)_0+r(p_2)_0}{r^2+k^2}},`$
$`(p_y)_+`$ $`=`$ $`{\displaystyle \frac{r^2(p_y)_0r(p_1)_0}{r^2+k^2}},`$
$`(p_1)_+`$ $`=`$ $`{\displaystyle \frac{rk^2(p_y)_0+k^2(p_1)_0}{r^2+k^2}},`$
$`(p_2)_+`$ $`=`$ $`{\displaystyle \frac{rk^2(p_x)_0+k^2(p_2)_0}{r^2+k^2}},`$
$`(p_3)_+`$ $`=`$ $`(p_3)_0.`$
Now, using the relation between the momenta and the quasi-velocities

$$p_x=\dot{x},p_y=\dot{y},p_1=k^2\omega _x,p_2=k^2\omega _y,p_3=k^2\omega _z,$$

we deduce that

$`\dot{x}_+`$ $`=`$ $`{\displaystyle \frac{r^2\dot{x}_0+rk^2(\omega _y)_0}{r^2+k^2}},`$
$`\dot{y}_+`$ $`=`$ $`{\displaystyle \frac{r^2\dot{y}_0rk^2(\omega _x)_0}{r^2+k^2}},`$
$`(\omega _x)_+`$ $`=`$ $`{\displaystyle \frac{r\dot{y}_0+k^2(\omega _x)_0}{r^2+k^2}},`$ (24)
$`(\omega _y)_+`$ $`=`$ $`{\displaystyle \frac{r\dot{x}_0+k^2(\omega _y)_0}{r^2+k^2}},`$
$`(\omega _z)_+`$ $`=`$ $`(\omega _z)_0.`$
Finally, integrating equations (5.1) at time $`\overline{t}=x_0/\dot{x}_0`$ with initial conditions given by (5.1) we obtain that if $`t>\overline{t}`$

$`x(t)`$ $`=`$ $`{\displaystyle \frac{r^2\dot{x}_0+rk^2(\omega _y)_0}{r^2+k^2}}(t\overline{t}),`$
$`y(t)`$ $`=`$ $`{\displaystyle \frac{r^2\dot{y}_0rk^2(\omega _x)_0}{r^2+k^2}}(t\overline{t})+\dot{y}_0\overline{t}+y_0,`$
$`\omega _x(t)`$ $`=`$ $`{\displaystyle \frac{r\dot{y}_0+k^2(\omega _x)_0}{r^2+k^2}},`$ (25)
$`\omega _y(t)`$ $`=`$ $`{\displaystyle \frac{r\dot{x}_0+k^2(\omega _y)_0}{r^2+k^2}},`$
$`\omega _z(t)`$ $`=`$ $`(\omega _z)_0.`$
### 5.2 Particle with constraint

Let us consider the motion of a particle of mass $`1`$ in $`\text{}^3`$ subjected to the following constraint

$$\varphi =(y^2x^2z)\dot{x}+(zy^2xy)\dot{y}+x\dot{z}=0.$$

In addition, let us assume that there is a central force system centered at the point $`(0,0,1)`$ with force field given by

$$F=xdxydy+(1z)dz.$$

Then, the Lagrangian function of the particle is

$$L=TV=\frac{1}{2}(\dot{x}^2+\dot{y}^2+\dot{z}^2+x^2+y^2+z^22z),$$

and the constraint defines a generalized differentiable codistribution $`D`$, whose singular set is $`S=\{(x,y,z):x=0,z=y^2\}`$.

On $`R`$, the regular set of $`D`$, the dynamics can be computed following the standard symplectic procedure  to obtain $`\mathrm{\Gamma }_{L,D}=\mathrm{\Gamma }_L+\lambda Z`$, where

$`\mathrm{\Gamma }_L`$ $`=`$ $`\dot{x}{\displaystyle \frac{}{x}}+\dot{y}{\displaystyle \frac{}{y}}+\dot{z}{\displaystyle \frac{}{z}}x{\displaystyle \frac{}{\dot{x}}}y{\displaystyle \frac{}{\dot{y}}}(z1){\displaystyle \frac{}{\dot{z}}},`$
$`Z`$ $`=`$ $`\left((y^2x^2z){\displaystyle \frac{}{\dot{x}}}+(zy^2xy){\displaystyle \frac{}{\dot{y}}}+x{\displaystyle \frac{}{\dot{z}}}\right),`$
and $`\lambda `$ is given by

$$\lambda =\frac{\mathrm{\Gamma }_L(\varphi )}{Z(\varphi )}=\frac{2x\dot{x}^2+y\dot{y}\dot{x}2y\dot{y}^2x\dot{y}^2+\dot{y}\dot{z}+x^3+y^3yz+x}{(y^2x^2z)^2+(zy^2xy)^2+x^2}.$$
(26)
Consequently, the motion equations on $`R`$ are

$`m\ddot{x}+x`$ $`=`$ $`\lambda (y^2x^2z),`$
$`m\ddot{y}+y`$ $`=`$ $`\lambda (zy^2xy),`$ (27)
$`m\ddot{z}+z1`$ $`=`$ $`\lambda x,`$
with the constraint equation $`\varphi =0`$.

From the discussions of , we know that in this case there is an integral surface, $`C`$, of the constraint $`\varphi `$, that is, a surface on which all motions satisfy the constraint. This surface is

$$C=\{(x,y,x):zx^2y^2+xy=0\}.$$

Note that $`SC`$. Therefore, if a motion takes place on the cone-like surface $`C`$, it is confined to stay on this critical surface, unless it reaches a singular point. In this case, the space of allowable motions is suddenly increased (in fact, $`T\text{}^3`$), and the motion can “escape” from $`C`$. In addition, this proves that the unique way to pass from one point of the exterior of the $`C`$ to the interior, or viceversa, is through the singular set $`S`$.

In particular, we are interested in knowing

1. Is there any trajectory satisfying equations (5.2) which passes through the singular set?
2. if so, which are the possible momentum jumps due to the changes in the rank of the codistribution $`D`$?

So far, we do not know an answer for the question of the existence of a motion of (5.2) crossing $`S`$. It seems that on approaching a singular point, the constraint force can become increasingly higher (26). Consequently, this force possibly “disarranges” the approaching of the motion to $`S`$. Numerical simulations are quite useless in this task, because of the special nature of the problem: the hard restriction given by the fact that a motion crossing the cone-like surface $`C`$ must do it through the singular part $`S`$. Indeed, the numerical simulation performed in  crosses the surface $`C`$ through points which are not in $`S`$.

Concerning the second question, let us suppose that there is a trajectory of the dynamical system (5.2), $`q(t)=(x(t),y(t),z(t))`$, that passes through a singular point at time $`t_0`$, i.e. $`x(t_0)=0`$ and $`z(t_0)=y^2(t_0)`$. The rank of the codistribution $`D`$ at the immediately preceding and posterior points is 1, meanwhile at $`q(t_0)`$ it is 0 (Case 3). So, a possible jump of the momentum can be induced by the change in the rank of $`D`$.

A direct computation shows that the projector $`𝒫`$ depends explicitly on the base point $`qQ`$. Equivalently, we have that $`D_{q(t_0)}^+`$ depends strongly on the trajectory $`q(t)`$. In fact, taking two curves $`q_1(t)`$, $`q_2(t)`$ passing through $`q(t_0)`$ at time $`t_0`$, and satisfying $`x_1(t)z_1(t)y_1^2(t)`$ and $`z_2(t)y_2^2(t)x_2(t)`$ when $`tt_0^+`$ respectively, one can easily see that $`D_{q_1(t_0)}^+D_{q_2(t_0)}^+`$ (the expression $`f(t)g(t)`$ when $`tt_0^+`$ means that $`lim_{tt_0^+}f(t)/g(t)=0`$).

Consequently, we are not able to give an answer to question (ii) (in case the first one was true) unless we assume some additional information: for example, that the balance between $`x(t)`$ and $`z(t)y^2(t)`$ is the same for $`tt_0^{}`$ and $`tt_0^+`$. In such a case, $`D_{q(t_0)}^{}=D_{q(t_0)}^+`$ and we would conclude that there is no jump. In mechanical phenomenae of the type sliding-rolling, as the ones studied in Section 5.1, this kind of “indeterminacy” will not occur in general.

In spite of the fact that the most natural thing in this case seems to be to think that there is no jump of momenta, a mathematical explanation of it is still to be found.

## Acknowledgements

This work was partially supported by grant DGICYT (Spain) PB97-1257. J. Cortés and S. Martínez wish to thank Spanish Ministerio de Educación y Cultura for FPU and FPI grants, respectively. The authors wish to thank F. Cantrijn for helpful comments and suggestions. We also wish to thank the referees for their valuable criticism which contributed to improve this paper.

Jorge Cortés † Manuel de León ‡ Sonia Martínez §

Instituto de Matemáticas y Física Fundamental
Consejo Superior de Investigaciones Científicas
Serrano 123, 28006 Madrid, SPAIN
e-mail: †j.cortes@imaff.cfmac.csic.es ‡mdeleon@imaff.cfmac.csic.es
§s.martinez@imaff.cfmac.csic.es

David Martín de Diego

Departamento de Economía Aplicada (Matemáticas)
Facultad de CC. Económicas y Empresariales
Universidad de Valladolid
Avda. Valle Esgueva 6, 47011 Valladolid, SPAIN
e-mail: dmartin@esgueva.eco.uva.es


text:
# An algorithm for cleaning extended interferences in radio astronomical records

## 1 Introduction

In the course of astrophysical observations one often encounters interferences that corrupt the arriving useful signal or the statistical noise estimates in observation records. This is especially pronounced in records of radio astronomical observations. In the processing of data, a number of noises are recorded that require special cleaning, which is frequently done manually, for instance, in the standard reduction system FADPS at RATAN-600 (Verkhodanov et al., 1993).

To perform the tasks of record cleaning from discontinuous jamming, extensive use is made of technical means at the detector level (see Berlin et al., 1995). However, when working with the archival data, one has to come across different kinds of noises that disturb the initial record. Such noises may be ordinary discontinuous jamming, radar signals, and also various jumps caused by the instability of ADC in data recording.

Recent work concerned with the development of the archive ODA (Kononov et al., 1998; Kononov, Pavlov, 1999; Kononov et al., 1999) of RATAN-600 observational data coming from the broad-band radiometers made it possible to solve the problems which demand the use of old data. Among such problems one may isolate, for instance, investigation of the background radiation (Parijskij, Korol’kov, 1986) and search for variability of radio astronomical objects. The processing of these data has naturally called for automatic search and removal of interferences that corrupt the signal. Some work over control of noises by programme means based on nonparametric evaluation of the mean was done at RATAN-600 several years ago (Erukhimov et al., 1990; Chernenkov, 1996; Shergin et al., 1997). One should also mention the programme development based on a gradient approach of removing discontinuous jamming which was initiated by Bursov (1987). The work over the creation of the cleaning procedure described below was started by Pavlov in 1997 (Pavlov, 1998) under the supervision of the former of the authors of this paper, when the problem of search and removal of solitary interference from the radar was resolved.

It is proposed in this paper to develop an algorithm capable of removing non-single interferences on a record. The algorithm employs both the robust (immune to influence) estimates of the mean on the basis of non-parametric rank algorithms (Huber, 1981) and the procedures based on extreme statistics, namely, minimax approach (Verkhodanov, Gorokhov, 1995), which enable fast evaluation of statistics changes in a given interval.

Below is presented a detailed description of removal of two classes of interferences: 1) radar and discontinuous jamming; 2) jumps in records.

## 2 Algorithm of radar interference cleaning

The problems presented by the signal from a radar (Fig. 1) are that it is not a $`\delta `$-function, that is, it cannot be detected with the aid of the “gradient” approach without confusion with a real source, and it is quasiperiodic, i.e. it is “fuzzed up” in the region of spatial frequencies, which hampers its removal on the Fourier plane. Besides, because the signal is extended, the clipping (namely, removal of all points exceeding a specified level with subsequent replacement by pixels according to given rules) of its spatial spectrum will result in corruption of the signal. For this reason the authors decided to use statistical techniques to detect and remove interferences produced by the radar. The developed procedure permitted getting rid of simple discontinuous jamming. In principle, within the frame of FADPS one can establish a procedure of removing interferences from radars on the basis of the already available modules (Verkhodanov et al., 1992). For instance, such an algorithm can be synthesized from the modules of low noise subtraction, subtraction of the signal of Gaussian shape with a size larger than the specified one, clipping the record at a given level, and taking the sum of the record made with the subtracted low-frequency noises and detected sources. However, because of the multiplicity of Gaussians being computed and also the repetition of vector operations of addition/subtraction increasing the calculation time of the script operation, and due to the ambiguities arising in superposition of a useful signal and a radar signal, the authors made a decision to abandon this approach.

Prior to describing of the developed algorithm we will explain some terms being employed in this paper. By the window is meant the interval on the abscissa axis (for our data this is a time interval measured in seconds) inside of which a signal search is made. By sliding the window along the record we mean the sequential filling of the array of pixels in the appropriate interval when the window is shifted by 1 pixel along the record. By the robust evaluation of the mean is implied a simple median of distribution of pixel values in the given interval, while the robust dispersion is the dispersion estimate made via the median absolute deviation (MAD) under the assumption of Gaussian noises (in this case if $`\mu `$ is MAD, then $`\sigma `$=$`\mu `$/0.674 in accordance with the definition of second-order moments for the normal distribution).

When developing the algorithm, we used the archival data (Kononov et al., 1999) written in the form of the FITS-like F format (Verkhodanov et al., 1993) of the data processing system FADPS at RATAN-600.

The algorithm eliminating interferences from the radar comprises several steps:

1. Copying of data from the F file to the working array and evaluation of the robust dispersion $`\sigma `$ of the total record.
2. Search for minimum and maximum values in the window sliding along the one-dimensional array. The window size is chosen equal to the beam of the radio telescope. As the window is sliding, an array of maximum peak-to-peak difference is filled in the current interval. The maximum scatter is computed as the difference of the maximum and minimum values of the pixels (max–min) in the current limits along the axis of shift.
3. When looking for the maximum among the values of differences, which exceeds the level $`n\sigma `$, where $`n`$ is the assigned level (by default 3), analysis is made to search for the maximum locations by way of simple comparing the array values.
4. In the region where the local maximum is found the behaviour of the statistics of the surrounding pixels is determined in case the signal does not fall fully within the window. With this aim in view, the current position of the interval with respect to the signal being detected is found: right and left of the maximum, by means of robust averaging of neighbouring points. Depending on the current vicinity (right or left) of the observed signal, the location of the interval in which the parameters will be estimated is chosen. The real maximum must always fall within the given interval.
5. Once the window position has been found, statistical estimate of principal characteristics of the detected signal is made. For this purpose the following formulae (Fomalont, 1989) are used:

$$I=\underset{t}{}I_t$$

$$X=\frac{1}{I}\underset{t}{}X_tI_t$$
$`(6)`$
$$B=\sqrt{\frac{1}{I}(X_tX)^2I_t},t=1,r;$$

where $`I_t`$ is the intensity value of the $`t`$-th point in the current interval, $`X_t`$ is the location of this point in the record, $`X`$ is the coordinate of the source centre of gravity, $`I`$ is the integral intensity of a possible sought-for signal, $`B`$ is the object size, $`r`$ is the size of the surveyed interval in pixels.
6. Comparison of the found parameters with the specified ones to search for an interfering signal. If the determined amplitude is smaller than the specified one or the size (duration) of the signal exceeds the given one (i.e. a real source is detected), the window is then shifted by one pixel and a return to the point of calculation of (max–min) statistics is accomplished.
7. If the parameters correspond to the expected properties of an interference, the boundaries of the signal to be removed are defined.
8. In the procedure of replacement of the points of the record within the limits found, provision is made for variants of replacing the generated noise with characteristic properties (dispersion and mean) of the noise process of the present record or ordinary regression by two estimated means at the edges of the record bordering on the current interval. To generate a noise Gaussian process, the standard library function of the language “C” rand() is used, which outputs uniformly distributed values. For conversion to the Gaussian noise, we applied a simple procedure of averaging with a shifted mean and normalizing:

$$N_i=\sigma \sqrt{12m}\frac{2\underset{i}{\overset{m}{}}Rm}{2m},$$

where $`R`$ is the vector of random numbers of the uniform distribution of the function rand() returned by the generator, $`m`$ is the number of the points involved in averaging of uniformly distributed numbers. The larger $`m`$, the better the distribution corresponds to normal. We employed $`m`$=32 in our procedure. The value of the current system time is used as the origin of the generator.
9. If there is one more maximum in a specified interval, then return to the point of search for the location of the maximum in the current window (see above) is performed.
10. If the end of the record is not reached, a shift is made by the size of the interval of cleaning along the axis.
11. When the record is finished the data are output to the F file.

The schematic of the radar interference removal algorithm described above is displayed in Fig.2.

The results of operation of the algorithm presented in Fig.3.

Fig.4 shows the power Fourier spectra of the two records displayed in Fig.3: prior to cleaning and after it. To highlight the difference, one and the same curve of the low-frequency background is subtracted from both records. The difference between the spectra is caused by elimination of the discontinuous jamming and quasiperiodic radar signal.

## 3 Algorithms of removing “jumps”

Apart from interferences treated by the described algorithm, there is a class of interferences not infrequently found in archival records: jumps of the mean or “steps” (Fig.6). We have developed a procedure of search for the jump on the basis of analysis of the mean over four neighbouring sliding intervals.

The algorithm is described as follows:

1. Reading data from the F file to the working array and estimation of the robust dispersion $`\sigma `$ of the total record.
2. Filling the array $`M_n`$, where every element is an estimate of the mean in the sliding window of a specified size. It is desirable that the size of the window should be chosen not shorter than 5 pixels.
3. Location of the interval in which a jump occurred. Here we analyse the array of the mean values estimated by a robust technique with the aim of revealing the position of the abrupt increase of the mean ($`>n\sigma `$, by default $`n`$=3) value (gradient approach) by way of comparing the values of four elements of the array $`M`$. The jump is considered to be detected in the $`i`$-th interval provided the following conditions are satisfied at the same time:

$$|_{i2}M_{i1}|<\sigma ,$$

$$|M_{i1}M_{i+1}|<n\sigma ,$$

$$|M_{i+1}M_{i+2}|<\sigma .$$

Thus the scheme of coincidence of the mean in two neighbouring intervals at the edges and of uneven transition in the central window is realized.
4. Revision of the jump location (for instance, on the $`k`$-th pixel) by means of comparison.
5. Decrease of the values of all pixels following the $`k`$-th pixel by the value of the jump amplitude.
6. Decrease of all values in the array of the means in the sliding window, which follow the current value of the mean.
7. If there is no limitation on the length of the record, a shift along the record and jump to the condition of the check point are performed. In other words, the cycle is completed.
8. Recording to the F file.

In Fig.5 is shown a schematic diagram of the above-described algorithm of searching for and removing the jump in the records. The results of performance of the algorithm are demonstrated in Fig.6.

## 4 Conclusions

Through analysing statistical properties in the current interval of the record the developed algorithms find and remove interferences of different kinds quickly and effectively.

The procedure is involved as a command “clip” and included in the standard data processing system FADPS at RATAN-600.

Provision is made for possible extension of the algorithm to two-dimensional data and its application to removal of cosmic particle signals on CCD images.

* The authors owe debts of gratitude to N.F. Vojkhanskaya for valuable comments in the course of preparing the manuscript.