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 y = state vx = state vy = state #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
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 vx = state #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.