SIMPLE PENDULUM : Solving Second Oder ODE with Python

Aim:

Create Simulation of Simple pendulum With MATLAB as a solution for second-order ODE

Objectives:

• Write a program that solves ODE which represents the equation of motion of a simple pendulum with damping

Given Details

1. Simulate the motion between 0-20 sec
2. for angular displacement = 0
3. angular velocity=3 rad/sec
4. at time t = 0
5. In the above equation,
6. g = gravity in m/s²
7. L = length of the pendulum in m
8. m = mass of the ball in kg
9. b=damping coefficient

Given Dimensions

1. L = 1 metre
2. M = 1 kg
3. b = 0.05
4. g = 9.81 m/s²

Theory

A heavy body suspended by a light inextensible thread is called a simple pendulum. This is impracticable, so in practice, a small metal sphere is suspended from rigid support using a cotton or silk thread. The point from which the bob is suspended is called the point of suspension. The center of gravity of the bob is called the centre of oscillations. The distance between the point of suspension and the center of oscillations is called the length of the pendulum.

The weight of the bob acts vertically downwards and the tension in the string acts vertically upwards. These two forces are equal and opposite and the bob has an equilibrium position in between. The bob displaced by the θ angle from the equilibrium position.

PYTHON PROGRAM For Animation of Simple Pendulum

"""
Simple Pendulum : Second Oder ODE
by Slomo Thomas
"""
#modules Imported
import math
import numpy as np 
import matplotlib.pyplot as plt
from scipy.integrate import odeint


def motion(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	
# damping co efficient
b = 0.05 
# gravity
g = 9.81 
# length
l = 1    
# mass hanging 
m = 1  

# initial condition
#angular displacement and angular velocity
theta_0 = [0,3] 

# time points
t = np.linspace(0,20,250)

# solve	ode
theta = odeint(motion,theta_0,t,args = (b,g,l,m))

# animation script
ct = 1
x0 = 0
y0 = 0

for i in theta[:,0]:
	x1 = 1*math.sin(i)
	y1 = -1*math.cos(i)

	filename = '%05d.png' %ct
	ct = ct + 1

	plt.figure()
	plt.plot([-1,1],[0,0],"green",linewidth=15)
	plt.plot([x0,x1],[y0,y1],"red")
	plt.xlim([-1.5,1.5])
	plt.ylim([-1.5,1.0])
	plt.plot(x1,y1,'o',markersize =20)
	plt.title('Oscilating Pendulum')
	plt.savefig(filename)
	

Explanation of the script

Steps involved in the program and the explanation as follows:

• Import math = This module helps in making use of certain functions like sin, cos and pi etc...
• Import matplotlib.pyplot = This module is used to create plots using the code.
• from scipy.Integrate import odeint = Scipy module let us differentiation and integration functions.
• Import numpy as np = numpy module lets us use some math functions like 'linspace'
• def = it is a defines the function or it calls for the function.
• A function called 'model ' is created to define the ode equations.
• Simple pendulum second order ODE is converted into first order ODE to make things simple.
• Inputs and intial conditions are provides
• Initial conditions are displacement and velocity.
• timestep size is provided
• theta is created to call the function. we have to be careful call the function because the range of errors.
• For the plot between the diaplacement, velocity vs time we have used 1st program
• Second program is used for the animation purpose.
• Each point has to be created in animation to tell the program where the string and the ball are at that moment of time.
• Python created random filename which will become diffcult to stich all the plot as one gif file. So, we have created at count 'ct' which increase by 1 along with timestep and this count can be used in naming the file.

Odeint Script

"""
Ode
@author: slomo
"""

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

def model(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	
# damping co efficient
b = 0.05 
# gravity
g = 9.81 
# length
l = 1    
# mass hanging 
m = 1  

# initial condition
#angular displacement and angular velocity
theta_0 = [0,3] 

# time points
t = np.linspace(0,20,250)

# solve	
theta = odeint(model,theta_0,t,args = (b,g,l,m))

plt.plot(t,theta[:,0],'b-',label=r'$\frac{dtheta_1}{dt}=\theta2$')
plt.plot(t,theta[:,1],'r--',label=r'$\frac{dtheta_2}{dt}=\frac{b}{m}\theta2-\frac{g}{L}sin\theta_1$')
plt.xlabel('Time')
plt.ylabel('Displacement')
plt.legend(loc='best')
plt.show()

Animation Video

video link

CONCLUSION

Using PYTHON, Simulation of a Simple Pendulum corresponding to the solution of a Second-order ODE is created

REFERENCE

[1] WIKIPEDIA

[2] MATLAB DOCUMENTATION

[3] NCERT PHYSICS TEXT

[4] GOOGLE IMAGES