Week 3 - Solving second order ODEs

Aim: To solve the ODE of Simple Pendulum Motion and animate the simulation of a pendulum using Python.

Theory: A Simple pendulum consists of a bob or a point mass attached to a string(of negligible mass) suspended to a fixed point. In ideal conditions, it follows a SHM i.e. when given a small displacement from its mean position(or equilibrium position), a restoring force(here its weight i.e. mgsin) accelerates the bob in the opposite direction of its initial displacement i.e. towards the mean position. At an extreme position, its potential energy gets converted to kinetic energy which gives it sufficient velocity to return to its mean position. But at equilibrium, all its potential energy is being converted to kinetic energy and thus has a maximum velocity which further displaces the bob from its mean position but this time in the opposite direction and can be easily understood as shown below.

The forces acting on pendulum in case of ideal conditions is shown below

Now, this motion is only in ideal conditions where we neglect the effect of air drag(i.e. air friction). in the case of air drag which acts as damping force the velocity reduces if not given constant force and this motion can be represented by the following equation which is derived using Newton's 2nd law,

`(d^2theta)/dt^2 + (b/m)(d theta)/dt + g/lsintheta =0`

where,

g = gravity in m/s2,

L = length of the pendulum in m,

m = mass of the ball in kg,

b=damping coefficient.

Now, in order to solve the ODE in Python we first break it into two single order equations,

let,

`theta = theta 1`

`(d theta)/dt = (d theta1)/dt = theta2`

Inserting the values from above into the equation, we get

`(d theta2)/dt + (b/m)theta2 + g/l sintheta1 =0`

thus the two ODEs are as follows,

Displacement `theta 2`

Velocity `(d theta2)/dt`

The code for simulating the above ODE using the given information is as follows

L=1 metre,

m=1 kg,

b=0.05.

g=9.81 m/s2.

import numpy as np
from scipy.integrate import odeint

import matplotlib.pyplot as plt
import math

# function that returns dz/dt
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

b=0.05
g =9.81
l=1
m=1

#initial points
theta_0= [0,3]

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

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

#plot results
plt.plot(t,theta[:,0], 'b-',label=r'$\frac{d\theta_1}{dt}=\theta2$')
plt.plot(t,theta[:,1], 'r--',label=r'$\frac{d\theta2}{dt}=-\frac{b}{m}\theta_2-\frac{g}{L}sin\theta_1$')
plt.ylabel('Plot')
plt.xlabel('Time')
plt.legend(loc='best')
plt.show()


#for creating the images of different position of simple pendulum 
ct=1
y=theta[:,0]

for d in y:
	x1=0;y1=0

	x2=-l*math.sin(d)
	y2=-l*math.cos(d)

	filename = 'test%05d.png' % ct
	ct=ct+1
		
	#for saving the figure in the folder 
	plt.figure()
	plt.plot([-1.5,1.5],[0,0],linewidth=8)
	plt.plot([x1,x2],[y1,y2],linewidth=3)
	plt.plot(x2,y2,'o',markersize=15)
	plt.xlim([-1.5,1.5])
	plt.ylim([-1.5,1.5])
	plt.savefig(filename)

Result & Explanation

The plot generated from the above simulation is shown below

From the figure, one can see that the initial displacement value at t=0 is 0 and the velocity is 3m/s

for damping coefficient =0.05 and time range of 20 seconds, the graph is damping towards rest.

For creating the simulation video, the 'for' loop is used to generate different positions of the pendulum which is later saved as a .png file and using the following command in the command prompt, the gif is generated.

The link to the simulation video is,

Simple Pendulum Simulation Using Python

Conclusion

The code for generating a simulation of ODE for Simple Pendulum motion is completed as per the objective.