Week 3 - Solving second order ODEs

AIM: WAP to simulate the transient behaviour of a simple pendulum and to create an animation of it's motion.

OBJECTIVE:

1. Create an animation of the pendulum motion

2. Simulate the transient behaviour of a simple pendulum

GOVERNING EQUATION:

Ordinary differential equation for oscillating motion of pendulum is,

`(d^2theta)/dt^2 +b/m +g/l*sintheta=0`

b=damping coeffecient

l=length of pendulum in metre

m=mass of ball in kilogram

g=gravitational acceleration in `m/s^2`

INPUTS:

time span=0-20 sec

angular displacement=0 rad

angular velocity=3 rad/sec

l=1 metre

m=1 kg

b=0.005

g=9.81 `m/s^2`

PYTHON PROGRAM:

To solve the pendulum equation, we need to break the second order equation into two first order ode equations. This is done by following the below mentioned steps.

The above mentioned equations are formed in the below mentioned python code. [θ1;θ2][θ1;θ2] is considered as theta vector.

PYTHON CODE:

#different modules like math, scipy, numpy and matplotlib are imported 
import math
from scipy.integrate import odeint
import matplotlib.pyplot as plt
import numpy

def model(theta,t,b,g,l,m): #function definition
	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

#inputs
b=0.05
g=9.81
l=1
m=1
t=numpy.linspace(0,20,150) #an array of time intervals are formed
theta_0=[0,3] #initial displacement and velocity
theta=odeint(model,theta_0,t,args=(b,g,l,m)) #ode solver is used to solve the equation by providing the inputs, time span and differential equations

#first figure is plotted
plt.figure(1)
plt.plot(t,theta[:,0],'r-',label='Displacement')
plt.plot(t,theta[:,1],'b--',label='Velocity')
plt.xlabel('Time')
plt.xlabel('Velocity and Displacement')
plt.legend()
plt.show()

#hinged point of pendulum
x0=0
y0=0
ct = 0 #counter
for i in theta[:,0]: #for loop is initiated for knowing the dynamic position of the pendulum at each and every instant
        #Determining the place of the pendulum at each time step
	x1 = (1*math.sin(i))
	y1 = (-1*math.cos(i))
        
        #Plotting the pendulum and saving at each time step
	filename = str(ct) + '.png'
	ct = ct+1	
	plt.figure()
	plt.plot([-1,1],[0,0], 'r', linewidth=3)
	plt.plot([x0, x1], [y0, y1])
	plt.plot(x1,y1, 'o', markersize = 20)
	plt.xlim([-1.5, 1.5])
	plt.ylim([-1.5, 0.5])
	
	plt.savefig(filename)

In the code all the inputs are provided like initial displacement, time span, initial velocity. The odeint solver is used to solve the pendulum equation and the results are stored in results vector. It has angular displacement and angular velocity solutions at different time intervals and a graph is plotted. Also, the pendulum motion is simulated.

Animation of the Pendulum:

I have created the animation by capturing the frames and then combining them into GIF using imgflip website.