Simple Pendulum

The simple pendulum consists of a mass $m$, a length $L$, and angle $\theta$ measured with respect to the vertical downward direction. The diagram below details the coordinates and parameters.

It's easy to use Newton's law to calculate the force components, but it's also easy to use Lagrangians, and this will warm you up for when we have to do the double pendulum.

We take the coordinate origin $(0,0)$ to be the point where the pivot is connected to the support. Then from the diagram, the cartesian coordinates of the mass $m$ are given by:

$(x,y) = [L\sin\theta,-L\cos\theta]$

The derivatives, which will be needed later, are given by

$[\dot x,\dot y] = [L\dot\theta\cos\theta,L\dot\theta\sin\theta]$

The total kinetic energy $KE$ and gravitational potential energy $PE$ is:

$KE = \frac{1}{2} m (\dot x^2 + \dot y^2) = \frac{1}{2}mL^2\dot{\theta^2}$
$PE = mgy = -mgL\cos\theta$

The motion is totally described by the angle $\theta$. This gives us a Langrangian:

$\mathcal{L} = \frac{1}{2} mL^2\dot{\theta^2} + mgL\cos\theta$

To calculate the motion, we then solve the differential equation

$\frac{d}{dt}\frac{\partial \mathcal{L}}{\partial \dot\theta} - \frac{\partial \mathcal{L}}{\partial\theta} = 0$

to get:

(1) $\ddot\theta = -\omega_0^2\sin\theta$

where $\omega_0 \equiv g/L$ is the natural resonant frequency of oscillation.

The constants of the motion are: $\theta(0) = \theta_0$ and $\dot\theta(0) = 0$.

This equation cannot be solved analytically, so we will have to employ numerical techniques to be able to plot the motion. The techniques we will consider are the Euler method, the modified Euler-Cromer, and the Runge-Kutta RK2 (2nd order). They are pretty much based on the same concepts, the latter just being more accurate.

As described in the chapter on numerical integration, these 3 techniques are all based on the definition of derivatives: that $\frac{dx}{dt}$ = limit$(\delta t\to 0)$ of $\frac{x(t+\delta t)-x(t)}{\delta t}$. If the step size in a numerical integration is small enough, we can unravel this formula and use:

$x(t+\delta t) = x(t) + \delta t\cdot \dot x(t)$

This formula is the basis of the Euler approximation: you start at $x(t)$, and use the slope $\dot x(t)$ (slope at $t$) to find the value of $x$ at $t+\delta t$: $x(t+\delta t)$. The difference in the 3 numerical integration techniques has to do propagate the curve from $x(t)$ to $x(t+\delta t)$, and whether you actually use: the slope at the beginning of the interval, $\dot x(t)$; the slope at the end of the interval, $\dot x(t+\delta t)$; or something in between, like $\dot x(t+\frac{1}{2}\delta t)$.

In this particular example, we start with the 2nd derivative which is an explicit function not of $t$, but of $\theta$ itself:

$\ddot\theta = f(\theta) = -\omega_0^2\theta$ and of course $\theta =\theta(t)$.

So the Euler approximation would give, as a perscription for numerical integration:

$\dot\theta(t+\delta t) = \dot\theta(t) + \delta t\cdot f(\theta(t)) = \dot\theta(t) - \delta t\cdot \omega_0^2\theta$
$\theta(t+\delta t) = \theta(t) + \delta t\cdot\dot\theta(t)$

In words: we use the value of the second derivative $\ddot\theta$ evaluated at $t$ to propagate $\dot\theta$ from $t$ to $t+\delta t$, but use the value of $\dot \theta$ at the beginning of the interval ($\dot \theta(t)$) to propagate $\theta$ from $t$ to $t+\delta t$.

Another way to understand the Euler approximation is through the calculus of Taylor expansion:

$\theta(t+\delta t) = \theta(t) + \frac{d\theta}{dt}\large | \normalsize _t$

That is, we keep the first term in the Taylor expansion and evaluate the derivative at time $t$.

The Euler-Cromer evaluates the derivative at the end of the interval:

$\theta(t+\delta t) = \theta(t) + \frac{d\theta}{dt}\large | \normalsize _{t+\delta t}$

and the Runge-Kutta (RK2) technique evalutes it at the halfway point: $\theta(t+\delta t) = \theta(t) + \frac{d\theta}{dt}\large | \normalsize _{t+\frac{1}{2}\delta t}$

In the case of the simple pendulum (and remember we are not making the small angle approximation that $\sin\theta\sim\theta$), we have a second derivative that is a function of the "position", and we need to make this into a 1st order differential equation. To do this, we can write two equations:

$\omega = \dot\theta$
$\dot\omega = -\omega_0^2\sin(\theta)=f(\theta)$

and expand around some value. so instead of doing a Taylor expansion in $t$, we do a Taylor expansion

The modified (aka "implicit") Euler technique (aka Euler-Cromer) differs in that it uses the value of $\dot \theta$ at the end of the interval ($\theta (t+\delta t)$) to propagate $\theta$ from $t$ to $t+\delta t$:

$\dot\theta(t+\delta t) = \dot\theta(t) + \delta t\cdot f(\theta(t)) = \dot\theta(t) - \delta t\cdot \omega_0^2\theta$
$\theta(t+\delta t) = \theta(t) + \delta t\cdot\dot\theta(t+\delta t)$

And finally, the 2nd order Runge-Kutta (RK2) method would be:

$\dot\theta(t+\delta t) = \dot\theta(t) + \delta t\cdot f(\theta(t)+\frac{1}{2}k_1) = \dot\theta(t) -\delta t\omega_0^2\sin(\theta(t)+\frac{1}{2}k_1)$
$k_1 = \delta t\cdot f(\theta(t)) = -\delta \omega_0^2 \sin(\theta(t))$
$\theta(t+\delta t) = \theta(t) + \delta t\cdot\dot\theta(t+\delta t)$

In the simulation below, the first box below the pendulum shows the energy (somewhat arbitrary units), where the yellow line is the initial energy and the curve shows the energy using the values for $\dot\theta$ and $\theta$. Deviations between the curve and the yellow line give you an idea of the accuracy of the simulation.

The 2nd plot is the "phase space" plot of $\omega\equiv\dot\theta$ vs $\theta$, and the bottom plot is just $\theta(t)$ vs $t$. The open circle (and the oscillation path) is what we would get if we used the small angle approximation $\sin\theta\sim\theta$.

Euler Euler-Cromer Runge-Kutta
Energy (arbitrary units):