PYTHON Program to solve second order ODE of simple pendulum

Aim:- To write a program that solves the ODE which represents the equation of motion of a simple pendulum with damping.

Theory:- A simple pendulum consists of a mass m hanging from a string of length L and fixed at a pivot point P. When displaced to an initial angle and released, the pendulum will swing back and forth with periodic motion.

We will come across many differential equations of different orders in engineering which can we be easily solved using programming. Similarly, here we solve the second order differential equation which represents the motion of simple pendulum which is shown below.

`(d^2theta) /dt^2+b/m (d theta)/dt+g/L sintheta=0`

Where,

`theta =`angular displacement `t=`time

`b=`damping coefficient `m=`mass of bob

`g=`acceleration due to gravity `L=`length of string

Code:-

'''
Solving second order ODE
'''
#importing libraries
import math
import matplotlib.pyplot as plt
import numpy as np
from scipy.integrate import odeint

#function for converting second order to first order
def model(theta,t,b,g,l,m):

	theta1=theta[0]
    theta2=theta[1]
	dthetha1_dt=theta2
	dthetha2_dt= -(b/m)*theta2 - (g/l)*math.sin(theta1)
	dtheta_dt = [dthetha1_dt, dthetha2_dt]
	return dtheta_dt

#inputs
b=0.1       #damping coefficient
g=9.81      #acceleration due to gravity
l=1         #length of string
m=1         #mass of bob

#initial conditions (angular displacement, angular velocity)
theta_0=[0,3]

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

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

#plot results
plt.figure()
#plot of angular displacement vs time
plt.plot(t,theta[:,0],'b-',label='angular displacement')
#plot of angular velocity vs time
plt.plot(t,theta[:,1],'r:',label='angular velocity')
plt.legend(loc='best')

#adding labels
plt.xlabel('Time (s)')
plt.ylabel('Angular displacement & Angular velocity')
plt.show()

#simulation of pendulum
ct=1    #counter variable

#for loop for finding position of bob at each interval
for i in (theta[:,0]):
    #fixed point of pendulum string
    x0=0
    y0=0

    #free end/bob coordinates
    x1=l*math.sin(i)
    y1=-l*math.cos(i)
    filename=str(ct)+'.png'

    #plotting
    plt.figure()
    #plotting line as fixed support
    plt.plot([-1, 1],[0,0],lw='5',color='k')
    #plotting string and bob
    plt.plot([x0,x1],[y0,y1],lw='2', color='b')
    plt.plot(x1,y1,'o',markersize=15,markerfacecolor='r')

    #title and axis limits
    plt.title('Simple Pendulum')
    plt.xlim([-1.25,1.25])
    plt.ylim([-1.25,0.25])

    #saving figures
    plt.savefig(filename)

    ct=ct+1 #incrementing counter variable

Explanation:-

• The libraries are imported.
• Function which converts second order into first order ODE is defined which returns an array consisting of angular displacement and angular velocity.
• Inputs, initial conditions are defined. And time is defined using linspace() command between 0 s to 20 s with 150 intervals.
• Then odeint function is called to solve the ODE.
• And using plot command angular displacement and angular velocity w.r.t. time are plotted with proper labels and legend and it's shown using show() command.

Code for the simulation of pendulum

• First counter variable is defined ast ct. Then we use for loop to find bob postion for different values of angular displacement.
• Coordinates of fixed end and free end are defined. And filename is defined.
• Pendulum is plotted with fixed support, string and bob. And title and axis limits are defined.
• These figures are saved in the same folder as that of the program with respective file name. And counter variable is increamented and the process repeats.
• These figures are later stiched together to create animation.

Result:-

Saved figures which are converted to animation.

Figure 1

From the figure we can see that the amplitude of both angular displacement and angular velocity decreases with time as its a damped motion. Blue line represnt the angular displacement and red dotted line represents the angular velocity of the pendulum.

Animation link:- Youtube link