Solving second order ODE using PYTHON

AIM: To write a program in python to solve ODE and to create the simple pendulum animation

OBJECTIVES OF THE PROJECT: The main objective of this project is to write a program in Python that will simulate the pendulum motion.

GOVERNING EQUATIONS: Basic math concepts like the resolution of vectors and locating coordinates the specific points in the plane.

NEWTON'S 2nd LAW

`F=ma`

F = force in newtons

m = mass of the body in kg

a =accleration in m/s2

Second order differential equation of simple pendulum is:

`(d^2θ)/dt^2 +b/(ml)(dθ)/dt +g/l sin θ =0`

• To convert this into first order ODE replace
• `θ= θ_2`
• `(d θ_1)/dt = θ_2`
• So second order ODE Equation will transform into
• `(dθ_2)/dt^2 +b/(ml)θ _2 +g/l sinθ _1 =0`

where g = gravity in m/s2,

L = length of the pendulum in m,

m = mass of the ball in kg,

t = time in seconds

θ = angular displacement rad/sec,

b = damping coefficient.

PROGRAMMING LANGUAGE USED: PYTHON

create a new file with an extension .py as ODE.py

import math
import numpy as np
import matplotlib.pyplot as plt
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*L)*theta2 - (g/L)*math.sin(theta1)

	dtheta_dt=[dtheta1_dt,dtheta2_dt]
	return dtheta_dt

#parameters
#damping coefficient
b=0.02
#acceleration due to gravity
g=9.81
#length of the pendulum
L=1
#mass of the pendulum bob
m=0.1

#defining initial condition for angular displacement and velocity
theta_0=[math.pi/2,0]

#defining the time
t=np.linspace(0,20,150)

#SOLVING THE ODE
theta=odeint(model,theta_0,t,args=(b,g,L,m))	

#plotting angular dispacement and velocity with time
plt.figure(1)
plt.plot(t,theta[:,0],'b-',label=r'$(dispacement)\frac{d\theta_1}{dt}=\theta_2$')
plt.plot(t,theta[:,1],'r--',label=r'$(velocity) \frac{d\theta_2}{dt}=-\frac{b}{m}\frac{theta_2}{L}-\frac{g}{L}sin\theta_1$')
plt.ylabel('plot')
plt.xlabel('time')
plt.legend(loc='best')
plt.show()

# ANIMATION OF SIMPLE PENSULUM 
#conssider a fixed point 
x0=0
y0=0

#extaracting the angulat displacement values
ang_s=theta[:,0]
ct=1
for alpha in ang_s:
	x1=L*math.sin(alpha)
	y1=-L*math.cos(alpha)
	#creating a new file every time 
	filename = 'simple%05d.png' % ct
	ct = ct + 1	
	#joining the fixed poinnt and moving bob coordinates
	plt.figure()
	plt.plot([x0,x1],[y0,y1],'k-',LineWidth=2)
	
	plt.axis([-2,2,-1.5,0.3])
	#creating pendulum bob
	plt.plot(x1,y1,'bo',MarkerSize=24)
	plt.plot([-1,1],[0,0],'y-',LineWidth=11)
	#referance line
	plt.plot([0,0],[0,-15],'r--')
	plt.title('ANIMATION OF SIMPLE PENDULUM')
	#saving the figures after plotting them
	plt.savefig(filename)

EXPLANATION TO THE CODE:

• Create a new file and import the modules like math, numpy and matplotlib.
• Write a function called model which takes the input arguments b,g,l,m,theta,t and returns the dtheta_dt.
• define the initial parameters with variable names as per your choice.
• solve the ordinary differential equation using the odeint which is imported from scipy.integrate.
• theta is a solution od ODE containing two-column vector
• first column is angular displacement and the second column is angular velocity.
• plot the angular displacement and angular velocity against the time.
• to create the simple pendulum animation create the plots for different angular displacements.
• store the plots in separate files and use image magick to create the animation( magick *.png output.gif).
• Minor errors were made , to avoid indentation errors use spacing at required places and make sure you are using either tab spacing or spacebar spacing don't use both in the same loop to avoid the error.

RESULTS:

the plot of angular displacement, angular velocity against time.

the Lync for an animation of the pendulum is given below

https://youtu.be/LcyOfWYE9HQ