Week 3 - Solving second order ODEs

Objective :

To write a program that motion of a simple pendulum with damping

Governing equation

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.

Inputs to equation

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.

Python program for ODE

import numpy as np
from scipy.integrate import odeint
import matplotlib.pyplot as plt
import math as mt
def model(theta,t,b,g,l,m):
	theta1=theta[0]
	theta2=theta[1]
	dtheta1_dt = theta2
	dtheta2_dt = -((b/m)*theta2)-((g/l)*mt.sin(theta1))
	dtheta_dt = [dtheta1_dt , dtheta2_dt]
	return dtheta_dt

#inputs 
b = .05
g = 9.81
l = 1
m = 1

#intial conditions
theta_0 = (0,mt.pi/2)
#time
t = np.linspace(0,20,100)
#solve ODE
theta = odeint(model,theta_0,t,args=(b,g,l,m))

#plot results
plt.plot(t,theta[:,0],'b-',label=r'$\frac{d\theta_1}{dt}= \theta_2$')
plt.plot(t,theta[:,1],'r--',label=r'$\frac{d\theta_2}{dt}=-\frac{b}{m}\theta_2-\frac{g}{l}sin\theta_1$')
plt.ylabel('Plot')
plt.xlabel('Time')
plt.legend(loc='best')
plt.show()


#Animation of Simple pendulum


x0=0    #this is fixed point at the wall
y0=0

ct=1
for ang in theta[:,0]:  #angle

	x1=l*mt.sin(ang)
	y1=-l*mt.cos(ang)

	filename='a'+str(ct)+'.png'
	ct=ct+1

	plt.figure()
	plt.plot([-1,1],[0,0],'red',linewidth=10)
	plt.plot([x0,x1],[y0,y1],'black')
	plt.plot(x1,y1,'o',markersize=20)
	plt.xlim(-1.5,1.5)
	plt.ylim(-1.5,1.5)
	plt.savefig(filename)

Output

Angular Displacement and angular velocity variation with time

Bob Oscillation animation

Output files are animated and uploaded in youtube.com

https://youtu.be/n55YdJZN9GU