Week 3 - Solving second order ODEs
Youtube link for Solving second order ODEs : https://youtu.be/xfnLSloWSlU
''' A program for solving second order ODEs and simulate its motion'''
import math
import matplotlib.pyplot as plt
import numpy as np
from scipy.integrate import odeint
def ode_function(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
# Inputs
# damping coefficient
b = 0.05
# acceleration due to gravity(m/s^2)
g = 9.81
# length of string (m)
l = 1
# mass of bob (kg)
m = 1
# Initial Conditions
theta_0 = [0,3]
# time points
t_span = np.linspace(0,20,100)
# Solve ODE
theta = odeint(ode_function,theta_0,t_span,args = (b,g,l,m))
# Plotting the graph
plt.plot(t_span,theta[:,0],'b-',label=r'$frac{dtheta1}{dt}=theta2$')
plt.plot(t_span,theta[:,1],'r--',label=r'$frac{dtheta2}{dt}=-frac{b}{m}theta2-frac{g}{l}sintheta1$')
plt.xlabel('Time (second)')
plt.ylabel('Plot')
plt.title('Damped Harmonic motion of a Simple Pendulum')
plt.legend(loc='best')
plt.grid()
plt.show()
# Iteration using for loop to create animation of motion of pendulum
loop_counter =1
tmp = theta[:,0]
for i in tmp:
x_start = 0
y_start = 0
x_end = l*math.sin(i)
y_end = -l*math.cos(i)
# Plotting for animation
plt.figure()
# Plotting String of pendulum
plt.plot([x_start,x_end],[y_start,y_end],linewidth=3)
# Plotting Bob of pendulum
plt.plot(x_end,y_end,'o',markersize=20,color='r')
# Plotting the rigid support through which string of pendulum is attached
plt.plot([-0.5,0.5],[0,0],linewidth=8,color='k')
# Plotting hatched vertical line for refernce
plt.plot([0,0],[0,-1.5],'--',color='y')
plt.xlim([-1.25,1.25])
plt.ylim([-2,1])
plt.title('Animation of transient behaviour of a Simple Pendulum')
filename = 'test%05d.png' % loop_counter
plt.savefig(filename)
loop_counter = loop_counter + 1