Week 3 - Solving second order ODEs

The motion of a simple pendulum is described by a second-order differential equation derived from Newton’s laws or Lagrangian mechanics. The angular position θ(t)\theta(t)θ(t) and angular velocity ω(t)\omega(t)ω(t) can be defined using the following function:

b is the damping coefficient (viscous friction),

m is the mass of the pendulum bob,

L is the length of the pendulum,

g is the acceleration due to gravity,

ω=dθdt\omega = \frac{d\theta}{dt}ω=dtdθ​ is the angular velocity.

L=1 metre,

m=1 kg,

b=0.05.

g=9.81 m/s2.

Simulate the motion between 0-20 sec, for angular displacement=0,angular velocity=3 rad/sec at time t=0.

When you do this you will get the position of the pendulum as a function of time.

This function numerically solves the damped pendulum equation using the Runge-Kutta method:

python

import numpy as np

import matplotlib.pyplot as plt

from scipy.integrate import solve_ivp

def damped_pendulum_equations(t, y, g, L, b, m):

theta, omega = y # y[0] = theta, y[1] = omega

dtheta_dt = omega

domega_dt = - (g / L) * np.sin(theta) - (b / (m * L**2)) * omega # Damping term included

return [dtheta_dt, domega_dt]

def solve_damped_pendulum(theta0, omega0, L, g, b, m, t_max, dt):

t_eval = np.arange(0, t_max, dt)

sol = solve_ivp(damped_pendulum_equations, [0, t_max], [theta0, omega0], t_eval=t_eval, args=(g, L, b, m))

return sol.t, sol.y[0], sol.y[1] # Returns time, angular position, and velocity

# Example parameters

L = 1.0 # Length of pendulum (meters)

g = 9.81 # Acceleration due to gravity (m/s²)

b = 0.05 # Damping coefficient (kg m²/s)

m = 1.0 # Mass of the pendulum bob (kg)

theta0 = np.radians(30) # Initial angle in radians

omega0 = 0.0 # Initial angular velocity

t_max = 10 # Simulation time (seconds)

dt = 0.01 # Time step

# Solve and plot results

t, theta, omega = solve_damped_pendulum(theta0, omega0, L, g, b, m, t_max, dt)

plt.figure(figsize=(10,4))

plt.subplot(1,2,1)

plt.plot(t, theta, label="Angular Position (rad)")

plt.xlabel("Time (s)")

plt.ylabel("Theta (rad)")

plt.legend()

plt.grid()

plt.subplot(1,2,2)

plt.plot(t, omega, label="Angular Velocity (rad/s)", color='r')

plt.xlabel("Time (s)")

plt.ylabel("Omega (rad/s)")

plt.legend()

plt.grid()

plt.tight_layout()

plt.show()