Module 3 - Project: pendulum solutions

Hide code cell content
import numpy as np
import matplotlib.pyplot as plt
from scipy.optimize import fsolve
plt.style.use('fivethirtyeight')

Module 3 - Project: pendulum solutions#

In this notebook, you will compare two solutions for the motion of a swinging pendulum

  1. linearized solution for \(\ddot{\theta}(t) = -\frac{g}{L}\theta\)

  2. numerical solution for \(\ddot{\theta}(t) = -\frac{g}{L}\sin\theta\)

You can choose your own pendulum length \(L\) for the analysis. Try tying a string to an object and recording the swinging of the pendulum. You can estimate the initial angle and its period of oscillation.

1. Linear solution#

The linear solution is the first order Taylor series approximation for the equation of motion

\(\ddot{\theta}(t) = -\frac{g}{L}\sin\theta \approx -\frac{g}{L}\theta\)

The solution to this simple harmonic oscillator function is

\(\theta(t) = \theta_0 \cos\omega t + \frac{\dot{\theta}_0}{\omega}\sin\omega t\)

where \(\theta_0\) is the initial angle, \(\dot{\theta}_0\) is the initial angular velocity [in rad/s], and \(\omega=\sqrt{\frac{g}{L}}\).

Here, you can plot the solution for

  • \(L = 1~m\)

  • \(g = 9.81~m/s^2\)

  • \(\theta_0 = \frac{\pi}{3} = [60^o]\)

  • \(\dot{\theta}_0 = 0~rad/s\)

g = 9.81 # m/s/s
L = 1 # m
w = np.sqrt(g/L) # rad/s
t = np.linspace(0, 4*np.pi/w) # 2 time periods of motion
theta0 = np.pi/3
dtheta0 = 0

theta = theta0*np.cos(w*t) + dtheta0/w*np.sin(w*t)

plt.plot(t, theta*180/np.pi)
plt.xlabel('time (s)')
plt.ylabel(r'$\theta(t)$ [degrees]')
Text(0, 0.5, '$\\theta(t)$ [degrees]')
../_images/05df57437f10d3cf321303827a21e153985c3e1b0b1306a2cbee6b901f05404f.png

In the plot above, you should see harmonic motion. If it starts from rest at \(60^o\) the pendulum will swing to \(-60^o\) in \(\frac{T}{2} =\frac{\pi}{\omega}\) seconds.

2. Numerical solution#

In the numerical solution, there is no need to linearize the equation of motion. You do need to create a state variable as such

\(\mathbf{x} = [x_1,~x_2] = [\theta,~\dot{\theta}]\)

now, you can describe 2 first order equations in a function as such

  1. \(\dot{x}_1 = x_2\)

  2. \(\dot{x}_2 = \ddot{\theta} =-\frac{g}{L}\sin\theta\)

putting this in a Python function it becomes

g = 9.81
L = 1
def pendulum(t, x):
    '''pendulum equations of motion for theta and dtheta/dt
    arguments
    ---------
    t: current time
    x: current state variable [theta, dtheta/dt]
    outputs
    -------
    dx: current derivative of state variable [dtheta/dt, ddtheta/ddt]'''
    
    dx = np.zeros(len(x))
    dx[0] = x[1]
    dx[1] = -g/L*np.sin(x[0])
    return dx

This function pendulum defines the differential equation. Next, integrate with solve_ivp to find \(\theta(t)\) as such.

Note: Here, I am using some previously defined values:

  • theta0: the initial angle

  • dtheta0: the initial angular velocity

  • t: the timesteps from the previous linear analysis \(t = (0...2T)\)

Hide code cell source
from scipy.integrate import solve_ivp

sol = solve_ivp(pendulum, [0, t.max()], [theta0, 0], t_eval = t)
plt.plot(sol.t, sol.y[0]*180/np.pi)
plt.xlabel('time (s)')
plt.ylabel(r'$\theta(t)$ [degrees]')
Text(0, 0.5, '$\\theta(t)$ [degrees]')
../_images/c4f9cbf9dcd8f7f4b7a56bcafb80e159cf420e9bb14c181f26f7f3e2d31645a3.png

Comparing results#

At first glance, it would seem both solutions are close enough, but here you can plot them together and compare.

plt.plot(t, theta*180/np.pi)
plt.plot(sol.t, sol.y[0]*180/np.pi)
plt.xlabel('time (s)')
plt.ylabel(r'$\theta(t)$ [degrees]')
Text(0, 0.5, '$\\theta(t)$ [degrees]')
../_images/bbca107a1e01c732ac23b7b545d1f0a39cda665379ffe67045b82540efd53dae.png

You should notice a difference in time period between the two solutions, which means the natural frequencies are different.

You can also consider angular velocity and angular acceleration, below a set of 3 plots are shown for \(\theta(t),~\dot{\theta}(t),~and~\ddot{\theta}(t)\).

dtheta = -theta0*w*np.sin(w*t) + dtheta0*np.cos(w*t)
ddtheta = -theta0*w**2*np.cos(w*t) - dtheta0*w*np.sin(w*t)

ddtheta_numerical = -g/L*np.sin(sol.y[0]) # using equation of motion
plt.figure(figsize= (10, 11))
plt.subplot(311)
plt.plot(t, theta*180/np.pi)
plt.plot(sol.t, sol.y[0]*180/np.pi)

plt.ylabel(r'$\theta(t)$ [degrees]')

plt.subplot(312)
plt.plot(sol.t, sol.y[1])
plt.plot(t, dtheta)

plt.ylabel(r'$\dot{\theta}(t)$ (rad/s)')


plt.subplot(313)
plt.plot(t, ddtheta)
plt.plot(t, ddtheta_numerical)
plt.xlabel('time (s)')
plt.ylabel(r'$\ddot{\theta}(t)$ (rad/s$^2$)')
Text(0, 0.5, '$\\ddot{\\theta}(t)$ (rad/s$^2$)')
../_images/106c9bf1a1c3f14b0ad4dc9d18516d9bfe8f54d9433cb4b8cf28bfa223e7d078.png

Question prompt#

  • Where do you see the biggest difference between the linearized solution and the numerical solution?

  • What are ways you can quantify the error between then two results? Is one solution more accurate?

  • What other calculations can you compare between these two solutions?