Week 3 - Solving second order ODEs

Objective:

• To simulate the transient behavior of a simple pendulum using ODE
• To create an animation of its motion

Problem Statement:

1. To write a program in python that solves the ODE which corresponds to the motion of a simple pendulum with damping and to plot its angular displacement and angular velocity w.r.t. time.
2. To simulate the motion between 0-20 sec, for angular displacement = 0, angular velocity = 3 rad/sec at time = 0.

Theory:

A simple pendulum is defined to have a point mass, also known as the pendulum bob, which is suspended from a string of length L with negligible mass (Figure). Here, the only forces acting on the bob are the force of gravity (i.e., the weight of the bob) and tension from the string. The mass of the string is assumed to be negligible as compared to the mass of the bob.

credited: https://openstax.org/books/university-physics-volume-1/pages/15-4-pendulums

Consider the torque on the pendulum. The force providing the restoring torque is the component of the weight of the pendulum bob that acts along the arc length. The torque is the length of the string L times the component of the net force that is perpendicular to the radius of the arc. The minus sign indicates the torque acts in the opposite direction of the angular displacement:

`tau = -l(mgsintheta)`

`I*alpha = -l(mgsintheta)`

`I * (d^2theta)/(dt^2) = -l(mgsintheta)`

`ml^2(d^2theta)/(dt^2) = -l(mgsintheta)`

and by rearranging the above equation we have,

`(d^2theta)/(dt^2) = -g/lsintheta`

and if we consider the damper in the system the above equation changes slightly as

`(d^2theta)/dt^2 = -(b/m)(d theta)/dt - (g/l)sintheta`

Governing equations:

To code the ODE in python we do break the higher-order ODE and here it is second-order ODE into the first-order ODE as,

`=> theta = theta1`,

`=> (d theta)/dt = (d theta_1)/dt = theta_2`,

`=> (d^2theta)/dt = (d^2theta_1)/dt = (d theta_2)/dt`,

`=> (d^2theta)/dt + (b/m)(d theta)/dt + (g/l)sintheta = 0`,

`=> (d^2theta_1)/dt^2 + (b/m)(d^2theta_1)/dt^2 + (g/l)sintheta_1 = 0`,

`=> (d theta_2)/dt + (b/m)theta_2 + (g/l)sintheta_1 = 0`,

`therefore (d theta_1)/dt = theta_2`and `(d theta_2)/dt = -(b/m)theta_2 - (g/l)sintheta_1`,

and in the matrix form, we can write as,

`d/dx[[theta_1],[theta_2]] = [[theta_2],[-(b/m)theta_2 - (g/l)sintheta_1]]`,

Solution Procedure:

• Initially, I import the required modules to solve various functions required in this problem.
• Given the inputs i.e., b, l, g, m which is given in the question.
• Here now we provide an array of the initial conditions of theta for initial angular displacement and initial angular velocity.
• Again here we provide an array of time span to store 200 equally spaced intervals between time 0-20 seconds.
• Here I define the function named ode_function which include all the parameter that I mention above, which eventually to solve the second-order ODE.
• We call the ode_function by the name 'odeint' command and with all the proper syntax of the calling function and result obtained is shown below the code as to plot a graph which eventually shows the variation of the angular displacement and angular velocity w.r.t. time.
• For creating an animation of the harmonic motion of the pendulum, we used for loop and a ct variable is initiated, the value of angular displacement is stored in a separate variable called temp and iteration is performed inside the for loop.
• In this execution, one end of the pendulum is fixed and at another end, the bob is hung with the help of a string and move to the solution of ODE and is plotted for each time interval.

I do hereby mention all the required steps in my code commands and easily understandable.

I hereby attached my code below it.....

PYTHON CODE OF SECOND ORDER ODE AND FOR ANIMATION

# Importing the required modules
import numpy as np
import math
import matplotlib.pyplot as plt
from scipy.integrate import odeint

def ode_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 # damping coefficient
g = 9.81 # acceleration due to gravity
l = 1    # length of string
m = 1    # mass of bob

# Initial conditions
theta_0 = [0,3]
# Time points
t = np.linspace(0,20,200)
# Solve ODE 
theta = odeint(ode_function,theta_0,t,args = (b,g,l,m))

# Plotting the graph
plt.plot(t,theta[:,0],'b', label = r'$\frac{d\theta1}{dt} = \theta2$')
plt.plot(t,theta[:,1],'--r', label = r'$\frac{d\theta2}{dt}=-\frac{b}{m}\theta2-\frac{g}{l}sin\theta1$')
plt.ylabel('Plot',)
plt.xlabel('Time (s)')
plt.title('Damped Harmonic Motion of a Simple pendulum', color = 'maroon')
# To show a legend on the plot at lower right
plt.legend(loc ='lower right')
# To show a grid
plt.grid()
# Function to show the plot
plt.show()

#  Iteration using for loop to create animation of motion of pendulum

ct = 1
temp = theta[:,0]
for i in temp:
	x_start = 0
	y_start = 0
	x_end = l*math.sin(i)
	y_end =-l*math.cos(i)
	filename = 'test%0.5d.png' % ct
	ct = ct + 1

	plt.figure()
	plt.plot([x_start,x_end],[y_start,y_end], linewidth = 3)
	plt.plot(x_end,y_end, 'o', markersize = 20, color = 'red')
	plt.plot([-0.5,0.5],[0,0], linewidth = 10, color = 'blue')
	plt.plot([0,0],[0,-1.5], '--', color = 'green')
	# Setting the x-axis and y-axis range
	plt.xlim([-1.25,1.25])
	plt.ylim([-2,1])
	# Giving a title to my graph
	plt.title('Animation of transient behaviour of a simple Pendulum')
	# Function to save the fig
	plt.savefig(filename)
	

Errors:

I haven't encountered a major error while writing and solving the code but indentation error I do sometimes encounter.

Overcome Errors:

I do overcome those errors by carefully observing the name of error shown by the system in the console.

Results:

The above plot shown is the angular displacement and angular velocity w.r.t. time of the pendulum motion which is decreasing because of the damping effect which exponentially decreases the amplitude and velocity of motion.

For animation, I hereby generate 200 plots and were stitched together using ImageMagick Display to form animation. I hereby attached the link of the animation below...

https://youtu.be/-xx_UTWyqg8

Conclusion:

• Here we could understand that any higher-order ODE or here it is second-order ODE can be solved by converting them in the first-order ODE using python.
• Simulation of the transient behavior of simple pendulum and its effect when we consider damping on its motion is explained by this code.

Notations:

`theta =`angular displacement (position)

`t = time`(second)

`b =`damping coefficient `((N*s)/m)`

`g =`acceleration due to gravity `(m/(s^2))`

`l =`length of string `(m)`

`m =`mass of bob `(kg)`

`tau =`torque `(N*m)`

`I =`moment of inertia `((kg)/(m^2))`

`alpha =`angular acceleration `((rad)/(s^2))`