Solving second order ODE of movement of Simple pendulum in Python

Objective

To plot the angular velocity and angular displacement of a pendulum motion vs time. And to create an animation of motion of the pendulum.

Introduction

In Sublime using python, the motion of simple pendulum was animated by ploting various images which represent the position of pendulum at different times. And also for given initial angular displacement and angular velocity, the ordinary differential equation for motion of simple pendulum was solved using odeint solver.

Simple pendulum

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

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

The above equation is the second order ordinary differential equation which represents the motion of a simple pendulum.

where,

b = Damping coefficient

m = Mass of the ball in kg

L = Length of the string in m

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

But ODE solver does not solve a second order equation so this equation was must breakdown to first order equation and the equation is

`theta = theta_1`

`(d theta_1)/(dt) = theta_2` ------------------------------------------ Equation 1

`(d theta_2)/(dt) = -(b/m) theta_2 - (g/L) sin theta_1` ------------- Equation 2

Equation 1 represents the angular displacement of the pendulum and equation 2 represents the angular velocity of the pendulum.

Program

# Simplt pendulum 

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

# Function for solving differential equations
def pend(theta, t, b, g, l, m):

	theta_1 = theta[0]
	theta_2 = theta[1]
	dtheta1_dt = theta_2
	dtheta2_dt = -(b/m)*theta_2-(g/l)*math.sin(theta_1)
	dtheta_dt = [dtheta1_dt, dtheta2_dt]
	return(dtheta_dt)

# Input variables

# Damping coefficient
b = 0.05

# Mass of the pendulum ball in kg
m = 1

# Acceleration due to gravity in m/s^2
g = 9.81

# Length of the string in m
l = 1

# Initial conditions (Displacement, velocity and time)
theta = [math.radians(0),3]
t = np.linspace(0,20,500)

# Calling ODE function
result = odeint(pend, theta, t, args = (b,g,l,m))

