Python Program: Solving second order ODEs

Aim: Write a Python code to simulate motion of a pendulum with damping.

Objective:

• Solving second order ODE which represents motion of pendulum.
• Simulating the motion of pendulum.
• Creating an animation of its motion.

Theory:

Simple Pendulum:

It is said that Galileo measured the periods of a swinging chandelier in church by his pulse beats. He observed that the motion of the chandelier was periodic. The system is a kind of pendulum. This periodic motion is simple harmonic for small displacements from mean position.

Simple Harmonic Motion:

Simple harmonic motion is a special type of periodic motion where displacement is sinusoidal function of time. In SHM, the restoring force on the moving object is directly proportional to and opposite of the object's displacement vector. It results in an oscillation which in absence of drag or any other dissipation of energy, continues indefinitely.

Damped Simple Harmonic Motion:

We know that the motion of a simple pendulum swinging in air, dies out eventually. This happens because the air drag and the friction at the support oppose the motion of the pendulum and dissipate its energy gradually. The pendulum is said to execute damped oscillations. In damped oscillations, the energy of the system is dissipated continuously, but for small damping the oscillations remain approximately periodic.

Ordinary Differential Equations:

An ordinary differential equation (ODE) are derivatives which contain one or more dependant variables and an independant variable. The term ordinary is used in contrast with partial differentials (PDEs can have more than one independant variables). An ODE can be classified as linear, non-linear, homogenous, non-homogenous. Based on the order of equation, we have, first order, second order, 'nth' order (where n is a natural number) equation. In engineering, ODE is used to describe the transient behaviour of a system and our interest for this project lies in pendulum.

Dynamics of Pendulum:

For predicting the motion of pendulum, we need to understand the dynamics and the equation that governs its motion. In general, motion of a pendulum can be considered simple harmonic and hence is given by:

`(d^2theta)/(dt^2) + (g/l)sintheta = 0`

where,

• `theta` = angular displacement of pendulum bob
• `t` = time (in sec)
• `g` = acceleration due to gravity
• `l` = length of string

But in actual practice motion of pendulum is damped due to dissipative effects. Hence, we need to introduce correction terms:

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

where,

• `b` = damping coefficient
• `m` = mass of pendulum bob (in kg)

This is differential equation of second order and we will use this to simulate the pendulum. But we need to simplify it. We'll use matrix method to do so.

Solution of ODE by matrix method:

Let `theta` be the angular displacement, then angular velocity will be.

`theta` = `theta1`

`(d theta)/(dt) = theta2`

`(d^2theta)/(dt^2) = d/(dt)((d theta)/(dt)) = (d theta2)/(dt)`

Substituting the above in damped ODE for pendulum we get;

`(d theta2)/(dt) + b/m theta2 + g/l sintheta1 = 0`

Hence, we arrive at these first order ODEs

`:. (d theta1)/(dt) = theta2`

`:. (d theta2)/(dt) = -(b/m theta2 + g/l sintheta1)`

Which can be written in matrix from as,

`d/(dt) [[theta1],[theta2]] = [[theta2],[-(b/m theta2 + g/l sintheta1)]]`

where,

• `theta1` = Angular displacement of pendulum bob
• `theta2` = Angular velocity of pendulum bob

Program:

The above equations can be expressed as shown below to obtain a python code.

#Second Order ODE

#Imporing required python libraries
from numpy import *
from scipy.integrate import odeint
from matplotlib.pyplot import *

#Function to form ODE
def model(theta,t,b,g,l,m):
    theta1 = theta[0]
    theta2 = theta[1]
    dtheta1_dt = theta2
    dtheta2_dt = -(b/m)*theta2 - (g/l)*sin(theta1)
    dtheta_dt = (dtheta1_dt, dtheta2_dt)
    return dtheta_dt

#Inputs
print("Provide Inputs in SI Units\n")
b = float(input("Enter Damping Coefficient : "))
g = 9.81 
l = float(input("Enter Length of Pendulum : "))
m = float(input("Mass of Bob : "))
print("\n")

#Initial Condition
theta_0 = [0,(pi/4)]

#Time input
t = linspace(0,20,500)

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

#Plotting Angular Displacement and Angular Velocity w.r.t. Time
plot(t,theta[:,0],label='Angular Displacement')
plot(t,theta[:,1],label='Angular Velocity')
grid()
legend()
title('Graphs for Pendulum - Premkumar Iyer')
savefig("Graph.png")
show()

x_1 = []
y_1 = []

#Simulation of Pendulum Motion
ct = 1 #loop counter
for i in range(len(t)):
  x_0 = 0
  y_0 = 0

  a = theta[i,1]

  x_1.append(l*sin(a))
  y_1.append(-l*cos(a))

  #Filename
  #The output jpg files were explicitly renamed so that it complies with cygwin
  if(ct<10):
        filename = 'pendulum00' + str(ct) + '.png'
        
  if(ct>10):
        if(ct<100):
            filename = 'pendulum0' + str(ct) + '.png'
        
        if(ct>=100):
            filename = 'pendulum' + str(ct) + '.png'

  #Plotting the motion
  plot([-1.5,1.5],[0,0],"k--",lw=3)
  plot([x_0,x_1[i]],[y_0,y_1[i]],lw=5)
  plot(x_1[i],y_1[i],"ro",ms=20)
  xlim(-1.5,1.5)
  ylim(-1.5,0.5)
  grid()
  title('Motion of a Simple Pendulum With Damping - Premkumar Iyer')
  savefig(filename)
  show()

  ct = ct+1

Working of Code:

1. The code begins with importing the required libraries. We'll use numpy (for calculations), scipy.integrate (for solving DE), matplotlib.pyplot (for plotting graphs)

2. A function is defined to from the required ODE in the form of code. Then, inputs are given (rather prompted to the user as input command is used). The initial conditions are set up. linspace command produces a linearly spaced array and we use it to produce time intervals.

3. The solution of ode is obtained using the function odeint which was present in scipy.integrate library. The result is stored in theta. This gives us solution for each time interval.

4. Plotting of the graph of angular displacement and angular velocity w.r.t time is done by calling in values from theta.

5. Empty arrays of x_1 and y_1 are created which would store the position of pendulum bob. Simulation of pendulum motion is facilitated by for loop which is iterated as per the time span.

6. In the loop, initial location is defined, again theta is called in to determine position of bob with repect to x and y axis. The values are appended in x_1 and y_1 arrays. Nested if statement is used to determine filename of the series of graphs which would be obtained in later blocks of code. A loop counter was initialized earlier for this purpose.

7. plot command is extensively used for this, first we plot the support onto which the string is attached. The we plot the string by joining a line between the origin and the bob. The bob is also plotted by calling in values from the array.

8. To obtain a consistent plot, the plot window is limited by using xlim and ylim commands, te grid is also kept on. A suitable title is provided to the plots and are saved using savefig command. show command displays the figures in output screen.

Output:

Output cell:

Graph:

The animation:

Result & Conclusion:

• The angular displacement and angular veocity are sinusoidal for a pendulum. As the motion is damped, the amplitude of motion reduces gradually.
• The time taken for pendulum to stop is inversely proportional to the damping coefficinet.
• Higher the damping coefficient, earlier the pendulum seizes its motion and vice versa.

References:

• Python for Mechanical Engineers - Course provided by Skill Lync
• Google Colab - Google LLC
• JupyterLab - Project Jupyter
• Anaconda - Anaconda Inc.
• Supporing Material: Physics-II, Textbook for Class XI - NCERT
• Supporting Material: Higher Engineering Mathematics - BS Grewal
• Supporting Material: WikiMedia

Material: Python Code, Notebook and Other Resources can be accessed here