Week 3 - Solving second order ODEs

Aim:

Write a program to solve a differential equation of a simple pendulum. Plot the angular displacement and angular velocity with respect to time in the range of 20seconds. Animate the pendulum position accordingly.

Theory:

The transient behaviour of the simple pendulum with damping is given by the following equation.

`frac{d^2theta}{dt^2}+fracbm*frac{d theta}dt+fracgl*sintheta =0`

where,

b = Damping coefficient = 0.05

m = mass of the ball in Kg =1 Kg

g = gravity in m/s2 =9.81

l = length of the pendulum in m = 1m

Since the required outputs are angular displacement and angular velocities. We need to integrate the above equation to obtain the angular displacement, Say θ_t and angular velocity. Say ν_t

`theta_t = int(frac{d theta}dt)dt = inttheta_1dt`

Where θ₁ =`frac{d theta}dt`

`nu_t = int(frac{d^2theta}{dt^2})dt = intfrac{d theta_1}{dt}dt = int theta_2dt`

Where θ₂ =`frac{d^2theta}{dt^2}`

Substituting θ₁ and θ₂ in governing equation, we get

`theta_2=-b/m*(d theta)/(dt)-g/l*sintheta`

`theta_1 = (d theta)/(dt)`

Integrating above equations we get θ_t and ν_t

Given at t= 0sec,

angular displacement, `theta = 0`

angular velocity, `nu = 3` m/s

For the position of free end with respect to the fixed end,

`x =+l*sin theta`

`y = -l*costheta`

Code:

Assuming

`Delta theta=[theta_1,theta_2]` is an array of two differential quantities.

# Simple pendulum from integrating ordinay differential equations
# By manikanteshwara reddy

# Importing modules
import math
import numpy as np
import matplotlib.pyplot as plt
from scipy.integrate import odeint

# Defining a function for Delta theta.
def delta(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
# Input arguments
# Damping co-efficient
b = 0.05
# Acceleration due to gravity
g = 9.81
# Length
l = 1
# Mass
m = 1

# Starting condition of the pendulum i.e, Angle and Velocity
theta_0 = [0,3]
# Time inteval, taken for 20sec
t = np.linspace(0,20,200)

#Solution for deleta theta for the given initial conditions
theta = odeint(delta,theta_0,t,args = (b,g,l,m))
# print(theta[:,0]*(180/math.pi))
# Plotting delta theta with respect to time
plt.figure()
# Displacement
plt.plot(t,theta[:,0],'b-',label=r'$\frac{d\theta_1}{dt}=\theta_2$')
# Velocity
plt.plot(t,theta[:,1], 'r-',label=r'$\frac{d\theta_2}{dt}=-\frac{b}{m}\theta_2-\frac{g}{l}sin\theta_1$')
# For plotting the centre line (note:X and Y inputs are array of 200 values) >> plt.plot(t,np.linspace(0,0,200),'m-')
# OR a straight line with start and end points
plt.plot([0,20],[0,0],'k--')
plt.ylabel('plot')
plt.xlabel('time')
plt.legend(loc='best')
plt.savefig('Simple pendulum')

# Animation of simple pendilum
# Fixed end location (assumed w.r.t (0,0))
x0 = 1.5
y0 = 1.5
# Counter
ct = 1
# For loop to give seperate image for each location of the pendulum
for a in theta[:,0]:
	# Location of other end
	x1 = x0 + l*math.sin(a)
	y1 = y0 - l*math.cos(a)
	filename ='Simple pendulum%05d.png' %ct
	ct = ct +1
	plt.figure()
	# fixed end
	plt.plot(x0, y0,'r*')
	# Pendulum string
	plt.plot([x0,x1],[y0,y1],'g-')
	# Free end
	plt.plot(x1,y1,'o', markersize=10, markerfacecolor = 'cyan')
	plt.xlim([0,3])
	plt.ylim([0,3])
	plt.axis('off')
	plt.savefig(filename)

The above code is succesfully run as shown and finished in 62.2 seconds.

The outputs are the plots of displacement and velocity with respect to time. And from the ''for loop'', a series of images are obtained which are later stitched using image magick command.

Result:

Handling differential equations and understanding the use of arrays to arrange the values of different properties is learnt. We insert indexes into the given equations to divide the equation into fragments. The indexing method is found to be helpful to solve an array of two differential quantities while integrating them. Use of markers in the matplotlib module is learnt.