Week 3 - Solving second order ODEs
Program to Create the Animation of Pendulum Motion
Aim:-To create python program to create the animation of pendulum motion
Equation used:- The ordinary differential equation which describes the pendulum motion is given by,
`frac{d^2theta}{dt^2}+frac{b}{m}*frac{d theta}{dt}+frac{g}{l}*sintheta = 0`
`theta` :Angle between the vertical axis and the pendulum thread.
m:Mass of the bob.
l:Length of the thread.
g:Acceleration due to gravity.
b:Damping coefficient.
Objectives:-
* To create python program which gives the solution for the ode which describes the pendulum motion .
*To create the animation of the pendulum moment and hence analyze its motion.
Procedure:-
The above given equation describes the pendulum motion. The necessary initial conditions are given to find out the solution for this ODE. At the beginning of the program necessary modules like math, matplotlib. pyplot, numpy, scipy are imported. Each modules are required to perform certain calculations. In this program one of the function of scipy module is used for find out the solution for the ode which describes the pendulum motion. The variables for length of the pendulum thread, mass of the bob, damping coefficient,gravity are assigned with given values. The linspace function of numpy is used to generate the different time points between 0 to 20 seconds. The ODE is solved by the odeint function of scipy module which takes function as argument that generates the ode . So a function is created to generate the ODE for the pendulum motion which takes mass of the bob, thread length, damping coefficient, time points, initial angular displacement and angular velocity as the arguments and generates and returns ODE. The above given ODE is generated by combining three terms of the equation. The initial conditions and the ODE is passed to the odeint function which generates the solution for the given time point. All possible values of angles are stored in the thetas list. Since , for each time point two solution for angles are possible. So they are separated by using forloop and store them in two different lists. To generate the png file describing the pendulum position , the logic is used in such a way that one point or end of the thread is fixed to a support at d distance more than the sum of length of the string and bob diameter. The another point is computed by the algebra using corresponding value of angle between vertical reference line and the thread line. Each positions are plotted and the bob is created by using marker attribute of plot function using appropriate marker size. Each png files are named by using string formatting method. When all the png files are generated , an animation is created with the help of imagemagic software.
Conclusion:-
*Using the ODE we can easily analyze the motion of the pendulum.
*The motion of the pendulum seen using the animation or gif exactly matches the motion of the pendulum in reality.
######## Program to Create the Animation for Pendulum Motion ########
import numpy as np #Importing numpy module.
from scipy.integrate import odeint # from the scipy.integrate module importing odient function to find the solution for ode.
import math #Importing math module for mathematical calculations.
import matplotlib.pyplot as plt #Importing matplotlib.pyplot module as plt.
def ode_generator(initial_conditions,t,b,g,l,m):#Creating a function to generate the ode.
theta1 = initial_conditions[0] #Assigning the 1st condition to the theta1 variable(i.e,angular_displacement).
theta2 =initial_conditions[1] #Assigning the 2nd condition to the theta2 variable(i.e,angular_velocity).
dtheta1_dt = theta2 #considering theta2 as derivative of theta1 w.r.to t.
dtheta2_dt = (-b/m)*theta2 -(g/l)*math.sin(theta1)#Required form of ode.
dtheta_dt = [dtheta1_dt,dtheta2_dt]
return dtheta_dt #Returning the generated ode to the main function.
l = 1 #Length of the thread (in m).
g = 9.81 #Acceleration due to gravity (in m/s^2).
b = 0.05 #Dampping Coefficient.
m = 1 #Mass of ball (in kg).
d = 0.25#Diameter of the bod(in m).
#Initial conditions
angular_displacement = math.radians(0)#Assigning initial angular displacement as 0 (in rad).
angular_velocity = 3 #Assigning initial angular velocity as 3(in rad/sec).
initial_conditions = []#Declaring an empty list for accomadation of initial conditions.
#Appending given conditions to the initial_conditions list.
initial_conditions.append(angular_displacement)
initial_conditions.append(angular_velocity)
#Time points
t = np.linspace(0,20,300)
#Solving ordinary differential equation generated by ode_generator function
thetas = odeint(ode_generator,initial_conditions,t,args = (b,g,l,m))#Finding the solution for the ode using odeint function of the scipy.integrate module.
theta_set1=[]#Declaring the list to store the first set of possible solutions.
#Forloop to separate the two possible set of solution generated.
for theta in thetas:
theta_set1.append(theta[0])#Appends the first value of theta in the each array to the theta_set1 array.
#Position where the pendulum is connected to the support.
x0,y0 = l,l+0.5
#Initializing the image number.
imgnum = 1
#Forloop to create plots for each position of the pendulum.
for i in range(0,len(set1)):
x1 = x0 +(l*math.sin(theta_set1[i]))#Calculates the position of the bob along x-axis for corresponding angle.
y1 = y0 -(l*math.cos(theta_set1[i]))#Calculates the position of the bob along y-axis for corresponding angle.
plt.figure() #Function which creates new figure.
imgname = 'Img%05d.png' %imgnum #Naming format for the plot to be generated.
imgnum = imgnum + 1 #Incrementing the image number.
plt.plot([x0,x1],[y0,y1],color = 'k')#Plotting the length of the thread for different angles.
plt.plot([0,10],[y0,y0],linewidth=5,color = 'k')#Plotting the base support for the pendulum.
plt.plot(x1,y1,marker='o',markersize = 40,markeredgecolor = 'k',markerfacecolor = 'r')#plotting the bob using marker attribute of the plot function.
plt.xlim(0,2)#Limitting the x-axis length to 2m.
plt.ylim(0,1.52)#Limitting the y-axis length to 1.52m.
plt.savefig(imgname) #Saving plot as given fileformat.
#################################### END ###########################################