Solving 2nd Order ODE For A Simple Pendulum Using Python

Objective: To write a program that solves the following ODE. This ODE represents the equation of motion of a simple pendulum with damping

The motion of a simple pendulum is a basic classical example of simple harmonic motion, consisting of a small bob and a massless string.Movement of the pendulum depends on the Newtons second law. When this law is written down,a second order Ordinary Differential Equation is formed that describes the position of the endof the pendulum w.r.t time.

Reference: http://file.scirp.org/Html/7-1720775_73856.htm

The 2nd order ODE for a simple pendulum with damping is given below.

`(d^2 θ)/(dt^2 )+b/m (dθ)/dt+g/L sin⁡θ=0`.....(1)

The above ODE is of order 2.To solve it we convert it into two simultaneous ODE of order 1.

Let `θ=θ_1`

and `(dθ_1)/dt=θ_2`...........(2)

Therefore the equation 1 becomes,

`(dθ_2)/dt+b/m (dθ_1)/dt+g/L sinθ_1=0`

`⇒(dθ_2)/dt=-b/m (dθ_1)/dt-g/L sinθ_1`.........(3)

Here, `θ_1` is the angular displacement and `θ_2` is the angular velocity of the pendulum.

Procedure:

In python to solve two simultaneous ODEs following steps is taken.

1.Importing libraries.

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

numpy library to use linspace() command which returns evenly spaced samples.

odeint is imported from scipy.integrate library to solve differential equation.

math library for trignometric equation.

matplotlib library for plotting graphs.

2. Defining constant coefficients of equation 1.

b = Damping coefficient = 0.05

g = gravitational acceleration = 9.81 `ms^(-2)`

l = length of string = 1 m

m = mass of pendulum = 1 Kg

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

3. Initial conditions.

At t = 0 , angular displacement for the pendulum is 0 and angular velocity is 3 rad per second.These conditions are stored in an array.

#initial condition
theta_0 = [0,3]

4. Time span.

A time array of 20 sec with 240 intervals is created upon which the ODE will be solved. It is created using numpy module.

#time plot
t = np.linspace(0,20,240)

5. Defining a simultaneous differential function.

A function is created which represents the system of differential equation, i.e.equation 2 & 3.

#Defining the differential equations

def system(theta,t,b,g,l,m):
	theta1 = theta[0]
	theta2 = theta[1]
	dtheta1_dt = theta2
	dtheta2_dt = -(b/m)*theta2-g*math.sin(theta1)
	dtheta_dt=[dtheta1_dt,dtheta2_dt]

	return dtheta_dt

The inputs of this function are theta (empty array which stores solved answers), t (time array),b (damping coefficient),g (gravitational acceleration), l (length of string) and m (mass of pendulum).'theta1' is the angular displacement and 'theta2' is the angular velocity. This function returns an array which has 2 columns and numerous rows. The 1st column is angular displacement and the second column is angular velocity.

6. Solving the ODEs.

Using 'odeint' command the ODEs are solved. The 'odeint' is a function in the scipy module which integrates the differential equation upto a time step and after getting that solution it keeps on integrating it. The 'odeint' function is generic means it will choose the best solver based upon the provided system. For 'odeint' command to work properly we pass the differential system function, initial conditions, time span and the arguments.

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

we are keeping the solved answers in an array of 2 columns 'theta'.

7. Creating animation.

The 'theta' array contains the solved angular displacement and velocity. Using the 1st column the position of the pendulum is plotted and it is joined to a fixed point(10,10) to look like a swinging pendulum.The position of the swinging end of the pendulum is calculated using geometry.

#For animation creating plots
f=1
for i in range(0,240):
	filename = str(f)+\'.png\'
	f= f+1
	plt.figure()
	plt.plot([10,l*math.sin(theta[i,0])+10],[10,10-l*math.cos(theta[i,0])],marker=\"o\")
	plt.xlim([0,20])
	plt.ylim([0,20])

	
	plt.savefig(filename)

Each figure is saved using 'plt.savefig' command separately with name coresponding to the 'i' value in '.png' format.For creating animation of the pendulum saved plots are stiched together to form animation using external software.

Ploting time vs angular displacement graph and time vs angular velocity graph using the 1st and 2nd column of theta array.

plt.plot(t,theta[:,0],\'b-\')
plt.plot(t,theta[:,1],\'r--\')
plt.show()

Conclusion:

The 2nd order ODE for a simple pendulumn with damping is solved using python and the calculated angular displacement and velocity is plotted against time. Using the saved plots animation is created by using external software.

The YouTube link to the pendulum animation is given below:

https://www.youtube.com/watch?v=7_OK6ocEyQM

Full code:

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

#Defining the differential equations

def system(theta,t,b,g,l,m):
	theta1 = theta[0]
	theta2 = theta[1]
	dtheta1_dt = theta2
	dtheta2_dt = -(b/m)*theta2-g*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 plot
t = np.linspace(0,20,240)

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


#For animation creating plots
f=1
for i in range(0,240):
	filename = str(f)+\'.png\'
	f= f+1
	plt.figure()
	plt.plot([10,l*math.sin(theta[i,0])+10],[10,10-l*math.cos(theta[i,0])],marker=\"o\")
	plt.xlim([0,20])
	plt.ylim([0,20])

	
	plt.savefig(filename)
	
	
plt.plot(t,theta[:,0],\'b-\')
plt.plot(t,theta[:,1],\'r--\')
plt.show()