All Courses
All Courses
Courses by Software
Courses by Semester
Courses by Domain
Tool-focused Courses
Machine learning
POPULAR COURSES
Success Stories
Question: Write a program in Python that will simulate the pendulum motion, just like the one shown in the start of this challenge. Use, L=1 metre, m=1 kg, b=0.05. g=9.81 m/s2. Simulate the motion between 0-20 sec, for angular displacement=0, angular velocity=3 rad/sec at time t=0. When you do this you will get the position…
Sandesh Jogdand
updated on 14 Jul 2020
Question:
Write a program in Python that will simulate the pendulum motion, just like the one shown in the start of this challenge. Use,
L=1 metre,
m=1 kg,
b=0.05.
g=9.81 m/s2.
Simulate the motion between 0-20 sec, for angular displacement=0, angular velocity=3 rad/sec at time t=0.
When you do this you will get the position of the pendulum as a function of time. Use this information to animate the motion of the pendulum,
Objective: To obtain velocity and displacement with respect to time plots by solving a 2nd order ODE for a simple pendulum using a Python program and animate the motion of the pendulum using the information.
Solution:
Using the mass-spring-damper system method, a simple pendulum motion can be described using the equation of motion as:
d2dt2(θ)+(cm)⋅ddt(θ)+(gl)⋅sin(θ)=0
where,
θ = angular displacement in [radians]
c = damping co-efficient [Ns/m]
m = Mass of the pendulum in [Kg]
l = Length of the string in [meters]
g = acceleration due to gravity in [m/s^2]
As the equation has a highest degree of ‘2’, it is a 2nd order ODE equation and needs to be converted into 2 simple 1st order ODEs to sole it using programming. We do that by substituting 1st order derivative term with a variable ‘θ2’ such that,
θ = θ1,
To solve the ODE in python (integrate) we use ‘odeint’ command from the ‘scipy’ module. Here, we define the ODE as a function, and we solve it for certain initial conditions over an array of ‘t’ time steps which are determined by the program itself. We also have to define the constants/arguments which are a part of the function in specific order.
By integrating the ODEs we get an array/list with 2 columns and multiple rows as solution. The 1st column has the solution of 1st ODE that is ‘θ1’ or the angular displacement and the 2nd column has solution for 2nd ODE – ‘θ2’ = θ1’ = ‘ω’ or the angular velocity. Using the 1st solution, θ1, we can find the co-ordinates of the pendulum for each solution and by assuming the pendulum to be fixed at (0,0) we can plot it as a line for every solution and get the animation.
Python Program:
import math
import matplotlib.pyplot as plt
from scipy.integrate import odeint
import numpy as np
import matplotlib.path as mpath
def equation_of_motion_model(theta,t,c,g,l,m):
theta1 = theta[0]
theta2 = theta[1]
dtheta1_dt = theta2
dtheta2_dt = -(c/m) * theta2 - (g/l) * math.sin(theta1)
dtheta_dt = [dtheta1_dt, dtheta2_dt]
return dtheta_dt
#the motion of a pendulum is governed by following equation: d^2/dt^2(theta) + (c/m)*d/dt(theta) + g/l*sin(theta)
#Input variables:
c = 0.05 #damping co-efficient Ns/m
m = 1 #Mass of the pendulum in Kg
l = 1 #Length of the string in meters
g = 9.81 #acceleration due to gravity in m/s^2
#Initial conditions to start the simulation from: thetax=[theta_initial, omega_initial]
theta0 = [0,3]
#defining time window and time steps for calculation:
t = np.linspace(0,20,301)
#now, integrating the ODEs to get solution arrays to be plotted:
theta = odeint(equation_of_motion_model,theta0,t,args=(c,g,l,m))
#print(theta)
'''
plt.plot(t,theta[:,0],'b-',label='Angular displacement')
plt.plot(t,theta[:,1],'r--',label='Angular velocity')
plt.legend(loc='best')
plt.xlabel('Time in seconds')
plt.ylabel('Model outputs')
plt.grid()
plt.show()
'''
'''
x = equation_of_motion_model(thetax,c,m,g,l)
print(x)
'''
#pendulum animation:theta[0] represents angular displacement of the pendulum hence is to be plotted as a plot for each timestep.
#as length of the pendulum string is 1m, the angle pendulum makes with the vertical is theta,
x0 = [0]
y0 = l
x1 = theta[:,0]
x3 = []
y3 = []
ct = 1
for x2 in x1:
x = l * math.sin(x2)
y = -1 * (l * math.cos(x2))
sphere = mpath.Path.unit_circle()
plt.plot([0,x],[0,y],'g',marker = sphere, markersize = 6)
plt.plot([-0.75,0.75],[0,0],'black',linewidth = 7)
plt.plot([0,0],[-1.2,0.2],'b--')
plt.ylim(-1.25,0.05)
plt.xlim(-1.25,1.25)
plt.grid()
ct = ct + 1
filename = str(ct)+'.png'
plt.savefig(filename)
plt.show()
#x3.append(x)
#y3.append(y)
#print(x3)
#print(y3)
j = []
k = []
#plt.plot([x3],[y3],'r')
#plt.show()
Procedure:
Plots:
1. Example plot:
2. Challenge problem plot:
3. Flipbook animation:
4. Error faced:
Due to the incorrect order in which the arguments were input in the ‘odeint’ command, incorrect output plots were being generated. After going through the ‘odeint’ tutorials on various sites and scipy.integrate.odepack section on scipy.org as well through my IDE pycharm I realised there is specific order to it and after jumbling them around I found the right order.
Conclusion:
1. The method used in this example can be used for even higher ordered derivative functions to be solved using programming. The splitting of higher order ODEs into a bunch of 1st order ODEs is simpler to solve and with the help of 'def' and function method can be automated with relative ease.
2. This program with some changes can be used to describe any spring-mass-damper system and can be altered to output various results for change in lets say damping coefficient, or length, initial conditions, time period etc. for a comparative study in design phase.
Leave a comment
Thanks for choosing to leave a comment. Please keep in mind that all the comments are moderated as per our comment policy, and your email will not be published for privacy reasons. Please leave a personal & meaningful conversation.
Other comments...
To analyze and compare different cases with various contact types and different parameters as well as notch positions to compare the results
Objective: To analyze and compare different cases with various contact types and different parameters as well as notch positions and compare the results for the same boundary conditions. Given model is a crash tube with hollow square configuration meshed with all Quad QEPH 24 shell elements. Two halves of the tube are…
14 Oct 2023 01:21 PM IST
CHT analysis for a Graphics Card with heat single source
Project – 1 CHT Analysis on a Graphics Card. Objective: To perform Conjugate Heat transfer Analysis on a graphics card which consists 1.GPU which is the only heat source and is in contact with the 2.cooling fins which loose the heat to incoming air stream via convection, 3.other electronics components and a PCB on…
10 Jan 2021 08:30 AM IST
To analyse and compare results for different types of materials/material properties with the same boundary conditions.
Objective: To analyse and compare results for different types of materials/material properties with the same boundary conditions. #The results are discussed in ppt attched Case 1:Law 2-M2_Plas_Johns_Zeril – Elasto-plastic: In this case the material is considered to be an elasto-plastic with input parameters…
10 Jan 2021 08:27 AM IST
To perform side impact analysis on a partial BIW model as per the Euro NCAP standards for velocity, positioning and dimensions
Objective: To perform side impact analysis on a BIW model minimised as per the Euro NCAP standards for velocity, positioning and dimensions The Setup: A 10 Inches thick pole is used as a rigid wall, impact occurs with the BIW model travelling at 32Km/hr or 8.9mm/ms such that the horizontal centreline of the model…
10 Jan 2021 08:26 AM IST
Related Courses
0 Hours of Content
Skill-Lync offers industry relevant advanced engineering courses for engineering students by partnering with industry experts.
© 2025 Skill-Lync Inc. All Rights Reserved.