Week 3 - Solving second order ODEs

Solving second order ODE

OBJECTIVE:

To write a program to simulate the transient behaviour of a simple pendulum and to create an animation of it is motion.

In the above equation,

g = gravity in m/s2,

L = length of the pendulum in m,

m = mass of the ball in kg,

b=damping coefficient.

Now for solving the second order ODE, we will first convert this second order ODE into a number of first order ODEs. we will write a separate function where we will assume two variables theta1 and theta2, where theta2 is the first order differential of theta1 and theta2=theta. Now using the equation given we will make an equation for theta2 where, theta2=-b/m*theta1- g/l*sin(theta1)

(d^2 θ/dt^2) +(dθ/dt)*b/m+(g/l)*sin θ=0

θ=θ1

(d^2 θ1/dt^2 )+(dθ1/dt)*b/m+(g/l)*sinθ1=0

dθ1/dt= θ2

d^2 θ1/dt^2=dθ2/dt

θ2/dt+(b/m)*θ2+(g/l)*sinθ1=0

θ2/dt=-(b/m)*θ2-(g/l)*sinθ1

dθ1/dt=θ2

STEP BY STEP:

Initially we are importing required modules using 'import' command.

Function model is created with required arguments, inside which the term which are broken to solve second order ODE are calculated.

Inputs are provided for the function.

theta_0 ranges between 0 to 3 m/s^2.

odeint command is used to solve ode function.

Then results are plotted using plotting command.

For creating animation we open a for loop, ranging from 0 to length of t.

Images are saved using savefig command.

CODE FOR SIMPLE PENDULUM:

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

#function
def model(theta,t,b,g,m,l):
	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 required
b=0.05
g=9.81
l=1
m=1
theta_0=[0,3]
t=np.linspace(0,20,200)

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

#plotting result
plt.figure()
plt.plot(t,theta[:,0],'g',label=r'Angular displacement')
plt.plot(t,theta[:,1],'b',label=r'Angular velocity')
plt.xlabel('plot')
plt.ylabel('Time')
plt.legend(loc='best')
plt.show()

#Animation step
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-')
	plt.plot([x0,x1],[y0,y1])
	plt.plot(x1,y1,marker='o',markersize='10')
	plt.axis([3.5,6.5,3,6])
	plt.savefig(filename)
	ct=ct+1

OUTPUT:

https://drive.google.com/drive/folders/11x5ogJ8AJ4sbUDEm7pxWsLkoidkMOVxC