## Dynamical Systems and Constraints: Part 2 - Bead on a Wire

This post begins a multi-part exploration of the role of constraints in mechanical systems. It builds on the general framework discussed in the last post but does so by examining specific examples. This month the focus will be on the simplest example of a constraint: a bead sliding on a wire.

Despite its trivial appearance, this problem is surprisingly complex due to the presence of the constraint and, as will be shown below, provides significant obstacles to obtaining an accurate numerical solution.

The following figure shows a portion of the wire, the position, , and the velocity, , of the bead, and the two forces acting upon it; namely the force of gravity and the the normal force . With this information, it is, in principle, easy to solve for the motion of the bead via Newton's laws. There are several equivalent ways of carrying this out but I'll be focusing on the most geometric of these.

The position of a particle in a plane is given by

Since the bead must remain on the wire, only those values of satisfying are permissible and the position can now be written as

The velocity is then given by

The normal to the curve, , is obtained from the requirement that motion consistent with the constraint must always be perpendicular to the normal to the wire at the point in question:

Using the constraint equation and then taking the differential yields the component form

from which it is easy to read off the components of the normal (the sign being chosen so that is in opposition . Subsequent normalization gives

The forces on the bead are due to gravity

and the normal force

For notational convenience in what follows, define

Summing the forces and inserting them into Newton's second law leads to the following equations (eliminating from both sides)

Expressing Newton's equations in state-space form then yields

These equations are easy to code and couple to a standard solver. I chose to implement this system in Python and SciPy as follows:

def roller_coaster(state,t,geometry,g): #unpack the state x = state[0] y = state[1] vx = state[2] vy = state[3] #evaluate the geometry f,df,df2 = geometry(x) tau = 1.0/(1.0+df**2) #return the derivative return np.array([vx,vy,-g*tau*df,-g*tau*df**2])

where the function geometry can be any curve in the plane. For example, for a bead sliding down a parabolic-shaped wire , the corresponding function is

def parabola(x): return x**2, 2*x, 2

and the assignment

geometry = parabola

provides this shape to the solver.

The first case studied with this formalism was where the wire shape is a simple straight segment, parallel to the edge of an inclided plane and refered to as such for the sake simplicity. The resulting motion was consistent with the shape

and the energy

was very well conserved.

These good results failed to hold when the shape of the wire was made parabolic. The resulting motion stayed on the parabola for a short time but eventually departed from the constraint by dropping below

while simultaneously showing a large deviation

from the initial energy.

The situation was no different for other geometeries tried. The problem seems to lie in the fact that the normal force is unable to respond to numerical errors. As the constraint begins to be violated, the wire is unable to conjure additional force to enforce it. The Newton method, while attractive from its simplicity, is simply inadequate for working with all but the most trivial of constraints.

The usual way of addressing this problem is to eliminate on degree of dynamical freedom. Either or has to go and the conventional approaches seem to favor as the parameter that stays. Again, there are several ways to affect this change but in this case both are worth examining.

In the first method, Newton's equations are used to eliminate the in favor of by combining their equations of motion by first expressing

and then by inserting the left-hand side into the equation yielding

Solving this to give

and then using

and

The final equation, after substitution for has been made and has been isolated is

In the second method, the constraint is put into the Lagrangian before forming the Euler-Lagrange equations. The resulting Lagrangian becomes

The conjugate momentum is

and the resulting Euler-Lagrange equation is

Simplifying a bit yields

which happily is the same equation as obtained from Newton's laws.

The last step is to solve this differential equation. Rewriting this equation in state-space form gives a simple system to solve numerically:

Like the original Newtonian system, this one was easily coded up into Python.

def roller_coaster_Lagrange(state,t,geometry,g): #unpack the state x = state[0] vx = state[1] #evaluate the geometry f,df,df2 = geometry(x) tau = 1.0/(1.0+df**2) return np.array([vx,-df*tau*(g+df2*vx**2)])

The constraint is always honored since is obtained after the fact from .

In addition, the conservation of energy

shows that there wasn't any error in either the theory or the numerical implementation.

However, the key question is does the solution accurately capture the motion by correctly placing the bead at at the correct time. There is no easy way to answer this question. The plot of the component as a function of time

certainly suggests that the Lagrange implementation is correctly capturing not just the dynamics but also the kinematics. That said, a more complete justification will have to wait for another time and another post.