Week 3 - Solving second order ODEs

AIM: To write a Python program to simulate the transient behaviour of a simple pendulum and to create an animation of its motion.

GOVERNING EQUATION:

The equation of motion of a simple pendulum with damping is-

`(d^2theta)/(dt^2) + b/m (d theta)/(dt) + g/Lsin(theta) = 0`

where,

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

L = length of the pendulum in `m`

m = mass of the ball in `kg`

b = damping coefficient in `(Ns)/m`

OBJECTIVES:

1. To simulate the transient behaviour at a simple pendulum.

2. To plot the graph of angular displacement `theta` and angular velocity `v` against time `t`.

3. To create an animation of its motion.

PROCEDURE FOR SOLVING ODE:

The 2nd order differential equation is converted into two 1st order differential equations.

`theta_1 = theta`

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

`(d theta_2)/dt = (d^2theta_1)/dt^2 = (d^2theta)/dt^2 = -b/m theta_2 - g/Lsin(theta_1)`

GIVEN:

g = 9.81 `m/s^2`

L = 1 `m`

m = 1 `kg`

b = 0.05 `(Ns)/m`

`theta_ 0` = 0 radian

`v_0` = 3 `(rad)/s`

`t_0 = 0` sec

`t_(end) = 20` sec

PYTHON CODE:

# Simple pendulum
# Importing libraries
import numpy as np 
from scipy.integrate import odeint
import math
import matplotlib.pyplot as plt

# function 
def model(theta,t,b,g,l,m):
	theta1 = theta[0]
	theta2 = theta[1]
	dtheta1 = theta2
	dtheta2 = -(b/m)*theta2 - (g/l)*math.sin(theta1)
	dtheta = [dtheta1,dtheta2]
	return dtheta


# inputs

b = 0.05
g = 9.81
l = 1
m = 1

# Initial conditions
theta0 = [0, 3]

# Time points
t = np.linspace(0,10,200)

# solve ODE
theta = odeint(model,theta0,t, args = (b,g,l,m))

# Plotting

plt.figure(1)
plt.plot(t,theta[:,0],'b-', label=r'$\frac{d\theta_1}{dt} = \theta_2$')
plt.plot(t,theta[:,1],'r--', label=r'$\frac{d\theta_2}{dt} = -\frac{b}{m}\theta_2 - \frac{g}{L}sin\theta_1$')
plt.legend(loc='best')

# Axis Label
plt.xlabel('Time(s)')
plt.ylabel('Values')

# Title of the Plot
plt.title('Angular Displacement and Angular Velocity v/s Time')
plt.grid()
plt.show()

# For animation
displacement = theta[:,0]

# Defining counter variable
ct = 1

plt.figure()
for theta in displacement:

	filename = str(ct) + '.png'
	ct = ct + 1
	
	plt.figure()
	x = l*math.sin(theta)
	y = -l*math.cos(theta)
	plt.plot([0,x],[0,y],'r-',lw=3)
	plt.plot(x,y,'ko', markersize=20)

	# Axis Limit
	plt.xlim([-2,2])
	plt.ylim([-2,0.5])
		
	# Axis Label
	plt.xlabel('x-axis')
	plt.ylabel('y-axis')

	# Title of the Plot
	plt.title('Simple Pendulum')

	# Saving Plots
	plt.savefig(filename)

PROCEDURE FOR CREATING ANIMATION:

I have stored all figures in a folder. After that, I have used imgflip(https://imgflip.com/gif-maker) to generate the GIF.

GRAPH:

OUTPUT FOLDER:

ANIMATION: