Week 3 - Solving second order ODEs

Aim: To write programme to solve ODE of simple Pendulum and create animation of motion of simple pendulum`

Equation used to solve Simple Pendulum motion:

(d^2 θ)/〖dt〗^2 +b/m*dθ/dt+g/l sin⁡θ=0`

`θ=θ_1`

`(dθ_1)/dt=θ_2`

`(d^2 θ)/(dt^2 )=(d^2 θ_1)/(dt^2 )=d/dt ((dθ_1)/dt)=d/dt (θ_2 )`

`d/dt [θ_2 ]=[-(b/m) θ_2-(g/L)sin⁡θ_1 ]`

Introduction:

A simple pendulum consists of a mass m hanging from a string of length L and fixed at a pivot point P. When displaced to an initial angle and released, the pendulum will swing back and forth with periodic motion. mass on a spring application, this model problem is representative of a wide variety of problems described by second order linear differential equations.

Objectives:

1. To solve second order ODE equation
2. To create animation of simple pendulum motion

Solution procedure:

from scipy.integrate import odeint 
import matplotlib.pyplot as plt 
import math
import numpy as np
# defining model 
def model(theta,t,b,g,l,m): 
	theta1=theta[0] # from ODE solution 
	theta2=theta[1] # from ODE solution 
	dtheta1_dt=theta2 # Frist ODE 
	dtheta2_dt=(-(b/m)*theta2-(g/l)*math.sin(theta1)) # Second ODE equation 
	dtheta_dt=(dtheta1_dt,dtheta2_dt) # combining two ODE's
	return dtheta_dt # returning dtheta_dt value 
b=0.05 # damping coeffient 
g=9.81 # gravity constant 
l=1 # length of string 
m=1 # mass of the bob
# initial condition 
theta_0=(0,3) # initial condition for displacement 
# Time point 
t=np.linspace(0,20,150) # To create time array
# sove ODE
theta=odeint(model,theta_0,t,args=(b,g,l,m)) # intigrating ODE 
# plot result
plt.plot(t,theta[:,0],'b-',label=r'$\frac{d\theta1}{dt}=\theta2$') # plot of displacement 
plt.plot(t,theta[:,1],'r--',label=r'$\frac{d\theta2}{dt}=-\frac{b}{m}\theta2-\frac{g}{l}sin\theta1$')# plot of acceleration 
plt.ylabel('plot') 
plt.xlabel('time')
plt.legend(loc='best')
plt.show()

# animation 
x0=0 # initial condition 
y0=0
ct=1 # intitial figure number
for i in theta[:,0]: # for loop for displacment 
	x1=l*math.sin(i) 
	y1=-l*math.cos(i)
	filename=str(ct)+"test%05d.png"%ct # saving figure in png format 
	ct=ct+1 # figure number increment 
# plotting 
	plt.figure() 
	plt.plot([-0.5,0.5],[0,0],linewidth=2,color='r') # support line 
	plt.plot([x0,x1],[y0,y1],linewidth=2,color='b')# length of pendulum 
	plt.plot(x1,y1,'o',color='g')# bob 
	plt.xlim([-1.5,1.5]) # graph x dimension 
	plt.ylim([-1.5,1.5]) # graph y dimension 
	plt.savefig(filename) 

Conclusion:

Due to damping effect amplitude and displacement reduces over the time