Solving second order ODE using Python

Objective : To Write a program to simulate the transient behaviour of a simple pendulum and to create an animation of it's motion.

Theory :

A simple pendulum consists of a mass 'M' attached to the one end of a light inextensible string of length 'l' , and the other end of the string attached to a fixed point.

When the bob in equilibrium is subjected to an external force it displaces from its initial position and oscillates in a circular arc in the same vertical plane,

This motion can be described in mathematical form and a code can be written in Python to get the actual representation of motion.

In Engineering, ODE is used to describe the transient behavior of a system. A simple example is a pendulum

The way the pendulum moves depends on the Newtons second law. When this law is written down, we get a second order Ordinary Differential Equation that describes the position of the bob w.r.t time.

Mathematical Equation :

The following is the mathemaical equation :

In the above equation,

g = acceleration due to gravity in m/s2,

L = length of the pendulum in m,

m =mass of the ball in kg,

b=damping coefficient.

Simulating the motion between 0-20 sec, for angular displacement=0,angular velocity=3 rad/sec at time t=0.

CODE :

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


def model(theta,t,b,g,m,l):

	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 = 0.1


theta_0 = [0,3]
t= np.linspace(0,20,400)


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


#Plot results
plt.plot(t,theta[:,0],'g' ,label=r'displacement')                                             # g -> green 
plt.plot(t,theta[:,1],'b' ,label=r'velocity')                                                 # b -> blue
plt.xlabel('plot')
plt.ylabel('Time')
plt.legend()
plt.savefig('output.png')
plt.show()


#animation
ct = 1
for theta in theta[:,0] :
	x0 = 0
	y0 = 0                                                         # x0 y0 fixed end of pendulum
	x1 =  l*math.sin(theta)
	y1 = -l*math.cos(theta)                                        # using sin and cos function to locate the cordinates of mass
	plt.figure()
	plt.plot([-1,1],[0,0],"black",linewidth=5)                    # horizontal support
	plt.plot([x0,x1],[y0,y1],"red")                                # pendulum string by connecting 2 points
	plt.plot(x1,y1,'o',markersize=8)                              # x1 , y1 = center of mass / pendulum representation  !
	plt.xlim([-1,1])
	plt.ylim([-1.5,1])
	filename = str(ct) + '.png'
	plt.savefig(filename)
	ct = ct +1

plot :

Animation link :

https://youtu.be/JQrEcIn8MeU

Conclusion : from the graph it can be seen that the values of angular displacement and angular velocity keeps on decreasing with time , this effect is know as damping effect.