Week 3 - Solving second order ODEs

AIM-

Using the second ODE to describe the transient behaviour of a system of a simple

the pendulum on python scripting.

OBJECTIVE-

In Engineering, ODE is used to describe the transient behaviour of a system. A simple example is a pendulum

The way the pendulum moves depends on Newton's 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, you will be able to simulate any of these 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-

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

After doing this is how the program looks like

PYTHON PROGRAM-

# import library
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.figure(1)
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(2)
	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)

GRAPHICAL REPRESENTATION -

PENDULUM VISUALIZATION -