In this note, we visualize the dynamics of the the following ODE system:

\[\begin{cases} \dot{\theta} = \Omega \\ \dot{\Omega} = -\sin(\theta) \end{cases}\]

where \(\theta\in [-\pi,\pi]\) is the angle of displacement, \(\Omega = \frac{d\theta}{dt}\in\mathbb{R}\) is the angular velocity. We solve the system using different choices of initial conditions \((\theta(0),\Omega(0))^T = (\theta_0, \Omega_0)\).

Without damping, the mechanical energy should be conserved as:

\(E(t) = \frac12 \Omega(t)^2 - \cos(\theta(t)) = E_0, \forall t\in [0,T_m)\) where $T_m$ is the final time of observation.

Below we numerically solve the dynamical system from $t=0$ to $t=t_m$. We make use of scipy’s built-in ODE integrator which by default uses fourth-order explicit Runge-Kutta time-stepping. We briefly describe the numerical scheme, consider the general initial value problem:

\[\frac{d\mathbf{x}}{dt} = \mathbf{f}(t,\mathbf{x})\]

where \(\mathbf{x}\in\mathbb{R}^d\) is a $d$-dimensional vector and \(\mathbf{f}:\mathbb{R}^+ \times \mathbb{R}^d \rightarrow \mathbb{R}^d\) is a time-dependent vector field. In order to obtain a numerical solution, we need to rely on a discretized time grid \(T^{(0)}, T^{(1)},\cdots, T^{(N_t)}\). The Runge-Kutta method can be roughly considered as a weighted average of local slopes (taken at a number of well-chosen points on the grid) to approximate \(\frac{d\mathbf{x}}{dt}\). Assuming a uniform grid, \(\Delta t = T^{(n+1)}-T^{(n)}\),\(\mathbf{x}_n := \mathbf{x}(T^{(n)})\), the RK4 scheme for time step \(T^{(n)}\mapsto T^{(n+1)}\) is computed as follows:

\[k_1 = \mathbf{f}(T^{(n)},\mathbf{x}_n)\] \[k_2 = \mathbf{f}(T^{(n)}+\frac{\Delta t}{2}, \mathbf{x}_n + \Delta t \frac{k_1}{2})\] \[k_3 = \mathbf{f}(T^{(n)}+\frac{\Delta t}{2}, \mathbf{x}_n + \Delta t \frac{k_2}{2})\] \[k_4 = \mathbf{f}(T^{(n)}+\Delta t, \mathbf{x}_n + \Delta t {k_3})\] \[\mathbf{x}_{n+1} = \mathbf{x}_n + \frac{1}{6}\Delta t(k_1 + 2k_2 + 2k_3 + k_4)\]

The period \(P\) of pendulum oscillations can be found from our dynamical system. \(P\in \mathbb{R}^+\) is such that \(\theta(t+P) = \theta(t)\).

\(E(t) = E_0 = \frac12\Omega(t)^2 - \cos(\theta(t))\) solving for \(\Omega(t)\) yields:

\[\Omega(t) = \pm \sqrt{2E_0 + \cos(\theta)}\]

Consider the points where velocity is 0, substitute in \(\Omega = 0\): \(E_0 = -\cos(\theta)\)

solving for \(\theta\) we obtain:

\(\theta_m = \cos^{-1}(-E_0)\) $\theta_m$ would correspond to the highest point the pendulum can reach after releasing. Given \(E_0\), without loss of generality we consider a quarter-period in phase space, by symmetry, effectively \(\theta\in [0,\theta_m]\). Then:

\[\Omega = \frac{d\theta}{dt} = \sqrt{2E_0 + \cos(\theta)}\]

then:

\[\frac{dt}{d\theta} = \frac{1}{\sqrt{2E_0 + \cos(\theta)}}\] \[dt = \frac{1}{\sqrt{2E_0 + \cos(\theta)}}d\theta\] \[\frac{P}{4} = \int_0^{\theta_m}\frac{1}{\sqrt{2E_0 + \cos(\theta)}}d\theta\]

Alternatively, we may also solve the dynamical system numerically, and obatin the period \(P\) by observing the solution \(\theta(t)\). In our numerical evaluation below, we consider $P$ as a mapping of \(E_0\):

\[E_0 \mapsto \int_0^{\theta_m}\frac{1}{\sqrt{2E_0 + \cos(\theta)}}d\theta\]
  • Download Jupyter notebook here.