Week 3 - Solving second order ODEs

Objective:- In this project we will create a python program to solve the ODE representing the equation of motion of a simple pendulum with damping, Plot the graph of angular displacement and angular velocity with respect to time. And finally simulate the pendulum motion.

Solution Procedure:-

Lets take the ODE equation in case of simple pendulum-

In this equation

b = damping coefficient =0.05

m = mass of bob attached to string (kg) =1kg

L= length of string (in meter) = 1m

For obtaining the solution of 2nd order ODE, we have to convert it into two first order ODE equations.

T1= T[0]

T2 = T[1]

dT1_dt = T2

dT2_dt = -(b/m)*T2 -(g/l)*sin(T1)

dTh_dt = [dT1_dt, dT2_dt]

here T(0) = theta

T(1) = theta1

Python Script:-

Code for Plot the Graph between angular displacement and angular velocity with respect to time-

#Python code to Plot the Graph of angular displacement and angular velocity with respect to time

import numpy as np 
from scipy.integrate import odeint
import math
import matplotlib.pyplot as plt 

plt.rc('figure', max_open_warning = 0) #More than 20 figures have been opened Warning Removed

#funcion that returns dtheta/dt
def function(theta,t,b,g,l,m):
    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

b = 0.05 #damping Coefficient
g = 9.81
l = 1    #Lenth of string
m = 1    #Mass of bob

#initial conditions
theta_0 = [0,3]


#Time array
t = np.linspace(0,20,200)

#odeint function to solve ode

theta = odeint(function,theta_0,t,args = (b,g,l,m))

#plot the Graph of angular displacement and angular velocity with respect to time
plt.plot(t,theta[:,0],'r-.',label=r'Angular Disp')
plt.plot(t,theta[:,1],'g-',label=r'Angular Velocity')
plt.xlabel('Time')
plt.ylabel('Plot')
plt.legend()
plt.savefig('out.png')
plt.show()

Output:-

Code for Animation:

#Python code to create animation

import numpy as np 
from scipy.integrate import odeint
import math
import matplotlib.pyplot as plt 

plt.rc('figure', max_open_warning = 0) #More than 20 figures have been opened Warning Removed

#funcion that returns dtheta/dt
def function(theta,t,b,g,l,m):
    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

b = 0.05 #damping Coefficient
g = 9.81
l = 1    #Lenth of string
m = 1    #Mass of bob

#initial conditions
theta_0 = [0,3]


#Time array
t = np.linspace(0,20,200)

#odeint function to solve ode

theta = odeint(function,theta_0,t,args = (b,g,l,m))

ct = 1
x0 = 0
y0 = 0

for i in theta[:,0]:
    
    x1 = (1*math.sin(i))
    y1 = (-1*math.cos(i))            
    filename = str(ct)+'.png'
    ct = ct +1    
    plt.figure()
    plt.plot([-1,1],[0,0],"Gray",linewidth=10)
    plt.plot([x0,x1],[y0,y1],"Green")
    plt.plot(x1,y1,'o',markersize = 20)
    plt.xlim([-1.5,1.5])
    plt.ylim([-1.5,1.0])
    plt.savefig(filename)

Animation:-