Week 3 - Solving second order ODEs

OBJECTIVE:

To write a python program solve the second-order Ordinary differential equation for simple Pendulum and to simulate oscillation of a simple pendulum and to create an animation of its motion.

PROBLEM STATEMENT:

To write a program that solves the following Ordinary Differential Equation which represents the equation of motion of a simple pendulum with damping.

`(d^2θ)/dt^2+b/m(dθ)/dt+g/Lsintheta=0`

where,

θ = angular displacement in radians

b = Damping Coefficient

L = Length of Pendulum in m

m = mass of the ball in kg

g = acceleration due to gravity in `m/s^2`

THEORY:

Simple Pendulum is one in which a point mass is suspended in a string or rod of negligible mass pivoted at the free end along which the pendulum oscillated to and fro about the equilibrium position due to the influence of gravity. The simple pendulum traces the Simple Harmonic Motion where the restoring is directly proportional to the displacement and acts in the direction opposite to that of displacement.

PYTHON PROGRAM:

#Solving Second Order ODE for Simple Pendulum

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

#Defining Function
def model(theta, t, b, g, l, m):

	theta1 = theta[0]
	theta2 = theta[1]
	dtheta1_dt = theta2
	dtheta2_dt = -(b/m)*theta2 - (g/1)*math.sin(theta1)
	dtheta_dt  = [dtheta1_dt, dtheta2_dt]
	return dtheta_dt

b=0.05
g=9.81  #m/s
l=1     #metre
m=1     #kg

#initial condition
theta_0 = [0, 3]

#Time Points
t = np.linspace(0,20,200)

#Solve ODE
theta = odeint(model, theta_0, t, args = (b,g,l,m))

#Plot Results
plt.plot(t, theta[:,0], 'b-', label=r'$\frac{d\theta_1}{dt}=\theta2$')
plt.plot(t, theta[:,1], 'r--', label=r'$\frac{d\theta2}{dt}=-\frac{b}{m}\theta_2-\frac{g}{L}sin\theta_1$')
plt.ylabel("Plot")
plt.xlabel("Time")
plt.legend(loc="best")
plt.title('Damped Oscillation of Simple Pendulum')
plt.show()

ct =1

for Theta in theta[:,0]:

	x0 = 0
	y0 = 0

	x1 = 1*math.sin(Theta)
	y1 =-1*math.cos(Theta)

	filename = 'frame %03d.png' %ct
	ct = ct+1

	plt.figure()
	#Plotting of Horizontal String Support
	plt.plot([-0.5,0.5],[0,0],'c',linewidth=4)
	#Plotting of Pendulum String
	plt.plot([x0,x1],[y0,y1],'k')
	#Plotting of Bob
	plt.plot(x1,y1,'o', markersize=25)
	#X-axis
	plt.xlim([-2,2])
	#Yaxis
	plt.ylim([-2,2])
	#Heading
	plt.title('Simple Pendulum simulation')
	plt.savefig(filename)
	

SOLUTION PROCEDURE:

• The above program is used to simulate the harmonic motion of a simple pendulum based on the Ordinary differential equation.
• Necessary libraries are imported into the program
• The function named 'model' is defined in which the Second-order ODE and its parameters given above are defined.
• The input data and initial conditions like angular displacement, angular velocity, time span as give are typed in the sublime text editor.
• The graph of angular displacement and angular velocity is obtained by using an inbuilt function 'odeint'. Then the graph is plotted by entering the time and theta values.
• Based on the time span and the number of frames between the given time span, a for loop is executed with the help of counter 'ct' command.
• For tracing the motion of a simple pendulum frames of the pendulum motion is saved in definite order by using the counter variable.
• The frames are stitched using Image Magick to create the animation video.
• Then the coordinates for placing the pendulum set up is typed according to our convenience following the above syntax.

RESULTS AND DISCUSSIONS:

The plot of Angular displacement and Angular velocity Vs time is obtained for damped oscillation of a Simple Pendulum by substituting the given values. Also, the simulation of the simple pendulum was created and uploaded in youtube.

YOUTUBE LINK FOR THE SIMULATION OF SIMPLE PENDULUM:

https://www.youtube.com/watch?v=YyhIQckK9lc&t=6s

THANK YOU