Week 3 - Solving second order ODEs
AIM
To simulate the transient behaviour of a simple pendulum
INTRODUCTION
A pendulum is a body suspended from a fixed support so that it swings freely back and forth under the influence of gravity.
A simple pendulum consists of a mass m hanging from a string of length L and fixed at a pivot point . When displaced to an initial angle and released, the pendulum will swing back and forth with periodic motion.
The simple harmonic motion of a pendulum shown
d ^2 θ /dt^2 + (b/m)(dθ/dt) + (g/ L) sin θ = 0
second order ordinary differential equation - position w.r.t time
where g= gravity (m/s^2)
l=length of pendulum(m)
m=mass of the ball in (kg)
b=damping coefficient
Forces acting on the pendulum can be resolved into sin θ and cos θ
PROCEDURE
- Intially we are importing math and mathplotlib,numpy( for linspace() command),odeint.
- Matplotlib is a plotting library for the Python programming language
- In the function program second order differential equation is converted to first order equation
- dtheta_dt valuesare stored in the 2 columns (angular displacement and angular velocity)
- theta_0=[0;3]; is a 2D array initial condition for initial displacement and initial velocity in 2D array.
- Then we are initializing b, g , l , m as inputs.
- theta_0=[0;3]; is a 2D array initial condition for initial displacement and initial velocity in 2D array.
- t is the vector defining the beginning and end limits of integration .
- Differential equations are solved in Python with the Scipy.integrate using function ODEINT
- Plot command used to plot the links in the graphs ,here we are plotting between (angular displacement', 'angular velocity’) VS time
- plt.show() is used to show the plot.
- In the for loop we are defining the points (x0,y0) (x1,y1).
- Then are plotting the coordinates
- Finally using cygwin and image magick ,png file is converted to mpg file ( magick %d.png[1-150] movie.mpg )
PROGRAM
import numpy as np
from scipy.integrate import odeint
import matplotlib.pyplot as plt
import math
#function that returns dz/dt
def model(theta,t,b,q,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 #b=damping coefficient
g=9.81 #g= gravity (m/s^2)
l=1 #l=length of pendulum(m)
m=1 #m=mass of the ball in (kg)
#initial coditions
theta_0=[0,3] #initial conditions for initial displacement and initial velocity in 2D array.
#time plots
t=np.linspace(0,20,150) #we are integrating from t = 0 to t = 20, and want to take 150 time steps,
#solve ode
theta=odeint(model,theta_0,t,args=(b,g,l,m))
#plot results
plt.figure(1)
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.xlabel('time')
plt.ylabel('plot')
plt.legend(loc='best')
plt.show()
ct=1
for i in theta[:,0]:
x0=0
y0=0
x1= x0+l*math.sin(i)
y1= y0-l*math.cos(i)
filename=str(ct)+'.png'
ct=ct+1
plt.figure(2)
plt.plot([-0.5, 0.5],[0,0],linewidth='3',color='r')
plt.plot(x1,y1,'o',markeredgecolor='b',markerfacecolor='b',markersize='15')
plt.plot([x0 ,x1],[y0 ,y1],linewidth='5',color='g')
plt.xlim([-2,2])
plt.ylim([-2,2])
plt.title('simple pendulum')
plt.savefig(filename)
plt.close()
RESULTS
Figure(1) graph between ((angular displacement and angular velocity) vs time)
FIGURE(2)
ANIMATION OF SIMPLE PENDULUM MOTION
ERROR
Initially i couldnt able to get the perfect video for simple pendulum