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 the 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
- We first consider Ɵ = Ɵ₁
- Hence,`(dƟ) /(dt)=(dƟ1)/ (dt) = Ɵ2`
- We then find Ɵ2
- `(d^2Ɵ)/(dt^2)=(d^2Ɵ1)/(dt^2) =d/(dƟ).((dƟ1)/(dt))=d/(dt).Ɵ2`
- 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
- In this program, we will first import the math library, numpy, scipy, matplotlib.pyplot.
- We are using scipy library for integration purposes.
- Then we will use the function to dene the above function of Ɵ1 and Ɵ2 to calculate them.
- Then we will dene the inputs such as damping coefficient, length of the pendulum, acceleration due to gravity, the mass of the ball. We will also use the time in the form of line spacing in 150 divisions.
- Now for solving the ode we will use odeint.
- Now we will use plot function to plot the pendulum ball.
- The png files are saved and then using an 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:https://www.youtube.com/watch?v=SNMOYw8N6ko&feature=youtu.be