4.9: Example - the hanging pendulum

We have employed the simple harmonic oscillator as a useful example to demonstrate the different solvers encountered so far. The SHO is sometimes used as a linearized model of the motion of a hanging pendulum. However, a real hanging pendulum’s motion is governed by a nonlinear ODE. The geometry of the problem is shown in 4.14.

Figure 4.14: The real hanging pendulum.

The independent variable is time \(t\), the dependent variable (which we want to solve for) is the angle \(\theta\) made by the pendulum with respect to vertical. The pendulum’s length is \(l\) and the gravitational force constant is \(g\). The ODE is derived as follows:

1. We consider the forces acting on the pendulum when it is displaced from the vertical by an angle \(\theta\). In particular, gravity pulls the pendulum down with a constant vertical force \(m g\).
2. The vertical force may be resolved into two components: One in the \(\hat{r}\) direction – in line with the pendulum’s support, and one in the \(\hat{\theta}\) direction. Elementary trig says the force in the \(\hat{\theta}\) direction is \(F_{\theta} = -m g \sin(\theta)\). The negative sign appears since the force pulls the pendulum bob back to the equilibrium point at \(\theta = 0\).
3. Newton’s second law states \(F = m a\). In the case of the pendulum bob, \(F\) is the force acting in the \(\hat{\theta}\) direction, and \(a\) is the acceleration in the \(\hat{\theta}\) direction. Since the pendulum bob is attached to the rod of length \(l\), the angular acceleration is \(l d^2 \theta / d t^2\).
4. Using the above information, the equation governing the pendulum’s motion is

   \[\nonumber m l \frac{d^2 \theta}{d t^2} = -m g \sin(\theta)\]

   or

   \[ \frac{d^2 \theta}{d t^2} = -\frac{g}{l} \sin(\theta) \label{eq:4.3}\]

    

Observe that for small \(\theta\) we may approximate \(\sin(\theta) \approx \theta\) which gives back the SHO. But for large pendulum swings this approximation doesn’t work and we need to use the full, non-linear expression for the restoring force on the pendulum bob. In this case an analytic solution for \ref{eq:4.3} exists – the analytic solution may be expressed using the Jacobi elliptic functions sn and cn. However, rather than playing around with higher transcendental functions, we’ll just use RK4 to compute the solutions numerically.

To proceed, we break the second order ODE \ref{eq:4.3} into two first order ODEs. We get

\[\begin{aligned} \nonumber \frac{du}{dt} &= v \\ \nonumber \frac{dv}{dt} &= -\frac{g}{l} \sin(u)\end{aligned} \label{eq:NonlinearPendEqTwoODEs}\]

Now we rearrange this system into a vector form for use with the Runge-Kutta implementation. We get

\[\begin{aligned} \begin{pmatrix} du/dt \\ dv/dt \end{pmatrix} = \begin{pmatrix} v \\ -(g/l) \sin(u) \end{pmatrix} \end{aligned} \label{eq:4.4}\]

This form is ready for implementation. The vectorized form \ref{eq:4.4} is implemented in a function called by the RK4 integrator itself, and a test wrapper is the main program run by the user. This is the usual software architecture previously shown in 2.2.

Running the RK4 implementation gives the plot shown in 4.15. This plot was obtained using \(g = 9.8\) m/sec\(^2\) (the value of \(g\) at the earth’s surface) and \(l = 0.5\)m, a reasonable value for the length of a pendulum. For initial conditions, the pendulum swing started with the bob pointing almost vertically up. Therefore, the pendulum takes a long time for the bob to accelerate. It then swings quickly around the downward hanging position and continues until pointing almost vertically up again. The pendulum then repeats this motion again and again, never stopping because the original equation does not include any dissipation (friction).

Figure 4.15: The hanging pendulum ODE solved using RK4. This plot shows the motion of the pendulum when the bob starts pointing almost vertically up.