Solving 2nd Order Ode - Simple Pendulum

Objectives:

• Simulate the transient behaviour of a simple pendulum and generate an animation for the same

Initial Conditions:

• Length of the pendulum = 1 m
• Mass of the ball = 1 Kg
• Damping coefficient = 0.05
• Gravity in m/s^2 = 9.81
• Simulation time = 0 - 20 seconds
• Initial position of mass = 0 Degrees
• Initial Velocity of mass = 3 m/s

Solving Second Order ODE in Python:

The transient behaviour of the pendulum is governed by the following second order ODE –

`(d^2 θ)/(dt^2 )+b/m dθ/dt+g/L sin⁡θ=0`

In order to solve the above second order ODE, we need to break it down into a set of first order ODE. To do so we first assume

`θ=θ_1`

On differentiating w.r.t time we get

`dθ/dt= (dθ_1)/dt=θ_2 ------------------------- (1)`

Again, on differentiating w.r.t time we get

`(d^2 θ)/(dt^2 )= (d^2 θ_1)/(dt^2 )= (dθ_2)/dt ---------------------------(2)`

Thus, we achieve two first order ODE’s,

• For Angular Displacement

`(dθ_1)/dt=θ_2`

• For Angular Velocity

`(dθ_2)/dt=-b/m*θ_2-g/L*sin⁡〖θ_1 〗`

Plot for Displacement and Velocity :

Code Methodology :

• Import required packages
• A function has been defined which represents the Original Second order ODE in the form of two single order ODE\'s as described above.
• The geometric inputs and constants have been defined
• To solve an ODE, it is required to provide some initial conditions, which have been set ( Intial position of mass is 0 Degrees and initial Velocity is 3 m/s)
• The time scale over which we want to study the behaviour of the pendulum has then been set.
• The variation of displacement and angular velocity w.r.t time has then been plotted
• In order to generate an animation, a counter variable has been defined initially. A for loop is then set up to take every value of theta from the angular position array and calculate the position of the centre of mass at each time point. The counter varable is then used to save a new image for each position of the pendulum. After every iteration the counter variable is updated and hence a new image is generated at each position

Link to Animation:-

# Program to simulate transient behaviour of simple pendulum
# Loash Suri

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

def model(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


# Defining constants
b=0.05
g=9.81
l=1
m=1


# Initial conditions
theta_0=[0,3]        # Setting initial position of mass at 0 degrees and setting initial velocity as 3 m/s

# Time Points
t=np.linspace(0,20,150)      # Defining time points for integration

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

# Plot the results
plt.plot(t,theta[:,0],\'b-\',label=r\'Displacement\')
plt.plot(t,theta[:,1],\'r--\',label=r\'Velocity\')
plt.legend(loc=\'best\')
plt.show()


# Animation plot
counter=1
# Fixed point
x1=0
y1=0
for theta in theta[:,0]:
	# Coordinates of Centre of mass
	x2=l*math.sin(theta)
	y2=-l*math.cos(theta)
	filename= str(counter)+ \'.png\'
	counter=counter+1
	plt.figure()
	# Drawing horizontal support of string
	plt.plot([-5,5],[0,0],\'b\',linewidth=5)
	# Drawing the circular mass
	plt.plot(x2,y2,\'o\', markersize=30)
	# Setting x axes range
	plt.xlim([-4,4])
	# Setting y axes range	
	plt.ylim([-4,2])
	
	# Drawing the pendulum string
	plt.plot([x2,x1],[y2,y1],\'r\')
	plt.savefig(filename)