Solving second order ODE by using Python

Aim: to write a program to solve second order ode and to simulate the pendulum motion by using python

Objective:

1. To write a program to solve second order (ODE) and to simulate the pendulum motion
2. To understand transient behavior of pendulum
3. To make an animation

Theoretical background:

an ordinary differential equation (ODE) is a differntial equations containing one or more functions of one independent variable and the derivative of those functions, In Engineering ODE are used to describe transient behavior of system example – Newton second law, simple pendulum

To solve second order ordinary differential equations (ODE) we need to convert second order ODE into two first order ODE, suppose if we have to solve third order ODE then we need to convert that into three first order ODE.

Let’s say given equation below is second order differential equation

`(d^2 θ)/(dt^2 )+b/m (dθ)/dt+g/L sin⁡θ=0`

To convert this equation into two first order ODE

Let’s say

`(dθ)/dt=(dθ_1)/dt=θ_2`

`(d^2 θ)/(dt^2 )=(d^2 θ_1)/(dt^2 )=d/dt ((dθ_1)/dt)=d/dt (θ_2 )`

So by substituting

`(dθ_2)/dt+b/m θ_2+g/L sin⁡θ=0`

`(dθ_2)/dt=-b/m θ_2+g/L sin⁡θ`

`(dθ_1)/dt=θ_2`

So we can represent those two first order ODE in matrix format or we can say solution vector

``d/dx [θ_1/θ_2 ]= [θ_2/(-b/m θ_2+g/L sin⁡θ )]`

In the above equation,

g = gravity in [ m/s²] ,

L = length of the pendulum in [m],

m = mass of the ball in [kg],

b = damping coefficient.

Intial input:

Length of the pendulum = 1 [meter]

Mass of the ball = 0.1 [kg]

Damping coefficient = 0.02

Acceleration due to gravity = 9.81 [m/s²]

Code:

import math
import matplotlib.pyplot as plt 
import numpy as np 
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)*theta2 - (g/l)*math.sin(theta1)
	dtheta_dt=[dtheta1_dt,dtheta2_dt]
	return dtheta_dt


b=0.02
g=9.81
l=1
m=0.1

#initial condition
theta_0=[0,3]
# time points
t = np.linspace(0,20,150)
# solve ODE
theta = odeint(model,theta_0,t,args = (b,g,l,m))

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

ct=1
x1= 0
y1= 0

n=150
for i in range (0,n):

	x2=x1+l*math.sin(theta[i,0])
	y2=y1-l*math.cos(theta[i,0])
	filename='pendulum%05d.png' % ct
	plt.figure()
	plt.plot([-1.5,1.5],[0,0],linewidth = 5)
	plt.plot([x1,x2],[y1,y2],linewidth = 4)
	plt.plot(x2,y2,marker='o',markersize='10')
	plt.xlim(-1.5,1.5)
	plt.ylim(-1.5,1.5)
	plt.savefig(filename)
	ct=ct+1

 

Procedure:

• Import math,matplotlib,numpy and scipy module
• define the function which convert 2nd order ODE to two first order ODE
• define input like damping coeifficient , accelration due to gravity,length , and mass of ball
• make the time point array over which ODE is going to intergrate
• define the initial condition like angular displacement and velocity
• we have used odeint function from scipy module to intergrate ode over time span
• plot the output of odeint and time
• make the animation

Procedure to make animation:

1. run main program we will get 150 png files .
2. open command prompt, you will see working directory.
3. to change directory use cd followed by directory of all png files.
4. once your directory changed, try checking files in that location by using dir and hit enter.
5. now use magick “*.png output.gif” command to get animation.

Error Encountered:

• you need to import scipy module to use odeint function
• Indentation Error

Animation:

https://drive.google.com/file/d/1cUcAkboTxF7BuGlqnX7LUlqLnsCXfswB/view?usp=sharing

Result:

• program to solve second order ode has been successfully executed.