Solving Second Order ODEs using Python

AIM

Simulation of the ,motion of a simple pendulum by solving the second order ODE and to create an animation of its motion.

OBJECTIVES

• To simulate the motion of the pendulum by solving the second order ODE for motion between 0-20 seconds for initial angular displacement = 0, and initial angular velocity =3 rad/s.
• To create an animation of the motion

THEORY

The ordinary differential equation for the motion of simple pendulum is given by :

where

θ = Angular displacement (in radians)

b = Damping coefficient (in Nm/s)

m =Mass of the ball (in kg)

g =Acceleration due to gravity (in m/s²)

l =Length of the string (in m)

Steps to solve the second order ODE

1.Consider θ=θ₁

2.Hence

3.Now we get two first order ODEs

This vector on solving gives the value of θ₁ and θ₂

INPUT DATA

Initial angular displacement(θ at t=0) = 0,

Initial angular velocity (dθ/dt at t=0) = 3 rad/sec,

Damping coefficient, b =0.05 Nm/s,

Mass of the ball, m =1 kg,

Length of the string, l =1 m,

Acceleration due to gravity, g =9.81 m/s²

Time period t =0 to t =20 sec

CODE

Programming Language used- PYTHON

#A program to solve the 2nd order ODE for motion of pendulum and simulate its motion
#Importing required modules
import math
import matplotlib.pyplot as plt
import numpy as np
from scipy.integrate import odeint

#Input parameters
#Damping Coefficient(in Ns/m)
b=0.05
#Acceleration due to gravity(in m/s^2)
g=9.81
#Length of the string(in m)
l=1
#Mass of the ball(in kg)
m=1

#Initial condition(angular displacement=0 and angular velocity=3 rad/s)
theta_0=[0,3]

#Time span
time_span=np.linspace(0,20,150)

def model(theta,time_span,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

#Solving the ODE
theta=odeint(model,theta_0,time_span,args=(b,g,l,m))

#Plotting the Angular displacement and Angular velocity
plt.figure()
plt.plot(time_span,theta[:,0],'b-')
plt.plot(time_span,theta[:,1],'r--')
plt.legend(['Angular displacement','Angular velocity'])
plt.xlabel('Time')
plt.ylabel('Angular displacement and Angular Velocity')
plt.title('Plot for Angular displacement and Angular velocity')
plt.show()

#creating animation
ang_dis=theta[:,0]
#counter
ct=1
for i in ang_dis:

	#starting coordinates
	x0=0
	y0=0
  
	#ending coordinates
	x1=l*math.sin(i)
	y1=-l*math.cos(i)

	#plotting displacement and velocity
	plt.figure()
	plt.plot([-0.5,0.5],[0,0],linewidth=5)
	plt.plot([x0,x1],[y0,y1],linewidth=3)
	plt.plot(x1,y1,'o',markersize=15)
	plt.xlim([-1.5,1.5])
	plt.ylim([-1.5,1.5])
	filename='pendulum%05d.png'%ct
	plt.savefig(filename)
	ct=ct+1

PROCEDURE

1. First the required modules(i.e math, matplotli.pyplot, numpy and scipy) are imported.
2. Then the input parameters(i.e. b,g,l,m) are declared.
3. The time array is defined from 0 to 20 sec using ‘linspace()’ command taking 150 values.
4. Then the function for the ODE is defined as a set of two first order ODEs.
5. The ‘odeint’ command is used on the function passing the required arguments to get the solution of the ODEs and storing it in ‘theta’.
6. Then the Angular Displacement and Angular Velocity were plotted against Time using the solution of the ODE.
7. A counter variable is initiated and set to 1
8. Then the displacement values are taken from the solution array ‘theta’ and using a for loop, the individual values of displacement are plotted.
9. These plots are saved and each time the counter iss increased by 1.
10. All the saved plots are stitched together and made into an animation file using ImageMagick.

RESULTS

Plot-

Animation of motion of Pendulum-

The animation can be viewed from the following link-

https://drive.google.com/file/d/1oj-SKrBQ6hXVVX-2mcHTFVgpt_2L1WHQ/view?usp=sharing

DISCUSSION

The ‘numpy’ module is used for ‘linspace()’ command. Similarly the ‘math’ and ‘matplotlib’ modules are used for trigonometric commands and plotting graphs respectively. Also the ‘scipy.integrate’ module is used for ‘odeint’ command to solve the ODE(integrate the ODE).

Syntax for ‘odeint’, y = odeint(model,y0,t)

where,

model: Function name that returns derivative values at requested y and t values as dy/dt = model(y,t)

y0: Initial conditions of the differential states

t: Time points at which the solution should be reported. Additional internal points are often calculated to maintain accuracy of the solution but are not reported.

The plot is made by taking the position of pendulum at each value in time.The position is taken from the solution of the ODE i.e. the angular displacements.

The Angular Displacement and Angular Velocity plot is obtained. The animation shows the position of pendulum at 150 instances in between 0 to 20 seconds.

ERRORS ENCOUNTERED

Error due to use of ‘linspace()‘ instead of ‘np.linspace()’

CONCLUSION

So, the second order ODE was solved and the animation of the motion of the simple pendulum was successfully created.