Solving Second Order ODE in Python

Aim: To write a program using Python as a scripting language to simulate the transient behavior of a simple pendulum and to create an animation of its motion.

Introduction:

In engineering, an ODE is used to describe the transient behaviour of a system. A simple example is a pendulum.

Newton's second law determines the movement of a pendulum. When this law is written down, we get a second order Ordinary Differential Equation which describes the position of the 'ball' w.r.t time.

Here, the position(theta) of the bob w.r.t time is determined by using the second-order ODE which is;

where, b= damping coefficient

g= gravity in m/s^2

m= mass of the ball in Kg

L= length of the pendulum in m

In order to solve the second-order ODE, we convert it into a number of first-order ODEs. We will then create a separate function where we will assume two variables; theta1 and theta2.

theta2 will be the first-order differential of theta1 and theat2= theta. Now, by using the equation given we will make a new relation for theta2, where theta2=-b/m*theta-g/l*sin(theta1).

CODE:

#Simulation of motion of Pendulum
#Importing the basic libraries/ modules
import math
import matplotlib.pyplot as plt
import numpy as np
from scipy.integrate import odeint

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

#Inputs 
b=0.05   #damping coefficient
g=9.81   #gravity in m/s^2
l=1      #length of the pendulum in m 
m=1    #mass of the ball in Kg

#Initial Conditions
theta_0=[0,3]

#Time Points 
t=np.linspace(0,20,150)

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

#Plotting
plt.plot(t,theta[:,0],'b--',label='Displacement')
plt.plot(t,theta[:,1],'r--',label='Velocity')
plt.xlabel('Time')
plt.ylabel('Plot')
plt.xlim([0,20])
plt.legend(loc='best')
plt.show()

#Animation
ct=1
x0=5
y0=5
for i in range(len(t)):
	x1=x0+l*math.sin(theta[i,0])
	y1=y0-l*math.cos(theta[i,0])
	filename='pendulum'+str(ct)+'.png'


	plt.figure()
	plt.plot([4,6],[5,5],'k-')                         #Horizontal support of pendulum
	plt.plot([x0,x1],[y0,y1])                          #String of pendulum
	plt.plot(x1,y1,marker='o',markersize='10')         #Bob tied at the end
	plt.axis([3.5,6.5,3,6])
	plt.savefig(filename)
	ct=ct+1

EXPLANATION OF THE CODE:

1. Importing all required modules and use ode solver from scipy module. odeint is used to solve the ode function.

2. Set the input parameters like mass, length etc.

3. Also, set the initial conditions like angle, time points. Simulate the motion for 0 - 20s for angular displacement = 0; angular velocity = 3 rad?sec at time t=0.

4. Now, calculate the velocities and the angle of displacement.

5. Plot the graphs and save all the figures. Here, plot command is made use for plotting the graphs had savefig command is used to save the images of the figures generated.

6. Now finally, create an animation using all the figures to depict the simulation of the simple pendulum. To create the animation we make use of for loop.

OUTPUT:

Successfully plotted a graph of Time vs Plot for a simple pendulum and also generated an animation depicting the motion of a simple pendulum.

ANIMATION: