Double Pendulum, more Iowa rides and the slug
As promised a couple of days ago, I went ahead and solved the double pendulum problem in classical mechanics. Definitely not as interesting as the Schwarzschild metric, but still pretty cool. I can never figure out how to make cool animations in matplotlib (any suggestions???) so I didn't make any. Regardless, the plot I'll present does show the crazy motion of the double pendulum (assuming I got all my algebra right!).
What's the motivation behind this problem? I remember sitting in Mechanics in my sophomore year wondering how the hell anyone could go about solving the double pendulum problem -- the equations of motion are super complicated! The method I employed to solve the two body problem and to solve for geodesics in the Schwarzschild metric still applies here. I set up the Lagrangian for the system, and apply the Euler-Lagrange equations, generating my equations of motion for the system. With these in hand, I replace my generalized coordinates and their respective time derivatives with new coordinates. This reduces the order of the system, but adds an extra equation per set of generalized coordinate. With a system of first order differential equations, I can apply any appropriate numerical method for first order systems. Forward Euler has a small stability region, so the best move, especially with unstable systems (like the double pendulum) is to use an implicit method, like backward Euler. The problem with backward Euler method is that you have to solve for the roots of a system of nonlinear equations at each time step. The best way to go about this is to use the corresponding forward Euler time step as the initial guess for the root. The SciPy optimization scipy.optimize.root method converge in a relatively small number of iterations (less than 5 or 10 iterations) when using the 'hybr' solver. Ultimately backward Euler is substantially faster (and more accurate) than forward Euler, because I can use a time step size about 3 orders of magnitude larger in backward Euler. This difference in time step (apparently) far outweighs the cost of solving the system of nonlinear equations at each time step. This is a testament to the power and speed of MINPACK. All I can say is that I'm deeply grateful that I live in 2015 as I don't have to write all these fortran routines myself.
Anyways, on to the double pendulum! I'm not going to solve it completely here, as I don't really want to type up the equations of motion. I'll set up the problem, and it'll become clear that moving from the double pendulum from the single pendulum is no small task. Imagine that we've got a single pendulum (just a weight attached to a massless rod that doesn't flex) hanging from some surface. The Lagrangian for this system is the following.
$$L = T-V \\
L = \frac{1}{2}m_{1}(l_{1}\dot{\theta_1})^2 - m_1gl_1(1-cos(\theta_1))$$
I've placed the origin at the place where the pendulum is attached to the surface, such that $x_1 = l_1sin(\theta_1)$ and $y-1 = l_1cos(\theta_1) $. I'm using the subscripts $m_1$, $\theta_1$ and $l_1$ because we're going to develop the double pendulum problem in a little bit. $g$ is the gravitation constant here on the surface of Earth, $l_1$ is the length of the rod of the pendulum, and $m_1$ is the mass of the object at the end of the pendulum. Notice that I could just as easily write the Lagrangian as follows:
$$L = \frac{1}{2}m_{1}(l_{1}\dot{\theta_1})^2 - m_1gl_1cos(\theta_1)$$
The equations of motion will be the same; when we apply the Euler-Lagrange equations
$$\frac{d}{dt}\left( \frac{\partial L}{\partial \dot{q}} \right) = \frac{\partial L}{\partial q} $$
the right side of the equation will be the same regardless of how we set up the Lagrangian. This will come in handy later. Now let's move on to the double pendulum. The coordinates of the second pendulum are as follows:
$$x_2 = l_1sin(\theta_1) + l_2sin(\theta_2) \\
y_2 = l_1cos(\theta_1) + l_2cos(\theta_2) $$
or
$$ x_2 = x_1 + l_2sin(\theta_2) \\
y_2 = y_1 + l_2cos(\theta_2) $$
Oh fudge. Now we have to differentiate with respect to time and then square these variables! We get some cancellations, but its still pretty ugly.
$$\dot{x_2} = l_1cos(\theta_1)\dot{\theta_1} + l_2cos(\theta_2)\dot{\theta_2} \\
\dot{y_2} = -l_1sin(\theta_1)\dot{\theta_1} - l_2sin(\theta_2)\dot{\theta_2} \\
\dot{x_2}^2 = (l_1cos(\theta_1)\dot{\theta_1})^2 + (l_2cos(\theta_2)\dot{\theta_2})^2 + 2l_1l_2cos(\theta_1)cos(\theta_2)\dot{\theta_1}\dot{\theta_2} \\
\dot{y_2}^2 = (l_1sin(\theta_1)\dot{\theta_1})^2 + (l_2sin(\theta_2)\dot{\theta_2})^2 + 2l_1l_2sin(\theta_1)sin(\theta_2)\dot{\theta_1}\dot{\theta_2} \\
\dot{x_2}^2 + \dot{y_2}^2 = (l_1\dot{\theta_1})^2 + (l_2\dot{\theta_2})^2 + 2l_1l_2\dot{\theta_1}\dot{\theta_2}(sin(\theta_1)sin(\theta_2) + \cos(\theta_1)\cos(\theta_2)) \\
\dot{x_2}^2 + \dot{y_2}^2 = (l_1\dot{\theta_1})^2 + (l_2\dot{\theta_2})^2 + 2l_1l_2\dot{\theta_1}\dot{\theta_2}cos(\theta_1 - \theta_2)$$
Notice that lil' trig substitution that allows us to get from the penultimate line to the last line. In the end, our Lagrangian is pretty nasty looking.
$$ L = T - V \\
L = \frac{1}{2}m_1(l_1\theta_1)^2 + \frac{1}{2}m_2((l_1\dot{\theta_1})^2 + (l_2\dot{\theta_2})^2 + 2l_1l_2\dot{\theta_1}\dot{\theta_2}cos(\theta_1 - \theta_2)) \\ - (m_1g(l_1 + l_2 - l_1cos(\theta_1)) + m_2g(l_1 + l_2 - l_1cos(\theta_1) - l_2cos(\theta_2)) $$
Just like before, I can drop all the freestanding $l_1$ and $l_2$ terms in the potential energy term.
$$L = \frac{1}{2}m_1(l_1\theta_1)^2 + \frac{1}{2}m_2((l_1\dot{\theta_1})^2 + (l_2\dot{\theta_2})^2 + 2l_1l_2\dot{\theta_1}\dot{\theta_2}cos(\theta_1 - \theta_2)) \\ + (m_1gl_1cos(\theta_1)) + m_2g( l_1cos(\theta_1) + l_2cos(\theta_2)) $$
This simplifies things a little bit. Now we have to apply the Euler-Lagrange equations to get the equations of motion. I'm not going through all the steps, but needless to say that it gets a little complicated, especially when you have to take the total derivative with respect to time of the left side of the Euler-Lagrange equation. It can be done however, and we can isolate $\ddot{\theta_1}$ and $\ddot{\theta_2}$ terms.
\begin{array}{rcl} \ddot{\theta_1} &=& \frac{4l_2^2m_2(\dot{\theta_1}sin(\Delta\theta)cos(\Delta\theta)\Delta\dot{\theta}) - 2l_1gm_2sin(\theta_2)cos(\Delta\theta) + \dot{\theta_2}sin(\Delta\theta)\Delta\dot{\theta} - gl_1sin(\theta_1)(m_1+m_2)}{l_1^2(m_1+m_2) + 4l_1^2m_2cos^2(\Delta\theta)} \\
\ddot{\theta_2} &=& \frac{2l_1}{l_2}(\ddot{\theta_1}cos(\Delta\theta) - \dot{theta_1}sin(\Delta\theta)\Delta\dot{\theta}) + \frac{g}{l_2}sin(\theta_2) \end{array}
where $\Delta\theta = \theta_1 - \theta_2$ and $\Delta\dot{\theta} = \dot{\theta_1} - \dot{\theta_2}$. Note that the equation for $\ddot{\theta_2}$ has a $\ddot{\theta_1}$ term. This is not a problem when actually solving the problem; we simply plug in the expression for $\ddot{\theta_1}$. This allows us to create a coupled system of four first order differential equations when we make the following substitutions. If you notice any errors in my algebra, please feel free to comment and let me know!
$$ \theta_a = \theta_1 ; \theta_b = \dot{\theta_1} \\
\theta_c = \theta_2 ; \theta_d = \dot{\theta_2} $$
Next, we solve with backwards Euler method. Below is a plot I made in matplotlib. The black line is the position of the first mass in time, and the red line is the position of the second mass as time progresses.
Below is a plot that shows how $\theta_1$ and $\theta_2$ change with time. The time units are arbitrary.
What's the motivation behind this problem? I remember sitting in Mechanics in my sophomore year wondering how the hell anyone could go about solving the double pendulum problem -- the equations of motion are super complicated! The method I employed to solve the two body problem and to solve for geodesics in the Schwarzschild metric still applies here. I set up the Lagrangian for the system, and apply the Euler-Lagrange equations, generating my equations of motion for the system. With these in hand, I replace my generalized coordinates and their respective time derivatives with new coordinates. This reduces the order of the system, but adds an extra equation per set of generalized coordinate. With a system of first order differential equations, I can apply any appropriate numerical method for first order systems. Forward Euler has a small stability region, so the best move, especially with unstable systems (like the double pendulum) is to use an implicit method, like backward Euler. The problem with backward Euler method is that you have to solve for the roots of a system of nonlinear equations at each time step. The best way to go about this is to use the corresponding forward Euler time step as the initial guess for the root. The SciPy optimization scipy.optimize.root method converge in a relatively small number of iterations (less than 5 or 10 iterations) when using the 'hybr' solver. Ultimately backward Euler is substantially faster (and more accurate) than forward Euler, because I can use a time step size about 3 orders of magnitude larger in backward Euler. This difference in time step (apparently) far outweighs the cost of solving the system of nonlinear equations at each time step. This is a testament to the power and speed of MINPACK. All I can say is that I'm deeply grateful that I live in 2015 as I don't have to write all these fortran routines myself.
Anyways, on to the double pendulum! I'm not going to solve it completely here, as I don't really want to type up the equations of motion. I'll set up the problem, and it'll become clear that moving from the double pendulum from the single pendulum is no small task. Imagine that we've got a single pendulum (just a weight attached to a massless rod that doesn't flex) hanging from some surface. The Lagrangian for this system is the following.
$$L = T-V \\
L = \frac{1}{2}m_{1}(l_{1}\dot{\theta_1})^2 - m_1gl_1(1-cos(\theta_1))$$
I've placed the origin at the place where the pendulum is attached to the surface, such that $x_1 = l_1sin(\theta_1)$ and $y-1 = l_1cos(\theta_1) $. I'm using the subscripts $m_1$, $\theta_1$ and $l_1$ because we're going to develop the double pendulum problem in a little bit. $g$ is the gravitation constant here on the surface of Earth, $l_1$ is the length of the rod of the pendulum, and $m_1$ is the mass of the object at the end of the pendulum. Notice that I could just as easily write the Lagrangian as follows:
$$L = \frac{1}{2}m_{1}(l_{1}\dot{\theta_1})^2 - m_1gl_1cos(\theta_1)$$
The equations of motion will be the same; when we apply the Euler-Lagrange equations
$$\frac{d}{dt}\left( \frac{\partial L}{\partial \dot{q}} \right) = \frac{\partial L}{\partial q} $$
the right side of the equation will be the same regardless of how we set up the Lagrangian. This will come in handy later. Now let's move on to the double pendulum. The coordinates of the second pendulum are as follows:
$$x_2 = l_1sin(\theta_1) + l_2sin(\theta_2) \\
y_2 = l_1cos(\theta_1) + l_2cos(\theta_2) $$
or
$$ x_2 = x_1 + l_2sin(\theta_2) \\
y_2 = y_1 + l_2cos(\theta_2) $$
Oh fudge. Now we have to differentiate with respect to time and then square these variables! We get some cancellations, but its still pretty ugly.
$$\dot{x_2} = l_1cos(\theta_1)\dot{\theta_1} + l_2cos(\theta_2)\dot{\theta_2} \\
\dot{y_2} = -l_1sin(\theta_1)\dot{\theta_1} - l_2sin(\theta_2)\dot{\theta_2} \\
\dot{x_2}^2 = (l_1cos(\theta_1)\dot{\theta_1})^2 + (l_2cos(\theta_2)\dot{\theta_2})^2 + 2l_1l_2cos(\theta_1)cos(\theta_2)\dot{\theta_1}\dot{\theta_2} \\
\dot{y_2}^2 = (l_1sin(\theta_1)\dot{\theta_1})^2 + (l_2sin(\theta_2)\dot{\theta_2})^2 + 2l_1l_2sin(\theta_1)sin(\theta_2)\dot{\theta_1}\dot{\theta_2} \\
\dot{x_2}^2 + \dot{y_2}^2 = (l_1\dot{\theta_1})^2 + (l_2\dot{\theta_2})^2 + 2l_1l_2\dot{\theta_1}\dot{\theta_2}(sin(\theta_1)sin(\theta_2) + \cos(\theta_1)\cos(\theta_2)) \\
\dot{x_2}^2 + \dot{y_2}^2 = (l_1\dot{\theta_1})^2 + (l_2\dot{\theta_2})^2 + 2l_1l_2\dot{\theta_1}\dot{\theta_2}cos(\theta_1 - \theta_2)$$
Notice that lil' trig substitution that allows us to get from the penultimate line to the last line. In the end, our Lagrangian is pretty nasty looking.
$$ L = T - V \\
L = \frac{1}{2}m_1(l_1\theta_1)^2 + \frac{1}{2}m_2((l_1\dot{\theta_1})^2 + (l_2\dot{\theta_2})^2 + 2l_1l_2\dot{\theta_1}\dot{\theta_2}cos(\theta_1 - \theta_2)) \\ - (m_1g(l_1 + l_2 - l_1cos(\theta_1)) + m_2g(l_1 + l_2 - l_1cos(\theta_1) - l_2cos(\theta_2)) $$
Just like before, I can drop all the freestanding $l_1$ and $l_2$ terms in the potential energy term.
$$L = \frac{1}{2}m_1(l_1\theta_1)^2 + \frac{1}{2}m_2((l_1\dot{\theta_1})^2 + (l_2\dot{\theta_2})^2 + 2l_1l_2\dot{\theta_1}\dot{\theta_2}cos(\theta_1 - \theta_2)) \\ + (m_1gl_1cos(\theta_1)) + m_2g( l_1cos(\theta_1) + l_2cos(\theta_2)) $$
This simplifies things a little bit. Now we have to apply the Euler-Lagrange equations to get the equations of motion. I'm not going through all the steps, but needless to say that it gets a little complicated, especially when you have to take the total derivative with respect to time of the left side of the Euler-Lagrange equation. It can be done however, and we can isolate $\ddot{\theta_1}$ and $\ddot{\theta_2}$ terms.
\begin{array}{rcl} \ddot{\theta_1} &=& \frac{4l_2^2m_2(\dot{\theta_1}sin(\Delta\theta)cos(\Delta\theta)\Delta\dot{\theta}) - 2l_1gm_2sin(\theta_2)cos(\Delta\theta) + \dot{\theta_2}sin(\Delta\theta)\Delta\dot{\theta} - gl_1sin(\theta_1)(m_1+m_2)}{l_1^2(m_1+m_2) + 4l_1^2m_2cos^2(\Delta\theta)} \\
\ddot{\theta_2} &=& \frac{2l_1}{l_2}(\ddot{\theta_1}cos(\Delta\theta) - \dot{theta_1}sin(\Delta\theta)\Delta\dot{\theta}) + \frac{g}{l_2}sin(\theta_2) \end{array}
where $\Delta\theta = \theta_1 - \theta_2$ and $\Delta\dot{\theta} = \dot{\theta_1} - \dot{\theta_2}$. Note that the equation for $\ddot{\theta_2}$ has a $\ddot{\theta_1}$ term. This is not a problem when actually solving the problem; we simply plug in the expression for $\ddot{\theta_1}$. This allows us to create a coupled system of four first order differential equations when we make the following substitutions. If you notice any errors in my algebra, please feel free to comment and let me know!
$$ \theta_a = \theta_1 ; \theta_b = \dot{\theta_1} \\
\theta_c = \theta_2 ; \theta_d = \dot{\theta_2} $$
Next, we solve with backwards Euler method. Below is a plot I made in matplotlib. The black line is the position of the first mass in time, and the red line is the position of the second mass as time progresses.
Below is a plot that shows how $\theta_1$ and $\theta_2$ change with time. The time units are arbitrary.
Pretty cool, huh? Moving away from talking about math/physics, Tom and I have been taking some good rides the last couple of days. Yesterday we took a short (~30 mile) ride, but it was brutal! We rode out into a 20 mph head wind! As I've come to kind of take for granted here in Iowa, the weather was perfect -- 80 degrees and sun. I'm getting a rather unpleasant farmer's tan. It doesn't seem to matter how much sunscreen I apply, I just keep getting tanner. I think it's rides like these that I'll remember and cherish when I get older; spending time with my parents is such a rare treat anymore. I'm excited for Tom to come out and chill with me in Nevada/Utah in a week or so.
Here's something I've been thinking about lately. In the imperial system we talk about our weight, while in the metric system we talk about our mass. We use pounds and kilograms kind of interchangeably, even though they are fundamentally different. A pound is a unit of force, and the kilogram is a unit of mass! The slug is the imperial unit of mass. A pound is the amount of force required to accelerate a slug of mass by $1 \frac{ft}{s^2}$. In imperial units my mass is about 5.6 slugs.