Solving second order ODE for Simple Pendulum under Damping Condition

Objective:

To write a program that solves the ODE that represents the equation of motion of a simple pendulum with damping. And also, to simulate the transient behavior of a simple pendulum and to create an animation of its motion.

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

Problem Statement:

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 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.

A simple pendulum consists of a mass m hanging from a string of length L and fixed at a pivot point. When displaced to an initial angle and released, the pendulum will swing back and forth with the periodic motion. By applying Newton's second law for rotational systems, the equation of motion for the pendulum may be obtained as:

`(d^2theta)/dt^2 +g/L *sin theta=0`

The primary forces acting on the bob are the gravitational force that makes it move in the first place and the force exerted by the string to keep it moving along a circular path. Also, there may be a damping force from friction at the pivot or air resistance or both.

Assuming that the damping force is proportional to velocity. Since arc length and central angle are themselves proportional (with proportionality constant. The negative sign is because the damping force has to be opposite the direction of motion. When this term is included in the model, the equation becomes:

`(d^2theta)/dt^2 = - b/m *(d theta)/dt - g/L *sin theta`

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.

Given values for the parameters:

1. L=1 metre,
2. m=1 kg,
3. b=0.05.
4. g=9.81 m/s2.
5. Time period to be considered 0-20 sec,
6. angular displacement=0, angular velocity=3 rad/sec at time t=0.

Solution Procedure:

Steps involved in Coding:

1. Import module numpy as np, from scipy. integrate import odeint, matplotlib. pyplot as plt and math module.
2. Defining of the function that returns the value of theta differentiated w.r.t time
3. Defining the constant values of l, g, b, and m as per the given condition.
4. Declaring the initial conditions of theta from 0 to 3
5. Using linspace from numpy module for creating a numeric sequence. Syntax of np. linspace (start point, last point, number of values to be taken in between the start and endpoint)
6. Solving the hence formed ODE using odeint command.
7. Plotting result by using plt command, thereby forming the required graphs
8. Naming the filename and sequential stitching of the plots formed for the creation of animation

Code:

"""
Title: Solving second degree Ordinary Differential Eqaution for a Simple Pendulum under damping conditions
Submitted by M K Pragya
___________________________________________________________________________________________________________________________________________________________________
"""
#importing the required modules

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

#definition of the function that returns dtheta/dt

def mode1(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


# Defining the Constants
b=0.05   			# Damping Coefficient (higher the coefficient , faster the oscillations die out)
g=9.81				# Acceleration due to Gravity
l=1				    # Length of the String of the pendulum in metre(m)
m=1					# Mass of the ball of the Pendulum in kg

#initial conditions
theta_0=[0,3]

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

# Solving ODE for the simple pendulum
theta= odeint(mode1, theta_0, t, args=(b,g,l,m))

# Plotting results

fig=plt.figure()
fig.patch.set_facecolor('lemonchiffon')																									# setting the background for the plot screen
plt.plot(t, theta[:,0], 'b-', linewidth= 2, label=r'$frac{dtheta_1}{dt}=theta2$ (Displacement)')										# displacement
plt.plot(t, theta[:,1], 'r--', linewidth=2 ,label=r'$frac{dtheta2}{dt}=-frac{b}{m}theta_2-frac{g}{L}sintheta_!$ (Velocity)')		# velocity
plt.ylabel('Plot')																														# Y- Axis Title
plt.xlabel('Time')																														# X- Axis Title
plt.legend(loc='best')																													# Printing the Legend
plt.savefig('output.png')
plt.show()

#_____________________________________________________________________________________________________________________________________________________________________

# Animation

# Defining Constants
ct=1																			# Defining the counter variable for plotting the positions

# Initial Conditions
x0=0
y0=0

# For Plotting the positions of the Pendulum

for i in theta[:,0]:

	# Using trigonometry to find coordinates of subsequent points

 	x1=(l*math.sin(i))
 	y1=(-l*math.cos(i))

 	# For generating plots

 	filename='test%05d.png'% ct  												# Saves the required plot by ordered name
 	ct=ct+1
 	fig=plt.figure()
 	fig.patch.set_facecolor('khaki')											# Setting the background for the plot screen

 	# Creating the Pendulum

 	plt.plot([-1.0,1], [0,0],, linewidth=25)					# Surface to which Pendulum is attached
 	plt.plot([x0,x1],[y0,y1],color="darkblue", linewidth=8)						# String of the Pendulum
 	plt.plot(x0,y0,marker='8', markerfacecolor='black', markersize=15)			# Hinge point of the pendulum from the surface
 	plt.plot(x1,y1,marker='o', markerfacecolor='crimson', markersize=20)		# Ball of the Pendulum

 	# Formatting of the plot

 	plt.title('Motion of a Simple Pendulum under Damping Condition', fontname="Times New Roman",fontsize=15,fontweight='bold')
 	plt.xlabel('X-Axis',fontname="Times New Roman",fontsize=8)					 # For forming the static image
 	plt.ylabel('Y-Axis', fontname="Times New Roman",fontsize=8)

 	plt.xlim([-1.5,1.5])
 	plt.ylim([-1.5,1])

 	# For creating the Animation

 	plt.axis('off')																
 	plt.savefig(filename)

Errors Incurred:

Following are the errors incurred and the way the errors were overcome:

Results:

Plots are formed for Angular Displacement, Angular Velocity with respect to Time.

Image of the pendulum hence formed:

Link to the animation:

https://youtu.be/bJCrKGGAFqc

Discussion:

In a vacuum with zero air resistance, a pendulum will continue to oscillate indefinitely with a constant amplitude. However, the amplitude of a simple pendulum oscillating in the air continuously decreases as its mechanical energy is gradually lost due to air resistance.