Week 3 - Solving second order ODEs

AIM: To find the motion of a simple pendulum by solving the governing equation of motion.

GOVERNING EQUATIONS:

1. In general, the simple pendulum follows a simple harmonic motion. When the bob of a pendulum is subjected to motion in presence of air resistance, damping occurs. The equation of motion of such a system is given by d2θdt2+bmdθdx+glsinθ=0d2θdt2+bmdθdx+glsinθ=0
   Where
   θθ- Angular displacement (in radians)
   b- Damping factor (N.s/m)
   m- Mass of the bob suspended (kg)
   l- Length of the string (m)
   
   
   If this second order differential equation is solved, the motion of the pendulum can be traced out.

OBJECTIVE:

• To solve the second order ODE by writing a script in python
• Plot the angular displacement vs time
• Plot the angular velocity vs time
• Make an animation of the simple pendulum under given conditions

OPERATING CONDITION OF PENDULUM:

• Length= 1 m
• Mass attached= 0.1 kg
• Damping coefficient= 0.05 Ns/m
• Acceleration due to gravity= 9.81 ms2ms2
• Initial displacement= 0 m
• Initial angular velocity= 3 m/s

MODIFICATIONS TO GOVERNING EQUATION:

We have a second order ODE as the equation of motion. In order to solve this using the ode45() function in MATLAB we have to perform a few modifications to the second order ODE.

The second order ODE has to be converted to a first order ODE and solved using ode45().

To convert the ODE, we make a few assumptions and name the following:

• θ=θ1θ=θ1
• dθ1dt=θ2dθ1dt=θ2

With these assumptions the second order ODE can be written in a matrix form as:

IMAGE

The column vector containing θ1andθ2θ1andθ2 are the angular displacement and angular velocity respectively.

SYNTAX OF ode45():

The ode45 function returns an array of the dependent and independent variable. As input parameters, ode45 requires a user defined function that describes the differential equation, a time range over which the ODE is to be solved and the initial conditions of the problem.

[ (Results vector) ]= ode45(@user_defined_function, span vector, initial conditions)

CODE:

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


#defining the function for ODE
def mode1(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

#Initializing important parameters
#damping coefficient
b=0.05

#accelaration due to gravity
g=9.81

#lenght of string
l=1

#mass of bob
m=0.1

#initial condition
theta_0 = [math.pi/6,0]

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

counter=1 #for animation

#solving ODE by function call
theta = odeint(mode1,theta_0,t,args = (b,g,l,m))

plt.figure(1)
plt.plot(t,theta[:,0],'b-',label = r'$frac{dtheta_1}{dt}=theta2$')
plt.plot(t,theta[:,1],'r--',label = r'$frac{dtheta_2}{dt}=frac{b}{m}theta_2-frac{g}{l}sintheta_1$')
plt.ylabel('Plot')
plt.xlabel('time')
plt.legend(loc='best')

plt.figure(2)
plt.plot([-1,1],[0,0],lw=2)
for i in range(len(theta)):
    
    filename= str(counter) + '.png'
    counter=counter+1
    plt.figure()
    #plotting the string
    ln = plt.plot([0, np.sin(theta[i, 0])], [0, -np.cos(theta[i, 0])], color='k', lw=2)
    #plotting the bob
    bob = plt.plot(np.sin(theta[i, 0]), -np.cos(theta[i, 0]),'o',markersize=20,color='b')
    plt.xlim(-1,1)
    plt.ylim(-1.5,0)
    plt.savefig(filename)

EXPLAINATION OF CODE:

• All inputs have to be declared and assigned
• A vector containing the time over which the ODE is to be solved has to be created
• Initial conditions have to be declared and assigned
• Write the function to solve ODE and return the solution vector
• The function written for ode45() has to be called with correct arguments
• Assign the coordinates of fixed end as x0 and y0
• In a for loop, determine the position of the moving end point (x1 and y1)
• Plot the results using plot()
• Create the support by drawing a line of width 15 and black color
• Create the rope by plotting between fixed and moving coordinates
• Create a ball at the moving end using 'markersize' and 'markerfacecolor' command
• Create a pivot at the fixed end using 'markersize' and 'markerfacecolor' command
• Write the instantaneous plot into an array named M
• Create an animation of the result

RESULT:

The code has been executed and the following animation has been created using Adobe premiere pro

https://youtu.be/Vfz0Gft2WnM