Week 3 - Solving second order ODEs

Simulation of Simple pendulum by Solving second order ODE using Python

Aim: Creating a program to solve second order ODE using Python

Objective:

• Create a program to solve the second order ODE.
• Plot the results.
• Use these results to simulate a simple pendulum.

Theory:

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 Newton’s second law. When this law is written down, a second order Ordinary Differential Equation is obtained that describes the position of the ‘ball’ w.r.t time.

In this program following ODE is solved. Some substitutions are made to solve the equation easily.

Next, variables such as b (Damping coefficient), g (Gravity), l (Length of pendulum) and m (Mass of pendulum), angular displacement & velocity and time are specified. A function is created to solve the terms in ODE. Actual ODE is solved by using ‘odeint’ command which solves complex equations in Python.

Then, the obtained values are plotted. For simulation of pendulum, coordinates are defined and the values are plotted. For creating animation external software is used.

Program:

RESULT:

# A program to solve the second order ode and simulate the motion of pendulum

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

# Input data

# Damping coefficient
b = 0.05   
# Gravity in m/s^2
g = 9.81             
# Length of pendulum in m
l = 1    
# Mass of pendulum in kg
m = 1	



# Angular displacement
theta0 = [0,3] 
# Time        		
t = np.linspace (0,20,150)    

# Function
def graph(theta,t):

	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

# Solving ODE
theta = odeint (graph,theta0,t)

# Plotting graph of ODE
plt.plot(t,theta[:,0],label = 'Displacement')			
plt.plot(t,theta[:,1],label = 'Velocity')
plt.xlabel('Time in sec')
plt.ylabel('Angular displacement and Velocity')
plt.legend()
plt.savefig('graph.png' + '.png')
plt.show()

# Count			
ct = 0

# Loop for pendulum animation
for i in theta[:,0]:

	# Defining coefficients			
	X1 = l*math.sin(i)
	Y1 = -l*math.cos(i)

	# Specifying file name
	filename = 'test%04d.png' % ct

	# Increasing count
	ct = ct+1

	# Plotting pendulum graph
	plt.figure()
	plt.plot([0,X1], [0,Y1])
	plt.plot([-0.5,0.5], [0,0], linewidth = 5)
	plt.plot(X1,Y1, marker='o', markersize=20)
	plt.xlim([-1,1])
	plt.ylim([-1.5,0.1])
	plt.title('Motion of simple pendulum')
	plt.savefig(filename)
	
# End

Angular displacement and velocity vs Time

Motion of simple pendulum

Conclusion:

It is observed that by increasing the value of angular displacement, swing of the pendulum is increased. The motion of the pendulum decreases w.r.t time, this happens due to the damping coefficient.