Week 3 - Solving second order ODEs

"""
Pendulum 
Solving second order ODE
"""
import math 
import numpy as np
from scipy.integrate import odeint
import matplotlib.pyplot as plt

#function called for returning dteta_dt
def model(teta,t,b,g,l,m):
	#equation of 2nd order is breaked down to single ode writen below
	teta1=teta[0]
	teta2=teta[1]
	dteta1_dt=teta2
	dteta2_dt=-(b/m)*teta2 - (g/l)*math.sin(teta1)
	dteta_dt=[dteta1_dt,dteta2_dt]
	return dteta_dt
#inputs 
b=0.02
g=9.81
l=1
m=0.1

#intial condition
teta_0=[math.pi/2,0]

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

#solving ODE
teta = odeint(model, teta_0,t, args = (b,g,l,m))
#Call odeint to generate the solution. To pass the parameters b,g,l and m to model,
# we give them to odeint using the args argument

#plot 
plt.plot(t,teta[:,0],'b-',label='displacement')
plt.plot(t,teta[:,1],'r--',label='velocity')
plt.ylabel('plot')
plt.xlabel('Time')
plt.legend(loc='best')
plt.show()

#program for motion of pendulum

#one end is fixed
x0=0
y0=0
m=teta[:,0]
ct=1

for r in m:
	#other end of rod has a ball which is in motion and its coordinates are x1,y1

	x1=(l*math.sin(r))
	y1=-(l*math.cos(r))
	filename='RM%d.png'%ct #to save all the images in each itteration

	plt.figure()
	plt.plot([x0,x1], [y0,y1],'b-')
	plt.plot(x1,y1,'o')
	plt.xlim([-2,2])
	plt.ylim([-2,2])
	plt.show()
	plt.savefig(filename)
	ct=ct+1
	

The report is attached below along with the animation