Solving a Second Order ODE for the Damped Oscillations of a Simple Pendulum

In mechanics and physics, simple harmonic motion is a special type of periodic motion of oscillation where the restoring force is directly proportional to the displacement and acts in the direction opposite to that of displacement and no other forces are involved.

In real oscillations, friction, or damping, slows the motion of the system. Due to frictional force, the velocity decreases in proportion to the acting frictional force. These oscillations are technically termed as Damped Oscillations. From Newton’s Second Law, on simplification, a second order differential equation is obtained as a function of displacement and angular velocity.

A simple pendulum in real conditions when made to oscillate, tends to dampen out after some time. The motion of the simple pendulum can be studied by a second order differential equation.

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

Here,

`theta` = angular displacement in radians

'b' = damping factor

`m` = mass of the pendulum

`g` = acceleration due to gravity

`l` = length of the pendulum

The above shown 2nd order differential equation gives the position of the pendulum bob at any instant of time with great accuracy. MATLAB can be used to calculate the oscillation of the pendulum, but it does not have a function to solve a 2nd order differential equation directly. So a 2nd order D.E. can be manipulated to form two 1st order D.E.'s and be solved. Matlab uses the function 'ode45' to solve a 1st order D.E.

Conversion of 2^nd order DE into two 1^st order DE.

Here, `theta` is the angular displacement. So, let `theta` be `theta_1`.

The differential of angular displacement is angular velocity. So,

`(d theta_1)/dt = theta_2`

This is the first differential equation.

The second one is obtained by differentiating `theta_2` with respect to time and we get,

`(d theta_2)/dt = -(b/m)*theta_2 - (g/l)*sintheta_1`

This is the second differential equation.

Both the 1st order differential equation can be represented as a solution matrix.

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

In the matrix form `[[theta_1],[theta_2]]` is the solution matrix, which can be addressed by a single letter, `theta`.

Therefore, `[[theta_1],[theta_2]]` = `theta`, in which,

`theta_1` = the angular position of the pendulum at various instances in time.

`theta_2` = the corresponding angular velocity at the same instances in time.

The function 'odeint' solves the input and returns 2 sets of values. The first column gives all the angular displacement `theta_1` in radians and the second column has angular velocity `theta_2`.

PART 1: Defining the function 'model(theta,t,b,g,l,m)'

• The libraries ‘numpy’, ‘scipy.integrate’, ‘math’, ‘matplotlib.pyplot’ are imported at the beginning in a fresh script

• A function is defined with return value being 'dtheta_dt' as shown

def model(theta,t,b,g,l,m)

• The name of the function is 'model', and the input values are theta,t,b,g,l,m.

• In the body of the function, `theta_2` and `(d theta_2)/dt` are calculated and are assigned as elements of the matrix 'dtheta_dt'.

• This matrix is returned to the main program.

Part 2: Writing the main program

• After the function body, initialize the parameters of the coefficients like damping factor, acceleration due to gravity, length of the pendulum and mass of the bob. And dependent variables like time interval, initial values of position and initial velocity.

• Call the function 'odeint' and assign the output to a variable. In calling the function the main dependent variables are specified. And the function 'model(theta,t,b,g,l,m)', time interval and the initial position and velocity are given as inputs. The function 'model(theta,t,b,g,l,m)' returns a matrix of `theta_2` and `(d theta_2)/dt`.

• The function 'odeint' accepts the values from the inputs and gives two outputs in the form of a matrix calculated by integration. A matrix of angular displacement `theta_1` and angular velocity `theta_2`.

• From the obtained values, a plot of displacement and velocity vs time is made which shows the rate of change of the quantities. The curve can be varied by changing the damping coefficient.

• Now, a counter is assigned to 1. A 'for' loop is opened to plot every position of the pendulum for the specified time interval. A pointer 'i' is assigned to the displacement column of the output matrix. So, whenever the 'for' loop increments the value of the pointer, the successive displacement value of the pendulum is obtained.

• Using trigonometry, the exact position of the pendulum bob is calculated and stored in temporary variables as coordinates.

• A plot of [X,Y] is done and a marker is assigned as per requirements to represent a bob of the pendulum. This plot is stored in an array of the counter, and then the counter is incremented and the loop is closed.

The Python program

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

# Function that returns dz/dt
def model(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

b = 0.05 # Damping factor
g = 9.81 # Acceleration due to gravity
l = 1 # Length of the pendulum
m = 1 # Mass of the pendulum

theta_0 = [math.pi/2,0] # Initial position in radians and velocity

t = np.linspace(0,10,1000)

theta = odeint(model,theta_0,t,args = (b,g,l,m))

plt.figure(1)   
plt.plot(t,theta[:,0],'b-',label=r'$\frac{d\theta_1}{dt}=\theta2$')
plt.plot(t,theta[:,1],'r--',label=r'$\frac{d\theta2}{dt}=-\frac{b}{m}\theta_2-\frac{g}{L}sin\theta_1$')
plt.ylabel('Plot')
plt.xlabel('Time')
plt.legend(loc='best')
plt.show()

x = []
y = []

for i in range(0,len(theta)):
    x0 = 0
    y0 = 0

    x.append(math.cos(theta[i,0]-math.pi/2))
    y.append(math.sin(theta[i,0]-math.pi/2))

    filename = '%03d.png' % i
    
    plt.figure(2)
    plt.clf()
    plt.plot([x0,x[i]],[y0,y[i]],'b-')
    plt.plot(x[i],y[i],'ro',markersize=20)
    plt.xlim([-1.5,1.5])
    plt.ylim([-1.5,0.5])
    plt.savefig(filename)
    

Code used in CYGWIN to create a '.mp4' file through 'imagemagick'

$ magick *.png Simple_pendulum.mp4

OUTPUT Parameters v/s time graph

The Damped Oscillations of a simple pendulum is shown in the link below

https://youtu.be/01WDlwDS3O4

ERRORS

Error: The syntax of the for loop is wrong.

Correction: There should be a colon, ':' at the end of 'for' statement.

Error: Syntax is wrong at assigning a 3 digit number for a filename.

Correction: '%03d.png' % i

This is the correct syntax.