Solving second order ODEs

Aim:

Write a program to simulate the transient behaviour of a simple pendulum and to create an animation of its motion.

Objective:

• Stimulate the motion between 0-20 seconds for angular displacement and angular velocity =3 rad/s at time t=0.
• Solve second order ODE of a simple pendulum with damping using odeint function in scipy library in python.
• Create an animation of the simple pendulum

Theory:

Consider a simple pendulum having mass ‘m’ attached to a string of length ‘l’ and having a damping coefficient ‘b’ is acting under the acceleration due to gravity ‘g’.

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

in the given equation g=9.81m/s^2, l=1m, m=1 kg, b=0.05 N-s/m

Steps to solve the problem:

• This will be solved by using odeint function of scipy library.

Code:

# Solving 2nd order differential equation
# By: Anubrat Padhee

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

# Inputs
# Damping coefficient unit-Ns/m
b=0.05;
# Acceleration due to gravity unit-m/s^2
g=9.81;
# Length of the string unit=m
l=1;
# Mass of the ball unit=kg
m=1;

# Initial conditions
# First one represents angular displacement and second one represents angular velocity
theta_0=[0,3]

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

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

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

# Plotting results
plt.figure()
plt.plot(t,theta[:,0])
plt.plot(t,theta[:,1])
plt.xlabel("Time")
plt.ylabel("angular displacement & angular velocity")
plt.legend(["angular displacement","angular velocity"])
plt.show()

# Creating animation of pendulum
ang_dis=theta[:,0]
# Counter variable
ct=1
for ang in ang_dis:
	# Starting coordinates
	x0=0
	y0=0
	# Ending coordinates
	x1=l*math.sin(ang)
	y1=-l*math.cos(ang)

	# Plotting
	plt.figure()
	plt.plot([-1.5,1.5],[0,0],linewidth=8)
	plt.plot([x0,x1],[y0,y1],linewidth=3)
	plt.plot(x1,y1,'o',markersize=15)
	plt.xlim([-1,1])
	plt.ylim([-1,1])
	filename=str(ct)+".png"
	plt.savefig(filename)
	ct=ct+1

Outputs:

Conclusion:

So, using the ordinary differential equation of simple pendulum and assuming angular velocity and displacement, an animation of simple pendulum was created.

Animation file:

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