Week 3 - Solving second order ODEs

SOLVING SECOND ORDER EQUATION USING PYTHON

AIM

Using second ODE to describe the transient behaviour of a system of simple pendulum on python scripting.

OBJECTIVE

In Engineering, ODE is used to describe the transient behavior of a system. A simple example is a pendulum

The way the pendulum moves depends on the Newtons second law. When this law is written down, we get a second order Ordinary Differential Equation that describes the position of the "ball" w.r.t time.

Similarly, there are hundreds of ODEs that describe key phenomena and if you have the ability to solve these equations then you will be able to simulate any of theses systems.

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

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.

Write a program in Python that will simulate the pendulum motion, just like the one shown in the start of this challenge.

use,

L=1 metre,

m=1 kg,

b=0.05.

g=9.81 m/s2.

Simulate the motion between 0-20 sec, for angular displacement=0,angular velocity=3 rad/sec at time t=0.

PROCEDURE

• integrate module is imported for ODE equation by the syntax from scipy.integrate import odeint.
• Then the variables and arguments are defined using function def mode1(theta,t,b,g,l,m):
• Then the second ODE equation is breakdown into two single ODE, the first equation will give displacement and the second equation will give velocity.
• Then the initial condition is defined by theta_0 = [ 0,3]
• To obtain the smooth curve the points are divide into 150 points between 0 to 20 sec by t = np.linspace(0,20,150)
• Then the function odeint are stored as theta.
• The displacement values are extracted and stored in null array THETA3 using append command.
• Then THETHA4 for loop is used to extract value from THETA3.

PROGRAM

# stimulation of pendulum
import math
import matplotlib.pyplot as plt 
import numpy as np 
from scipy.integrate import odeint

# Function
def mode1(theta,t,b,g,l,m):
	theta_1 = theta[0]
	theta_2 = theta[1]
	dtheta1_dt = theta_2
	dtheta2_dt = -(b/m) * theta_2 - ((g/l) * math.sin(theta_1))
	dtheta_dt  = [dtheta1_dt , dtheta2_dt]
	return dtheta_dt 
 # input

b = 0.05
g = 9.81
l = 1
m = 1

# initial condition 
theta_0 = [ 0,3]

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

# solve ODE
theta =odeint(mode1,theta_0,t,args=(b,g,l,m))
a = theta[:,0] 
# plot
plt.figure(1)
plt.plot (t, theta [:,0], 'b--',label=r'$\frac{d\theta_1}{dt}=\theta_2 $')
plt.plot (t, theta [:,1], 'r--',label=r'$\frac{d\theta_2}{dt}=-1\frac{b}{m}\theta_2-\frac{g}{L}sin\theta_1 $')
plt.xlabel('time')
plt.ylabel('plot')
plt.legend([0,12])
plt.legend (loc='best')
plt.show()


#Animation
THETA3 = []
for theta3 in theta[:,0]:
	THETA3.append(theta3)

ct = 1
for THETA4 in THETA3:
	x0 = 0
	y0 = 0
	x1 = 1*math.sin(THETA4)
	y1 =-1*math.cos(THETA4)
	filename='pendulum%05d.png'%ct

	ct=ct+1
	plt.figure(2)
	plt.plot([-0.2,0.2],[0,0])
	plt.plot([x0,x1],[y0,y1])
	plt.plot(x1,y1,'-o')
	plt.xlim([-1.5,1.5])
	plt.ylim([-1.5,1])
	plt.title('Motion of Pendulum')
	plt.savefig(filename)
	plt.show()

ERROR

Invalid Syntax error for scipy.integrate.

Used syntax as from scipy.integrate import odeint

OUTPUT

TIME PLOT

ANIMATION