2nd Order ODE Simulation in Python.

2nd Order ODE Simulation

AIM: To write a program in Python that will simulate the pendulum motion.

OBJECTIVE:

1. To solve the second order ODE in Python and plot graph for Velocity and displacment.
2. To simulate the motion of simple pendulum and show the animation

ASSUMPTIONS:

1. The value of acceleration due to gravity is (g) 9.81 m/s^2
2. The length of the pendulum is 1 meter.
3. The damping coefficient is 0.05.
4. The mass of the pendulum bob is 1 kg.

THEORY:

The motion of the pendulum can be easilty represented by the following differential equation:

`(d^2 θ)/dt^2 +b/m dθ/dt + g/L sin⁡θ = 0`

where:

`theta`= angular displacement of the bob

b = damping coefficient

g = accelration due to gravity

m = mass of the pendulum bob

L = length of the pendulum string

t = time period

`(dθ)/dt`=angular velocity

`(d^2θ)/(dt^2)` = angular acceleration

PYTHON CODE:

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

#function
def model(theta,t,b,g,l,m):
	theta1=theta[0]
	theta2=theta[1]
	dtheta1_dt= theta2
	dtheta2_dt= -(b/m)*theta2-(g/1)*math.sin(theta1)
	dtheta_dt=[dtheta1_dt,dtheta2_dt]
	return dtheta_dt

b=0.05     #damping factor
g=9.81    
l=1        #length of string
m=1        #Mass of the bob

theta_0=[0,2]    #initial conditions where angular displacement is 0 and velocity is 2m/s


t=np.linspace(0,25,125)          #line plots

theta= odeint(model,theta_0,t,args=(b,g,l,m))        # solving the ode from function by odeint function

abhi.figure()    #open new figure

abhi.plot(t,theta[:,0],'b-',label=r'$\frac{d\theta_1}{dt}=\theta_2$')     #plotting time and angular displacement


#plotting velocity and time
abhi.plot(t,theta[:,1],'r--',label=r'$\frac{d\theta_2}{dt}=\frac{b}{m}\theta_2-\frac{g}{l}sin\theta_1$')
abhi.ylabel('Plot')
abhi.xlabel('Time')
abhi.legend(loc='best')
abhi.show()


val=theta[:,0]
count=1                    #counter to count filename
for i in range(0,len(val)):
	abhi.figure()
	x1=0.5-l*np.sin(val[i])
	y1=0.5-l*np.cos(val[i])
	abhi.plot([0.5,x1],[0.5,y1])
	abhi.plot(x1,y1,'o',markersize=12)                  #bob
	abhi.plot([-1,2],[0.5,0.5],linewidth=5)             #string
	abhi.xlim([-1,2])                                   #x and y limits
	abhi.ylim([-1,2])
	filename=str(count)+ '.png'                        #saving file as .png
	abhi.savefig(filename)
	count=count+1                                       #increment filename by 1

Explanation:

1.Set physics parameters of the pendulum like mass, length, gravity and damping coefficient.
2.Set the initial conditions and timespan for which you want to run the simulation.
3.To use the ode solver, use odeint sovler from scipy module with main inputs as (t and theta) while you use model function (the one we created earlier) with inputs as (t,theta,l,b,g,m) and t and theta_0 as the initial conditions.
4.These results are stored in an array of theta
5.Since results contain 2 components, angle displacement and velocity, to plot both the graphs, use the (:,0) and (:,1) property of Python and plot those graphs.
6.To create animation, Use a for loop so that we can get X and Y co-ordinates using trigonometry, Here, 0.5 and 0.5 are the hinge points. Since they are (0.5,0.5) we have subtracted the co-ordinates of the bob from them so that proper rod of the pendulum can be seen.
7.Use the ct counter to make up for the file name using str(ct) and a '.png' file extension.
8.Use a image stitcher to stich all the images to create an animation.

CONCLUSION:

1. Python can be effectively used to solve 2nd order ODE's
2. The graph of plot versus time can be easily obtained.

For the animation of the pendulum is given below: