Simple Pendulum Simulation by solving 2nd Order ODE - Python

Simple Pendulum Simulation by solving second order ODE . 

.                                            

• The objective of this program is to simulate a simple pendulum by solving second order ODE into two first order ODE's. The way the pendulum moves depends on the Newtons second law. When this law is written down, we get a second order Ordinary Differential Equation that describes the position of the ball w.r.t time. 
• The governing differential equation representing the motion of simple pendulum with damping is given as                                  

•              
•  where , `Theta` = Displacement 

         b    = damping Co-efficient 

         g     = gravity m/s^2

         L     = length of string in (m)

         m    = mass of ball (kg)

• Upon solving the given Second order ode , we convert it into two first order ode's , which are  :

              `(d(theta1))/dt = theta2`         &     `(d(theta2))/dt = - b/m* theta2 - g/l sin*theta1`

Python Program : 

'''
Transient behaviour of simple pendulum by solving 2nd order ode

Code by AADIL SHAIKH
'''


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

# Function 
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
g = 9.81
l = 1
m = 1

#initial condition
theta_0 = [0,3]

# time points
t = np.linspace(0,20,300)

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

disp_data = theta[:,0]

# plot results 
plt.figure(1)
plt.subplot(2,1,1)
plt.plot(t,theta[:,0],'b-',label=r'$\frac{d\theta_1}{dt}=\theta2$')
plt.ylabel('Angular displacement',fontsize = 'large' , fontweight = 'bold')
plt.xlabel('Time',fontsize = 'large' , fontweight = 'bold')
plt.legend(loc='best')
plt.title('Angular_disp vs Time & Angular_vel vs Time',fontsize = 'large' , fontweight = 'bold')

plt.subplot(2,1,2)
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('Angular Velocity (rad/s)',fontsize = 'large' , fontweight = 'bold')
plt.xlabel('Time',fontsize = 'large' , fontweight = 'bold')
plt.legend(loc='best')
#plt.savefig('Angular.png')
plt.show()


ct = 1
for i in range(len(t)) :
	x0 = 0
	y0 = 0
	x1 = l*math.sin(theta[i,0])
	y1 = -l*math.cos(theta[i,0])

	filename = 'Pendu%04d.png' % ct
	ct = ct + 1
	
	plt.figure()
	plt.plot([x0,x1],[y0,y1],'red',linewidth=2)
	plt.plot([-1,1],[0,0],'g',linewidth = 4)
	plt.plot(x0,y0,'o',markersize =10,markerfacecolor = 'k')
	plt.plot(x1,y1,'o',markersize =30,markerfacecolor = 'm')
	plt.title('Pendulum system')
	plt.xlim([-1.2,1.2])
	plt.ylim([-1.5,0.5])
	plt.savefig(filename)

# https://youtu.be/V7UH3Ho2xZs 
#              
#                  THE END                                          #

Explanation : 

1. Import numpy, scipy for odeint solver, matlabplotlib for plots , math. 
2. We define the function next , where we write the two first order ode's mentioned in the start of this report. 
3. We combine both the ode's with one common output dtheta_dt . Which is the value we pass in theta where we call the function while solving with odeint.
4. All the inputs are defined as per conditions after the function .
5. Initial conditions are to simulate the motion between 0-20 seconds which is defined using linspace command, and angular dispalcement as 0 and angular velocity to be 3 rad/s
6. Initial conditions [0,3] called in the odeint pass that angular displacement is to be 0, and ang velocity to be 3 when t = 0. 
7. Then odeint is solved and arguements which were defined above are passed here with the initial conditions just described above.
8. Then plotting begins.  
9. Two subplots are made : one is Angular dispalcement vs time , and angular velocity vs time 
10. One data is of the angular displacement while the other is of angular velocity, and theyre called respectively and plotted.
11. Appropriate labels,legends and titles are given.
12. The equations are labeled in the plot itself.

Plot :      

1. Next is creating a Pendulum and Simulating the behaviour of this equation in it and plotting the animation. 
2. For this we used for loop, After defining the length of the string as 1, the co-ordinates are defined by the moment of the string according to sin and cos angles , as shown in the schematic picture in the beginning of report. 
3. The string is fixed at the pivot and moves at the x1, y1 location with a pendulum attached on it 
4. Pendulum ball is created using markersize . 
5. appropriate color, x-y axis limits and titles are defined.
6. A rigid base is created at the top with color green.
7. The for loops runs through the co-ordinates and saves the images in the filename giving each file a unique name, thus we end up with 300 images. 
8. Image magick is used to stitch the images together and form a mp4 video file which is uploaded in youtube. 
9. magic "Pendu%04d.png"[1-300] pendulum.mp4   command was used to stitch the images and animation is created simply. 

Animation : 

#                                                         THE END                                                                  #