Week 3 - Solving second order ODEs

SOLVING SECOND ORDER EQUATION USING PYTHON.

AIM
Using second ODE to describe the transient behaviour of a system of simple pendulum on python scripting.

OBJECTIVE

In Engineering, ODE is used to describe the transient behavior of a system. A simple example is a pendulum
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.
Similarly, there are hundreds of ODEs that describe key phenomena and if you have the ability to solve these equations then
you will be able to simulate any of theses systems.
Your objective is to write a program that solves the following ODE. This ODE represents the equation of motion of a simple
pendulum with damping
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.
Write a program in Python that will simulate the pendulum motion, just like the one shown in the start of this challenge.
use,
L=1 metre,
m=1 kg,
b=0.05.
g=9.81 m/s2.
Simulate the motion between 0-20 sec, for angular displacement=0,angular velocity=3 rad/sec at time t=0.

PROCEDURE

• integrate module is imported for ODE equation by the syntax from scipy.integrate import odeint.
  
• Then the variables and arguments are defined using function def mode1(theta,t,b,g,l,m):
  
• Then the second ODE equation is breakdown into two single ODE, the first equation will give displacement and the
  second equation will give velocity.
  
• Then the initial condition is defined by theta_0 = [ 0,3]
  
• To obtain the smooth curve the points are divide into 150 points between 0 to 20 sec by t = np.linspace(0,20,150)
  
• Then the function odeint are stored as theta.
  
• The displacement values are extracted and stored in null array THETA3 using append command.
  
• Then THETHA4 for loop is used to extract value from THETA3.

PROGRAM

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


#Damping coefficient
b = 0.02
#Gravitational acceleration
g = 9.81
# length
l = 1
#Mass 
m = 0.1

#Function that return dz/dt
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

#Initial Condition
theta_0 =[math.pi/2,0]
#Time Point
t = np.linspace(0,10,150)
#Solve ODE

theta = odeint(model,theta_0,t,args = (b,g,l,m))
print(theta)
plt.plot(t,theta[:,0])
plt.plot(t,theta[:,1])


plt.xlabel('Time')
plt.ylabel('Displacement/Velocty')
plt.legend(['Displacement','Velocity'])
plt.show()




#Pendulum Bob visualization 
ct = 1




# Starting point
x0 = 0
y0 = 0

L1 =1

for t1 in theta[:,0]:
	X1 = L1*math.sin(t1)
	Y1 = -L1*math.cos(t1)
	plt.figure()
	plt.plot([x0,X1],[y0,Y1])
	plt.plot(X1,Y1,'o',markersize = 15,markerfacecolor = 'red')
	plt.xlim([-1.5,1.5])
	plt.ylim([-1.5,1.5])
	file='test%05d.png'%ct 
	ct = ct+1
	plt.savefig(file)

THE GRAPHICAL REPRESENTATION

PENDULUM BOB VISUALIZATION