Realization of Nonholonomic Constraints and Singular Perturbation Theory for Plane Dumbbells

TL;DR

This paper investigates the dynamics of plane dumbbells under nonholonomic constraints realized through viscous friction, revealing complex long-term behaviors and challenges in applying singular perturbation theory.

Contribution

It provides a detailed analysis of the relation between viscous friction and ideal constraints, highlighting the limitations of classical asymptotic methods in this context.

Findings

Long-term behaviors of frictional and constrained systems can differ significantly.

Matching asymptotics require careful choice of time scales and cannot always handle secular terms.

A reduction procedure for plane dumbbells is developed, sometimes leading to analytic solutions.

Abstract

We study the dynamics of pairs of connected masses in the plane, when nonholonomic (knife-edge) constraints are realized by forces of viscous friction, in particular its relation to constrained dynamics, and its approximation by the method of matching asymptotics of singular perturbation theory when the mass to friction ratio is taken as the small parameter. It turns out that long term behaviors of the frictional and constrained systems may differ dramatically no matter how small the perturbation is, and when this happens is not determined by any transparent feature of the equations of motion. The choice of effective time scales for matching asymptotics is also subtle and non-obvious, and secular terms appearing in them can not be dealt with by the classical methods. Our analysis is based on comparison to analytic solutions, and we present a reduction procedure for plane dumbbells that leads to them in some cases.

Equations

$$R_i:=\frac{c}{2}\frac{|{\dot{r}}_i\cdot (r_1-r_2){|}^2}{|r_1-r_2{|}^2}.$$

$$R=\frac{c}{2}|{\text{pr}}_{r_1-r_2}(\dot{r_1}-\dot{r_2}){|}^2=\frac{c}{2}\frac{|(\dot{r_1}-\dot{r_2})\cdot (r_1-r_2){|}^2}{|r_1-r_2{|}^2}.$$

$$\frac{d}{dt}\frac{\partial L}{\partial {\dot{r}}_i}-\frac{\partial L}{\partial r_i}+\frac{\partial R}{\partial {\dot{r}}_i}=0.$$

$$\frac{\partial }{\partial {\dot{r}}_i}\frac{m}{2}(|{\dot{r}}_1{|}^2+|{\dot{r}}_2{|}^2)=m{\dot{r}}_i$$

$$\frac{\partial }{\partial {\dot{r}}_i}\frac{c}{2}|{\text{pr}}_{r_1-r_2}({\dot{r}}_i){|}^2=c{\text{pr}}_{r_1-r_2}({\dot{r}}_i)=c\frac{{\dot{r}}_i\cdot (r_1-r_2)}{|r_1-r_2{|}^2}(r_1-r_2),$$

$$\delta\big{(}(\dot{r_{1}}-\dot{r_{2}})\cdot(r_{1}-r_{2})\big{)}=(r_{1}-r_{2})\cdot(\delta\dot{r_{1}}-\delta\dot{r_{2}})\\ =(r_{1}-r_{2})\cdot\delta\dot{r_{1}}-(r_{1}-r_{2})\cdot\delta\dot{r_{2}}\,.$$

