Week 3 - Solving second order ODEs
import 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] #displacement
theta2 = theta[1] #velocity
dtheta1_dt = theta2
dtheta2_dt = -(b/m)*theta2-(g/l)*math.sin(theta1)
dtheta_dt = [dtheta1_dt, dtheta2_dt]
return dtheta_dt
#constants
b=0.05
g=9.81
m=1
l=1
#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 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\theta_2}{dt}=-\frac{b}{m}\theta_2-\frac{g}{L}sin\theta_1$')
plt.ylabel('plot')
plt.xlabel('time')
plt.legend(loc='best')
plt.show()