Simulation of a Simple Pendulum by Python

Objective: Write a program to simulate the transient behavior of a simple pendulum and to create an animation of its motion

Introduction: A pendulum is a weight suspended from a pivot so that it can swing freely.When a pendulum is displaced sideways from its resting, equilibrium position, it is subject to a restoring force due to gravity that will accelerate it back toward the equilibrium position. When released, the restoring force acting on the pendulum's mass causes it to oscillate about the equilibrium position, swinging back and forth.

Differential Equation of Motion of a Simple Pendulum:

CODE:

# Importing libraries

import math
import numpy as np
from scipy.integrate import solve_ivp
import matplotlib.pyplot as plt
import os

# Initial and end values
st = 0          # Start time (s)
et = 20.4       # End time (s)
ts = 0.1        # Time step (s)
g = 9.81        # Acceleration due to gravity (m/s^2)
L = 1           # Length of pendulum (m)
b = 0.5         # Damping factor (kg/s)
m = 1           # Mass of bob (kg)


# 1st order equations to solve in a function


"""
 theta1 is angular displacement at current time instant
 theta2 is angular velocity at current time instant
 dtheta2_dt is angular acceleration at current time instant
 dtheta1_dt is rate of change of angular displacement at current time instant i.e. same as theta2 
"""

def sim_pen_eq(t,theta):
		dtheta2_dt = (-b/m)*theta[1] + (-g/L)*np.sin(theta[0])
		dtheta1_dt = theta[1]
		return [dtheta1_dt, dtheta2_dt]



# main

theta1_ini = 0                 # Initial angular displacement (rad)
theta2_ini = 3                 # Initial angular velocity (rad/s)
theta_ini = [theta1_ini, theta2_ini]
t_span = [st,et+ts]
t = np.arange(st,et+ts,ts)
sim_points = len(t)
l = np.arange(0,sim_points,1)

theta12 = solve_ivp(sim_pen_eq, t_span, theta_ini, t_eval = t)
theta1 = theta12.y[0,:]
theta2 = theta12.y[1,:]
plt.plot(t,theta1,label='Angular Displacement (rad)')
plt.plot(t,theta2,label='Angular velocity (rad/s)')
plt.xlabel('Time(s)')
plt.ylabel('Angular Disp.(rad) and Angular Vel.(rad/s)')
plt.legend()
plt.show()



# Simulation

x = L*np.sin(theta1)
y = -L*np.cos(theta1)


for point in l:
		plt.figure()
		plt.plot(x[point],y[point],'bo',markersize=20)
		plt.plot([0,x[point]], [0,y[point]])
		plt.xlim(-L-0.5,L+0.5)
		plt.ylim(-L-0.5,L+0.5)
		plt.xlabel('x-direction')
		plt.ylabel('y-direction')
		filenumber = point
		filenumber=format(filenumber,"05")
		filename="image{}.png".format(filenumber)
		plt.savefig(filename)
		plt.close()


os.system("ffmpeg -f image2 -r 20 -i image%05d.png -vcodec mpeg4 -y movie.avi")	

EXPLANATION:

After solving the diffrential equation, the code will generate multiple frames of position of pendulum corresponding to each time instant in the array 't'. Each of these frames or figures will be saved with the name format as 'image00000.png','image00001.png','image00002.png' and so on. Then, the code will call the os and pass a comand line to the os to stich the images together to make a movie. Note that the definition of this comand line is in ffmpeg library. The figures and the movie will be saved in the same directory where the code file is saved

OUTPUT:

LINK:

https://youtu.be/fcGbmzE7-Y8