Week 3 - Solving second order ODEs

OBJECTIVE:

The core objective of the project is to simulate the simple harmonic motion of the simple pendulum by solving the second order ODE which defines the movement using Python. Once simulated, an equivalent animation of the same is created and checked for the correct movement.

PROBLEM STATEMENT:

The second order ODE that defines the SHM of the pendulum is as below:

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

where,

g = gravity in m/s2,

L = length of the pendulum in m,

m = mass of the ball in kg,

b = damping coefficient.

So, with this initial parameters, the ODE has to be solved for a time of 0 to 20 seconds with the initial condition at time = 0,

Angular displacement = 0 degrees

Angular velocity = 3 rad/s

PROCEDURE:

As a first step, the essential libraries needed for the problem are imported. Then the input parameters are specified as follows.

1. b = 0.05
2. m = 1 kg
3. g = 9.81 m/s^2
4. L = 1 m

Then the timesteps are created using the 'linspace' command from 0 to 20 seconds with a total number of 300 timesteps. Then the initial conditions at t = 0 is specified.

In order to solve the second order ODE, it is converted into two separate first order ODE's with the help of a user defined function. Then to solve them, an inbuilt function from the 'scipy' module is used. The function used is 'odeint' which integrates the system of two first order ODE's.

Once they are solved, then a plot showing the variation of angular displacement and angular velocity is created for checking the SHM behaviour. After that, a 'for' loop is created to calculate the position of the pendulum at different timesteps. With that position information, a series of plots is created.

Once all the plots are generated, they are stitched together to form the animation using a program called Imagemagick.

Program snippet:

# Simulating Simple pendulum motion by solving the ODE

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

# Defining the function to calculate the angular displacement and velocity at different time steps
def pendulum(theta_0,t,b,m,g,l):

	theta1 = theta_0[0]  # Angular displacement
	theta2 = theta_0[1]  # Angular velocity in rad/s

	dtheta1_dt = theta2
	dtheta2_dt = -(b/m)*theta2 - (g/l)*math.sin(theta1)

	theta = [dtheta1_dt,dtheta2_dt]   # Solution matrix

	# Returning the value as the output of the function
	return theta


b = 0.05  # Damping coefficient
m = 1     # Mass of the body in kg
g = 9.81  # Acceleration due to gravity in m/s^2
l = 1     # Length of the string in m

# Providing the timestep info
t = np.linspace(0,20,300)

# Initial conditions at time = 0
theta_0 = [0, 3]

# Calling the user defined function 
theta = odeint(pendulum,theta_0,t,args = (b,m,g,l))

# Plotting the displacement and velocity againt time
plt.figure(4)
plt.plot(t,theta[:,0],'r:',label='Angular displacement')
plt.plot(t,theta[:,1],'b-',label='Angular velocity')
plt.xlabel('TIME')
plt.ylabel('MOTION')
plt.title('SHM OF SIMPLE PENDULUM')
plt.legend()
plt.savefig("SHM motion.png")

# Creating the action of simple pendulum using plots
# Fixed point
x0 = 0
y0 = 0

ct = 1

# Loop to calculate the position of pendulum at each time step
for i in range(300):

	x1 = l*math.sin(theta[i,0])
	y1 = -l*math.cos(theta[i,0])

	plt.figure()
	plt.axis([-2,2,-2,2])
	plt.plot([-0.5,0.5],[0.07,0.07],'b-', linewidth = 10)
	plt.plot([x0,x1], [y0,y1],'r', linewidth = 4)
	plt.plot(x1,y1,marker ='o',markersize = 15,markerfacecolor = 'g')
	plt.savefig('pendulum%05d.png'% ct)
	ct = ct + 1

Encountered errors and debugging:

Error 1:

While defining the function, only 5 arguments are there instead of 6. The following error is deugged by correcting the number of arguments in the function definition in the beginning.

Error 2 :

The syntax for saving the figure was wrong. It is corrected by removing the "," in Line:57.

Error 3 :

While defining the 'for loop', an error occured because of wrong range for 'i', since 't' is not an integer scalar array. It is corrected by changing to 300.

Error 4:

Instead of '[]', '()' were used in calculating the coordinates of the pendulum in Line:58.While accessing a list, only '[]' should be used.

Error 5 :

An error was made by providing a wrong syntax for creating marker. So it is changed to " marker = 'o' ".

RESULTS:

Output:

Plot 1:

Inference:

From the above plot, it is evident that the ODE solved using Python are producing logical results. It clearly depicts the simple harmonic motion of the pendulum as the amplitude of both angular displacement and angular velocity are decreasing gradually as the time passes.

Animation of the pendulum:

Inference:

As seen from previous plot, the motion of the pendulum is in simple harmonic behaviour. It can be seen from the movement of pendulum is being slowed down for each swing.

CONCLUSION:

The animation created and the plot produced by solving the ODE are in good agreement. This indicates that this python program is in agreement with the experimental results.