Week 3 - Solving second order ODEs

Objective: To solve the second-order ODE equation and create an animated gif video of a simple pendulum.

Solution:

The second order ode represents the motion of the simple pendulum and the equation is as follows:

`(d^2 theta)/dt + b/m (dtheta)/dt + g/l sin theta = 0` ....... (1)

Python uses scipy integration to solve ODE but it solves only simple single-order ODEs

Equation (1) can be modified as follows:

let `theta = theta_1` and

`(d theta_1)/dt = theta_2` ------- (2)

`(d theta_2)/dt = -b/m theta_2 - g/l sin(theta_1)` ---- (3)

Procedure:

import math to solve mathematical functions

import numpy to create and use arrays elements

import odeint from scipy.integration to solve ODE

import os to get the folder details and save figures

import matplotlib.pyplot to create the graph figures

Create a function model with equations (2) and (3) to solve eqn(1) and get the values for theta_1 and theta_2 which are angular displacements and angular velocities.

Pass arguments

b bulk coffeficient = 0.01

g acc. due to gravity = 9.8 m/s^2

l length = 1 m

m mass of bob = 0.1 kg

After plotting and saving the plot figures in os,

the figures are created as gif video using imageio package, which is must be installed to access.

The code is as follows in detail:

# Code to solve the motion of simple pendulum equation and create an animated gif file
import numpy as np
import math
from scipy.integrate import odeint
import matplotlib.pyplot as plt
import imageio
import os  

# function that returns dtheta/dt
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
#parameters
b= 0.02
g= 9.81
l= 1
m= 0.1

#initial condition
theta_0 = (0,3)

# time points

t = np.linspace(0, 20, 150)

#solve ode
theta = odeint(model, theta_0, t, args =(b, g, l, m))


#plot results
plt.plot(t, theta[:, 0], 'b-', label=r'$\frac{d\theta_1}{dt} = \theta2$')
plt.plot(t, theta[:,1], 'r--', label=r'$\frac{d\theta_2}{dt}= -\frac{b}{m}\theta_2-\frac{g}{t}sin\theta_1$')
plt.ylabel('plot')
plt.xlabel('time')
plt.legend(loc='best')
plt.show()


ct=1     
filenames=[]                     
for Theta in theta[:,0]:
	
	#Co-ordinates of Fixed Point of String
	x0 = 0
	y0 = 0
	#Co-ordinates of End Point of String
	x1 = l*math.sin(Theta)
	y1 = -l*math.cos(Theta)

	filename = str(ct) + '.png'
	filenames.append(filename)
	ct = ct + 1

	plt.figure()
	# Plotting Origin, string and ball of pendulum
	plt.plot([-0.5,0.5],[0,0],'g',linewidth=5)
	plt.plot([x0,x1],[y0,y1],'r')
	plt.plot(x1,y1,'o',markersize=25)
	#Limiting X-axis and Y-axis
	plt.xlim([-1.5,1.5])
	plt.ylim([-1.5,1.5])
	#Title
	plt.title('Animation of Simple Pendulum')
	#To save the file name
	plt.savefig(filename)
  
#imageio module to generate animated gif video using saved files
with imageio.get_writer('Pendulum.mp4', mode='I', fps=10) as writer:
    for filename in filenames:
        image = imageio.imread(filename)
        writer.append_data(image)
for filename in set(filenames):
    os.remove(filename)		#removes the png files saved

Conclusion:

The simple pendulum motion equation, a second-order ODE has been solved and a animated video describing motion of simple pendulum with angular velocity 3 rad/sec has been created using python.

The youtube video links:

https://youtu.be/M6zi_dNPpDc