Simulation Of Simple Pendulum By Using Second Order ODE In Python
Aim:
To write and execute the program code for creating the animation of Simple Pendulum.
Description:
Transient behaviour is the study of the observing one or more variables,about,how it is changing with respect to time.Here, such transient behaviour of the Simple Pendulum is to be studied with the help of the Ordinary Differential Equation(ODE).Ordinary Differential Equation is the differential equation which contains one or more fuctions of one independent variable and the derivatives of those functions.
For Example,
is the Ordinary Differential Equation.Here, \'x\' is the independent variable,\'\' is the function of the independent variable and \'dx/dt\' is the derivatives of the function.Since the highest numbered derivative in this equation is one, it is called the First Order Ordinary Differential Equation.
ODEs are based on the fact that, If we know the function and possibly, some of the derivatives at a particular point, then by this information and together with the differential equation, we can determine that function over its entire domain.So that, the ODEs are highly recommended in physics.Generally, the motion of a body is described by its position and velocity as the time value varies. Newton\'s laws allow these variables to be expressed dynamically as a differential equation for the unknown position of the body as a function of time.Here, the transient behaviour of the Simple Pendulum is observed by Newton\'s Second Law using Second Order ODE.
The motion of the Simple Pendulum with damper is represented in ODE as,
where,
m - mass of pendulum bob ,kg
L - Length of pendulum ,m
b - damping coefficient ,Ns/m
g - gravitational acceleration ,m/sec^2
In order to solve the second order (or) any higher-order differential equations, it is recommended to convert it to a system of first order differntial equations and proceed to solve.This conversion is in the picture.
Coding:
#Program for Animation of Simple Pendulum
import math
import matplotlib.pyplot as plt
from scipy.integrate import odeint
import numpy as np
#Function Definition
def simple_pendulum(theta,t,b,g,m,l):
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
b = 0.05 #Damping Coefficient
g = 9.81 #Gravity
m = 1 #Mass of the bob
l = 1 #length of the bob
#Initial Conditions
theta_0=[0,3]
#Time
time = np.linspace(0,20,100)
#Solving Differential Equation
results_ode = odeint(simple_pendulum,theta_0,time,args =(b,g,m,l))
ct=1
for i in results_ode[:,0]:
x1 = 0
y1 = 0
x2 =l*math.sin(i)
y2 = -l*math.cos(i)
plt.figure()
filename = str(ct) + \'.png\'
ct = ct + 1
plt.hlines(0,-0.5,0.5)
plt.plot((x1,x2),(y1,y2))
plt.plot(x2,y2, linestyle=\'\', marker =\'o\',markersize= 20,zorder=3)
plt.xlim([-1.5,1.5])
plt.ylim([-1.5,1.5])
plt.savefig(filename)
#Plotting
plt.figure()
plt.plot(time,results_ode[:,0],\'b-\',label = r\'$\\frac{d\\theta_1}{dt}=\\theta2$\')
plt.plot(time,results_ode[:,1],\'r--\',label = r\'$\\frac{d\\theta2}{dt}=-\\frac{c}{m}\\theta_2-\\frac{g}{l}sin\\theta_1$\')
plt.xlabel(\'time\')
plt.ylabel(\'plot\')
plt.legend(loc =\'best\')
plt.show()
Coding Explaination:
import math
import matplotlib.pyplot as plt
from scipy.integrate import odeint
import numpy as np
#Function Definition
def simple_pendulum(theta,t,b,g,m,l):
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
The libraries required for this program are imported.Especially, in order to solve ordinary differential equation, \'scipy\' module of \'odeint\' solver is imported.Then the Function named \'simple_pendulum\' is defined to solve those converted simple ordinary differential equations.The calculated values are stored in the array \'dtheta_dt\' and \'return\' command is used to return those values when the function is called in the main program.
#Inputs
b = 0.05 #Damping Coefficient
g = 9.81 #Gravity
m = 1 #Mass of the bob
l = 1 #length of the bob
#Initial Conditions
theta_0=[0,3]
#Time
time = np.linspace(0,20,100)The parameters of simple pendulum such as length of the pendulum, mass of the bob, damping coefficient for damper are provided as inputs.Also the boundary conditions for angular displacement and angular velocity and the time interval between which the animation to be happen also defined as 0 to 20sec of 100 values is specified using \'np.linspace()\' command.
#Solving Differential Equation
results_ode = odeint(simple_pendulum,theta_0,time,args =(b,g,m,l))The function named \'simple_pendulum\' is called to solve the ODE using command \'odeint\' and storing the results in an array \'results_ode\'.
ct=1
for i in results_ode[:,0]:
x1 = 0
y1 = 0
x2 =l*math.sin(i)
y2 = -l*math.cos(i)
plt.figure()
filename = str(ct) + \'.png\'
ct = ct + 1
plt.hlines(0,-0.5,0.5)
plt.plot((x1,x2),(y1,y2))
plt.plot(x2,y2, linestyle=\'\', marker =\'o\',markersize= 20,zorder=3)
plt.xlim([-1.5,1.5])
plt.ylim([-1.5,1.5])
plt.savefig(filename)
Here, the endpoint of the pendulum is fixed and \'for loop\' is used to plot the positions of the bob with respect to time with the calcuated values in the \'results_ode\' array.The pendulum bob is created using \'markers\' and \'markersize\' and assigning it as a \'text\' providing \'zorder\' as \'3\' .The \'savefig()\' command is used to get each and every plot and save images as \'.png\' format from the array \'filename\'
#Plotting
plt.figure()
plt.plot(time,results_ode[:,0],\'b-\',label = r\'$\\frac{d\\theta_1}{dt}=\\theta2$\')
plt.plot(time,results_ode[:,1],\'r--\',label = r\'$\\frac{d\\theta2}{dt}=-\\frac{c}{m}\\theta_2-\\frac{g}{l}sin\\theta_1$\')
plt.xlabel(\'time\')
plt.ylabel(\'plot\')
plt.legend(loc =\'best\')
plt.show()
Here, the values of angular displacement and angular velocity is plotted with time and the variation is shown and the respective equations also labelled in the plot.
Error:
(i)Specified the angle in \'degrees\' in function definition which is wrong. In Python, it should be in radians and the error shows as,
Output:
The Youtube Link for Animation is linked here: