SOLVING SECOND ORDER ORDINARY DIFFERENTIAL EQUATION USING PYTHON

AIM

To solve the ODE of the simple pendulum and to plot the curve of various parameter by using PYTHON coding.

1. Time vs angular displacement
2. Time vs angular velocity

To simulate the motion of a pendulum for time t= 0 to 20 second.

To make the animation of the motion of pendulum.

GOVERNING EQUATIONS

`(d^2 ϴ)/(dt^2 ) + b/m (dϴ)/dt + g/l sin⁡ϴ= 0`

In the above equation,

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.

Given data

g = 9.81 `m/s^2`

L = 1 m,

m = 1 kg,

b = 0.05,

Angular displacement = 0 m and angular velocity = 3 rad/s at time t= 0 second.

BODY OF THE CONTENT

A simple pendulum consists of a mass attached at the end of the string and the other end of the string is fixed as shown in fig

The very first step to solve the above ODE is to convert the second order differential equation into first order differential equation

theta=theta1
theta1= dtheta/dt

dtheta1/dt = theta2

dtheta2/dt = -(b/m)*theta2 - (g/l)*sin(theta1)

The program is completed into two parts. First part which is function to solve the ordinary differential equation and second part is for plotting and making the animation of the simple pendulum. Both two part are written in the same program.

STEPS OF PRAGRAM

Before starting program make sure that scipy and numpy library is installed on your laptop.

1. Provide topic of the code which is under comment

2. Calling of library which will used for solving ODE, plotting and calculating i.e. math, matplotlib.pyplot, scipy and numpy library.

3. Making of function from line 11 to 26.

In line 12, Function takes the value of ϴ, time, damping coefficient, acceleration due to gravity, length, mass.

Line 15 and 16 is defining of ϴ1 and ϴ2.

In line 19 and 20 governing equation is written

In line 23 combined equation is written for governing equation.

In line 26, return is used to give the values to the function.

4. From line 29 to 39 input of the simple pendulum is provided.

In line 30 to 33 input parameter of pendulum where b, g, l and m are damping coefficient, acceleration due to gravity, length of the string, mass of the bob respectively.

In line 36 theta_0 include input parameter of angular displacement and angular velocity.

In line 39, time is from 0 to 20 second which is divided in 200 equal parts.

5. Line 42 is used to solve Differential equation which take an input of model, theta_0, t and args =(b, g, l, m) are function, angular displacement: angular velocity, time and pendulum parameter.

6. Line 45 and 46 is used to separate the angular displacement and angular velocity which is obtain after solving differential equation.

7. From line 50 to 62 having a code for plotting of the curve.

Line 51 is used to define figure.

Line 53 is used to plot the curve between time and angular displacement.

In line 55 is used to plot the curve between time and angular velocity.

Whereas at the end of line 53 and 55 defining of legend command is given.

In line 57 and 58 labelling of X and Y axis of the plot.

In line 60, legend command is used to display legend on the plot.

In line 62, show command is used to display the plot.

8. From line 65 to 96 code used for getting images of angular displacement and time for making animation of the pendulum arrangement.

In line 66 ct is used for count start time.

In line 68 for loop is used to run the program for different values of angular displacement.

9. From line 70 to 76, command is used to define the coordinate i.e. x0, y0, and x1, y1.

10. Line 79 is used to define the name of the image which to be save for making of animation where animation of pendulum is the name of the image and %03d is used to give numbering of the image in 3 digit.

11. From line 81 to 90 command is used to plot the image of the pendulum arrangement

In line 82, defining of figure is done.

In line 84, code is written to plot the fixed wall having a linewidth of 5.

In line 86, code is written for plotting of x0,y0 and x1,y1.

In line 88, code is written for making bob of the pendulum having shape is circle and size is 5.

In line 90, code is used to limit the X and Y axis.

12. Line 93 is used to save the image which is stored in line 79 named as filename.

Line 96 is used to increase the count time.

PYTHON PROGRAM CODE

"""
ORDINARY DIFFERENTIAL EQUATION TO SOLVE FOR SIMPLE PENDULUM
"""

#importing of libraries such as numpy, scipy.integrate, matplotlib.pyplot, math
import numpy as np 
from scipy.integrate import odeint
import matplotlib.pyplot as plt
import math

#Function that return dz/dt
def model(theta,t,b,g,l,m):

	#defining theta 1 and theta 2
	theta1 = theta[0]
	theta2 = theta[1]

	#defining the dtheta1/dt and dtheta2/dt
	dtheta1_dt = theta2
	dtheta2_dt = -(b/m)*theta2 - (g/l)*math.sin(theta1)
	
	#defining dtheta/dt
	dtheta_dt = [dtheta1_dt,dtheta2_dt]

	#giving value of dtheta/dt to function
	return dtheta_dt

	
#Input of the pendulum arrangement
b = 0.05						#Damping coefficient
g = 9.81						#Acceleration due to gravity
l = 1 							#Length of the string
m = 1							#Mass of the pendulum

#initial condition [angular_displacement,angular_velocity]
theta_0 = [0,3]

#Time points starts from 0 to 10 second
t = np.linspace(0,20,200)

#Solve ODE where model is function,theta_0 is initial input, t is time and args is other argument
theta = odeint(model,theta_0,t,args = (b,g,l,m))

#separating of angular displacement
angular_displacement= theta[:,0]
#separating of angular velocity
angular_velocity = theta[:,1]


#plot result which is stored in figure 1
plt.figure(1)
#olot of t vs angular displacement
plt.plot(t,theta[:,0],'b-',label=r'$\frac{d\theta_1}{dt}=\theta2$')
#plot of t vs angular 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$')
#labelling of plot
plt.ylabel('plot')
plt.xlabel('time')
#giving legend to plot
plt.legend(loc='best')
#command to show the plot
plt.show()


#ct is the count time
ct=1
#for loop fo making animation
for i in angular_displacement:

	#initial condition of string
	x0 = 1
	y0 = 1

	#final condition of the string
	x1=  l-l*math.sin(i)
	y1 = l-l*math.cos(i)

	#file name of the image to be saved in the directory 
	filename = 'animation of pendulam %03d.png' %ct
		
	#plotting of the figure
	plt.figure()
	#plot of the wall having linewidth of 5
	plt.plot([0.5,1.5],[1,1],linewidth=5)
	#ploting of the string
	plt.plot([x0,x1],[y0,y1])
	#plot of the pendulum
	plt.plot(x1,y1,'o',markersize=15)
	#defining the axis of the plot
	plt.axis([-0.5,2.5,-0.5,1.5])
	
	#saveing of the image in the directory
	plt.savefig(filename)

	#increament of the count time
	ct= ct+1

ERROR OBTAIN WHILE PROGRAMMING

1. Because of giving number in line 82 i.e. figure(2) the images are not getting correct as shown in figure of the directory of the computer

Instead of figure(2) write figure() in order to get correct image of the pendulum.

CONCLUSION

The plot obtains after solving the ODE is

1. Time vs angular displacement and Time vs angular velocity

2. Image obtain for the animation of the motion of pendulum.

YOUTUBE LINK