"""
Program for simulating the motion of simple pendulum
"""
import numpy as np
from scipy.integrate import odeint
import matplotlib.pyplot as plt
import math
# Defining inputs
m = 1 # mass of the ball(kg)
L = 1 #Length of the pendulum (m)
b = 0.05 #Damping Coefficient
g = 9.81 # Acceleration due to gravity(m/s2)
# Initial conditions i.e theta=0 and omega=3 rad/s at t=0
y0 = [0,3]
# Time interval for simulation of simple pendulum
t_interval = np.linspace(0,20,500)
#Function for ordinary differential equation
def model(y0,t_interval,m,L,b,g):
# Angular displacement
theta = y0[0]
# Angular velocity
omega = y0[1]
# Ordinary differential equation
dtheta_dt = omega
domega_dt = -((b/m)*omega + (g/L)*(theta))
dTHETA_dt = [dtheta_dt, domega_dt]
return(dTHETA_dt)
# Solving the first order differential equation using odeint function from scipy module.
solution = odeint(model,y0,t_interval,args = (m,L,b,g))
Angular_Displacement = []
Angular_Velocity = []
for i in range(0, len(solution)): #len(solution) means number of elements in the solution array
# Angular displacement
angular_displacement = solution[i,0]
Angular_Displacement.append(angular_displacement)
# Angular velocity
angular_velocity = solution[i,1]
Angular_Velocity.append(angular_velocity)
# Plotting a plot of variation of angular displacement and angular velocity with respect to time
plt.figure(2)
plt.plot(t_interval, Angular_Displacement, color = 'r', linewidth = 2 , label = 'Angular Displacement')
plt.plot(t_interval, Angular_Velocity, color = 'b', linewidth = 2, label = 'Angular_Velocity')
plt.legend()
plt.grid()
plt.xlabel('time')
plt.ylabel('Angular velocity')
plt.savefig('Plot of variation of angular velocity and angular displacement with respect to time')
plt.show()
# Creating an animation of the motion of simple pendulum
# Counter variable
ct = 1
# For loop for creating different positions of simple pendulum
for j in range(0,len(t_interval)):
Theta = Angular_Displacement[j]
Omega = Angular_Velocity[j]
# Defining co-ordinates of fixed end of string of pendulum
X0 = 0
Y0 = 0
# Defining co-ordinates of free end of string of pendulum
X1 = L*(math.sin(Theta))
Y1 = -L*(math.cos(Theta))
# Plot for animation
plt.figure(2) # Creating an empty figure window
plt.plot([-1,1],[0,0], color ='r', linewidth = 3) # Plot for creating support
plt.plot([X0,X1],[Y0,Y1],color = 'k', linewidth = 2) # Plot for creating string of the pendulum
plt.plot(X1,Y1, marker = 'o',markerfacecolor = 'm', markeredgecolor = 'k',markersize =15)
# Inserting text in a plot
plt.text(-1, 0.3,'Angular Displacement =' + str(Theta) + 'rad', color ='r')
plt.text(-1, 0.2,'Angular Velocity=' + str(Omega) + 'rad/s', color = 'b')
# Giving title to a graph
plt.title('Animation of motion of simple pendulum', fontsize = 15, fontweight = 'bold')
plt.grid()
plt.savefig('Animation of simple pendulum %05d.png' %ct)
plt.show()
filename = str(ct) + '.png'
ct = ct + 1
# Explaination of the code
# Importing the required basic modules
1) import.math = This module is required for operating basic trigonometric functions.
2) import numpy as np = This is used to execute the linspace command to create the
arrays.
3) from scipy.integrate import odeint = This module is used to execute odeint function
from scipy module which is used to solve the
first order ordinary differential equations.
4) import matlab.pyplot as plt = Used for importing the plotting module to obtain the
plots.
Output: All the required module gets imported.
# Defining the required inputs
m = 1 # mass of the ball(kg)
L = 1 #Length of the pendulum (m)
b = 0.05 #Damping Coefficient
g = 9.81 # Acceleration due to gravity(m/s2)
y0 = [0,3] # y0 represents the initial conditions. 0 represents angular displacement
and 3 represents angular velocity.
t_interval = np.linspace(0,20,500) # The integration is done over time 0-20 secs over
500 points. It is presented as a row vector. In
other time interval of 20 secs is divided over
500 points.
Output = The required inputs get defined to represent the angular displacement and
angular velocity over a time interval.
# Defining the required function for solving ODE
def model(y0,t_interval,m,L,b,g):
# Here model is the function name which contains the required inputs like initial conditions, time interval, mass, length, damping co efficient and gravity.
# Angular displacement
theta = y0[0] # 0 represents the displacement
# Angular velocity
omega = y0[1] # 1 represents the velocity
# Ordinary differential equation
dtheta_dt = omega # Differentiation of theta (Angular Displacement gives
us the angular velocity (omega)
domega_dt = -((b/m)*omega + (g/L)*(theta)) # Solving the equation for domega_dt
gives the following term.
dTHETA_dt = [dtheta_dt, domega_dt] # Both the above terms are stored in
a matrix dTHETA_dt.Both the terms
are integrated over a time interval
using the ode_int function.
return(dTHETA_dt) # If we define a function we need to call the function
by return command.
Output: The function gets defined by input arguments like initial conditions, time interval, mass, length, damping coefficient, gravity and with return command the function gets called.
solution = odeint(model,y0,t_interval,args = (m,L,b,g))
# The solution variable solves the ordinary differential equations using the function odeint.
Angular_Displacement = [] # Creation of empty arrays for angular displacement Angular_Velocity = [] # Creation of empty arrays for angular velocity
for i in range(0, len(solution)):
# Angular displacement
angular_displacement = solution[i,0]
Angular_Displacement.append(angular_displacement)
# Angular velocity
angular_velocity = solution[i,1]
Angular_Velocity.append(angular_velocity)
# The above for loop code is used to append the angular displacement and angular velocity values in the empty arrays. It is done for the range of the values for elements in the solution array.
plt.figure(2)
plt.plot(t_interval, Angular_Displacement, color = 'r', linewidth = 2 , label = 'Angular Displacement')
plt.plot(t_interval, Angular_Velocity, color = 'b', linewidth = 2, label = 'Angular_Velocity')
plt.legend()
plt.grid()
plt.xlabel('time')
plt.ylabel('Angular velocity')
plt.savefig('Plot of variation of angular velocity and angular displacement with respect to time')
plt.show()
Output: The above commands are used to plot the angular displacement and angular velocity with respect to time interval.
Final step: To create the animation of simulation pendulum for different positions of theta and plotting the results of angular displacement and angular velocity accordinly.
Couldn't create the animation as only images are created. I will require help in this regard.
# Counter variable
ct = 1
# For loop for creating different positions of simple pendulum
for j in range(0,len(t_interval)):
Theta = Angular_Displacement[j]
Omega = Angular_Velocity[j]
# Defining co-ordinates of fixed end of string of pendulum
X0 = 0
Y0 = 0
# Defining co-ordinates of free end of string of pendulum
X1 = L*(math.sin(Theta))
Y1 = -L*(math.cos(Theta))
# Plot for animation
plt.figure(2) # Creating an empty figure window
plt.plot([-1,1],[0,0], color ='r', linewidth = 3) # Plot for creating support
plt.plot([X0,X1],[Y0,Y1],color = 'k', linewidth = 2) # Plot for creating string of the pendulum
plt.plot(X1,Y1, marker = 'o',markerfacecolor = 'm', markeredgecolor = 'k',markersize =15)
# Inserting text in a plot
plt.text(-1, 0.5,'Angular Displacement =' + str(Theta) + 'rad', color ='r')
plt.text(-1, 0.3,'Angular Velocity=' + str(Omega) + 'rad/s', color = 'b')
# Giving title to a graph
plt.title('Animation of motion of simple pendulum', fontsize = 15, fontweight = 'bold')
plt.grid()
plt.axis([-1.5,1.5,-1.5,1])
plt.savefig('Animation of simple pendulum %05d.png' %ct)
ct = ct + 1
plt.plot()