Week 3 - Solving second order ODEs

# Motion of Pendulum Challenge ODE

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

#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

#Input
l = 1
m = 0.1
g = 9.81
b = 0.05

#Initial condition
theta_0 = [0,3]

# Time points
t = np.linspace(0,20,100)

# Solve ODE
theta = odeint(model,theta_0,t,args = (b,g,l,m))

#plot
plt.plot(t,theta[:,0],'b--',label = r'{dtheta1}/{dt} = theta2 ')
plt.plot(t,theta[:,1],'r--',label = r'{dtheta2}/{dt} = -{b}/{m}theta2 -{g}/{l}sin/theta1')
plt.xlabel('Time')
plt.ylabel('Plot')
plt.xlim([0,15])
plt.legend(loc = 'best')
plt.show()