The ideal pendulum
Numerical solutions to motion of an ideal pendulum can generate useful dynamical system datasets for testing a bunch of proof-of-concept research ideas. Below, I show solutions to an undamped and damped ideal pendulum coded in Python.
import numpy as np
from scipy.integrate import solve_ivp
import matplotlib.pyplot as plt
The state vector for the ideal pendulum is composed of the angle, \(\theta(t)\), that the pendulum makes with the vertical direction, and the angular velocity, \(\omega(t)\), both being functions of time, \(t\). Hence, in the code you will see:
θ, ω = state_vector
Undamped pendulum
The system of ordinary differential equations (ODEs) that describes motion of an undamped ideal pendulum is as following:
\(\begin{equation}\begin{cases} \frac{d \theta(t)}{dt} = \omega(t) \\ \frac{d \omega(t)}{dt} = - \frac{g}{L} \sin(\theta(t)) \end{cases}\end{equation}\)
where \(g\) is the gravitational acceleration and \(L\) is the length of the pendulum.
We first define a function that generates the right-hand-side (RHS) vector for this system of ODEs:
def undamped_pendulum(t,
state_vector,
g=9.8,
L=1.0):
θ, ω = state_vector
return [ω, - g/L * np.sin(θ)]
Next, we create a generic function that will solve an initial value problem (IVP) given the time vector, the RHS vector, the initial condition, \(\theta(t=0) = \theta_0\) and \(\omega(t=0) = \omega_0\), and any other parameters:
def generate_pendulum_trajectory(θ_0,
ω_0,
T=10.0,
n_points=200,
b=0.1,
g=9.8,
L=1.0,
damped=False,
method='BDF'):
time_vector = np.linspace(0, T, n_points)
if damped:
sol = solve_ivp(damped_pendulum,
[0, T],
[θ_0, ω_0],
method=method,
t_eval=time_vector,
args=(b,),
rtol=1e-9,
atol=1e-9)
else:
sol = solve_ivp(undamped_pendulum,
[0, T],
[θ_0, ω_0],
method=method,
t_eval=time_vector,
args=(g,L,),
rtol=1e-9,
atol=1e-9)
return sol.t, sol.y.T
Now you can simply solve the undamped pendulum for a specific choice of parameters:
time, trajectory = generate_pendulum_trajectory(θ_0=0.1,
ω_0=0.0,
T=20.0,
n_points=500,
b=0.1,
g=9.8,
L=9.8,
damped=False,
method='BDF')
To test my code, I create the figures below which are my reproductions of figures shown in section 1.2.3 in Prof. Scheinerman’s textbook Invitation to Dynamical Systems. The parameters’ values used in the textbook are \(g = 9.8 \frac{m}{s^2}\) and \(L = 9.8 m\) for simplicity.
The motion of an undamped pendulum with \(\theta_0 = 0.1\) and \(\omega_0 = 0.0\).
The motion of an undamped pendulum with \(\theta_0 = 3.0\) and \(\omega_0 = 0.0\).
The motion of an undamped pendulum with \(\theta_0 = 0.0\) and \(\omega_0 = 2.0\). Interestingly, here I wasn’t quite able to get Prof. Scheinerman’s figure and the solution strongly depends on the type of solver used (I’ve used BDF).
Damped pendulum
The system of ODEs that describes motion of a damped ideal pendulum is as following:
\(\begin{equation}\begin{cases} \frac{d \theta(t)}{dt} = \omega(t) \\ \frac{d \omega(t)}{dt} = - \frac{g}{L} \sin(\theta(t)) - b \cdot \omega(t) \end{cases}\end{equation}\)
where \(b\) is the damping parameter.
We define a function that generates the RHS vector for this system of ODEs:
def damped_pendulum(t,
state_vector,
b=0.1,
g=9.8,
L=1.0):
θ, ω = state_vector
return [ω, - g/L * np.sin(θ) - b * ω]
And we solve it with the same function generate_pendulum_trajectory() as before:
time, trajectory = generate_pendulum_trajectory(θ_0=3.0,
ω_0=0.0,
T=30.0,
n_points=500,
b=0.1,
g=9.8,
L=9.8,
damped=True,
method='BDF')
Here’s one example of the trajectory of the damped pendulum, where you can observe the angle of the pendulum decaying with time:
The motion of a damped pendulum with \(\theta_0 = 3.0\) and \(\omega_0 = 0.0\).