$$\{\begin{array}{l}m{\ddot{r}}_1+{\lambda }_1(r_1-r_2)=0\\ m{\ddot{r}}_2+{\lambda }_2(r_1-r_2)=0.\end{array}$$

$$\{\begin{array}{l}m{\ddot{r}}_{12}+{\lambda }_{12}{r}_{12}=0\\ m\ddot{r}+\lambda r_{12}=0\\ {\dot{r}}_{12}\cdot r_{12}=0,\dot{r}\cdot r_{12}=0.\end{array}$$

$${\dot{r}}_{12}=i\dot{\theta }\rho e^{i\theta },\text{ and }{\ddot{r}}_{12}=(-{\dot{\theta }}^2+i\ddot{\theta })\rho e^{i\theta }.$$

$$m(-{\dot{\theta }}^2+i\ddot{\theta })\rho e^{i\theta }+{\lambda }_{12}\rho e^{i\theta }=0.$$

$${\dot{r}}_{12}=i\dot{\theta }\rho e^{i\theta },{\ddot{r}}_{12}=(-{\dot{\theta }}^2+i\ddot{\theta })\rho e^{i\theta },$$

$$r(t)=r(0)+\int_{0}^{t}i\beta\rho e^{i(\theta_{0}+\omega t)}dt=r(0)+\frac{\beta\rho}{\omega}e^{i(\theta_{0}+\omega t)}\Big{|}_{0}^{t}\\ =r(0)+\beta\rho e^{i\theta_{0}}\,\frac{e^{i\omega t}-1}{\omega}=r(0)-\frac{\beta\rho}{\omega}e^{i\theta_{0}}+\frac{\beta\rho}{\omega}e^{i(\theta_{0}+\omega t)}$$

$$r_{1,2}=r\pm r_{12}=r(0)-\dot{r}(0)/\omega +(\beta /\omega \pm 1)\rho e^{i({\theta }_0+\omega t)},$$

$$F=\frac{c}{2}\left(\frac{|\dot{r_1}\cdot (r_1-r_2){|}^2}{|r_1-r_2{|}^2}+\frac{|\dot{r_2}\cdot (r_1-r_2){|}^2}{|r_1-r_2{|}^2}\right),$$

$$\{\begin{array}{l}m{\ddot{r}}_1+\lambda (r_1-r_2)+c\frac{{\dot{r}}_1\cdot (r_1-r_2)}{|r_1-r_2{|}^2}(r_1-r_2)=0\\ m{\ddot{r}}_2-\lambda (r_1-r_2)+c\frac{{\dot{r}}_2\cdot (r_1-r_2)}{|r_1-r_2{|}^2}(r_1-r_2)=0.\end{array}$$

$$\{\begin{array}{l}m{\ddot{r}}_{12}+2\lambda r_{12}=0\\ m\ddot{r}+c\frac{\dot{r}\cdot r_{12}}{|r_{12}{|}^2}{r}_{12}=0.\end{array}$$

$$\dot{r}\cdot r_{12}=\text{Re}[(\alpha +i\beta )\rho e^{i\theta }\overline{\rho e^{i\theta }}]=\alpha {\rho }^2,\ddot{r}=[(\dot{\alpha }-\omega \beta )+i(\dot{\beta }+\omega \alpha )]\rho e^{i\theta },$$

$$\{\begin{array}{l}\dot{\beta }+\omega \alpha =0\\ m(\dot{\alpha }-\omega \beta )+c\alpha =0.\end{array}$$

$$\{\begin{array}{l}\dot{\beta }=-\omega \alpha \\ \epsilon \dot{\alpha }=-\alpha +\epsilon \omega \beta .\end{array}$$

$$\{\begin{array}{l}{\dot{\beta }}^0+\epsilon {\dot{\beta }}^1+\cdots =-\omega ({\alpha }^0+\epsilon {\alpha }^1+\dots )\\ \epsilon {\dot{\alpha }}^0+{\epsilon }^2{\dot{\alpha }}^1+\cdots =-{\alpha }^0-\epsilon {\alpha }^1-\cdots +\omega (\epsilon {\beta }^0+{\epsilon }^2{\beta }^1+\dots ).\end{array}$$

$$\{\begin{array}{l}{\dot{\beta }}^1=-\omega {\alpha }^1\\ {\dot{\alpha }}^0=-{\alpha }^1+\omega {\beta }^0.\end{array}$$

$$\{\begin{array}{l}\alpha =\epsilon \omega B_0+O({\epsilon }^2)\\ \beta =B_0+\epsilon (-{\omega }^2{B}_{0}t+B_1)+O({\epsilon }^2).\end{array}$$

$$\{\begin{array}{l}{\beta }^{\prime }=-\epsilon \omega \alpha \\ {\alpha }^{\prime }=-\alpha +\epsilon \omega \beta .\end{array}$$

$$\{\begin{array}{l}{\beta }^{0\prime }=0\\ {\alpha }^{0\prime }=-{\alpha }^0\end{array}\{\begin{array}{l}{\beta }^{1\prime }=-\omega {\alpha }^0\\ {\alpha }^{1\prime }=-{\alpha }^1+\omega {\beta }^0.\end{array}$$

$$\{\begin{array}{l}{\beta }^0=\beta (0)=:{\beta }_0\\ {\alpha }^0=\alpha (0)e^{-\tau }=:{\alpha }_0{e}^{-\tau }.\end{array}$$

$$\{\begin{array}{l}{\alpha }^1={\beta }_0\omega (1-e^{-\tau })\\ {\beta }^1=-{\alpha }_0\omega (1-e^{-\tau }).\end{array}$$

$$\begin{cases}\alpha=\alpha_{0}e^{-t/\varepsilon}+\varepsilon\omega\beta_{0}(1-e^{-t/\varepsilon})+O(\varepsilon^{2})\\ \beta=\beta_{0}-\varepsilon\omega\big{(}\beta_{0}\omega t+\alpha_{0}\omega(1-e^{-t/\varepsilon})\big{)}+O(\varepsilon^{2})\,.\end{cases}$$

$$\begin{gathered} \alpha =(1-4{\epsilon }^2{\omega }^2{)}^{-1/2}(-\left(\frac{1-\sqrt{1-4{\epsilon }^2{\omega }^2}}{2}{\alpha }_0+\epsilon \omega {\beta }_0\right)e^{-\frac{1-\sqrt{1-4{\epsilon }^2{\omega }^2}}{2\epsilon }t} \\ +\left(\frac{1-\sqrt{1+4{\epsilon }^2{\omega }^2}}{2}{\alpha }_0-\epsilon \omega {\beta }_0\right)e^{-\frac{1+\sqrt{1-4{\epsilon }^2{\omega }^2}}{2\epsilon }t}); \end{gathered}$$

$$\begin{gathered} \beta =(1-4{\epsilon }^2{\omega }^2{)}^{-1/2}(\left(\frac{1+\sqrt{1-4{\epsilon }^2{\omega }^2}}{2}{\beta }_0-\epsilon \omega {\alpha }_0\right)e^{-\frac{1-\sqrt{1-4{\epsilon }^2{\omega }^2}}{2\epsilon }t} \\ +\left(-\frac{1-\sqrt{1-4{\epsilon }^2{\omega }^2}}{2}{\beta }_0+\epsilon \omega {\alpha }_0\right)e^{-\frac{1+\sqrt{1-4{\epsilon }^2{\omega }^2}}{2\epsilon }t}). \end{gathered}$$

$$\alpha$$

$$\beta$$

Full text

Keywords: linear velocity constraints, knife-edge, Chaplygin sleigh, viscous friction, small parameter, multiple time scales, fast and slow motion, slow manifold, matching asymptotics, drift dynamics

Introduction

Nonholonomic mechanics is experiencing something of a revival in recent decades, see [3] for a historical perspective, and [18, 19] for some of the new themes. In this paper we look at approximating nonholonomic constraints by anisotropic viscous friction for plane ”dumbbells”, pairs of masses joined by a weightless connecting mechanism of some sort [11, 2.1]. They are traditional toy models for exploring various mechanical effects and behaviors, see e.g. [4, 11, 23], and a natural starting point for gaining insight into approximations of nonholonomic systems.

The idea of de-idealizing nonholonomic constraints by ”realizing” (replacing) them with forces of viscous friction goes at least as far back as Carathéodory [6]. Carathéodory considered the case of the Chaplygin sleigh and concluded that motions of the skidding sleigh, with knife-edge constraints replaced by viscous friction forces, did not converge to the constrained motions. After the development of perturbation theory for systems of ordinary differential equations in the 1950-60s Fufaev [10] analyzed the skidding sleigh in terms of slow and fast motions, and showed that convergence does in fact take place for $t>0$; this is presented in his book with Neimark [20, IV.3]. Both the Chaplygin and the Carathéodory-Fufaev sleighs are equivalent to dumbbells of the type we consider.

In 1981 Brendelev [5] and Karapetian [14], see also [1, 1.6], [17], generalized Fufaev’s analysis to general systems with linear velocity constraints. If one realizes them by viscous forces then in the limit of infinite friction frictional motions converge to constrained motions for positive time. The idea is that the viscous system quickly evolves towards the ”slow manifold” of the constrained system (fast motion), and then remains in its vicinity (slow motion). More recently, Eldering [8] clarified the nature of convergence to nonholonomic dynamics by interpreting the results of Brendelev and Karapetian in terms of the geometric theory of singular perturbations, which goes back to Fenichel [9]. Also recently, Deppler et al. considered realization of constraints by more general viscoelastic forces, which may provide a closer approximation of the actual physics [7].

But analysis of concrete examples in light of what the general theory implies for them is hard to come by, in particular it turns out that the convergence involved does not guarantee the kind of approximating behavior one might expect. The general results provide convergence on a finite time interval, but the relations between fast and slow motions, estimation of their time scales, and long term behavior remain largely out of the picture. A quantitative account of transient effects occurring in such realizations was developed in the theory of singular perturbations by O’Malley and Vasil’eva, see e.g. [21, Ch. 8], by using perturbative expansions and two time scales. In applications their approach came to be known as the method of matching asymptotics. Unfortunately, unlike applications to celestial mechanics and electrical engineering [13, 21], biology [12] and chemistry [22], applications to nonholonomic dynamics appear to be understudied and underappreciated in the literature. We will show that the choice of appropriate time scales is non-trivial even in the simplest cases, complicating application of matching asymptotics, and the long term behaviors of nonholonomic systems and their realizations may differ. Moreover, we encounter interesting situations, where secular terms appear in perturbative expansions of equations that are non-oscillatory, so the standard methods for eliminating them, such as averaging [13, 3.2], [21, Ch.11], do not apply.

In this paper we first present a general approach for reducing the equations of motion of plane dumbbells to a form that simplifies analysis of their solutions (Section 1). This involves representing positions of the masses as complex numbers, and manipulating equations in complex form before specific generalized coordinates are selected to take advantage of their symmetries and integrals of motion. We then focus on two examples, which we call the ”cart sleigh” (Sections 2-3) and the ”double spear” (Section 4), where the equations of motion turn out to be completely integrable, i.e. one can find analytic solutions. This allows us to compare approximate solutions obtained by the method of matching asymptotics to the exact ones, and assess their range of applicability precisely. This leads to a better understanding of the nature of convergence of motions in frictional realizations to the limiting motions of nonholonomic systems.

The two examples we chose display two in some sense opposite behaviors in approximating nonholonomic systems in the limit of infinite friction. Whereas the skidding double spear ”shadows” (approximates uniformly in time) its nonholonomic limit, the long term behavior of the skidding cart sleigh is completely different for arbitrarily large friction, despite the convergence of motions on finite time intervals. The latter effect seems to be underappreciated in the standard descriptions in terms of fast/slow motions, and limits the sense in which such descriptions can be ”trusted”.

Our analysis is based on comparisons to analytic solutions. Discriminating between these two cases based on the equations of motion directly, and finding effective methods for approximating long term behavior of frictional solutions emerge as interesting open problems. Our conclusions are summarized in Section 5.

1 Holonomic and nonholonmic dumbbells

By a dumbbell we understand a pair of masses joined by a weightless connecting mechanism. The simplest examples are rigid rods, telescoping rods with or without damping, or springs with zero or non-zero equilibrium lengths. The connection is meant to physically implement forces between the masses acting along the segment connecting them. The masses themselves can be thought of as small balls free to move without friction along the plane, or as small balls mounted on knife-edges that constrain directions of their instantaneous velocities. Since we are interested in realization of such constraints strict constraints may be replaced with forces, e.g. of viscous friction, that approximate them. The knife-edges produce non-holonomic constraints linear in velocities, and we will only consider cases where they are either perpendicular or parallel to the connecting segments, although our approach generalizes to cases where they are attached at some other fixed angle. We only consider dumbbell motion in the plane. Schematic diagrams for various types of dumbbells are presented on Fig.1.

We start by describing how the equations of motion are derived in vector form. Let $r_{1},r_{2}$ be the vectors of positions, and $m_{1}$, $m_{2}$ be the masses of a dumbbell, then the kinetic energy is $T=\frac{1}{2}m_{1}|\dot{r}_{1}|^{2}+\frac{1}{2}m_{2}|\dot{r}_{2}|^{2}$. For simplicity we assume that the masses are always equal, $m_{1}=m_{2}=m$, so $T=\frac{m}{2}(|\dot{r}_{1}|^{2}+|\dot{r}_{2}|^{2})$. Perpendicular knife-edge is represented by the constraint $\dot{r}_{i}\perp(r_{1}-r_{2})$, or in terms of the dot product $\dot{r}_{i}\cdot(r_{1}-r_{2})=0$. For the parallel knife edge $\dot{r}_{i}\parallel(r_{1}-r_{2})$. To give a dot product form it is convenient to introduce a linear transformation $J=\begin{pmatrix}0&-1\\ 1&0\end{pmatrix}$ that rotates every vector by $90^{\circ}$ counterclockwise, note that $J^{*}=J^{-1}=-J$. The constraint can then be represented by $\dot{r}_{i}\cdot J(r_{1}-r_{2})=0$. The rigid rod constraint is holonomic, $|r_{1}-r_{2}|=\text{const}$, but it is often convenient to work with its differential consequence obtained by differentiating $(r_{1}-r_{2})\cdot(r_{1}-r_{2})=\text{const}$, namely $(\dot{r_{1}}-\dot{r_{2}})\cdot(r_{1}-r_{2})=0$.

The forces of viscous friction realizing knife-edges are described by the Rayleigh dissipation function [11, 2.4]. For the knife-edge attached to the $i$-th mass perpendicular to the connecting rod the Rayleigh function is $R_{i}:=\frac{c}{2}|\text{pr}_{r_{1}-r_{2}}\dot{r}_{i}|^{2}$ , where $c>0$ is the friction coefficient, and $\text{pr}_{a}b$ denotes the orthogonal projection of $b$ to $a$. Since $\text{pr}_{a}b=\frac{b\cdot a}{|a|^{2}}a$ we have explicitly

[TABLE]

For parallel knife edges $r_{1}-r_{2}$ is replaced by $J(r_{1}-r_{2})$, note that $|J(r_{1}-r_{2})|=|r_{1}-r_{2}|$. Similarly, for damped rods the Rayleigh function is

[TABLE]

In the absence of constraints the equations of motion in vector form are given by

[TABLE]

where $\frac{\partial}{\partial a}$ stands for the partial gradient along the vector $a$, and $L:=T-U$ is the Lagrangian function. A straightforward computation shows that $\frac{\partial|Aa|^{2}}{\partial a}=2A^{*}Aa$, where $A$ is a linear transformation, and this along with the chain rule is enough to differentiate all the functions that we need in this paper. In particular,

[TABLE]

because $\text{pr}_{a}^{*}=\text{pr}_{a}^{2}=\text{pr}_{a}$ for orthogonal projections.

The constraints are handled according to d’Alembert’s principle of zero virtual work. For the constraints linear in velocities that we are considering the following rule (sometimes called Jourdain’s principle [2], [11, 2.4]) suffices: take the variation of the constraint with respect to velocities, multiply it by a Lagrange multiplier, and add the term in the dot product with $\delta\dot{r_{i}}$ to the $i$-th equation in (1). For example, to implement the differentiated rigid rod constraint we take the variation:

[TABLE]

Thus, we will add $\lambda(r_{1}-r_{2})$ to the equation for $r_{1}$, and $-\lambda(r_{1}-r_{2})$ to the one for $r_{2}$. This works even for non-linear constraints [2].

The reasons we prefer to write equations in vector form initially rather than go straight to the generalized coordinates, as e.g. in [4], are twofold. Vector equations provide geometric insight and manifest symmetries that are often lost once some specific choice of generalized coordinates is made. Second, it is not always immediately clear what generalized coordinates are most beneficial for analyzing or solving the system. We shall see that when working from vector equations such choices often ”suggest themselves”. This can be seen as an intuitive version of non-holonomic reduction without the abstract formalism and technicalities of the general case [15]. An important ingredient in this reduction process is the observation that vectors in the plane can be interpreted as complex numbers. Then applying $J=\begin{pmatrix}0&-1\\ 1&0\end{pmatrix}$ to a vector corresponds to multiplying it by the imaginary unit $i$, and the vector equations turn into scalar equations for complex numbers. The dot product also has a simple expression in terms of complex multiplication: $z\cdot w=\text{Re}[z\overline{w}]=\text{Re}[\overline{z}w]$. Complex numbers can be represented in Cartesian, exponential and polar forms, which provides a rich selection of real variables to choose from for generalized coordinates, and since complex numbers can be also multiplied and divided this selection is enriched by applying the idea to ratios of the original vectors and their combinations.

2 Cart sleigh

Let us illustrate the outline of the previous section with a simple example of a telescopic rod dumbbell with perpendicular knife-edges. One can think of this dumbbell as a toy model of a sleigh with short sharp rails placed on ice (or a pair of skates affixed to a rod). If the knife-edges are replaced with wheels one gets a common simple model of a two-wheeled cart [20, III.3]; we therefore call this dumbbell the cart sleigh.

The Lagrangian of the cart sleigh is just the kinetic energy $L=\frac{m}{2}(|\dot{r}_{1}|^{2}+|\dot{r}_{1}|^{2})$, the knife-edges are implemented by constraints $\dot{r}_{1}\cdot(r_{1}-r_{2})=0$ and $\dot{r}_{2}\cdot(r_{1}-r_{2})=0$. Even before writing the equations of motion subtracting the second constraint from the first gives us $(\dot{r_{1}}-\dot{r_{2}})\cdot(r_{1}-r_{2})=0$, which is none other than the differential form of the rigid rod constraint. In other words, in this case a telescopic rod behaves as if it was rigid! Now $\frac{d}{dt}\frac{\partial L}{\partial\dot{r}_{i}}=m\ddot{r}_{i}$, and applying Jourdain’s principle to each constraint we get the equations of motion

[TABLE]

Since we have $r_{1}-r_{2}$ in both equations it seems natural to introduce a new vector variable $r_{12}:=\frac{1}{2}(r_{1}-r_{2})$ (the reason for $\frac{1}{2}$ will become clear shortly). Subtracting we get a self-contained equation $m\ddot{r}_{12}+(\lambda_{1}-\lambda_{2})r_{12}=0$ with the transformed constraint $\dot{r}_{12}\cdot r_{12}=0$. We now need a complementary variable, for which $r:=\frac{1}{2}(r_{1}+r_{2})$ is a natural candidate. The $\frac{1}{2}$ coefficients ensure that both $r_{1},r_{2}$ can then be recovered without fractions: $r_{1}=r+r_{12}$, $r_{2}=r-r_{12}$, and $r$ has the physical interpretation of being the center of mass of the dumbbell (and the midpoint of the rod). Adding the equations and the constraints we get $m\ddot{r}+(\lambda_{1}+\lambda_{2})r_{12}=0$ and $\dot{r}\cdot r_{12}=0$, which directly tells us that the velocity of the center is always perpendicular to the rod. Since the Lagrange multipliers are also unknown we might as well set $\lambda_{12}:=\lambda_{1}-\lambda_{2}$, $\lambda:=\lambda_{1}+\lambda_{2}$, which results in the system:

[TABLE]

Some information about the motion can be extracted directly from this vector system. Taking dot product of the first equation with $\dot{r}_{12}$, and the second one with $\dot{r}$ we get right away that $\ddot{r}_{12}\cdot\dot{r}_{12}=\frac{d}{dt}\frac{1}{2}|\dot{r}_{12}|^{2}=0$, so $|\dot{r}_{12}|=\text{const}$, and similarly $|\dot{r}|=\text{const}$, i.e. both the center, and each mass relative to the center (which is what $\pm r_{12}$ represent) move with constant speeds.

Now it is time for complex notation. Since $|r_{12}|=\rho=\text{const}$ it is natural to choose exponential form for $r_{12}:=\rho e^{i\theta}$, then

[TABLE]

With this notation the first equation in (3) is

[TABLE]

As $\lambda_{12}$ is real-valued we get from separating real and imaginary parts

[TABLE]

Since $\theta$ enters only through $\dot{\theta}$ it is natural to set $\omega:=\dot{\theta}$, so that $\lambda_{12}=m\omega^{2}$ and $\dot{\omega}=0$, i.e. $\omega=\text{const}$. The dumbbell rotates with constant angular velocity (this is expected from the constant speed and length of $r_{12}$).

We also notice that $r$ enters the equations only through $\dot{r}$, $|\dot{r}|=\text{const}$, and $\dot{r}\perp r_{12}$. This means that $\dot{r}$ is a constant (real!) multiple of $ir_{12}$, $\dot{r}=i\beta r_{12}$ $\ddot{r}=i\beta\dot{r}_{12}=-\beta\omega\rho e^{i\theta}$ reducing the second equation to $-m\beta\omega\rho e^{i\theta}+\lambda\rho e^{i\theta}=0$, so $\lambda=\beta m\omega$. The values of $\beta$ and $\omega$ can be found from initial conditions. We now have ${r}_{12}=\rho e^{i(\theta_{0}+\omega t)}$, and $\dot{r}=i\beta\rho e^{i(\theta_{0}+\omega t)}$, which yields

[TABLE]

This means that the center of mass uniformly rotates along a circle with the center $r(0)-\beta\rho e^{i\theta_{0}}/\omega=r(0)-\dot{r}(0)/\omega$ of radius $\dot{r}(0)/\omega$, while the masses uniformly rotate along a circle centered at the moving point of the radius. Their resulting motion is therefore epicyclic (as in geocentric astronomical models), with $r$ tracing the deferent, and $r_{12}$ the epicycle. Due to the equal frequencies and phases this is a very special case of epicyclic motion however. For the masses the resulting motion is

[TABLE]

i.e. the resulting trajectory is itself a circle with the same center as the deferent, and the radius smaller/larger by the radius of the epicycle. Thus, the center of mass and both masses move along concentric circles with constant angular velocities, Fig. 3. The above discussion applies to the case $\omega\neq 0$, the $\omega=0$ case can be obtained by taking the limit $\omega\to 0$ in (4), which describes uniform rectilinear motion in the direction of the initial push.

Of course, for the cart sleigh one can predict the final answer using physical intuition, but things usually do not work out so simply and neatly. Still, in a number of cases the reduction procedure outlined above leads to an analytic solution of dumbbell equations, and even when it does not at least some qualitative features and integrals of motion can be found as a matter of course. Moreover, when constraints are realized by large viscous forces, which is our point of interest, equations are reduced to a form where asymptotic methods of perturbation theory can be applied fruitfully by encapsulating friction coefficients into small parameters, with nonholonomic limits representing the unperturbed motion. The viscous motions then reveal some surprises already in the case of the cart sleigh, as we show next.

3 The sleigh skids

Consider what happens if we kick the cart sleigh in a way that produces initial velocities inconsistent with the constraints? The formal answer is that we can not, literally. The equations of motion for a constrained system with initial conditions not conforming to the constraints have no solutions. At least no classical solutions, but it is also unclear what ”non-classical solutions” might mean here. This is not to say that the sleigh simply can not start moving, that would still be a solution, the trivial one, while no solution means that the theory literally predicts nothing in this case. But under the usual intuition of mechanical idealizations a ”kick” amounts to an ”infinite” impulse force imparting finite instantaneous velocity, and it can act in any direction. So ideally speaking we ”can” in fact kick the cart sleigh in a non-conforming way, and it ”should” do something. One could suggest that the ”infinite” reaction force in response to the ”kick” cancels its component normal to the constraints, so it instantly turns initial conditions into conforming ones. But does this intuition reflect the behavior under what the ”kick” is supposed to approximate?

In this section we will investigate what happens if the constraint idealization is relaxed in the spirit of Carathéodory and Fufaev, i.e. if skidding is allowed. We realize the constraints with large viscous friction, and discuss if it makes sense to talk about ”solutions” to non-holonomic equations of motion (as limits) when the initial values do not conform to the constraints. After deriving the equations of motion we first apply the popular method of matching asymptotics to it, which is a standard method in singular perturbation theory [13, 3.2], [21, Ch.8], but is rarely if ever used in the context of non-holonomic dynamics. Part of our goal is to call attention to it, and to highlight the advantages and the challenges of its application. Then we derive the analytic solution and compare it to the matched asymptotic expansion, which illustrates both its strengths and limitations.

In the case of the cart sleigh it made no difference whether we allowed the connecting rod to telescope or made it rigid, but when the knife-edges are allowed to skid it does make a difference. We pick the rigid rod realization, which is simpler. The new elements are the viscous forces, which replace the knife-edge constraints and are given by the Rayleigh function

[TABLE]

and the explicit rigid rod constraint in the differential form $(\dot{r_{1}}-\dot{r_{2}})\cdot(r_{1}-r_{2})=0$. The vector equations of motion are obtained as before (see (2) for taking the derivative of the Rayleigh function):

[TABLE]

After the familiar substitution $r_{12}:=\frac{1}{2}(r_{1}-r_{2})$, $r:=\frac{1}{2}(r_{1}+r_{2})$, adding and subtracting the equations while taking into account that $\dot{r}_{12}\cdot r_{12}=0$ results in

[TABLE]

Interpreting $r,r_{12}$ as complex numbers we set $r_{12}=\rho e^{i\theta}$ with $\rho=\text{const}$, and the first equation tells us that $\omega:=\dot{\theta}=\text{const}$ and $\lambda=m\omega^{2}/2$ as in the constrained case. The first equation can now be solved explicitly, and gives $r_{12}(t)=r_{12}(0)e^{i\theta}=\rho e^{i(\theta_{0}+\omega t)}$. Unlike before, however, we do not have that $\dot{r}\perp r_{12}$ or that $|\dot{r}|=\text{const}$. Nonetheless, since we expect $\dot{r}=i\beta r_{12}$ in the limit let us set $\dot{r}=zr_{12}$, where $z:=\alpha+i\beta$ is now a variable, and complex valued. Then

[TABLE]

so after separating the real and imaginary parts the second equation reduces to a pair of equations for $\alpha,\beta$:

[TABLE]

Since we plan to investigate what happens when $c\to\infty$ it makes sense to introduce a small parameter $\varepsilon:=m/c$, and rewrite the system in the Tikhonov form [21, Ch.8]

[TABLE]

This is the reduced system we will analyze. The system is linear and can be solved exactly, but as our goal is general insight let us disregard that for the moment, and apply a method that would work for non-linear systems as well. Setting $\varepsilon=0$ (i.e. $c=\infty$) gives $\alpha=0$ and $\dot{\beta}=0$, i.e. $\beta=\text{const}$: this is the cart sleigh considered in the previous section, as expected. To get the next order of approximation in $\varepsilon$ we proceed as follows. Set $\alpha=\alpha^{0}+\varepsilon\alpha^{1}+\dots$, $\beta=\beta^{0}+\varepsilon\beta^{1}+\dots$, so upon substitution into (5) we have

[TABLE]

We then equate the terms in each order of $\varepsilon$.

The zero order terms simply reproduce what we found above: $\alpha^{0}=0$, $\dot{\beta}^{0}=0$, and in the first order we get

[TABLE]

Let $\beta^{0}=B_{0}=\text{const}$, then $\alpha^{1}=\omega B_{0}$ and $\dot{\beta}^{1}=-\omega^{2}B_{0}$, so $\dot{\beta}^{1}=-\omega^{2}B_{0}t+B_{1}$, where $B_{1}$ is another constant. Thus, to the first order in $\varepsilon$ we find

[TABLE]

We could try to find $B_{0},B_{1}$ by using initial values for $\alpha,\beta$, but that would be premature. The theory implies that this regular expansion is only valid for ”slow motions”, into which the system settles after a short transient period, so the integration constants are not directly related to the initial values.

To find the transient ”fast motions” the standard approach is to introduce stretched time $\tau:=t/\varepsilon$, then $\frac{d}{dt}=\frac{1}{\varepsilon}\frac{d}{d\tau}$. Denoting derivatives with respect to $\tau$ by ′ we transform (5) into

[TABLE]

Expanding as in (6) produces to the first order in $\varepsilon$:

[TABLE]

Thus, $\beta^{0}=\text{const}$, $\alpha^{0}=A_{0}e^{-\tau}$, and since these are fast motions we are now justified to find the integration constants from the initial values, i.e.

[TABLE]

We now turn to matching. The idea is that right after the initial kick the system undergoes the fast motion, and its asymptotic values of $\alpha,\beta$ become the initial values for the subsequent slow motion. This finally allows us to connect the constants $B_{0},B_{1}$ to the initial values.

Note that for the slow motions $\alpha^{0}=0$, which matches with the limit of the fast motion at temporal infinity, $\lim_{\tau\to\infty}\alpha_{0}\,e^{-\tau}=0$. The same type of matching for $\beta^{0}$, which remains constant in the zero order during the fast motion gives $B_{0}=\beta_{0}=\beta(0)$. Since the initial values are accounted for in the zero order, in the first order we should set $\alpha^{1}(0)=\beta^{1}(0)=0$. Solving (9) we then have

[TABLE]

When $\tau\to\infty$ this should match first order terms in (7) with $t\to 0$. For $\alpha^{1}$ the match is automatic, and for $\beta^{1}$ we determine $B_{1}=-\alpha_{0}\omega$. Combining (7), (10) and (11) we get what is called the two timing or matched perturbative expansion to the first order in $\varepsilon$:

[TABLE]

The approach we outlined is the method of matching asymptotics mathematically justified by O’Malley and Vasil’eva in the case of singular perturbations. The expansion is valid as asymptotic expansion on a finite time interval $[0,T]$ with $T\sim O(1)$, see [21, 8.3]. The last restriction on the length of time interval is not a formality as we will see.

We can now attempt to use fast motion to define ”generalized solutions” for constrained systems with non-conforming initial values. They should be the classical solutions with initial values replaced by the values ”at the end” of the fast motion, that is by their limits as $\tau\to\infty$. In our example this amounts to replacing $(\alpha_{0},\beta_{0})$ with $(0,\beta_{0})$. This is equivalent to the intuition of reaction forces instantly canceling the component of the initial ”kick” orthogonal to the constraints. But it turns out that if we do so we can not count on the frictional solutions staying close to the constrained solutions beyond $T\sim O(1)$.

To see this let us consider system (5) from a different perspective. The phase portrait of the system in the $\alpha$-$\beta$ plane along with a sample trajectory is shown on Fig.5. As we already remarked, (5) is linear and can be solved exactly. Looking for solutions of the form $\beta=e^{\xi t}$ we find $\alpha=-\frac{\xi}{\omega}\,e^{\xi t}$, and the characteristic equation for $\xi$, namely $\varepsilon\xi^{2}+\xi+\omega^{2}=0$. Its two roots are $\xi_{1,2}=-\frac{1\pm\sqrt{1-4\varepsilon^{2}\omega^{2}}}{2\varepsilon}$, both real and negative for $\varepsilon\omega<1/2$. The general solution is of the form $\beta=b_{1}e^{\xi_{1}t}+b_{2}e^{\xi_{2}t}$, and enforcing initial values one finds:

[TABLE]

[TABLE]

The formulas are cumbersome, but one can see by inspection that for any $\varepsilon>0$ the exponents are strictly negative, so $\alpha(t),\beta(t)\xrightarrow[t\to\infty]{}0$. Intuitively, this is apparent from the phase portrait, where the origin is a globally attractive equilibrium. Note that the constrained cart sleigh is represented by the $\beta$-axis, with constrained dynamics consisting of fixed points. Recall that $\dot{r}=(\alpha+i\beta)r_{12}$, and $r_{12}(t)=\rho e^{i\theta}$, which means that $\dot{r}(t)\xrightarrow[t\to\infty]{}0$. Unlike for the constrained sleigh of the previous section the skidding sleigh’s center of mass tends to rest asymptotically no matter how large the viscous friction is, and despite the fact that the friction forces quickly become ”almost” orthogonal to $\dot{r}$, so should do ”almost” no work against the overall motion of the sleigh. Thus, the ”generalized solutions” to the constrained equations will drift further and further apart from the solutions for $\varepsilon>0$ as time goes on despite the latter getting closer and closer to the slow manifold of the constrained solutions $\alpha=0$. This is because they are approaching it by approaching the origin, while for the constrained motion $\beta$ remains constant, and so does not approach [math].

As one can see from Fig. 6 for a typical trajectory masses move along winding in spirals, and the center of mass asymptotically approaches a point, around which the masses spin with angular velocity $\omega$. This should be compared to Fig. 3 depicting the motion of a constrained sleigh, where the masses and the center of mass are moving along concentric circles. Note that even the skidding sleigh does not come to a stop despite losing energy to the constraints. Instead, the system ”finds” an asymptotic motion where constraints no longer do any work, and the residual energy is conserved. In this case it can be attributed to the idealization of the ”infinitesimal length” knife-edges. Generally speaking, the energy is conserved in constrained motions, so as long as there exist non-trivial ones we can expect non-trivial asymptotic motions that conserve energy.

An issue with matching asymptotics is already apparent from the second equation in (12), which includes a term proportional to $t$ in the first order of $\varepsilon$. If we take this term at face value then for large $t$ our $\beta$, and hence $\dot{r}$, would become arbitrarily large as $t\to\infty$. Terms of this nature first appeared in celestial mechanics when approximating the motion of the planets, and came to be called ”secular terms”. Their presence indicates that (12) can not be relied upon for large $t$. At the end of 19th century Lindstedt and Poincare developed a method for ”eliminating” secular terms, i.e. constructing expansions that do not contain them, and therefore have a chance of holding uniformly in time. Unfortunately, their approach, and the method of averaging developed later, rely on detecting resonances and avoiding them in expansions, which is specific to equations that in their unperturbed form describe something close to periodic motions [13, 3.1], [21, 10.1]. But (5) is non-oscillatory for $\varepsilon=0$ or small $\varepsilon>0$, so the standard apparatus of averaging does not apply. In fact, this is a general feature one can expect from perturbative expansions for realization of non-holonomic constraints by forces of viscous friction: all the pain of secular terms, no benefit of averaging.

We can sidestep the difficulty in this case only because we have the exact analytic solution. To better understand long term behavior let us expand (13) to the first order in $\varepsilon$ using that $\sqrt{1-4\varepsilon^{2}\omega^{2}}=1-2\varepsilon^{2}\omega^{2}+O(\varepsilon^{4})$:

[TABLE]

This expansion, unlike (12), is not just $O(1)$, one can show that it converges to the exact solution for all $t>0$ when $\varepsilon\to 0$. The time scales manifestly present in (3) are the stretched time $t/\varepsilon$ of fast motion, and the compressed time $\varepsilon t$ of (very) slow motion. In the sense of matching asymptotics, which relies on the apparent form of the terms, the ordinary slow time scale $t$ is not present at all! When we ”forced” it into the expansion (12) by naively using matching asymptotics the result was a secular term that restricted its validity. This does not preclude (12) from holding for $t\sim O(1)$ of course, but its form is quite misleading the long term. The value of $\beta$ drifts to [math] on the very slow time scale of $\varepsilon t$, so on $t/\varepsilon$ and even $t$ scales it registers as staying constant at zero order. The secular term emerges in the first order to compensate for the drift.

Figure 7 displays graphs of the exact solution (13), its approximation (3), and the two timing perturbative expansion (12), both to the first order. As expected, the latter remains close to the analytic solution for small $t$, but then deviates from it significantly.

Note also that even the analytic solution (13) is only valid if $\varepsilon\omega<1/2$. When $\omega>1/2\varepsilon$, i.e. if the sleigh is made to spin fast enough initially, it is seen from the formula that the behavior changes qualitatively. The square root $\sqrt{1-4\varepsilon^{2}\omega^{2}}$ becomes imaginary, and $\beta$ undergoes damped oscillations with frequency $\sqrt{\omega^{2}-1/4\varepsilon^{2}}$. Their amplitudes are multiplied by $e^{-t/2\varepsilon}$ however, so they all but disappear after the fast motion. For such initial values both (12) and (3) are misleading even in the zero order, as they ”predict” constant $\beta^{0}$ during the fast motion, rather than damped oscillations. This is indicated by the presence of $\varepsilon\omega$ terms in the first order, which then will not be ”small”. This illustrates the issue discussed by Holmes [13, 3.2.4] in the context of multiple scales expansions: in a valid expansion its terms must be kept ”well-ordered” by magnitude.

Neither the correct time scales $t/\varepsilon$, $\varepsilon t$ nor the role of $\varepsilon\omega$ as the more adequate small parameter are apparent from the original system (5), which suggests that it is not the optimal form for studying the behavior of the skidding cart sleigh. In hindsight, the stretched time system (8) is a more attractive option since it displays the role of $\varepsilon\omega$ explicitly. We can do even better by multiplying the second equation by $-\varepsilon\omega$ and setting $\gamma:=-\varepsilon\omega\alpha$. Then the system becomes

[TABLE]

In this form the system is no longer singularly perturbed in the usual sense (it still is in a relevant technical sense [21, Ch.10]), it is apparent that the ”correct” small parameter is $\epsilon:=\varepsilon^{2}\omega^{2}$ rather than $\varepsilon$ or $\varepsilon\omega$, and the regular and compressed time scales $\tau:=t/\varepsilon$ and $\epsilon\tau:=\varepsilon\omega^{2}t$ are the ”correct” ones from (3). In particular, it is more effective in investigating the fast motion. Already in the zero order we have

[TABLE]

This is better than the zero order approximation in (12), and one may proceed with the matching aymptotics as we did above. However, in general the correct time scales can not be discerned from the system by such simple manipulations and more complicated techniques are required, e.g. treating different time scales as independent variables and solving partial differential equations involving them, see [13, 3.2.2], [21, 11.4].

4 Skidding double spear

In the cart sleigh both knife-edges were perpendicular to the connecting rod; now we will consider a ”double spear”, where they are both parallel. The constrained motion is restricted to a straight line and not very interesting, so we start directly with the skidding version, see Fig. 8.

The Rayleigh function now is:

[TABLE]

where $J=\begin{pmatrix}0&-1\\ 1&0\end{pmatrix}$ is the $90^{\circ}$ counterclockwise rotation matrix. The vector equations of motion are obtained as before (see (2) for taking the derivatives):

[TABLE]

where $\lambda$ is the Lagrange multiplier for the rigid rod constraint $(\dot{r_{1}}-\dot{r_{2}})\cdot(r_{1}-r_{2})=0$. In terms of the center of mass $r:=\frac{1}{2}(r_{1}+r_{2})$, and the radius vector $r_{12}:=\frac{1}{2}(r_{1}-r_{2})$ the system is

[TABLE]

As before we set $r_{12}=\rho e^{i\theta}$, which reduces the first equation to $m(-\dot{\theta}^{2}+i\ddot{\theta})+2\lambda+ic\dot{\theta}=0$, so with $\omega:=\dot{\theta}$ we have $\lambda=m\omega^{2}/2$ and $m\dot{\omega}+c\omega=0$. In contrast to the cart sleigh, the rotational angular velocity is not constant in this case.

Setting $\dot{r}=zr_{12}=(\alpha+i\beta)r_{12}$ reduces the second equation to

[TABLE]

and with $\varepsilon:=m/c$ we obtain

[TABLE]

The last equation decouples from the other two, and the solution to it is $\omega=\omega_{0}e^{-t/\varepsilon}$. We therefore have

[TABLE]

which means that $r_{12}(t)$ quickly settles into the limit value of $\rho e^{i(\theta_{0}+\varepsilon\omega_{0})}$. Moreover, (17) reduces to a non-autonomous system in two variables:

[TABLE]

The fast time scale $t/\varepsilon$ is now explicitly present in the coefficients, which makes the substitution $\tau:=t/\varepsilon$ all the more natural, so that with primes denoting $\tau$ derivatives

[TABLE]

With the benefit of past experience we set right away $\gamma:=\varepsilon\omega_{0}\beta$ and $\epsilon:=\varepsilon^{2}\omega_{0}^{2}$ to finally obtain

[TABLE]

The regular expansions to the first order are $\alpha=\alpha^{0}+\epsilon\alpha^{1}+\dots$, $\gamma=\gamma^{0}+\epsilon\gamma^{1}+\dots$, and upon the substitution into (21) we have

[TABLE]

In the zero order in $\epsilon$ this gives

[TABLE]

with solutions

[TABLE]

Here $\gamma_{0}:=\gamma(0)$ and $\alpha_{0}:=\alpha(0)$ are the initial values, which we substituted directly because $\tau$ is already the time scale of the fast motion. In the first order in $\epsilon$ we have

[TABLE]

and the initial values are $\alpha^{1}(0)=\gamma^{1}(0)=0$ since $\alpha_{0}$, $\gamma_{0}$ are already accounted for in the zero order. Therefore, the first order solutions are:

[TABLE]

Note the presence of secular terms, which are however suppressed by negative exponents this time. Taking $\tau\to\infty$ we obtain an approximation of the values ”after” the fast motion: $\gamma=0$ and

[TABLE]

The compressed time scale is $t_{1}=\epsilon\tau$ with $\frac{d}{d\tau}=\epsilon\frac{d}{dt_{1}}$. Still using primes for derivatives we have

[TABLE]

or substituting expansions:

[TABLE]

In the zero order we get $\gamma^{0}=0$, which matches the fast motion. In the first order

[TABLE]

hence $\gamma^{1}=-e^{-t_{1}/\epsilon}\alpha^{0}$ and $\alpha^{0\,\prime}=-e^{-2t_{1}/\epsilon}\alpha^{0}$, so $\alpha^{0}=ae^{\frac{\epsilon}{2}e^{-2t_{1}/\epsilon}}$ for some constant $a$. When $t_{1}=0$ we have $\alpha^{0}=ae^{\frac{\epsilon}{2}}=\alpha_{0}+\frac{\gamma_{0}}{2}$, and matching with the fast motion gives $\alpha=(\alpha_{0}+\frac{\gamma_{0}}{2})e^{-\frac{\epsilon}{2}(1-e^{-2t_{1}/\epsilon})}=(\alpha_{0}+\frac{\gamma_{0}}{2})e^{-\frac{\epsilon}{2}(1-e^{-2\tau})}$. As seen from the last expression no new time scale appears in the answer even when we try to ”force” it. This is in contrast to the case of the cart sleigh, and may suggest that the regular expansion in $\tau$ already works for all times uniformly. But even in this example it does not, albeit for a different reason than for the cart sleigh. To see why let us again compare to the analytic solution.

We start by transforming (21) into a single second order equation for $\alpha$:

[TABLE]

The substitution $x:=\epsilon^{1/2}\,e^{-\tau}$ transforms it further into $x\alpha_{xx}-\alpha_{x}+x\alpha=0$, which is a well-known equation in mathematical physics. It turns into one known even better after the substitution $\alpha:=xu$, namely into

[TABLE]

with $\nu=1$. Equation (28) is the equation of cylindrical (Bessel) functions [16]. Their appearance is somewhat surprising here because they relate to behavior in time, and have nothing to do with geometry. The two linearly independent solutions to (28) are the Bessel functions of the first and the second kind $J_{\nu}(x)$ and $Y_{\nu}(x)$, so the general solution can be written as $\alpha(x)=ax\,J_{1}(x)+bx\,Y_{1}(x)$ with $x=\epsilon^{1/2}\,e^{-\tau}$, and arbitrary constants $a,b$.

To solve for $a,b$ in terms of initial values we will use some standard properties of Bessel functions. Namely, the derivative formulas \big{(}x^{\nu}\,J_{\nu}(x)\big{)}^{\prime}=x^{\nu}\,J_{\nu-1}(x), \big{(}x^{\nu}\,Y_{\nu}(x)\big{)}^{\prime}=x^{\nu}\,Y_{\nu-1}(x), and the Wronskian identity:

[TABLE]

Note that $\frac{d}{d\tau}\alpha|_{\tau=0}=e^{-\tau}\gamma|_{\tau=0}=\gamma(0)=\gamma_{0}$, and $\frac{d}{d\tau}=-x\frac{d}{dx}$, therefore $\frac{d\alpha}{d\tau}(0)=-x\frac{d\alpha}{dx}|_{x=\epsilon^{1/2}}=-\epsilon^{1/2}\alpha_{x}(\epsilon^{1/2})$. We also have

[TABLE]

so the system for $a,b$ is

[TABLE]

Using Cramer’s rule and (29):

[TABLE]

Since $\epsilon=\varepsilon^{2}\omega^{2}_{0}$ and $\gamma=\varepsilon\omega_{0}\beta$ we have $\epsilon^{-1/2}\gamma_{0}=\text{sign}(\omega_{0})\beta_{0}$, and $\text{sign}(\omega_{0})\epsilon^{1/2}\gamma_{0}=\varepsilon\omega_{0}$. Thus,

[TABLE]

[TABLE]

Figure 9 shows the trajectories of the masses and of the center of mass for motion with the relatively large initial velocity to highlight the transient phase. As expected from (23) on the time scale of $t\sim\varepsilon$ after some spinning the spear settles into the nearly rectilinear motion characteristic of the constrained system. In other words, the frictional trajectory shadows a constrained one in this case.

Since formulas (30) are rather obscure we will make use of the following Taylor expansions at $x=0$:

[TABLE]

In (4) $C$ is the Euler constant. Clearly $J_{0},Y_{0}$ are even and $J_{1},Y_{1}$ are odd, so (30) will not change if all entries of $\epsilon^{1/2}$ are replace by $\varepsilon\omega_{0}$, and $\text{sign}(\omega_{0})$ is removed. Expanding into the powers of $\epsilon^{1/2}$ with this in mind we obtain

[TABLE]

Note that the logarithmic terms appearing in (4) cancel out in (4) to the order displayed. Taking into account that $\gamma=\varepsilon\omega_{0}\beta$, $\gamma_{0}=\varepsilon\omega_{0}\beta_{0}$ we also see that it matches the regular expansion obtained in (23) and (25) to this order.

Due to the exponential suppressing factor in $x=\epsilon^{1/2}e^{-\tau}$ the uniform convergence of (4) in time is implied by the mere local convergence of the Taylor series in $x$ at [math]. Keep in mind however that logarithmic terms do appear in higher orders, while they are by definition absent from the regular expansion, which therefore does not converge uniformly as mentioned earlier, despite the exponential suppression of secular terms. The challenge however is to extract such information from the equations of motion (5) and (17), rather than from the analytic solutions. On Fig.10 graphs of the regular perturbative, and the analytic expansions (to the second order) coincide, but both slightly but visibly deviate from the analytic solution for large times.

5 Conclusions

We showed that plane dumbbells provide useful toy models for analyzing viscous realizations of non-holonomic systems, both for their technical simplicity, and richness of displayed phenomena. Analysis of their motion was shown to be greatly aided by treating positions of the masses as complex numbers, which leads to an elementary version of symmetry reduction to suitable generalized coordinates. The reduced equations of motion are more tractable, and sometimes linear and/or can be solved analytically.

Even when an analytic solution is unavailable methods of singular perturbation theory can be used to analyze large viscosity approximations of non-holonomic motion quantitatively. This is because the ratio $m/c$ serves as a natural small parameter. A number of approximating methods have been developed for analyzing systems with a small parameter. Let us briefly survey some of them to provide context.

Perhaps the most straightforward method is the regular (Taylor) expansion of the solution in the powers of the parameter with iterative determinations of coefficients from the equations and initial conditions. Unfortunately, the resulting expansion works poorly in the long term, and not at all for singularly perturbed systems. A modification of it, known as the ”two timing” expansion, was theoretically justified by O’Malley and Vasil’eva in 1960s based on Tikhonov’s theory of singularly perturbed systems [21, 8.3]. We used Vasil’eva’s procedure for determining the coefficients of these expansions, known as ”matching asymptotics”. Although our examples are, in the end, linear, the generalization to non-linear systems is straightforward. The method is quite attractive for analyzing transient motions because it is straightforward to apply, and provides good accuracy. Nonetheless, the method’s utility is limited for a number of reasons: appearance of more than two time scales, especially in non-linear systems, secular terms, logarithmic behavior, etc. Moreover, equations of motion can be misleading, the time scales ”apparent” from them may not be the time scales that best represent the system’s behavior.

For oscillatory systems some refinements of the expansion methods aimed at long term approximations were developed, such as the Poincare-Lindstedt method and the method of averaging [13, 3.1], [21, 10.1]. Unfortunately, they are specific to oscillatory problems, as they are based on ”averaging out” oscillations when determining coefficients to eliminate the secular terms. Both methods can be seen as shortcuts for the general method of multiple time scales, where the latter are treated as independent variables to be determined by solving partial differential equations [13, 3.2.2], [21, 11.4]. The method is quite cumbersome even in linear, and multiply so in non-linear, cases.

One of the main reasons the matched perturbative expansions do not represent the motion on longer time scales is the appearance of the so-called secular terms, whose exponents are multiplied by polynomials. As we saw, in the examples with dumbbells, viscously perturbed systems are quite peculiar in this regard. Neither the unperturbed (constrained) solutions nor their perturbations are oscillatory, as in the paradigmatic examples of secular terms in celestial mechanics and electric engineering. We believe that the phenomenon is general and calls for a new approach to elimination of secular terms, specific to this context, which can serve as shortcuts to be solving partial differential equations for appropriate time scale.

In a different direction, it turned out that the attractive idea of using viscous realization to define ”generalized solutions” to constrained systems for initial values not conforming to the constraints, does not work in general. One problem is that the frictional trajectories do not always ”shadow”, i.e. approximate uniformly in time, any constrained trajectory, including those obtained by replacing non-conforming initial values by conforming ones ”at the end” of the fast motion. The overall picture that emerges is this. The theory implies that on a fixed time interval viscous motions converge to constrained motions when the damping coefficient $c$ goes to infinity [5, 14]. One can split the viscous motion into initial fast motion, during which the system quickly approaches the manifold of constrained motions, and subsequent slow motion in its neighborhood. One might expect that during the latter the viscous motion shadows some non-holonomic trajectory in the constrained manifold.

However, as we saw in the example of the skidding cart sleigh, just because a system stays close to the constrained manifold does not mean that it stays close to any particular trajectory in it, i.e. there may not be any solution to the constrained equations of motion that viscous solutions converge to. Indeed, one can see from (13) that for any $\varepsilon>0$ we have $\beta_{\varepsilon}(t)\xrightarrow[t\to\infty]{}0$, whereas in the limit $\beta=\beta_{0}=\text{const}$. In other words, no matter how small $\varepsilon$ is, for $t$ large enough the viscous trajectory will come apart from any constrained trajectory, despite remaining close to the constrained manifold at all times. In fact, the global behavior of the system with viscous friction is qualitatively different from that of the constrained system. In the former the sleigh asymptotically comes to a stop, whereas in the latter it keeps going forever. Therefore, in general there is no constrained solution that frictional solutions approach in the limit of infinite damping starting from non-conforming initial values. This phenomenon is called drift or creep dynamics [8].

This means that the fast/slow motion picture does not always tell the whole story of non-holonomic constraints realized by viscous friction. There is an additional complication in that the manifold of constrained motions is typically non-compact, e.g. for the skidding sleigh it is $\{(\alpha,\beta,\omega)\ |\ \alpha=0\}$ (at the level of reduced variables), and we generally expect a non-compact range of motions for constrained systems. Perturbation theory for non-compact invariant manifolds is developed in [8]. However, while such manifolds persist under perturbations, the same can not be said about global attractors which can change abruptly, and it is their structure that determines asymptotic behavior. And this can happen even when the configuration space is compact.

We expect the same conclusions to hold for the viscoelastic realization of constraints considered in [7], where convergence is also proved only for $t\sim O(1)$.

In the light of this let us look back at the classical example of the Chaplygin sleigh, which is equivalent to the dumbbell in Fig. 11. Carathéodory [6] reasoned heuristically that for arbitrarily large damping the trajectories of the skidding sleigh differ dramatically from the ideal ones. Fufaev [10], see also [20, IV.3], gave a fast/slow motion analysis of the situation explaining that the viscous motion does in fact approximate the constrained one shortly after $t=0$ (although not at $t=0$), which undermined Carathéodory’s reasoning. Fufaev did not, however, analyze long term behavior of the skidding sleigh other than to say that it stays close to the constrained manifold. As we saw for simpler dumbbells, the two motions may yet come apart. Could Carathéodory have been right, if only for a wrong reason?

Numerical simulations in [8] suggest that for velocities the answer is negative, and there is shadowing just as in the case of the double spear (Figure 8 in [8] shows that the asymptotic directions of motion differ, and so there is no shadowing for position variables). But the Carathéodory-Fufaev sleigh equations are more involved, and analytic solutions to them are not available. It would be interesting to find the answer definitively, and even more desirable to find a general approach to answering such questions without the recourse to analytic solutions.

Acknowledgement: The authors are grateful to Jaap Eldering for multiple suggestions and corrections during the preparation of the paper.