Solving second order ODEs using Python

OBJECTIVE:

To write a program that solves the ODE given below and simulate the results obtained to create an animation.

where θ = Angular displacement of the ball (in radians)

          t = time

          b = damping coefficient

          m = mass of the ball (in kg)

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

          L = Length of string  (in metres)

PROBLEM STATEMENT: 

An ODE is used to describe the transient behavior of a system. In mechanics and physics, the simple harmonic motion is a special type of periodic motion where the restoring force on the moving object is directly proportional to, and opposite of, the object's displacement vector.

In real time scenarios the air resistance slows down the velocity of the simple pendulum and for the same reason damping factor is introduced in the equation given below. Hence it is called as damped simple harmonic motion. (the oscillations are damped). From Newton's second law, an ODE is obtained for the motion of simple pendulum as shown in the figure below. The angular displacement and angular velocity are represented in the equation.

where θ = Angular displacement of the ball (in radians)

          t = time

          b = damping coefficient

          m = mass of the ball (in kg)

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

          L = Length of string  (in metres)

The problem statement is defined by the following:

1. To solve the above second order differential equation using python. This is done by breaking the second order ODE to two different first order ODE and then using "odeint" function present in scipy module of python.

2. Plotting the various angular displacements and angular velocities with respect to time after solving the differential equation.

3. To simulate the pendulum between 0-20 seconds for angular displacement = 0 and angular velocity = 3 rad/s at time t=0. This is done by creating images at different time intervals and converting it to gif image to get the desired simulation. The positions of the string are defined as shown in the figure below.

Fixed point of the string = (0,0)

Free end of the string = (l*sinθ, -l*cosθ) where l = length of the string (in metres)

SOLUTION PROCEDURE:

import matplotlib.pyplot as plt
from scipy.integrate import odeint 
import numpy as np 
import math
import os,shutil
import imageio as io

#To get the current working directory
parent_dir=os.getcwd()
#To create a directory
directory = "Pendulum Animation"

#Creating a folder in the defined working directory
path = os.path.join(parent_dir, directory) 

#Checking if path already exists
if not os.path.exists(path):
	os.mkdir(path) 

#If path exists then to delete all the files in the path
else:
	files_in_directory = os.listdir(path)
	#For getting the files in the directory
	filtered_files = [file for file in files_in_directory]
	for file in filtered_files:
		path_to_file = os.path.join(directory, file)
		#To delete the files present in the directory
		os.remove(path_to_file)

#Defining function for the ODE
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

#Damping coefficient
b=0.05
#Length of the pendulum (in metres)
l=1
#Acceleration due to gravity (in m/s^2)
g=9.81
#Mass the of bob (in kg)
m=1

#Initial Condititons
theta_0=[0,3]

#Line points for time interval 0-20 seconds
t = np.linspace(0,20,200)

#Solving the defined ODE
theta= odeint(model,theta_0,t, args = (b,g,l,m))

#Plotting the results
plt.plot(t,theta[:,0],'b-',label=r'$frac{dtheta_1}{dt}=theta2$,Angular Displacement')
plt.plot(t,theta[:,1],'r--',label=r'$frac{dtheta_2}{dt}=-frac{b}{m}theta2-frac{g}{L}sintheta_1$, Angular Velocity')
plt.ylabel('Angular Displacement, Angular Velocity')
plt.xlabel('Time')
plt.legend(loc='best')
plt.show()


#Defining empty arrays for storing the results
Angular_Displacement=[]
Angular_Velocity=[]

#Creating images array to store the images
images=[]

#Initialising counter variable for storing images
ct=1

#Creating array to store the different coordinate values
prev_x1=[]
prev_y1=[]
#Storing the values in the respective arrays
for i in range(0,len(theta)):
	displacement=theta[i,0]
	Angular_Displacement.append(displacement)
	velocity=theta[i,1]
	Angular_Velocity.append(velocity)

	#Coordinates for fixed end of pendulum
	x=0
	y=0

	#Defining the position and velocity 
	Position=Angular_Displacement[i]	
	Velocity=Angular_Velocity[i]
	
	
	#Coordinates for free end of pendulum
	x1=l*(math.sin(Position))
	y1=-l*(math.cos(Position))
	
	#Storing the coordinates in the defined array
	prev_x1.append(x1)
	prev_y1.append(y1)

for i in range(0,len(prev_x1)-1):


	filename=str(ct)+'.png'
	#incrementing the counter to get different file name for images
	ct=ct+1;

	#Coordinates for support of the pendulum
	plt.plot([-1,1],[0,0],'--',color='gray',linewidth=3.5)

	#Creating pivot point for the string of the pendulum
	plt.plot(x,y,'.',color='black')

	#Coordinates for defining the string of the pendulum
	plt.plot([x-0.015,prev_x1[i]-0.015],[y,prev_y1[i]],color = 'k',linewidth=0.5)
	plt.plot([x+0.015,prev_x1[i]+0.015],[y,prev_y1[i]],color = 'k',linewidth=0.5)

	
	#Defining the marker at the end of pendulum
	plt.plot(prev_x1[i],prev_y1[i],marker = 'o',markerfacecolor='m',markeredgecolor='k',markersize=15)

	#For marking the trajectory of the marker
	plt.plot([prev_x1[i-1],prev_x1[i],prev_x1[i+1]],[prev_y1[i-1],prev_y1[i],prev_y1[i+1]],'b--',linewidth=1.5)


	#Defining title for the animation
	plt.title("Motion of Simple Pendulum",fontsize=18)

	plt.grid()	#Coordinates for defining the string of the pendulum


	#Defining the grid by giving axes values
	plt.axis([-1.5,1.5,-1.5,0.5])

	#saving all the figures				
	plt.savefig(filename)	
	#appending the images to the defined array
	images.append(io.imread(filename))	

	#moving the array of images to the desired directory
	shutil.move(f'{parent_dir}/{filename}',f'{path}')	

	#Clearing the window screen after a run
	plt.clf()

#Converting images to gif and saving it in the desired directory
io.mimsave(f'{parent_dir}/Pendulum Animation/Animation.gif',images)	

Points to be Noted: 

1. The libraries import os,shutil and import imageio are used to create and work on files in directory and to create gif animation using the images.

2. A separate for loop is create to trace the trajectory of bob.

3. Also, in animation t=20 second is not visible since the loop will run only for the value upto t=19 second to avoid the indexing error as given in point 3 of error encountered heading.

4. The graph obtained after solving the ODE must be closed to get the animation file since the animation code is given next to the ODE solving block.

ERRORS ENCOUNTERED:

1. After running the code, the following image is obtained (for 20th second). Instead of capturing frame by frame, all the time intervals are captured in a single image. This is not a main error but it should be avoided since the animation captures these images.

The error is avoided by including the following line in the code. The plt.clf() clears the window everytime the image is plotted in the graph for every time interval.

The image obtained in 200th second is shown below.

2. The following error is obtained while running the code as given in figure below.

This error is obtained because the first plt.show() gets executed and plt window is opened. Only after closing that window the next set of code will executed for saving the images. Again, if we try to run the code the error is shown because the plot window is still opened and during last run some images get stored in the folder.

To avoid this error, the plot window should be closed following which the code executes perfectly as shown below.

3. The following error is obtained as shown in the figure below.

This error occured because for tracing the trajectory of the pendulum we have defined the values in the range of (0, length of results). For the iteration i+1 the value is not defined since it exceeds the defined list. 

To avoid this error, the loop is defined from (0, length of results-1) as shown in the figure below.

RESULTS AND DISCUSSION:

Therefore, a python code is created for solving the given second order ordinary differential equation for describing the motion of simple pendulum.

The Angular Displacement, Angular Velocity vs Time plot is obtained as shown in the figure below.

The animation for motion of pendulum is shown in figure below.

Youtube Link: https://youtu.be/pvUbRjPqbE4