Week 3 - Solving second order ODEs

"""
Problem solving of 2nd order differential equation for pendulum
"""
#The equation, d2theta/dt2+(b/m)dtheta/dt+(g/l)sintheta=0

#Now, if theta=theta1(taken) then dtheta/dt=dtheta1/dt=theta2
#Or, d2theta/dt2=d2theta1/dt2=d/dt(theta2)
#Now, d/dt(theta2)+(b/m)theta2+(g/l)sin(theta1)=0

import math
import matplotlib.pyplot as plt 
import numpy as np 
from scipy.integrate import odeint
#The function is 'model' based upon the ODE
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
# where, b=damping co-efficient
#m=mass , g=9.81m/s2, l=length
b=0.10
g=9.81
l=1.5
m=0.15
#initial condition where theta=0 and velocity=0, displacement=pi/2
theta_0=[0,5]


#Time period within which getting 100 points
t=np.linspace(0,10,100)
#Now solving the ODE of given model with co-efficients b,g,l,m within t time
theta=odeint(model,theta_0,t,args=(b,g,l,m))

#Plotting the results by  theta vs. time
#Solutin for displacement
plt.figure(1)
plt.plot(t,theta[:,0])
#Solution for velocity
plt.plot(t,theta[:,1])

d=theta[:,0]

ct=1
# Initial point
x0=0
y0=0
for i in d:
# End point
		x1=l*math.sin(i)
		y1=-l*math.cos(i)
		ct=ct+1
		filename=str(ct)+".png"
		plt.figure()
# For hanging support
		plt.plot([-2,2],[0,0])
# Position of string and support joint
		plt.plot([x0,x1],[y0,y1])
# "o" for sphere and markersize for sphere size
		plt.plot(x1,y1,'o',markersize=20)
# Plot chart axis size
		plt.xlim([-2,2])
		plt.ylim([-2,2])
		plt.ylabel('d')
		plt.xlabel('t')
		plt.savefig(filename)