Simple Pendulum#

Let consider a simple pendulum consists of a mass \(m\) attached to the end of a string of length \(L\)

simple-pendulum.jpg

Solve ODE#

The motion of a simple pendulum can be described by a second-order differential equation, but, we need to convert this into a system of first-order equations to solve it numerically.

The second-order differential equation for the angle \(\theta\) of the pendulum (with respect to the vertical) is given by $:nbsphinx-math:ddot{theta} + \frac{g}{L} \sin`(:nbsphinx-math:theta`) = 0 $

Here \(\ddot{\theta}\) represents the second derivative of theta with respect to time (angular acceleration). Where where \(g\) is the acceleration due to gravity, and \(L\) is the length of the pendulum. This equation arises from the balance of gravitational and rotational forces acting on the pendulum.

To convert this into a system of first-order differential equations, we introduce \(\omega\) as the angular velocity of the pendulum (\(\omega=\dot{\theta}\)).

The system of equations then becomes:

\(\frac{d\theta}{dt} = \omega\) and $:nbsphinx-math:frac{domega}{dt} = -\frac{g}{L} \sin`(:nbsphinx-math:theta`) $

We can now define a function to represent this system and then use odeint to solve it. Here’s how you would define the function in Python:

Import Modules#

[29]:
import numpy as np
from scipy.integrate import odeint
import plotly.graph_objects as go

Note

The functions create_line_trace, create_point_trace, create_arrow_trace, and and others were written in previous tutorials. Please include them in your notebook on top before starting to follow this tutorial. You can download it by clicking the Download icon on the Navigation Bar.

Define System Function#

[31]:
def simple_pendulum(y, t, g, L):
    theta, omega = y
    dtheta_dt = omega
    domega_dt = -(g / L) * np.sin(theta)
    return [dtheta_dt, domega_dt]

Constants and Conditions#

[32]:
# Constants
g = 9.81  # Acceleration due to gravity (m/s^2)
L = 1.0   # Length of the pendulum (m)
m = 5 # mass

# Initial conditions
theta0 = np.pi / 4  # Initial angle (45 degrees)
omega0 = 0.0        # Initial angular velocity
y0 = [theta0, omega0]

# Time points
t = np.linspace(0, 10, 500)

Solve ODE#

[33]:
# Solve ODE
solution = odeint(simple_pendulum, y0, t, args=(g, L))

# Extract the solutions
theta_vals, omega_vals = solution.T

# Compute angular acceleration
alpha_vals = -(g / L) * np.sin(theta_vals)

Plot the Solutions#

[34]:
# Create traces
trace1 = go.Scatter(x=t, y=theta_vals, mode='lines', name='Angular Displacement')
trace2 = go.Scatter(x=t, y=omega_vals, mode='lines', name='Angular Velocity')
trace3 = go.Scatter(x=t, y=alpha_vals, mode='lines', name='Angular Acceleration')

# Create a figure
fig = go.Figure()

# Add traces to the figure
fig.add_trace(trace1)
fig.add_trace(trace2)
fig.add_trace(trace3)

# Update layout for a better visualization
fig.update_layout(title='Simple pendulum Dynamics',
                  xaxis_title='Time (s)',
                  yaxis_title='Theta / Omega',
                  legend_title='Variable')

fig.show()

Animation#

[35]:
# Convert angular displacement to Cartesian coordinates for animation
x_vals = L * np.sin(theta_vals)
y_vals = -L * np.cos(theta_vals)

# Create the animation
curve_points = np.column_stack((x_vals, y_vals, np.zeros_like(x_vals)))
fig = create_particle_animation(curve_points, title="Simple Pendulum Animation")

# Show the animation
fig.show()

Energy#

Kinetic Energy is given by $E_k = \frac{1}{2}`mL\ :sup:`2:nbsphinx-math:omega``2 $ and Potential Energy is given by $E_p = mgh $ where \(h = L - L\cos(\theta)\)

[36]:
def kinetic_energy(omega, L):
    return 0.5 * m * L**2 * omega**2

def potential_energy(theta, g, L):
    return m * g * L * (1 - np.cos(theta))

# Calculate energies
KE = kinetic_energy(omega_vals, L)
PE = potential_energy(theta_vals, g, L)
TE = KE + PE  # Total Energy

# Plotting energies
energy_fig = go.Figure()
energy_fig.add_trace(go.Scatter(x=t, y=KE, mode='lines', name='Kinetic Energy'))
energy_fig.add_trace(go.Scatter(x=t, y=PE, mode='lines', name='Potential Energy'))
energy_fig.add_trace(go.Scatter(x=t, y=TE, mode='lines', name='Total Energy'))

energy_fig.update_layout(title='Energy of the Mass-Spring System',
                         xaxis_title='Time (s)',
                         yaxis_title='Energy',
                         legend_title='Type of Energy')

energy_fig.show()

Linear Momentum#

[37]:
P = m * L * omega_vals
fig = go.Figure()
fig.add_trace(go.Scatter(x=t, y=P, mode='lines',name='Linear Momentum'))
fig.update_layout(title='Linear Momentum Vs Time')
fig.show()