# Ploting Angular Displacement and Angular Velocity vs Time
plt.plot(t, result[:,0], \'r--\', label = r\'$\\frac{d\\theta_1}{dt}=\\theta_2$\'+ \'  Angular Displacement\')
plt.plot(t, result[:,1], \'b-\', label = r\'$\\frac{d\\theta_2}{dt}=\\frac{-b}{m}*\\theta_2-\\frac{g}{l}*sin\\theta_1$\' + \'  Angular Velocity\')
plt.xlabel(\'Time in sec\')
plt.legend()
plt.show()

# Program for creating pendulum movement animation
theta = np.linspace(-40,40,50)

# Fixed end coordinates and length of the string
x0 = 6
y0 = 7
l = 5

# End position of the string
x1 = []
y1 = []

# Variable for accessing angles values in reverse order
rev = 2*len(theta)

# Counter variable
forw = 1

# Loop for creating forward motion of the pendulumn
for i in theta:

	# Formulae for calculating the cordination of end points at various time
	x1.append(x0 + l*math.sin(math.radians(i)))
	y1.append(y0 - l*math.cos(math.radians(i)))
	

for j in range(0, len(theta)):
	
	# Ploting the position of pendulumn and saving the plot
	plt.figure(forw)
	plt.plot([5.5,6.5], [7,7], \'--\', [6, 6], [7, 2], \'k-.\', x1, y1, \'b--\')
	plt.plot([x0, x1[j]], [y0, y1[j]], x1[j], y1[j], \'or\', markersize = 15)
	plt.arrow(x1[j], y1[j], 0.75, 0, length_includes_head=True, head_width=0.075)
	plt.xlim(0, 13)
	plt.ylim(0, 10)
	plt.savefig(\'ODE img%04d\' %forw)
	forw = forw+1

	# Ploting the position of pendulumn and saving the plot			
	plt.figure(rev)
	plt.plot([5.5,6.5], [7,7], \'--\', [6, 6], [7, 2], \'k-.\', x1, y1, \'b--\')
	plt.plot([x0, x1[j]], [y0, y1[j]], x1[j], y1[j], \'or\', markersize = 15)
	plt.arrow(x1[j], y1[j], -0.75, 0, length_includes_head=True, head_width=0.075)
	plt.xlim(0, 13)
	plt.ylim(0, 10)
	plt.savefig(\'ODE img%04d\' %rev)
	rev = rev-1

Program explanation

First thing in Sublime was to save the file with .py extension along with a suitable name then only sublime compile that file with python otherwise sublime wouldn’t compile the program.

Next thing after saving the file in python was to import the necessary modules(In python a module is a file containing python code library) which were helpful for this program and the necessary modules were math, matplotlib.pyplot, numpy and odeint from scipy.integrate.

The mass of the ball, length of the string, damping coefficient, acceleration due to gravity and the initial conditions such as angular displacement, angular velocity and time and the values were

Mass of the ball = 1 kg

Length of the string = 1 m

Damping coefficient = 0.05

Acceleration due to gravity = 9.8 m/s^2

Initial angular displacement = 0 degree

Initial angular velocity = 3 rad/s

Time span = 20sec

After defining the initial variables, odeint solver called to solve ode equation of simple pendulum.

Syntax:

odeint(function, initial_condition_1, initial_condition_2, args = arguments)

Now, the solver returns the angular displacement and angular velocity of the pendulum at various time of given time span.

The output values returned by the odeint solver was saved in a variable defined as result as an array of two columns. First column contains the angular displacement values and the second column contains the angular velocity values.

Now the angular displacement and angular velocity were ploted with time in x-axis. Angular displacement was represented in red colour dashed line by using ‘r--’ condition and angular velocity was represented in blue colour continuous line by using ‘b-‘ condition. X-axis was labelled as ‘Time in sec’ using xlabel command. Since to two things were ploted legend was used to differentiated which is which.

In python legend() command used to show the legend but the name of the curve or line must defined in plot() command by using label condition.

Now to create an animation of motion of pendulum a new variable called theta was defined with contains the values of theta i.e., angular displacement (not called by odeint solver).

For this value of theta, the position of pendulum was calculated and ploted and the image saved with suitable and different name for each position.

The fixed end of the string fixed at coordinates of (6,7) and the length of the string was defined as 5units.

To calculate the fixed and of the pendulum basic trigonometry formula was used through for loop and the values were appended in arrays x1 and y1. And for each angle (theta) the position of pendulum was calculated, ploted and saved with name of ‘ODE img001.png’.

The limit of x-axis and y-axis was fixed to avoid changing of limit of axes using xlim and ylim command.

A pendulum oscillates back and forth, so an arrow included in to plot to represent the direction. Since an arrow included for same position of the ball of pendulum two images will be saved, one with right side arrow and another with left side arrow.

Two counter variables forw and rev were defined. ‘forw’ variable holds the value of 1 and increased by one unit i.e., 1,2,3,4,5…… for each iteration and rev holds the value 2xlen(theta) i.e., two times of length of theta array and the value of rev variable decreased by one unit i.e., 200,199,198,197,196… for each iteration.

‘forw’ variable represents the plots which had the right side arrow and ‘rev’ variable represents the plots which had left arrow.

For clear understanding lets look at the program, theta array holds 100 values therefore, rev variable holds 200 and as defined forw holds 1.

Assume all the end coordinates of string were calculated with respect to the theta values. Now for the first value of theta that is for degree -45 the position of pendulum ploted with a left arrow and right arrow. Now the plot which has the right side arrow was saved with name ‘ODE img0001.png’ and the plot which has the left side arrow was saved with name ‘ODE img0200.png’.

Now from the name you can saw that forw and rev used for saving image name.

The reason for creating two plots with two different sided arrow for the same position was, while oscillating first position of forward motion is the last position of reverse motion.

Therefore ‘ODE img0001.png’ was first image and ‘ODE img0200.png’ was last image of animation.

Like that ‘ODE img0100.png’ was last image of forward motion i.e., last image with right side arrow and ‘ODE img0101.png’ was the first image of reverse motion i.e, last image with left side arrow.

Like this, in for loop for each iteration two images will be saved, one is for forward motion and another is for reverse motion.

Result

The plots between time and angular velocity and angular displacement were done and an animation of motion of simple pendulum was also done.