Week 3 - Solving second order ODEs

AIM: To Write a program in Python that will simulate the pendulum motion and to create an animation of its motion.

OBJECTIVE:

• To write a program to solve the ODE and stimulate the transient behavior of simple pendulum.
• To create animation of the motion.

GOVERNING EQUATIONS:

• Equation of motion of a simple pendulum with damping

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

In the above equation,

g = gravity in m/s2,

L = length of the pendulum in m,

m = mass of the ball in kg,

b=damping coefficient.

`theta`=It is the corresponding angle with respect to the vertical.

• Here to solve the ode we are going to use the equations in first order,

`(dθ_1)/dt=θ_2` and `(d theta_2)/dt=-b/mtheta_2-g/lsintheta_1`

INPUT PARAMETERS:

Damping coefficient (b) =0.05

Gravitational constant (g) =9.81

Length of the pendulum (l) =1

Mass of the bob (m) =1

BODY OF THE CONTENT:

CODE FOR MOTION OF PENDULUM:

#Program to stimulate the pendulum motion
import math
import matplotlib.pyplot as plt
import numpy as np
from scipy.integrate import odeint

#defining function to give the ode equation
def ODE_func(theta,t):

	b=0.05     #damping co-efficient
	l=1        #length of pendulum in metre
	g=9.81     #Gravity in m/(s^2)
	m=1        #mass of ball  in kg

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

#initial conditiion
theta_0=[0,3]

#time span

t=np.linspace(0,20,200)

theta=odeint(ODE_func,theta_0,t)

plt.plot(theta[:,0],'g-',label= 'angular.displacement(rad)')
plt.plot(theta[:,1],'b-',label= 'angular velocity(rad/sec)')
plt.xlabel('Time')
plt.ylabel('theta')
plt.legend()
plt.title('Transient behaviour of simple pendulum')
plt.savefig("Transient behaviour plot.png")
plt.show()

#code for images of position of pendulum

ct=1

for i in theta[:,0]:
	x0=0
	y0=0
	l=1.5    #length of string

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

	filename = 'img%05d.png' %ct
	ct=ct+1
	plt.figure()
	plt.plot([-1,1],[0,0],'g-',linewidth=10)
	plt.plot([0,0],[0,-2],'g--')
	plt.plot([x0,x1],[y0,y1],'b-',linewidth=5)
	plt.plot(x1,y1,'ro',markersize=35)
	plt.xlim([-1.5,1.5])
	plt.ylim([-2.5,0.5])
	plt.savefig(filename)

EXPLANATION OF THE CODE:

1.)Explanation of code to plot the graph of transient behaviour of pendulum

• Firstlywe import the required math module,odeint from scipy module and the required libraries.
• Then we define the function ODE_func which includes values of damping coefficient(b),gravity(g),length of pendulum (l),mass of the ball (m) and this function gives the values of theta's i.e (dtheta1/dt) and (dtheta2/dt) and stores it in dtheta and returns dtheta whenever function is called.For solving the ode's we need a different library.
• Then we assign values for theta at initial level in matrix form because it cointains 2 values which are angular displacement and angular velocity.
• Then we define a time span and convert this span into equal interval by using linspace command.
• We will further use odeint command to solve the first order ode with the help of the ODE_func ,the theta matrix containing values of angular displacement and angular velocity and the timespan .
• Now we use the plot command to get graphs for angular displacement with respect to time and angular velocity with respect to time.
• We use the legend command to add the legends and the title command to give title to the plot.

2.) EXplanation of code to get the stimulation motion of pendulum

• Fisrtly we will use the counter that we help to keep adding up new figure as the for loop functions.
• Further we define the for loop in which we assign co-ordinates of string at fixed end (x0,y0)=(0,0).
• The we will calculate the co-ordinates of free end of the string(x1,y1)
• Now we start with plotting, here we first start we plotting the fixed line in which the string is pivoted.
• Further we use plot command to plot a line which helps to show equlibrium position of the ball
• Now we use plot to get the string of the pendulum.
• Further we use the plot command to get the ball and its position,so as the loop continues for different values of (x1,y1) we get the different position of the ball.
• Then we define the limits for the ball motion.
• Finally we use a command to save each and every graph that we get in to a folder.
• After we get the plots we use the imageMagick to stich the graphs and get the GIF of the motion of pendulum and then convert it to animatiion.

OUTPUTS:

• Plot of transient behaviour of simple pendulum:

• Link for the animation of pendulum motion

https://www.youtube.com/watch?v=ksj2Fhyo1e4&feature=youtu.be

CONCLUSION:

So using the Python program we were able to solve the ode for a simple pendulum and make an animation of its motion.