Python program to simulate the transient behaviour of a simple pendulum

AIM : To write a Python program to simulate the transient behavior of a simple pendulum and to create an animation of its motion.

BACKGROUND : In Engineering, ODE is used to describe the transient behavior of a system. A simple example is a pendulum. The way the pendulum moves depends on the Newton’s second law. When this law is written down, we get a second order Ordinary Differential Equation that describes the position of the "ball" w.r.t. time.

Similarly, there are hundreds of ODEs that describe key phenomena and if we have the ability to solve these equations then we will be able to simulate any of these systems.

Our objective is to write a program that solves the following ODE. This ODE represents the equation of motion of a simple pendulum with damping.

`frac{d^2theta}{dt^2}+frac{b}{m}frac{d theta}{dt}+frac{g}{L}sin theta=0`

g = gravity in m/s²,

L = length of the pendulum in m,

m = mass of the ball in kg,

b=damping coefficient.

INPUT VARIABLES :

L=1 metre, m=1 kg, b=0.05, g=9.81 m/s².The motion is simulated between 0-20 seconds.

INITIAL CONDITIONS :

At time t=0, angular displacement=0 ; angular velocity= 3 rad/sec.

CALCULATIONS : In order to solve the above ODE, the following calculation is done to get the solution through MATLAB. The basic idea involved here is to break down a single `n^{th}`order ODE into 'n' lower order ODEs and representing them in matrix form to obtain the solution vector.

let, `theta=theta1`

`frac{d theta}{dt}=frac{d theta1}{dt}=theta2` `

hence, `frac{d^2 theta}{dt^2}=frac{d^2 theta1}{dt^2}=frac{d}{dt}(frac{d theta1}{dt})=frac{d}{dt}theta2`

Using these expressions in the governing ODE,

`frac{d^2 theta1}{dt^2}+frac{b}{m}frac{d theta1}{dt}+frac{g}{L}sin theta1=0`

`frac{d theta2}{dt}+frac{b}{m}theta2+frac{g}{L}sin theta1=0`

hence, `frac{d theta2}{dt}=-frac{b}{m}theta2-frac{g}{L}sintheta1`

Representing in matrix form,

`frac{d}{dt}((theta1),(theta2))= [[theta2,0],[-frac{b}{m}theta2,-frac{g}{L}sin theta1]]`

SOLUTION ALGORITHM

1. The input variables will be entered as mentioned in the problem.
2. The time period will be split into 200 parts between limits of 0-20 seconds as mentioned in the problem.
3. The displacement and velocity as a function of time will be calculated in a function file and it will be called later.
4. First, the plotting will be done to represent displacement and velocity as a function of time with displacement and velocity taken along y-axis and time taken along x-axis.
5. Then an animation plot will be created to simulate the movement of pendulum.
6. The horizontal bar, vertical bar of the pendulum and the pendulum bob will be created.
7. The position of the pendulum at particular time will be captured and finally all the frames will be combined.
8. The output will be obtained in a video file to represent the motion in the form of an animated video.
9. The entire script will be written in Python using Sublime Text Editor.

PYTHON CODE :

"""
program to solve 2nd order ODE
the ode describes the transient behaviour of a simple pendulum

By Abhishek Arora

"""

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

def system(theta_init,time,b,g,l,m):

#angular displacement
	theta1= theta_init[0]
#angular velocity
	theta2= theta_init[1]

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

	dtheta_dt= [dtheta1_dt,dtheta2_dt]
	return dtheta_dt


#initial conditions

#damping coefficient value
b= 0.05

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

#length of the pendulum(m)
l= 1

#mass of the ball(kg)
m= 1

#initial conditions for angular displacement and angular velocity
theta_init= [0,3]

#time duration to integrate the ode is taken as 0 to 10 seconds
#it means we will be dipicting the motion of the pendulum in this range
time= np.linspace(0,20,200)

theta= odeint(system,theta_init,time,args= (b,g,l,m))

"""
odeint is an inbuilt solver to solve odes
first term in the bracket represents the function  that defines the equations
second term represents the initial conditions
third term represents the time array within which analysis is to be done
forth term are the input arguments
"""
#storign the results in 2 arrays
ang_disp=[]
ang_vel=[]

for i in range(0,len(theta)):
	result1= theta[i,0]
	ang_disp.append(result1)

	result2= theta[i,1]
	ang_vel.append(result2)

#plotting the results
plt.figure(1)
plt.plot(time,theta[:,0],'b-',label= 'angular displacement')
plt.plot(time,theta[:,1],'r-',label= 'angular velocity')
plt.legend()
plt.xlabel('time')
plt.ylabel('displacement and velocity')
plt.show()

#making an animation to simulate the movement
ct=1
for i in range(0,len(time)):
	disp= ang_disp[i]
	vel= ang_vel[i]

#coordinates for fixed end
	x0= 0
	y0= 0

#coordinates for free end
	x1=  l * (math.sin(disp))
	y1= -l * (math.cos(disp))
	

	

#plotting the pictures
	plt.figure(2)

#support
	plt.plot([-1,1],[0,0],'b')

#pendulum
	plt.plot([x0,x1],[y0,y1],'g')

#bob
	plt.plot(x1,y1,marker='o',color='r')

#providing the axis limits
	plt.axis([-1.5,2,-1.5,1])
	plt.title('simple pendulum animation')

	#saving the plots
	filename= str(ct) + '.png'
	plt.savefig(filename)
	ct= ct+1
	plt.clf()
	

OUTPUT : The following output is obtained after running the code :

* The animation of the output can be viewed from the following youtube link:

https://youtu.be/GOH_f6sRopw

CONCLUSION :

• The same model can be used to solve other order ODE and simulate any other system.
• Solving the ODEs is extremely helpful especially while doing vehicle body dynamics.

NOTE :

• ODE int is an inbuilt Python solver that is used to solve the ODE here.
• To create the animation video FFmpeg software is used.
• Command prompt with the following text can be used to create an animation file.