Week 3 - Solving second order ODEs

Objective:- To write a program that solves the ODE which represents the equatiion of motion of a simple pendulum and to create an animation of it's motion.

Governing equation:-

Theory:- In engineering ,ODE is used to describe th etransient behavior of a system.

eg:- Simple pendulum

The way the pendulum moves depends on newton's second law.When this law is written down, we get a second order differntial Equation that describes the position of 'ball' w.r.t. time.

Simillarly, there are hundreds of ODEs that describes key phenomenon and if you have the ability to solve thee equations then you will be able to simulate any of these systems.

Where,

g =gravity in m/s2,

L= length of pendulum in m,

m = mass of the ball in kg,

b = damping coefficient.

Data given:-

• L = 1 meter,
• m = 1kg,
• b = 0.05
• g= 9.81 m/s^2.

Simulation of Motion between 0-20 sec

• for angular displacement = 0
• angular velocity = 3 rad/sec
• at time t = 0.

Explanation of code:-

• First we have to import numpy module as np, scipy.integrate then import odeint for second ODE,math and matplotlib for plots.
• Define functions that returns the value of dtheta/dt.
• Then initialise the input values of b,g,l,& m.
• Give conditions of theta i.e. 0 and 3rad/sec.
• using linspace command we have given 200 motions in the interval of 0 to 20.
• Solved ode by using odeint command.
• plt.plot command is used to plot graph.
• For Loop is used to create an animation.
• 'marker' command is used to show pendulum'BOB'.
• magick *.png pooja.gif is used to collect all 100 images in order to play an animation in seperate file.

1. A program in python that will simulate the pendulum motion.

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


#function that returns dtheta/dt

def function(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


b = 0.05
g = 9.81
l = 1
m = 1

#initial condition
theta_0 = [0,3]


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


#solve ode
theta = odeint(function,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],'g--',label=r'$\frac{d\theta_2}{dt}=-\frac{b}{m}theta_2-\frac{g}{1}sin\theta_1$')
plt.xlabel('plot')
plt.ylabel('Time')
plt.legend()
plt.savefig('Output.png')
plt.show()

	
#Base points for Pendulum
x_0 = 0
y_0 = 0
	

ct = 1

#Pendulum animation
for i in theta[:,0]:


	
	
	x1 = l*math.sin(i)
	y1 = -(l*math.cos(i))
	filename = str(ct)+'.png'
	ct = ct+1
	plt.figure()
	plt.plot([-1.5,1.5],[0,0],linewidth = 6)
	plt.plot([x_0,x1],[y_0,y1],linewidth = 4)
	plt.plot(x1,y1,marker = 'o',markersize = 20)
	plt.xlim([-1.5,1.5])
	plt.ylim([-1.5,1.5])
	plt.savefig(filename)
	

2.Results:-

• Angular displacement, angular velocity w.r.t. Time:-

• Animation of Simple Pendulum can be seen below, through link:-

https://drive.google.com/file/d/1nozo6ni98MF-GQ2sqj65I8OlB3E66V7X/view?usp=sharing