Solving the 2nd order ODE for a pendulum using python programming

AIM:

To solve the second order ODE for a pendulum using python programming. 

PROBLEM STATEMENT:

1.To plot a graph between the variation of plot versus time.

2.To create the animation of the movement of the pendulum for various values of theta.

THEORY:

A pendulum is a weight suspended from a pivot so that it can swing easily. When the pendulum is displaced sideways from its resting position, it is subjected to a restoring force or a damping force due to the action of gravity that will accelerate it back to its equilibruim position. When the pendulum is released, the restoring force acting on the pendulums mass causes it to oscllate about its equilibrium position. This particular phenomena is being studied and analyzed using the Ordinary Differential Equation. The ordinary differential equation which is being derived for movement of the pendulum with respect to the various changes in the values of theta. The derived equation is as follows:

`frac{d^2theta}{dt^2} + frac{b}{m}cdotfrac{d theta}{dt} + frac{g}{l}sin(theta) = 0`

Where,

b = damping co-effiecient.

m = mass in kg.

g = acceleration due to gravity in m/s^2.

l = length of the pendulum in m.   

We are solving for the value of theta and also the rate at which the value of theta is changing from time to time with respect to velocity. In order to solve this, we form 2 ordinary differential equations which is needed respectively for displacement and velocity of the pendulum. The derivation goes as follows:

Let,

`theta = theta1`

Differentiating the following we get,

`frac{d theta}{dt} = frac{d theta1}{dt} = theta2`

Differentiating the equaton again with respect to t,

`frac{d^2theta}{dt^2} = frac{d^2theta1}{dt^2}`

From the above equation we get,

`frac{d^2theta1}{dt^2} = frac{d}{dt}(frac{d theta1}{dt}) = frac{d (theta2)}{dt}`

We must change the primary equation according to the following conditions given,

`frac{d^2theta1}{dt^2} + frac{b}{m}cdotfrac{d theta1}{dt} + frac{g}{l}cdot sin(theta1) = 0`_ _ _ _ _ _ _.1

Since, 

`frac{d theta1}{dt} = theta2`_ _ _ _ _ .2

`frac{d^2theta1}{dt^2} = frac{d theta2}{dt}`_ _  _ _ _ _ .3

From .1, .2, .3 the general equation changes as,

`frac{d theta2}{dt} = -frac{b}{m}cdottheta2 -frac{g}{l}cdot sin(theta1)`_ _ _ _ _.4

METHODOLOGY AND ALGORITHM:

1. To start with the program we are importing the necessary modules such as math, matplotlib, numpy, scipy.The math module is used for calculations and inserting the required math functions, the matplotlib module is used for plotting the final result in the form of graphs, the numpy module is used in this program for creating values which are evenly spaced sequence in a specified interval, the scipy module is used for intergrating the ordinary differential equations to get results.

2.The constant or inputs of the program are being defined such as 

b = damping co-efficient.

g = acceleration due to gravity.

l = length of the pendulum.

m = mass of the bob.

3.Since we are solving 2 Ordinary Differential Equation of first order. We have to set the initial conditions (i.e) Theta_0 = [math.pi/2,0]. The first ODE gives us the value of displacement(i.e) math.pi/2 gives us the displacement here and the second condition (0) is for the velocity. 

`theta2 = frac{d theta1}{dt}`

Since `theta1 =`Displacement.

