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 springstick 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_{1},
r_{2}. note that we use bold and overline to indicate a vector quantity.
other constants are:
 m_{1} = mass amount at springstick intersection
 m_{2} = 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_{1}, r_{2} 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_{1} = {r sin θ, −r cos θ}
r_{2} = {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_{1}, v_{2}
r_{1}' = v_{1} = {r θ' cos θ + r' sin θ,
−r' cos θ + r θ' sin θ }
r_{2}' = v_{2} =
{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
^{1}⁄_{2} m v^{2}. note that we use the vector dot product to square the velocity vector,
v^{2} = v · v. so we have
t = ^{1}⁄_{2} m_{1} v_{1} · v_{1} + ^{1}⁄_{2} m_{2} v_{2} · v_{2}
the potential energy is from the spring energy and the gravitational potential of the
two point masses.
v = ^{k}⁄_{2} (r − h)^{2} −
m_{2} g (r cos θ + s cos φ) −
m_{1} 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 =
0 =
0 =
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_{1} + m_{2}) (h − r) +
2 m_{1} (m_{1} + m_{2}) r θ'^{2} +
2 g m_{1} (m_{1} + m_{2}) cos θ +
2 s m_{1} m_{2} φ' ^{2} cos (θ − φ) −
k m_{2} (h − r) cos(2θ − 2φ)
)


2 m_{1} (m_{1} + m_{2})

θ'' =

1

(
k m_{2} (h − r) sin(2θ − 2φ) −
2 g m_{1} (m_{1} + m_{2}) sin θ −
4 m_{1} (m_{1} + m_{2}) r' θ' +
2 s m_{1} m_{2} φ' ^{2} sin(θ − φ)
)


2 m_{1} (m_{1} + m_{2}) r

φ'' =

k (h − r) sin(θ − φ)

s m_{1}

this is almost ready for the rungekutta 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 lefthand 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 nonlinear force. as the spring compressed to near zero (or to near its smallest possible length) the force would go to infinity.