Week 3 - Solving second order ODEs

Aim:  To write a Python program which solves the equation of motion of a simple pendulum with damping.

Objective:

• To learn how to solve an Ordinary Differential Equation (ODE) using Python
• To solve ODE of a simple pendulum that represents its equation of motion
• To simulate the pendulum motion
• Understand the importance behind motion study of a simple pendulum

Problem Statement :-

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

def model(initial,t,b,l,g,m):
	t1=initial[0]
	t2=initial[1]
	dt1=t2
	dt2=(-b*t2/m)-(g*math.sin(t1)/l)
	dt=[dt1,dt2]
	return dt

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

initial=[0,3]
t=n1.linspace(0,20,200)
disvelo=odeint(model,initial,t,args=(b,l,g,m))

plt.plot(t,disvelo[:,0],'-',label='displacement')
plt.plot(t,disvelo[:,1],'--',label='velocity')
plt.legend()
plt.xlabel('time')
plt.ylabel('displacementvelocity')
plt.show()

y=disvelo[:,0]
ct=1
for ang in y :
	x1=0;y1=0
	x2=l*math.sin(ang)
	y2=l*math.cos(ang)
	file='test%05d.png'%ct 
	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(file)
	ct=ct+1

Solution Procedure –

First of all, we import the libraries that are going to be used in this program. ‘math’ module for basic math commands, ‘matplotlib.pyplot’ module for plotting purposes, ‘numpy’ module for arrays and ‘odeint’ module from ‘scipy.integrate’ for integrating our second order ODE.

Before assigning the input values we write a function named ‘model’ with inputs as ‘theta’, ‘time’, ‘b’, ‘g’, ‘l’ and ‘m’. This function breaks our second order ODE into 2 first order ODEs and returns the value of dtheta_dt.

Then we assign the input values of damping coefficient (b), gravitational acceleration (g), length of pendulum (L) and mass of the bob (m).

After that we define a variable theta_0 which contains the initial angular displacement = pi/2 and initial angular velocity = 0. After that we assign the value of time span as 0 to 10 sec by the help of linspace command with 150 intervals.

Now to solve the ODE, we use the command ‘odeint’ which is going to integrate the ODEs and inside the command, we call the function ‘model’ with all the inputs and the arguments.

Then we use the plot command to generate the plot between displacement and time and velocity and time.

Now to create an animation of pendulum, first we define the initial values of the coordinates, which is going to be a fixed point. Then we used a ‘for’ loop with loop counter variable (i) with size equal to the length of time (t) array. Also count (ct) variable is also addressed to create separate file for each case.

Then we define the position of the end point of the string. To create the mass at one end of the pendulum, we use ‘o’ inside the plot command for end point of string, then we choose the color and choose the size by using ‘markersize’ command.

To name the file in a definite order, we create a variable named ‘filename’ with str of ct and extension as ‘.png’ and use ‘savefig’ command to save the files.

Now when we run the code, all 150 images are saved in the working folder. We stitched these images together to create the required GIF file.

 Result –

1. Displacement and velocity plot with respect to time-

2.  GIF of animation of Simple Pendulum –

3.Youtube video link :- https://youtu.be/eU_E4OkcLEw