Week 3 - Solving second order ODEs

Objective: To see a simulation of the simple pendulum by solving second-order ordinary differential equations using Python.

Theory: This ODE represents the equation of motion of a simple pendulum with damping

.............(1)

In the above equation,

g = gravity in m/s2,

L = length of the pendulum in m,

m = mass of the ball in kg,

b=damping coefficient.

Now, we cannot solve second order ODE in python so, for that, we convert this second order ODE into the first-order ODE which shown as follows:-

Let, say θ=θ1

Then equation(1) tends to

.............(2)

And, if θ1 change with respect to time we can say that

..............(3)

Also, we can say that

..............(4)

Now,by substituting equation(3) & (4) in equation (2),we get

................(5)

..............(6)

By, using equation (3) & (4) which is of first-order we will simulate simple pendulum by using Python.

Code:

In the following code, first, we have import all the necessary modules that are needed to solve ode after that we have allocated the useful parameter.

After then we have defined the ode function, in which we have created the matrix for dtheta_dt so, that we can calculate displacement and velocity. After that, we will call this function in the odeint line so that we can solve ode equation.

And then we plotted time versus displacement and time versus velocity in the same figure using subplot command.

Also, we have created the animation of the pendulum for the time span.

#Importing the modules which is required to solve ode
import numpy as np
from scipy.integrate import odeint
import matplotlib.pyplot as plt
import math

#Defining the ode function
def model(theta,t,b,m,l,g):

	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

#Input parameter;
t=np.linspace(0,20,150) #Simulation time
theta_0=[0,3]  #angular displacement and velocity
b=0.05 #damping coefficient
m=1 #mass of bob
l=1   #length of pendulum
g=9.81 #gravity

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

#Plotting angular displacement & angular velocity plot
plt.figure()
plt.plot(t,theta[:,0],color='red',linewidth=1)
plt.plot(t,theta[:,1],color='orange',linewidth=1)
plt.legend(['Angular displacement','Angular velocity'])
plt.show()


#Creating animatiom
ct=1

for i in theta[:,0]:
    x=-1*math.sin(i)
    y=-1*math.cos(i)
   
    filename='test%05d.png' % ct
    ct=ct+1
    plt.figure()
    plt.plot([-0.75,0.75],[-0.2,-0.2],color='black',linewidth=4) #allocting co-ordinates for fixed support
    plt.plot([0,x],[-0.2,y],color='red',linewidth=3)#allocting co-ordinates for pendulum string
    plt.plot(x,y,'o',markersize=13,Markerfacecolor='orange') #allocting co-ordinates for BOB
    plt.xlim([-1,1])
    plt.ylim([-1.5,0])
    plt.savefig(filename)

Output:-

Observation: By observing the above figure its quite clear that as time increases displacement and velocity decrease, because of the damping effect.

Simulation of the pendulum: