Nobody really cares about pendulums, but there are times when knowing a little about them might come in handy; If you doubt it, ask The Pit and the Pendulum's main character how he feels about the matter.
If we lived in the idealized world of college physics courses, life would be far easier for the poor man in the picture. The elder would have swung the pendulum from an angle no greater than 10 degrees from vertical, making it behave like an undamped harmonic oscillator. Undamped’ implies, among other things, the absence of air resistance.
The ideal scene of the previous paragraph can be mathematically modeled with the following differential equation. θ is the angle made by the pendulum from vertical at any given moment.
\[\frac{d^2 \theta}{dt^2} + \frac{g}{l} \sin{\theta} = 0\]
Eq 1.
If θ is always very small (θ < 10°), we could say that the sine of the angle is the angle itself (i.e. sin θ ≈ θ), this simplifies the equation to:
\[\frac{d^2 \theta}{dt^2} + \frac{g}{l} \theta = 0\]
Eq 2.
Equation 2 describes simple harmonic motion, with an analytical solution, but it only holds true for small angles (boring). We are interested in solving equation 1, and we're gonna use Python -equipped with some hardcore guns- to do it. In order to use Python as a mathematical tool you'll need to install and download Python as well as the SciPy, NumPy, and Matplotlib libraries.
SciPy is a library that holds lots of mathematical functions, classes, numerical methods, etc. Matplotlib is a library that graphs and plots stuff. When you install them, you're ready to go. If you're using Windows, go to file->all programs->Python 2.x->IDLE. That will open a Python code editor; start by copying this into it:
import scipy
import scipy.integrate
import numpy as np
import matplotlib as mpl
import matplotlib.pyplot as plt
That imports all the libraries and packets we need. scipy.integrate packet has a set of functions for numerical integration and differential equation solving. The class scipy.integrate.solve_ivp is similar to Matlab’s ode45, which is used for solving differential equations. scipy.integrate.solve_ivp only understands first order differential equations and equation 1, my friends, is a second order one. We must then express equation 1 as a system of first order differential equations:
\[ \begin{bmatrix} \text{theta} \\ \text{dthetadt} \end{bmatrix} = \begin{bmatrix} \theta \\ \dot{\theta} \end{bmatrix} \rightarrow \mathbf{dy2dt} \rightarrow \text{return: } \begin{bmatrix} \dot{\theta} \\ \ddot{\theta} \end{bmatrix} = \begin{bmatrix} \text{dthetadt} \\ \text{d2thetadt} \end{bmatrix} = \begin{bmatrix} \text{dthetadt} \\ -\frac{g}{l} \sin{\theta} \end{bmatrix} \]
As you can see, the system receives θ and dθ/dt, and outputs dθ/dt and d2θ/dt2. Taking that into account, we define this function:
# Define the pendulum model as a system of first-order differential equations
def dy2dt(t, y, L):
theta, dthetadt = y
d2thetadt2 = -(g/L) * np.sin(theta)
return [dthetadt, d2thetadt2]
Following the f(t,y,args * ) pattern, it shall be inferred by the reader that L (the length of the pendulum) is passed as an additional argument to the function. Having done that, the equation is ready to be solved with scipy.integrate.solve_ivp. Other parameters that this method will receive are the time range we want it to integrate through, the initial conditions as well as the integration method, in this case we'll be utilizing the Runge-Kutta 4,5 method.
# Constants
g = 9.8 # Gravity (m/s^2)
L = 1.0 # Length of the pendulum (1 meter)
# List of initial angles in degrees to drop the pendulum from
angles = [5, 10, 15, 30, 45, 60, 90]
# Custom colors for plotting each angle
colors = ["#f04859", "#ec5b5b", "#e6775f", "#df9763", "#dd9c63", "#d9b366", "#d4c669"]
# Create a new Matplotlib figure
fig = plt.figure()
# Iterate through each angle and its corresponding color
for j, c in zip(angles, colors):
# Define initial conditions: angle in radians, initial angular velocity
theta0 = np.radians(j) # Convert angle to radians
dthetadt0 = 0.0 # Initial velocity (0 degrees/s)
# Solve the system of ODEs using the Runge-Kutta 4,5 method (LSODA is default in latest SciPy)
sol = scipy.integrate.solve_ivp(
fun=lambda t, y: dy2dt(t, y, L),
t_span=(0, 10), # Time range (0 to 10 seconds)
y0=[theta0, dthetadt0], # Initial conditions
method='RK45', # Use Runge-Kutta 4,5 method
max_step=0.01 # Time step
)
# Extract results: time points and corresponding theta values (converted to degrees)
t = sol.t
thetas = np.degrees(sol.y[0]) # Convert theta from radians to degrees
# Plot theta(t) for the given initial angle
plt.plot(t, thetas, color=c)
# Add legend and labels to the graph
plt.legend([f'{angle} deg' for angle in angles], loc="upper right")
plt.xlabel('Time (s)')
plt.ylabel('Theta(t) (degrees)')
plt.title('Pendulum Motion for Various Initial Angles')
plt.show()The initial conditions of the problem are the angle θo (theta0) and ωo angular velocity (dthetadt0), which are the initial angle from which the pendulum is released and the initial velocity of the pendulum respectively. If it is just dropped without any impulse, then ωo = 0 [degrees/s].
We want to know the effect of dropping the same pendulum at different angles. That's why we start by iterating over different values for θo, starting from 5 to 90 degrees.
solve_ivp solves equation 1 for the time range from 0 to 10 seconds. The solution for the equestion are saved in sol.y (where y[0] is an array with all values of θ for the time range [0 to 10 second], while y[1] is an array containing the angular velocities of the pendulum within the same time range).
Finally, the different solutions are plotted using the Matplotlib library.
* You can download the code that renders this image from the downloads section.
As you can see, the period of the pendulum increases with the initial angle. Looking closer, you'll notice that for initial angles between 5 degrees and 15 degrees the period remains constant and the pendulum behaves like a simple harmonic oscillator with a period of 2 seconds approximately (for a 1 meter pendulum, based on equation 4), but when the angle is greater than 15 degrees this approximation doesn't work.
\[ T \approx 2\pi \sqrt{\frac{l}{g}} \]
Eq 4. Pendulum period approximation.
Equation 4 works only for small amplitudes. So, which is the general formula for the period of a pendulum regardless of its oscillating amplitude? This is a formula which is easily inferred applying the law of conservation of mechanical energy to a pendulum.
\[ T = \frac{4}{\sqrt{2}} \sqrt{\frac{l}{g}} \int_0^{\theta_o} \frac{d\theta}{\sqrt{\cos{\theta} - \cos{\theta_o}}} \]
Eq 5.
The integral from the right side of the equation has not a direct solution and, in fact, there are several ways to compute it. The scipy.integrate packet has the function .quad(f,a,b), which numerically integrates the function f over an interval [a,b]. f, however, must be provided following a different syntax than the used to pass a function to the ode class, as follows:
# Define function f and calculate its integral over the [0, theta0] interval f = lambda x: 1 / np.sqrt(np.cos(x) - np.cos(theta0)) F, erri = integrate.quad(f, 0, theta0) # Integrate numerically # Calculate the period (T) using the integral and other constants T = 4 * np.sqrt(l / g) * F
Using a similar procedure to that used to plot the solution of equation 1, we can also graph the period as a function of the initial angle θ₀:
The pink line shows how the period of a 1 meter pendulum changes as the initial angle θo increases. As we already know, when θo is small, the period of the pendulum remains constant and can be calculated with equation 4. When θo is greater than 10 degrees (approximately), the period starts rising significantly, following the model from equation 5.
The green line is the plot of the relation between equation 5 and equation 4. A way to interpret this behavior is considering equation 5 and equation 4 equal for small values of θo because
[Eq 5. / Eq 4] ≈ 1 for every θo < 10.
That's it! In a future post I'll talk about a pendulum damped model where air resistance is considered.
You can find the full code in the following Github repo: https://github.com/thecodebeatz/pendulum-simulation