Simulate the transient behaviour of a simple pendulum by solving a 2nd order differential equation and to create an animation of it s motion using Python

Title:

Simulate the transient behaviour of a simple pendulum by solving a 2nd order differential equation and to create an animation of it\'s motion using Python.

Introduction:

A simple pendulum consist of a bob of mass(m) attached at the end of a massless string which is fixed at one end(also called as origin), it performs a to and fro motion.

The motion of a simple pendulum is also called a simple harmonic motion (SHM).The real pendulum is the one whose oscillation die out with time.this effect is called as damping. The equation which governs this is:

`(d^2theta)/dt^2 + b/m((d theta)/dt) +g/L sintheta = 0`

this is a 2nd order differential equation because the first term has the power of 2.

`theta`= angular displacement

b = damping coefficient

m = mass of bob

g = accleration due to gravity

l = length of the pendulum

for solving the 2nd order ODE the above equation needs to be reduced to its simplest form

`(d^2theta)/dt^2 + b/m((d theta)/dt) +g/L sintheta = 0`

`theta = theta_1`

`(d theta)/dt = (d theta_1)/dt = theta_2`_________(1)

`(d^2 theta)/dt^2 = (d^2 theta_1)/dt^2 = d/dt((d theta_1)/dt) = (d theta_2)/dt`_____(2)

`(d^2theta_1)/dt^2 + b/m((d theta_1)/dt) +g/L sintheta_1 = 0`

`(d theta_2)/dt +b/m theta_2 +g/L sin theta_1 = 0`

the equation are define in coding as a function

Code Explanation:

• Importing basic libraries math, matplotlib (because intrest is to plot the result), numpy and from scipy importing \'odeint\'(it is the function which will integrate the ordinary diff equation)
• defining a function as model consisting theta, t, b,g,l,m
• then wrinting the code as explained above
• entering the given inputs(b = 0.05, g = 9.81, l = 1, m= 1)
• giving initial conditions to theta(0,3) as 0 is the angular displacement and 3 as the angular velocity
• then inputing the time values from 0 to 20 using numpy, and there will be 150 images for the position of the bob between 0 to 20 seconds.
• then solving the ODE using odeint function, where model represents the system, theta_0 and t represents the initial condition and time and inputing the arguments.
• than we will plot the graph where [:,0] states multiple rows and one column this represents the displacement values, and [:,1] represents the velocity values
• than taking the displacement values [:,0] as h, taking counter variable as 1 for counting the images generated.
• now taking x0, y0 as 0 as the fixed position of the pendulum.
• now using for loop for taking different values of theta at different position and time step and calculating the different positions of the bob
• now creating filenames by using ct as a string \'str\' and creating the images in .png format.
• now plotting the baseline where the fixed position of the bob is to be held
• now creating a circle using markersize as 25
• creating x and y limits
• creating title as\'pendulum motion\'
• now finally using plt.show() for displaying the plots created.

Python code

#Simulate the transient behaviour of a simple pendulum by solving a 2nd order differential equation and to create an animation of it\'s motion using Python.

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

#function that returns dz/dt
def model(theta,t,b,g,l,m):

	theta1 = theta[0]
	theta2 = theta[1]
	dtheta1_dt = theta2
	dtheta2_dt = -(b/m)*theta2 - (g/l)*math.sin(theta1)
	dtheta_dt = [dtheta1_dt, dtheta2_dt]
	return dtheta_dt

# inputs
b=0.05    #damping coefficient
g=9.81    #gravity
l=1       #length of the pendulum
m=1       #mass

# initial conditions
theta_0 = [0,3]   #[angular displacement,angular velocity]	

#time points
t = np.linspace(0,20,150)  # here we need to calculate the motion between 0 to 20 sec at 150 positions.

# solve ODE
theta = odeint(model,theta_0,t,args = (b,g,l,m))


#plot results
plt.plot(t,theta[:,0],\'b-\',label=r\'$\\frac{d\\theta_1}{dt}=\\theta2$\')
plt.plot(t,theta[:,1],\'r--\',label=r\'$\\frac{d\\theta2}{dt}=-\\frac{b}{m}\\theta_2-\\frac{g}{L}sin\\theta_1$\')
plt.ylabel(\'Plot\')
plt.xlabel(\'time\')
plt.legend(loc=\'best\')
plt.show()

#cerating images for different pendulum position
h = theta[:,0]
ct = 1

#initial fixed position of the pendulum
x0 = 0
y0 = 0

for i in h:
	x1 = x0 + (l*math.sin(i))
	y1 = y0 - (l*math.cos(i))
	ct = ct + 1
	filename = str(ct) + \'.png\'

	plt.figure()
	plt.plot([-1,1],[0,0],\'b-\',linewidth = 3)
	plt.plot([x0,x1],[y0,y1])
	plt.plot(x1,y1,\'o\',markersize = 25)
	plt.xlim([-2,2])
	plt.ylim([-2,1])
	plt.title(\'pendulum motions\')

	plt.savefig(filename)
	plt.show()

Result of the plot

errors faced while coding

this was the major error i faced, after the for loop while generating y1 values i initially used the + sign which gave me the plot of the pendulum on the upper side of the fixed point an image showing the error from the 150 images generated is shown below

video of motion of pendulum for different position and time step

the link for the video

https://youtu.be/UH-xFhyKm1A

<<<<<Thank You>>>>>