Week 3 - Solving second order ODEs

Aim/Objective: Solving second order ODE using Python

Formulae:

L=1 metre,

m=1 kg,

b=0.05.

g=9.81 m/s2.

Conditions/Parameters:

Simulate the motion between 0-20 sec, for angular displacement=0,angular velocity=3 rad/sec at time t=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.

Code:

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

#function
def model(theta,t,b,g,m,l):
	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

#inputs required
b=0.05
g=9.81
l=1
m=1

theta_0=[0,3]
t=np.linspace(0,20,200)

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

#plotting result
plt.figure()
plt.plot(t,theta[:,0],label=r'Angular displacement')
plt.plot(t,theta[:,1],label=r'Angular velocity')
plt.xlabel('plot')
plt.ylabel('Time')
plt.legend(loc='best')
plt.show()

#Animation step
ct=1
x0=5
y0=5
for i in range(len(t)):
	x1=x0+l*math.sin(theta[i,0])
	y1=y0-l*math.cos(theta[i,0])
	filename='pendulum'+str(ct)+'.png'
	


	plt.figure()
	plt.plot([4,6],[5,5],'k-')
	plt.plot([x0,x1],[y0,y1])
	plt.plot(x1,y1,marker='o',markersize='10')
	plt.axis([3.5,6.5,3,6])
	plt.savefig(filename)
	ct=ct+1

Explanation:

• Iimporting required modules using 'import' command.
• Function model is created with required arguments, inside which the term which are broken to solve second order ODE are calculated.
• Inputs are provided for the function.
• theta_0 ranges between 0 to 3 m/s^2.
• odeint command is used to solve ode function.
• Then results are plotted using plotting command.
• For creating animation we open a for loop, ranging from 0 to length of t.
• Using savefig command images are saved.

Graph: