Solving second order ODE's using PYTHON Programming

AIM:

1. To write a program in PYTHON that solves the Second order ODE corresponding to motion of a simple pendulum and plot its Angular Displacement and Angular Velocity.
   
   
2. Create an animation of the motion of pendulum by simulating the motion between 0-20sec with angular displacement = 0, angular velocity = 3rad/sec at time t = 0.

THEORY:

Simple Pendulum : 

A point mass attached to a light inextensible string and suspended from fixed support is called a simple pendulum. The vertical line passing through the fixed support is the mean position of a simple pendulum. 

Pendulum is a weight suspended from a pivot so that it can swing freely. When a pendulum is displaced sideways from its resting, equilibrium position, it is subject to a restoring force due to gravity that will accelerate it back toward the equilibrium position. When released, the restoring force acting on the pendulum's mass causes it to oscillate about the equilibrium position, swinging back and forth. The time for one complete cycle, a left swing and a right swing, is called the period. The period depends on the length of the pendulum and also to a slight degree on the amplitude, the width of the pendulum's swing.

PYTHON PROGRAM:

#Program to solve second order ODE's and simulate its pendulum motion

import matplotlib.pyplot as plt 	#Library for plotting functions
import math							#Library for math functions
from scipy.integrate import odeint	#For 'odeint' function
import numpy as np					#For 'linspace' function


#Function Definition
def model(theta,t,b,g,l,m):
		theta1=theta[0] #ODE 1
		theta2=theta[1]	#ODE 2
		dtheta1_dt=theta2
		dtheta2_dt=-(b/m)*theta2-(g/l)*math.sin(theta1)
		dtheta_dt=[dtheta1_dt,dtheta2_dt]
		return dtheta_dt


#Input Parameters
b=0.05	#Damping Coefficient
g=9.81	#Acceleration due to Gravity
l=1 	#Length of Pendulum
m=1 	#Mass of Pendulum

#Initial Condition
theta_0=[0,3] #Angular Displacement and Velocity
#Time point
t=np.linspace(0,20,300) #20 sec motion

#Solving ODE using 'odeint' function
theta=odeint(model,theta_0,t,args=(b,g,l,m))

#Plotting points
plt.plot(t,theta[:,0],'r-',label=r'$frac{dtheta_1}{dt}=theta_2$')
plt.plot(t,theta[:,1],'b--',label=r'$frac{dtheta_2}{dt}=-frac{b}{m}.theta_2-frac{g}{l}.sin(theta_1$)')
plt.xlabel('Displacement')
plt.ylabel('Time')
plt.title('Damped Harmonic motion')
plt.legend(loc='best')
plt.show()

#Animation of Simple pendulum
#Initial Position
xstart=0 
ystart=0
ct=1

for i in theta[:,0]: #File
		xend=l*math.sin(i)
		yend=-l*math.cos(i)

		plt.figure()

		#Base of Pendulum
		plt.plot([-1,1],[0,0],'black',linewidth=8)

		#String of Pendulum
		plt.plot([xstart,xend],[ystart,yend],'black')

		#Bob of Pendulum
		plt.plot(xend,yend,'o',markersize=20,color='r')

		#Graph Limits of x nad y axis
		plt.xlim(-1.5,1.5)
		plt.ylim(-1.5,1.5)

		#Saving the Frames
		filename=str(ct)+'.png'
		plt.savefig(filename)
		ct=ct+1

RESULT:

Damped Harmonic Motion:

Simple Pendulum Oscillation :