Week 3 - Solving second order ODEs

OBJECTIVE:

By using ODE plot graph and create animation file for simple pendulum.

PROCEDURE:

Plot graph for angular displacement and angular velocity

1. Import required library

Import numpy,scipy.integrate,matplotlib.pyplot and math library.

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

2. Diffrential equation

The pendulum moves depend on the newtons second law. This ball describe position of ball w.r.t. time.

Below equation shows motion of simple pendulum with damping.

Reduce above differential equation to single degree deferential equation.

dƟ1/dt=Ɵ2…………(1)

dƟ2/dt=(-b/m)*Ɵ2 – (g/l)*sin(Ɵ)

#function that return dz/dt
def model(theta,t,b,g,l,m):
	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

3. Input parameter

b=0.05………..Damping coefficient

g=9.81………..Acceleration due to gravity

l=1………………Length of string (m)

m=1…………….Mass of bob(kg)

b=0.05	#damping coefficient
g=9.81	#acceleration due to gravity(m/s^2)
l=1		#length of string(m)
m=1		#mass of bob(kg)

4. Initial Condition

Angular displacement=0 and angular velocity=3 rad/sec.

theta_0=[0,3]

5. Time points

Used motion between 0-20 sec with 500 spcing.

t=np.linspace(0,20,500)

6. Solve ODE

ODE solved by using odeint command.

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

7. Plot result

Plot graph between time vs angular displacement and time vs angular velocity.

#plot results
plt.plot(t,theta[:,0],label='Angular Displacement')
plt.plot(t,theta[:,1],label='Angular Velocity')
plt.legend()
plt.xlabel('time')
plt.ylabel('Angular Displacement and Angular Velocity')
plt.savefig('plot')

ANIMATION PLOT

1. Input Parameter

Use initial coordinate of (x,y) at (0,0) and counter variable ct=1.

#animation plot
#Input parameter
ct=1	#counter variable
x0=0	#initial x coordinate of bob
y0=0	#initial y coordinate of bob

2. Tracing path of bob

Use for loop to trace angular position of bob. For this tracing, plot (x,y) co-ordinate at each iteration of angular displacement.

Counter variable create to give filename.

Large no of images crated due to for loop, to avoid this issue use if loop to create image up to 54 .

for i in theta[:,0]:
	x1=(l*math.sin(i))
	y1=(-l*math.cos(i))

	filename='Pendulum' + str(ct) +'.png'
	ct=ct+1
	if ct==56:
		break

3. Plot Result and create animation file:

By using plt.plot( ) command plot position of bob at different angular position by using (x,y) coordinate. Create animation file from images.

plt.figure()
plt.plot([-1,1],[0,0],"black",linewidth=5)
plt.plot([x0,x1],[y0,y1],"red")
plt.plot(x1,y1,'o',markersize=15)
plt.xlim([-2,2])
plt.ylim([-2,2])
plt.savefig(filename)

RESULT:

Final code:

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

#function that return dz/dt
def model(theta,t,b,g,l,m):
	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


#input
b=0.05	#damping coefficient
g=9.81	#acceleration due to gravity(m/s^2)
l=1		#length of string(m)
m=1		#mass of bob(kg)

#initial condition
theta_0=[0,3]

#time points
t=np.linspace(0,20,500)

#solve ODE
theta=odeint(model,theta_0,t,args=(b,g,l,m))

#plot results
plt.plot(t,theta[:,0],'-','red',label='Angular Displacement')
plt.plot(t,theta[:,1],'--','blue',label='Angular Velocity')
plt.savefig('plot')


#animation plot
#Input parameter
ct=1	#counter variable
x0=0	#initial x coordinate of bob
y0=0	#initial y coordinate of bob

for i in theta[:,0]:
	x1=(l*math.sin(i))
	y1=(-l*math.cos(i))

	filename='Pendulum' + str(ct) +'.png'
	ct=ct+1
	if ct==56:
		break

	plt.figure()
	plt.plot([-1,1],[0,0],"black",linewidth=5)
	plt.plot([x0,x1],[y0,y1],"red")
	plt.plot(x1,y1,'o',markersize=15)
	plt.xlim([-2,2])
	plt.ylim([-2,2])
	plt.savefig(filename)

Angular displacement and angular velocity plot:

Animation: