MyPhysicsLab – Dangling Stick

This simulation shows a dangling stick which is a massless rigid stick with a point mass on each end. One end of the stick is attached to a spring, and gravity acts.
Click "show controls" and then you can change parameters in the simulation such as mass, gravity, and damping. You can drag either end of the stick with your mouse to change the starting position. The "reset" button puts the simulation in a motionless equilibrium. If you don't see the simulation try instructions for enabling Java. Scroll down to see the math!

Can you find the bug in this simulation? Play with it for a while and see if you can determine when the bug occurs. Perhaps you can even have a theory for why it occurs. The answer is found below.

Physics

dangling stick variables

The equations of motion are derived by the Lagrangian method. The three independent variables are:

r = length of spring
θ = angle of spring (0 = vertical down)
φ = angle of stick (0 = vertical down)

We use subscript 1 for the mass at the spring-stick intersection, and subscript 2 for the mass at the free end of the stick. The cartesian (x, y) position of the point masses is represented by R₁, R₂. Note that we use bold and overline to indicate a vector quantity.

Other constants are:

m₁ = mass amount at spring-stick intersection
m₂ = mass amount at free end of stick
S = length of stick
k = spring stiffness
h = spring rest length
g = gravity constant

Kinematics of Dangling Stick

First we calculate the kinematics. That is, we find expressions for the positions R₁, R₂ of the two masses in terms of the independent variables r, θ, φ. Then we differentiate to get the velocities. Note that kinematics is purely an exercise in geometry, there is no information about the forces involved at all. R₁ = {r sin θ, −r cos θ} R₂ = {r sin θ + S sin φ, −r cos θ − S cos φ)} The two components of the positions correspond to {x, y} respectively. Next differentiate to get the velocities of the masses, v₁, v₂ R₁' = v₁ = {r θ' cos θ + r' sin θ, −r' cos θ + r θ' sin θ } R₂' = v₂ = {r θ' cos θ + r' sin θ + S φ' cos φ, −r' cos θ + r θ' sin θ + S φ' sin φ}

Energy and the Lagrangian

Next we get separate expressions for the kinetic energy T and potential energy V of the system. The difference of the two is the Lagrangian L = T − V.

Kinetic energy is given by ¹⁄₂ m v². Note that we use the vector dot product to square the velocity vector, v² = v · v. So we have T = ¹⁄₂ m₁ v₁ · v₁ + ¹⁄₂ m₂ v₂ · v₂ The potential energy is from the spring energy and the gravitational potential of the two point masses. V = ^k⁄₂ (r − h)² − m₂ g (r cos θ + S cos φ) − m₁ g r cos θ

Equations of Motion

To find the equations of motion, we evaluate the Lagrangian equation once for each of the three independent variables. 0 = $\frac{d}{dt} \left( \frac{\partial L}{\partial r'} \right) - \frac{\partial L}{\partial r}$ 0 = $\frac{d}{dt} \left( \frac{\partial L}{\partial \theta'} \right) - \frac{\partial L}{\partial \theta}$ 0 = $\frac{d}{dt} \left( \frac{\partial L}{\partial \phi'} \right) - \frac{\partial L}{\partial \phi}$ The three Lagrangian equations can now be solved to find expressions for the second derivatives r'', θ'', φ''. After using a computer algebra program such as Mathematica we find the following equations of motion.

r'' =	1	( k (2 m₁ + m₂) (h − r) + 2 m₁ (m₁ + m₂) r θ'² + 2 g m₁ (m₁ + m₂) cos θ + 2 S m₁ m₂ φ' ² cos (θ − φ) − k m₂ (h − r) cos(2θ − 2φ) )
	2 m₁ (m₁ + m₂)

θ'' =	1	( k m₂ (h − r) sin(2θ − 2φ) − 2 g m₁ (m₁ + m₂) sin θ − 4 m₁ (m₁ + m₂) r' θ' + 2 S m₁ m₂ φ' ² sin(θ − φ) )
	2 m₁ (m₁ + m₂) r

φ'' =	k (h − r) sin(θ − φ)
	S m₁

This is almost ready for the Runge-Kutta algorithm to numerically solve the equations. All that is needed is to change these 3 second order equations to 6 first order equations by defining separate variables for r', θ', φ'. Let the new variables be respectively v_r, v_θ, v_φ. Then the left-hand side of the above equations of motion become v_r', v_θ', v_φ' and the three additional equations are simply r' = v_r
θ' = v_θ
φ' = v_φ

A Bug

There is a problem in the equations of motion when the length of the spring r goes to zero, because then one of the second derivatives (the one for θ'' ) goes to infinity. If you play with the simulation long enough you will likely see this happen.

Here is my idea for why it blows up: The rate that the angle of the spring changes depends on the torque (twist) applied, which is partially dependent on the length of the spring. When the length goes to zero, the torque goes to infinity.

In the real world, springs never reach zero length. To eliminate this bug, we would need to model the spring with a non-linear force. As the spring compressed to near zero (or to near its smallest possible length) the force would go to infinity.