Rate of change of displacement  = velocity`

4.A time aray is being created using the numpy module.Hence, a variable t is used which goeas as follows:

t = np.linspace(0,10,150)

The meaning of this state is that time vares from 0 to 10 seconds and will provide 150 outputs for these amount of time. Linespace is a command used for creating an evenly spaced sequence in a specified interval.

5.Solving the ODE: Uisng the odeint command we can solve the ordinary differential equation. To run this command we require the scipy module. This functions man purpose is to integrate a system of ODE. In order to integrae, we must define our function, initial conditions, time and arguments which are the constants or inputs of the program.

6. Defining the function(model):

#Defining the 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

Line no .4 and .5 in the function represent the first and second ODE respectively. We are assigning dtheta_dt, dtheta2_dt to a common variable dtheta_dt. We are doing this, so that we combine theta1 and theta2 into theta and we are returning the function using return dtheta_dt.

The odeint command is used. Its purpose is to take the initial conditons given and it is going to integrate it to the given time that is specified in the time steps. The output of odeint which solves for model,theta_0,t,args = (b,g,l,m) are stored in a variable theta. The odeint command changes with respect to the inputs given and it is quite generic and can be used for various range of ODE.

7.Plotting the result: For plotting the result we use the matplotlib module. The code for plotting goes as follows:

#plot results

plt.plot(t,theta[:,0],'b-',label = r'$frac{dtheta_1}{dt} = theta2$')
plt.plot(t,theta[:,1],'r--',label = r'$frac{dtheta_2}{dt} = -frac{b}{m}theta_2-frac{g}{t}sintheta_1$')
plt.xlabel('TIME')
plt.ylabel('PLOT')
plt.legend( loc = 'best')
plt.savefig('2ND ORDER ODE')
plt.show()

In the first line we are plotting time against theta[:,0]. This line represents the first order ODE related to displacement. The meaning of the solution vector theta is that the x co-ordinate of theta represents the displacement and the ':' operator means it take any number of values for the displacement and for the y co-ordinate of the solution vector theta, the velocity is '0' at the beginning. There are specific representations given as 'b-' which is a blue colour line that represents the plot. The legend of that plot is being written preceeding it.

Then the next line of code is similar to the first one but the difference is that it represents the first order ODE related to velocity. The meaning of the solution vector theta is that the displacement value and it is a solution vector which contains 2 columns and multiple rows. The first column is giving us the solution for the first ODE. Theta[ :,0] gives us the displacement values whereas theta[:,1] gives us the velocity value. We are giving the plot 

xlabel = 'TIME' and ylabel = 'PLOT'. Then we provide the location of the legend using the plt.legend command which selects the best position to place the legend using the legend command. We use the plt.savefig() command to save the command with the respective name assigned to it.

8.Code for animation of the pendulum: 

#Animation of the motion of pendulum

#Base of the pendulum
x0 = 0
y0 = 0

#Using the counter variable to append in each loop
counter = 1

#Using the for loop to generate range of values for theta

for i in theta[:,0]:

#Resolving the components respectively

	x1=l*math.sin(i)
	y1=-l*math.cos(i)

	filename = str(counter) + '.png'
	counter = counter + 1

	plt.figure()
	plt.plot([-1,1],[0,0],'red',linewidth = 10)
	plt.plot([x0,x1],[y0,y1],'black')
	plt.plot(x1,y1,'o',markersize=20)
	plt.xlim(-1.5,1.5)
	plt.ylim(-1.5,1.5)
	plt.savefig(filename)
	plt.plot()

1. At first we assign the base of the pendulum at x0 = 0 and y0 = 0.

2. Then a counter variable is used in order to append the value of counter on each iteration that takes place. Here there are 150 iterations that take place due to the command given in the np.linspace(0,10,150).

3. We are using a for loop for generating various values for theta using the variable 'i'. The components of the end point of the pendulum have been resolved as 

`x1 = lcdotmathcdot sin(i)`

`y1 = -lcdotmathcdotcos(i)`

 The reason for the negative sign is because in the plot if both the x and y component lie in the first quadrant, then the pendulum tends to swing about in the first quadrant itself. It is shown in the picture below:

If we would have used both the x and y component in th first quadrant then there will be no role for the accleration due to gravity. Hence, the plot is changed as follows with a negative sign for the y component and it is also the end point of the pendulum.On fixing the plot we get the result as follows:

 4. Then we use a variable called filename to save the output obtained in the form of images using the command: filename = str(counter) + '.png' . Using this command we can save the filename as the number with .png extension.

5. The plt.figure() command is used in order to save the output obtained on each iteration. The counter variable is added on each iteration.

6.Using the plt.plot() function, we are providing the length of the pivot which is represented with a linewidth = 10. Then the length of the pendulum is given using the co-ordinates (x0,y0) and (x1,y1) from the pivot point, it is assigned the colour black.

7.The marker size is used to represent the bob of the pendulum and it is specified at the end of (x1,y1), it is represented by 'o' with a size of 20.

8.The limits of are plots are provided at (-1.5,1.5) using the plt.xlim() function.

9.The file is saved using the plt.savefig(filename) function. To visualize the plot we use plt.show().

PYTHON CODE:

#Soving the 2nd order ODE of a pendullum using python programming

#Importing required libraries
import math
import matplotlib.pyplot as plt 
from scipy.integrate import odeint
import numpy as np 

#Defining the 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

#Input Parameters

b = 0.02 #(Damping Co-efficient)
g = 9.81 #(Acceleration due to gravity,Unit = m/s^2)
l = 1    #(Length of the pendulum,Unit = m)
m = 0.1  #(Mass of the bob,Unit = kg)

#Initial condition
theta_0 = [0,math.pi/2]

#Time points

t = np.linspace(0,10,150)

#solve ODE

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

#plot results

plt.plot(t,theta[:,0],'b-',label = r'$frac{dtheta_1}{dt} = theta2$')
plt.plot(t,theta[:,1],'r--',label = r'$frac{dtheta_2}{dt} = -frac{b}{m}theta_2-frac{g}{t}sintheta_1$')
plt.xlabel('TIME')
plt.ylabel('PLOT')
plt.legend( loc = 'best')
plt.savefig('2ND ORDER ODE')
plt.show()

#Animation of the motion of pendulum

#Base of the pendulum
x0 = 0
y0 = 0

#Using the counter variable to append in each loop
counter = 1

#Using the for loop to generate range of values for theta

for i in theta[:,0]:

#Resolving the components respectively

	x1=l*math.sin(i)
	y1=l*math.cos(i)

	filename = str(counter) + '.png'
	counter = counter + 1

	plt.figure()
	plt.plot([-1,1],[0,0],'red',linewidth = 10)
	plt.plot([x0,x1],[y0,y1],'black')
	plt.plot(x1,y1,'o',markersize=20)
	plt.xlim(-1.5,1.5)
	plt.ylim(-1.5,1.5)
	plt.savefig(filename)
	plt.plot()

OUTPUT GRAPH:

ANIMATION OF THE PLOT:

OUTPUT GENERATED FOR THE ANIMATION:

We have obtained 150 plots for the given experiment.

LINK TO THE ANIMATION OF THE PLOT:

https://www.youtube.com/watch?v=jSMDpkS0Z9k

RESULT & CONCLUSION:

From the given report we can conclude that the graph is plotted for change of theta with respect to the time and the animation of the movement of the pendulum for various values of theta is obtained.

REFERENCES:

1. https://en.wikipedia.org/wiki/Pendulum#:~:text=A%20pendulum%20is%20a%20weight,back%20toward%20the%20equilibrium%20position.