Solving second order ODEs (Python)

1.OBJECTIVE: - PROGRAM TO SIMULATE THE TRANSIENT BEHAVIOR OF A SIMPLE PENDULUM AND TO CREATE AN ANIMATION OF IT'S MOTION.

2.THEORY: -

In Engineering, ODE is used to describe the transient behavior of a system.

The way the pendulum moves depends on the Newton’s 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 these systems.

Our 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,

Gravity (g) = 9.81 m/s^2

Length of the pendulum (L) = 1 m

Mass of the ball (m) = 1 Kg

Damping coefficient (b) 0.05

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

3.PROGRAM:-

To solve 2^nd order ODE in Python we transform the 2^nd order equation into two 1^st order ODE

Using ode45() we solve these equation in Python.

3.1 INPUTS:-

• Gravity (g) = 9.81 m/s^2
• Length of the pendulum (L) = 1 m
• Mass of the ball (m) = 1 Kg
• Damping coefficient (b) 0.05
• Time Span = 0 to 20 Sec
• Initial condition (at t = 0 ) -

Angular Displacement = 0

Angular velocity = 3 m/s

3.2 PYTHON CODE:-

# PROGRAM TO SIMULATE THE TRANSIENT BEHAVIOR OF A SIMPLE PENDULUM AND TO CREATE AN ANIMATION OF IT'S MOTION
import math
import matplotlib.pyplot as plt 
import numpy as np 
from scipy.integrate import odeint

#Funtion (that returns dtheta/dt)
def function(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

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

#Initial condition
theta_0 = [0,3]

#Time points
t = np.linspace(0,20,200)

#SOLVE ODE
theta = odeint(function,theta_0,t,args = (b,g,l,m))

#PLOTTING
plt.plot(t,theta[:,0],'g-',label=r'Angular Displecement (rad)')
plt.plot(t,theta[:,1],'b--',label=r'Angular Velocity (rad/s)')
plt.title('Angular displacement & Angular velocity vs Time graph')
plt.xlabel('plot')
plt.ylabel('Time')
plt.legend()
plt.savefig('Output.png')
plt.show()


#ANIMATION
ct=1
x0 = 0
y0 = 0

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

	filename = 'test%05d.png' % ct
	ct = ct +1

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

3.3 EXPLANATION OF THE PYTHON CODE:-

1. Import module numpy as np, scipy.inedrate import odeint for second ODE, math and for plots matplotlib
2. Define a function named 'function' which returns the dtheta/dt value.
3. Input the given values (b,g,m,l).
4. Defining the initial condition from theta i.e. 0 to 3 rad/s.
5. Use np.linspace for creating numeric sequence of time points.
6. Solve ODE by using the odeint command.
7. Plot results by using plt.plot command.
8. For creating the animation run a for loop and save all images (created in every for loop) using savefig command.
9. Using those saved images create a animation in an another software.

3.4 OUTPUT:-

• Angular displacement & Angular velocity vs Time graph –

• Simulation video of the pendulum –

Google Drive Link -https://drive.google.com/drive/folders/19LmVFmp8drlEi7GRXVEFCaGzFu96Hztc?usp=sharing

4.CONCLUSION:-

From the graph and the video we can see that the velocity is reduced with time. This happens because of drag force due to the air resistance on the bob of the pendulum which is the cause of its damping,