Week 3 - Solving second order ODEs

# ODE solver Pendulum Problem
import numpy as np
from scipy.integrate import odeint
import matplotlib.pyplot as plt
import math


def model(theta,t,b,g,l,m): # function definition
	theta1=theta[0]
	theta2=theta[1]
	dtheta1_dt=theta2 # Second order diffierential equation simplified 
	dtheta2_dt=-(b/m)*theta2-(g/l)*math.sin(theta1) # Differential of theta 2 equal speed
	dtheta_dt=[dtheta1_dt,dtheta2_dt] # solution matrix
	return (dtheta_dt)


b=0.02 # Damping coefficient
g=9.81 # gravity
l=1 # rod length
m=0.1 # ball mass in kg
# initial condition
theta_0=[0,3] # At t=0 speed = 3 rad/sec
# time
t=np.linspace(0,20,150) # the time span from zero to 20 seconds with 150 steps
# 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}=\theta2$')
plt.plot(t,theta[:,1],'r--',label=r'$\frac{d\theta2}{dt}=-\frac{b}{m}\theta_2-\frac{g}{l}sin\theta_1$')
plt.xlabel('Plot')
plt.ylabel("time")
plt.legend(loc='best')
plt.show()

#Simple Pendulam Animation
ct = 1
x0 = 0
y0 = 0
for i in theta[:,0]:
	x1 = 1*math.sin(i)
	y1 = -1*math.cos(i)
	filename = '%shm.png' %ct
	ct = ct + 1

	plt.figure()
	plt.plot([-1,1],[0,0],"blue",linewidth=15)
	plt.plot([x0,x1],[y0,y1],"green")
	plt.xlim([-1.5,1.5])
	plt.ylim([-1.5,1.0])
	plt.plot(x1,y1,'o',markersize =30)
	plt.title('SHM Pendulum')
	plt.savefig(filename)