Week 3 - Solving second order ODEs

Aim: A program to solve the equation of 2nd order differential equation

The equation for damping is as follow :

where In the above equation,

g = gravity in m/s2,

L = length of the pendulum in m,

m = mass of the ball in kg,

b=damping coefficient.

and 

L=1 meter,

m=1 kg,

b=0.05.

g=9.81 m/s2.

Code explanation:

• create a function to calculate the dtheta_dt which contains the values of dtheta1_dt and  dtheta2_dt
• provide input of  g,length ,mass ,b
• create the angle and time values for the pendulum
• now using the odeint syntex we solve the equation
• plot the values

Code:

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

#function for crating the
def model(theta,t,b,g,m,l):
	theta1 = theta[0]
	theta2 = theta[1]
	dtheta1_dt = theta2
	dtheta2_dt = (-b/m)*theta2-(g/l)*math.sin(math.radians(theta1))
	dtheta_dt = [dtheta1_dt,dtheta2_dt]
	return  dtheta_dt

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

theta_0 = [0,60]
t = np.linspace(0,80,500) #time span

#solving ode
theta = odeint(model,theta_0,t,args = (b,g,l,m))

#plotting
plt.plot(t,theta[:,0],'g',label=r'Displacement')
plt.plot(t,theta[:,1],'b',label=r'Velocity')
plt.xlabel('time')
plt.ylabel('Displacement')
plt.savefig('plot.png')
plt.show()


#creating pngs for animation
ct = 1
for t in theta[:,1]:
	x0 = 0
	y0 = 0
	x1 = (l*math.sin(math.radians(t)))
	y1 = -(l*math.cos(math.radians(t)))

	image = str(ct)+'.png'
	ct = ct+1

	#plotting
	plt.figure()
	plt.plot([-1,1],[0,0],color='red',linewidth=4)
	plt.plot([x0,x1],[y0,y1])
	plt.plot(x1,y1,'o',markersize=15)
	plt.xlim([-1.5,1.5])
	plt.ylim([-1.5,1.5])
	plt.savefig(image)

Results:

Animation video : 

Conclusion:

Successfully created the graph of displacement vs time by solving the second-order differential equation.