Week 3 - Solving second order ODEs

Aim:  To observe the transient nature of the Simple Pendulum using Python.

Theory:  In this problem we consider a Simple Pendulum of mass ‘m’ attached to a string having length ‘l’ and the other end being fixed. ‘b’ is the Damping Coefficient and ‘g’ is acceleration due to gravity.

The mathematical equation for the motion of Simple Pendulum is given by an Ordinary Differential Equation;

`(d^2 Ɵ)/(dt^2 )+b/m*(dƟ)/dt+g/l*sin⁡θ`

Procedure: Steps to solve the problem

1. We first consider Ɵ = Ɵ₁

Hence,`dƟ/dt=(dƟ1)/dt= Ɵ2`

2. We then find Ɵ2

`(d^2 Ɵ)/(dt^2 )=(d^2 Ɵ1)/(dt^2 )=d/dt* ((dƟ1)/dt)=d/dt*Ɵ2`

3. We then store it in a 2 columns in which the first column will consist of angular displacement values while the second columns of angular velocity values.

4. In this program we will first import the math library, numpy, scipy, matplotlib.pyplot.

5. We are using scipy library for the integration purpose.

6. Then we will use the function to define the above function of Ɵ1 and Ɵ2 to calculate them.

7. Then we will define the inputs such as damping coefficient, length of pendulum, acceleration due to gravity, mass of ball. We will also use the time in the form of linspacing in 150 divisions.

8. Now for solving the ode we will use odeint.

9. Now we will use plot function to plot the pendulum ball.

10. The png files are saved and then using online video converter we will convert them to video format.

Program Code:

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

# function
def model(T,t,b,l,g,m):
	T1= T[0]
	T2 = T[1]
	dT1_dt = T2
	dT2_dt = -(b/m)*T2 -(g/l)*math.sin(T1)
	dTh_dt = [dT1_dt, dT2_dt]
	return dTh_dt


# inputs
# damping co-efficient (Ns/m)
b = 0.05 

# length of pendulum (m)
l = 1 

# gravity (m/s^2)
g = 9.81 

# mass of ball (kg)
m = 1 

# initial conditions
t_x = [0,3]

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

# solve ode
T = odeint(model,t_x,t,args=(b,l,g,m))

# plot
X = T[:,0]
ct = 1

for theta in X: #Using for loop to calculate the values over a range
	x1 = 0
	y1 = 0
	x2 = l*math.sin(theta)
	y2 = -(l*math.cos(theta))

	#Saving the file as png
	file = 'test%05d.png'%ct
	

	#PLotting the Data
	plt.figure()
	plt.plot([-1.25,1.25],[0,0],linewidth = 4)
	plt.plot([x1,x2],[y1,y2],linewidth = 3)
	plt.plot(x2,y2,'o',markersize =10)
	plt.xlim([-1.25,1.25])
	plt.ylim([-1.25,1.25])
	plt.savefig(file)
	plt.show()
	ct = ct+1

Results and Plots:

Sample image Simple Pendulum:

Animation Video: