Pendulum motion by solving ODE in python

Aim of the project :To write a program to simulate the transient behaviour of a simple pendulum and to create an animation of it\'s motion.

This was the second order ODE given to us for damped oscillatory motion of a simple pendulum.

In the above equation,

g = gravity in m/s2,

L = length of the pendulum in m,

m = mass of the ball in kg,

b=damping coefficient.

Now for solving the second order ODE, we will first convert this second order ODE into a number of first order ODEs.

We will write a separate function where we will assume two variables theta1 and theta2, where theta2 is the first order differential of theta1 and theta1 = theta. Now , using the equationn given we will make an eqauation for theta2 where

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

def model(theta,t_span,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

After having written the function, we will use in-built function odeint (which was imported using scipy.integrate library) to solving the differentail equation to get values of theta (both theta1 and theta2) at different time instances.We will get 500 values of theta1 and theta2 where theta1 will be first column of the solution and theta2 will be second column of the solution.

NOw , we will begin with the animation of the pendulum motion. For this, we will plot displacement points within the time period of 20seconds. We will keep one end of the string fixed and other end shall move according to the solution of the ODE i.e. column 1 of the theta solution. WE will use a simple for loop for execution. I used the origin as (5,5) for plotting animation.

for i in range(len(t_span)):

	x = x0 + l*math.sin(theta[i,0])
	y = y0 - l*math.cos(theta[i,0])
    
	filename = \'pendulum\'+str(ct) + \'.png\'


	plt.figure()
	plt.plot([4,6],[5,5],\'k-\')
	plt.plot([x0,x],[y0,y])
	plt.plot(x,y,marker = \'o\', markersize = \'5\')
	plt.axis([3,8,3,8])
	plt.savefig(filename)
	ct+=1

All the generated images were stored and animation was created using imagemagick.

Code :

import math
import matplotlib.pyplot as plt
from numpy import *
from scipy.integrate import odeint

l = 1   # length in mtr
b = 0.05  # damping coefficient
m = 1  # mass in kg
g = 9.81

# initial conditions

theta_0 = [0, 3]


#  timespan
t_span = linspace(0, 20, 300)

def model(theta,t_span,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


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

plt.figure()
plt.plot(t_span,theta[:,0],label =\'angular displacement\')
plt.plot(t_span,theta[:,1],\'r--\',label = \'angular velocity\')
plt.legend()
plt.show()

#creating animation
x0 =5
y0 =5

ct = 1
for i in range(len(t_span)):

	x = x0 + l*math.sin(theta[i,0])
	y = y0 - l*math.cos(theta[i,0])
    
	filename = \'pendulum\'+str(ct) + \'.png\'


	plt.figure()
	plt.plot([4,6],[5,5],\'k-\')
	plt.plot([x0,x],[y0,y])
	plt.plot(x,y,marker = \'o\', markersize = \'5\')
	plt.axis([3,8,3,8])
	plt.savefig(filename)
	ct+=1

Outputs:

A total of 300 images were used for creating a movie using imagemagik and following code line in command prompt.

ffmpeg  -i pendulum%d.png movie.avi

Youtube link for animation:

https://www.youtube.com/watch?v=qOAdhLPh5-E&feature=youtu.be

<iframe width=\"560\" height=\"315\" src=\"https://www.youtube.com/embed/qOAdhLPh5-E\" frameborder=\"0\" allow=\"accelerometer; autoplay; encrypted-media; gyroscope; picture-in-picture\" allowfullscreen></iframe>