Week 3 - Solving second order ODEs

mport numpy as np
from scipy.integrate import odeint
import matplotlib.pyplot as plt
import math

# function that returns dz/dt
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

b = 0.05
g = 9.81
l = 1
m = 1

# initial condition
theta_0 = [math.pi/2,0]

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

# 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}\theta2-frac{g}{l}sin\theta_1$')
plt.ylabel('plot')
plt.xlabel('time')
plt.legend(loc='best')
plt.show()