Week 3 - Solving second order ODEs
Aim:
To solve the second-order ODE function using Python which displays the motion of the simple pendulum and creates an animation of the motion.
Research and Objective:
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 the periodic motion. By applying Newton's second law for rotational systems, the equation of motion for the pendulum may be obtained
`(d^2 theta)/(dt^2) + g/L sin theta = 0`
Governing equations:
In this work a second order ordinary differential equation that describes the position of a ball with respect to time which has damping co-efficient
`(d^2 theta)/(dt^2) + b/m (d theta)/(dt)+ g/L sin theta = 0`
In the above equations,
g = gravity in `m/s^2`
L = length of the pendulum in m,
m = mass of the ball in kg,
b = damping co-efficient
given parameters
Length (L) = 1 m,
Mass (m) = 1 kg,
Damping co-efficient (b) = 0
Accelaration due to gravity (g) 9.81 `m/s^2`
For the time period = 0-20 sec
Angular dispalcement = 0
Angular velocity = 3 rad/sec at time t = 0
Body of the Objective:
import math
import matplotlib.pyplot as plt
import numpy asnp
from scipy.integrate import odeint # scinetifc equations
# function which represents the ODE system to define ordinary differential equations
def model(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
b = 0.05 # damping co-efficient
g = 9.81 # in m/s^2
l = 1 # length of the pendulum
m = 1 # mass in Kg
#initial condition for dispalcment, theta and velocity
theta_0 = [math.pi/3,0]
#time points time array linspace 10 sec
t = np.linspace(0,20,100) # numpy module 0 to 20 seconds 100 points
#solve ODE
theta = odeint(model,theta_0,t,args = (b,g,l,m)) # Odeint function that will intigrate a system of ordinary differenctal equations, arguments
#plot results
plt.plot(t,theta[:,0],'b-',label=r'$frac{dtheta_1}{dt}=theta2$') # column first ODE rep displacemnt values
plt.plot(t,theta[:,1],'r--',label=r'$frac{dtheta2}{dt}=-frac{b}{m}theta_2-frac{g}{L}sintheta_1$') # velocitiy
plt.ylabel('Plot')
plt.xlabel('time')
plt.legend(loc='best')
plt.show()
#Animation of simple pendulum at each angle and displacment
x = theta[0,1]
y = theta[:,0]
ct = 1
for angle in y:
x1 = 0; y1 = 0
x2 = l*math.sin(angle)
y2 = -l*math.cos(angle)
plt.figure()
plt.plot([-2.0,2.0],[0,0],'y',linewidth = 5)
plt.plot([x1,x2],[y1,y2],linewidth = 3)
plt.plot(x2,y2,'o',markersize = 20, markerfacecolor ='w')
plt.plot(x1,y1,'o',markersize = 10, markerfacecolor = 'r')
plt.xlim([-2.0,2.0])
plt.ylim([-2.0,0.5])
filename = str(ct) + '.png'
ct = ct +1
plt.savefig(filename)
from PIL import Image
import glob
#Create the frames
frames=[]
imgs=glob.glob("*.png") # to include all the .png files
import os
#natsort.natsorted(frames, reverse=True)
print(type(frames))
print(imgs)
print(frames)
image_folder = '/Users/Nagamani/Desktop/Skill-Lync/Py for Mech'
images = [image for image in os.listdir(image_folder) if image.endswith(".png")]
images = sorted(images, key=lambda string: int(string.split("_")[-1][:-4]))
print(images)
for i in images:
new_frame=Image.open(i) # open all the png files one by one
frames.append(new_frame) # append all the png files to an empty array 7 naming image object 'new_frame'
print(type(frames))
#save into a GIF file that loops forever
frames[0].save('Simple Pendulum.gif', format='GIF', append_images=frames[1:], save_all=True, duration=60, loop=0) # loop 0 means the loop goes on forever
# if you want the loop to go on for one time you will have to equate the loop to 1
Conclusions:
Please check the youtube link for the simple pendulum motion.