Week 3 - Solving second order ODEs
import numpy as np
from scipy.integrate import odeint
import matplotlib.pyplot as plt
import math
# functionthat 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.02
g = 9.81
l = 1
m = 0.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,ArgSpec = (b,g,l,m))
# plot results
plt.plot(t,theta[:,0],'b-',label=r'sfrac{dtheta_1}{dt}=theta2s')
plt.plot(t,theta[:,1],'r--',label=r'sfrac{dtheta2}{dt}=-frac{b}{m}theta_2-frac{g}{L}sintheta_1s')
plt.ylable('plot')
plt.xlabel('time')
plt.legend(Loc='best')
plt.show()