Solving second order Differential Equations using Python

AIM:  To simulate the transient behavior of a simple pendulum and to create an animation of it's motion.

OBJECTIVE : 

1. To write a program in Python that solves the second-order ODE corresponding to the motion of a simple pendulum and to plot its angular displacement and angular velocity wrt time.
2. To create the animation of the motion of pendulum by simulating the motion between 0-20 sec, for angular displacement=0, angular velocity=3 rad/sec at time t=0.

THEORY:  A simple pendulum is defined to have an object that has a small mass, also known as the pendulum bob, which is suspended from a light wire or string, such as shown in Figure 1.

                            Figure 1

BODY OF THE CONTENT :

1. At first math, matplotlib, numpy and scipy modules were imported for their various functions.
2. Initial inputs such as b,g,l,m are provided according to data in the question.
3. An array of initial theta was created to define initial angular displacement and initial angular velocity.
4. An array of time span was created to store 100 equally spaced intervals between time 0 to 20 seconds.
5. To solve the second order ODE, a function named ode_function which stores the above matrix data in it.
6. The given second order ODE was solved using 'odeint' command followed by proper syntax of the calling function and the result obtained was used to plot the graph of variation of angular displacement and angular velocity wrt time.
7. Now for creating the animation of obtained harmonic motion of pendulum, for loop was used and a loop counter variables was initiated,the values of angular position is stored in a separate variable called 'tmp' and is iterated in for loop. In this execution the one end of the string of pendulum was fixed and the other end moves according to the solution of the ODE solved and was plotted for each time interval of the time span.

PYTHON CODE FOR SOLVING ODEs AND  Its ANIMATION :-

''' Code for solving second order ODE's  and simulate its motion'''
import math
import matplotlib.pyplot as plt
import numpy as np
from scipy.integrate import odeint


def ode_function(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
# Inputs for the simple pendulum
# Damping coefficient
b = 0.05
# Acceleration due to gravity(m/s^2)
g = 9.81  
# Length of string (m)                                                               
l = 2  
# Mass of bob (kg)                                                                  
m = 2      

# Initial Conditions
theta_0 = [0,3] 
# time points
t_span = np.linspace(0,20,150)

# Solving Second order differential equation with scipy plugin

theta = odeint(ode_function,theta_0,t_span,args = (b,g,l,m))

# Plotting the graph for the Damped harmonic motion of a pendulum 
plt.plot(t_span,theta[:,0],'b-',label=r'$frac{dtheta1}{dt}=theta2$')
plt.plot(t_span,theta[:,1],'r--',label=r'$frac{dtheta2}{dt}=-frac{b}{m}theta2-frac{g}{l}sintheta1$')
plt.xlabel('Time (second)')
plt.ylabel('Plot')
plt.title('Harmonic motion of a Damped Simple Pendulum')
plt.legend(loc='best')
plt.grid()
plt.show()

# Iteration using for loop to create animation. Hence to obtain image at different position of motion of pendulum

loop_counter =1
tmp = theta[:,0]
for i in tmp:
	x_start = 0                                                           
	y_start = 0
	x_end =  l*math.sin(i)                                          
	y_end = -l*math.cos(i)
	# Plotting for animation
	plt.figure()
	# Plotting String of pendulum
	plt.plot([x_start,x_end],[y_start,y_end],linewidth=3)
	# Plotting Bob of pendulum
	plt.plot(x_end,y_end,'o',markersize=20,color='r')
	# Plotting the rigid support through which string of pendulum is attached
	plt.plot([-0.5,0.5],[0,0],linewidth=8,color='k')
	# Plotting hatched vertical line for refernce
	plt.plot([0,0],[0,-1.5],'--',color='y')
	plt.xlim([-2,2])
	plt.ylim([-3,2])
	plt.title('Animation of transient behaviour of a Simple Pendulum')
	filename = 'test%05d.png' % loop_counter
	plt.savefig(filename)
	loop_counter = loop_counter + 1

ERRORS : 

1. missed to provide loop_counter value 

Rectified by :-

RESULTS :

• The following plot was generated for the given program -

• The above plot of angular displacement and angular velocity shows the pendulum's motion to be gradually decreasing. This is due to damping effect which exponentially decreases the amplitude and the velocity of motion.

• For animation 150 plots were generated and were stitched together to form animation.The link of the animation of simulation-

       YOUTUBE LINK FOR SIMULATION

CONCLUSION :  In this way we can understand that any second order or any higher order ODE can easily be solved using Python, and simulation of transient behaviour of simple pendulum was made using it to explain the effect of damping on it's motion.