Week 3 - Solving second order ODEs

Python program to simulate the motion of a simple pendulum.

Aim: To plot a graph of velocity and displacement versus time using ode45 solver, and hence to create an animation of the simulation of a simple pendulum.

Theory:In mathematics, an ordinary differential equation (ODE) is a differential equation containing one or more functions of one independent variable and the derivatives of those functions.A Second order differential equation is an equation involving the unknown function y, its derivatives y' and y'', and the variable x.

Governing Equation:

`(d^2θ)/(dt^2)+b/m⋅(dθ)/(dt)+g/l⋅sin(θ)=0`

where

• `theta` - the angular dispacement (in degree)
• b - damping coefficient
• m - mass of the bob (in kgs)
• l - length of the string (in m)
• g - acceleration due to gravity (in `m/s^2`)

Code:

 
#program to simulate the movement of a simple pendulum
import math                       #importing required libraries
import matplotlib.pyplot as plt 
import scipy
import numpy as np
from  scipy.integrate import odeint

def ode_func(theta,t_span,b,l,m,g):           #defining the function

	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)



b=0.05;                           #Damping coefficient                                   
l=1;                              #length of the string
m=1;                              #mass of the bob
g=9.81;                            #acceleration due to gravity

#initial conditions
theta_0=[0,5];                    #determines the initial displacement and velocity
t_span=np.linspace(0,20,300);        #time span of 20 seconds



theta = odeint(ode_func,theta_0,t_span,args=(b,l,m,g)) #calling  odeint function
                           
 
plt.plot(t_span,theta[:,0],'-b',label=r'$\frac{d\theta_1}{dt}=\theta_2$')                #plot of displacement vs time
plt.plot(t_span,theta[:,1],'--r',label=r'$\frac{d\theta_2}{dt}=-\frac{b}{m}\theta_2-\frac{g}{l}sin\theta_1$') #plot of velocity vs time
plt.legend(loc='best')  #legend
plt.title(' Angular Displacement and Velocity vs time graph')  #establishing the title
plt.xlim(0,20)                  #x and y limits
plt.ylim(-10,10)                 #axis specification
plt.xlabel('time')               #labelling the x axis
plt.ylabel('Angular displacement and velocity')    #labelling the y axis
plt.show()   
       


#animation

plt.figure() 
ct=1                               #counter variable
for dis in (theta[:,0]):

	                              #value of theta of i^th row and first column  #determination of the coordinates
	x0=0
	y0=0
	    
	x1=3*math.sin(dis)
	y1=-3*math.cos(dis)
	  
	plt.figure()
	plt.plot([-1,1],[0,0],lw='2',color='r')             #plotting the slab
	plt.plot([x0,x1],[y0,y1],lw='2',color='b')        #plotting the string
	plt.plot(x1,y1,'o',markersize=10,markerfacecolor='r')     #plotting the bob using markers command                                                       #axis specification 
	plt.title('simulation of pendulum') 
	plt.xlim(-5,5)
	plt.ylim(-5,1)
	
	filename=str(ct)+'.png'                        #saving the file with unique file names
	plt.savefig(filename)
	ct=ct+1                                         #incrementing the counter variable
	 
	 
	 
	  
	  

Code Explanation:

• The code is started withimporting all the required libraries
• The values of the constants in the differential equation are determined

L=1 metre,

m=1 kg,

b=0.05.

g=9.81 m/s2.

• The initial conditions such as initial displacement and initial velocity is specified as variable named [theta_0].
• The span of time is determined using linspace command which generates 300 values between 0 to 20.
• The odeint function is used inorder to solve the second order differential equation.
• The displacements and the velocities for a span of 20 seconds are stored in the form of array inside theta.
• Using plot command the graph of Angular displcement and Angular velocity versus time is drawn.
• The xlabel and ylabel is specified
• The title of the graph is also determined using title command
• Inorder to distinguish between the two graphs legend command is used. The application of this command is reflected at the top right corner of the graph.
• All these plots are drawn in figure 1.
• The next figure is mainly used for the animation pupose.
• The counter variable is specified and the for loop begins for the purpose of the animation.
• Each value of displacement is specified using a variable dis throughout the length of the displacement values for a span of 20 seconds.
• The coordinates of the string is also specified.
• The bob is indicated using the markers command and its color is selected using markerfacecolor command.
• The title of each plot is kept same so that the same gets displayed in the animation.
• Each plot is saved with unique name , using the filename variable and it saves each plot with unique numbers.
• Finally, the counter variable is incremented and the for loop gets executed again.

Explanation of the function code:

• The function is named as ode_func
• An array of dtheta/dt has to be formed which contains theta1 and theta2 as the values of theta
• theta1 will be the first value of theta
• theta2 is the second value of theta
• dtheta_dt is a function of t,theta,b,l,m,g
• The values of dtheta1_dt and dtheta2_dt are determined using above formulas
• These two values will be stored in the array [dtheta_dt]
• Later the value is returned back to the function, and it will be used whenever the function will be called throughout the program.
• The same values are taken while plotting the graph.

Results and plots:

The graph is plot between the Angular velocity and Angular displacements versus time.It is found that the red dotted lines represent the Angular dispacement and the blue solid line represents the angular velocities of the bob.

The above image shows all the images of the pendulum with different angles,which will be stiched together to form the animation.

The link of the gif of the pendulum simulation is pasted below:

https://drive.google.com/file/d/176M--Qb6lIb5Zf4-Yje3PiyLHSgXsv6R/view?usp=sharing

Reference:

The following site is used to make the gif of the pendulum movement

https://ezgif.com/maker

Conclusion:

Python is a easier software which helps us to understand and execute a problem easily in the form of a program.This also help us to create animations easily.