Week 3 - Solving second order ODEs

AIM : To write a code and solve the ODE which represents the equation of motion of the simple pendulum on PYTHON.

OBJECTIVE : The objective is to write a program that solves the following ODE. This ODE represents the equation of motion of a simple pendulum with damping

`(d^2theta)/(dT^2) + (b)/(m)(d theta)/dt + (g)/(l) sintheta = 0`

EXPLANATON : A simple pendulum consists of a ball (point-mass) m hanging from a (massless) string of length L and fixed at a pivot point P. When displaced to an initial angle and released, the pendulum will swing back and forth with periodic motion. By applying Newton's second law for rotational systems, the equation of motion for the pendulum may be obtained

τ=Iα⇒ `−mgsinθL=(mL^2)(d^2θ)/(dt^2)`

and rearranged as

`(d^2θ)/(dt^2)+((g)/(L))sinθ=0` (without damping force)

`(d^2theta)/(dT^2) + (b)/(m)(d theta)/dt + (g)/(l) sintheta` (with damping force)

• The fundamental idea solving ODE is analytically and computational is when we have higher ODE ,we break down into set of first order ODE

`(d^2theta)/(dT^2) + (b)/(m)(d theta)/dt + (g)/(l) sintheta = 0`

Where θ = θ₁

`(d theta)/(dt^2) = (d theta1)/(dt) = θ2` ------------------------(1)

`(d^2 theta)/(d t^2) = (d^2 theta 1)/(d t^2) = (d)/(dt) ((d theta 1)/(d t)) = (d)/(dt) (theta 2)`----------------------(2)

`(d^2theta)/(dt^2) + (b)/(m)(d theta)/dt + (g)/(l) sintheta` = 0

`(d^2theta 1)/(dt^2) + (b)/(m)(d theta 1)/dt + (g)/(l) sintheta 1` = 0

`(d theta 2)/(dt) + (b)/(m)( theta 2) + (g)/(l) sintheta 1`= 0

• We have created two ODE’S

1. `(d theta 1)/(d t) = theta 2`
2. `(d theta 2)/(dt) = - (b)/(m)( theta 2) - (g)/(l) sintheta 1`

This is breaking down into 2^nd order ODE into series of 1^st order ODE

`(d)/(dt) ((theta1),(theta 2)) = ((theta 2),(- (b)/(m)( theta 2) - (g)/(l) sintheta 1))`

INPUTS : This ODE represents the equation of motion of simple pendulum with damping

`(d^2theta)/(dT^2) + (b)/(m)(d theta)/dt + (g)/(l)sintheta = 0`

In above equation,

g = gravity in m/s2

L = length of the pendulum in m,

m = mass of the ball in kg,

b = damping coefficient.

• Inputs of the matlab program
• L = 1 metre,
• m = 1 kg,
• b = 0.05.
• g = 9.81 m/s2
• motion between 0-20sec,
• angular displacement = 0
• angular velocity = 3 rad/sec at time t = 0

PROGRAM :

#solving 2nd order ODE 
import numpy as np 
from scipy.integrate import odeint
import matplotlib.pyplot as plt 
import math

#function that returns dtheta/dt
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    #damping co-eff
g=9.81    #gravity
l=1       #length of pendulum
m=1       #mass of the pendulum
#initial conditions
theta_0=[0,3]
#timeplots
t=np.linspace(0,20,200)

#solving ode
theta=odeint(model,theta_0,t,args=(b,g,l,m))
#plotinng results
plt.figure()
plt.plot(t,theta[:,0],'g-',label=r'$\frac{d\theta_1}{dt}=\theta 2 $')
plt.plot(t,theta[:,1],'b--',label=r'$\frac{d\theta2}{dt}=-\frac{b}{m}\theta_2-\frac{g}{l}sin\theta_1 $')
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],"black",linewidth=10)
	plt.plot([x0,x1],[y0,y1])
	plt.plot(x1,y1,marker='o',markersize='20')
	plt.axis([3.5,6.5,3,6])
	plt.savefig(filename)
	ct=ct+1

PROCEDURE :

• Import module numpy as np ,scipy.intedrate import odeint for second ODE ,math and for plots matplotlib.plt
• Define function that returns the value of dtheta/dt
• initalize the input values of b,g,m,l
• define initial condition form theta i.e 0 and 3 rad/sec
• Use np.linspace for creating numeric sequence, np.linspace(start point,last point,number of values take in between start and end point).
• solve ODE by use of odeint commond.
• plot result by using plt.plot commond to plot graph
• For animation use loop and name the filename
• Use plt.plot command to create every single point and after that stiches all the pics and create a
• gif

ERRORS : Dimension error

OUTPUT :

PLOTS :

https://youtu.be/1EdDEaejD7w

CONCLUSION : The above graph shows the angular displacement and angular velocity and it changes over the time due to the effect of the damping, the values of reduce angular displacement and angular velocity with time as it is clearly visible in above plot