Python Program To Solve The Equation Of Motion

PYTHON PROGRAM TO SOLVE THE EQUATION OF MOTION OF A SIMPLE PENDULUM WITH DAMPING

Objective:

To write a Python program that would solve the equation of motion of a simple pendulum with damping and simulate the pendulum motion.

Simple Pendulum:

Any weight which is suspended from a pivot and can swing freely when displaced from resting position can be described as a pendulum.

A simple pendulum is an idealized mechanical model comprising of a point mass (also called as bob) suspended by a massless and inextensible string.

When a simple pendulum is pulled to one side and released, it oscillates about its equilibrium position in simple harmonic motion.

After the pendulum bob is released from certain angle θ, the path traced by it is an arc of a circle courtesy of restoring force. Restoring force always acts towards the equilibrium position, which in our case is the lowest point of the trajectory.

The way the pendulum moves depends on the Newtons second law. When this law is written down, we get a second order Ordinary Differential Equation that describes the position of the "ball" w.r.t time.

In the above equation,

g = gravity in m/s2,

L = length of the pendulum in m,

m = mass of the ball in kg,

b=damping coefficient.

Procedure:

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

d²θ/dt² + (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.

1. Libraries used:

The following libraries are imported to execute the program.

• matplotlib

Matplotlib is one of the most popular Python packages used for data visualization.

It is a cross-platform library for making 2D plots from data in arrays.

Here in the program,it is used to plot the motion of pendulum with respect to time.

• math

In python a number of mathematical operations can be performed with ease by importing a module named “math” which defines various functions which makes our tasks easier.

Here,math module is used for performing certain trigonometric operations.

• numpy

NumPy is a general-purpose array-processing package. It provides a high-performance multidimensional array object, and tools for working with the arrays.

Here,numpy is used with linspace that returns evenly spaced values within a given interval of time of oscillation of pendulum.(i.e. number of elements are returned).

• scipy

SciPy is a python library that is useful in solving many mathematical equations and algorithms.

The scipy.integrate sub-package provides several integration techniques including solving an ordinary differential equation integrator.

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.

2. Variable initialization and declaration:

The following variables are initialized as inputs,

L=1 metre

m=1 kg

b=0.05

g=9.81 m/s2

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.

4. Time points:

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

5. Defining a new function:

A new function is created with the name simple pendulum.

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 ODE:

For 'odeint' command to work properly we pass the differential simple pendulum function, initial conditions, time span and the arguments.

The solved answers are stored in an array of 2 columns 'theta'.

Output:

A plot is created between angular displacement versus time and angular velocity versus time.

In this plot, blue curve indicates the angular displacement and red dotted curve indicates the angular velocity of the pendulum motion.

As the damping effects are present in the system with a damping coefficient of 3,the amplitude of motion of the pendulum tends to decrease with respect to time.

8.Animation:

Two empty arrays are created to add the values of angular displacements and angular velocities solved in the problem using 'append' command.

The positions of pendulum at fixed end and free end are defined using coordinates.

The position of bob is also defined at the free end of string.

The plots are defined inside a 'for' loop to create an image then and there while the 'for' loop runs for each time with a saved name in sequence.

The created images are then stitched together and converted to animation.

Program:

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

# inputs
m = 1            # (mass of ball in kg)
b = 0.05         # (damping coefficient)
g = 9.81         # (acceleration due to gravity in m/s^2)
L = 1            # (length of string in metres)


def simple_pendulum(theta, t):

    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


# initial conditions
theta_0 = [0, 3]            # (first element represents the angular displacement in radians and second element represents the angular velocity in rad/sec)

# time points
t = np.linspace(0, 20, 300)

# solving ode
theta_results = odeint(simple_pendulum, theta_0, t)
Angular_Displacement = []
Angular_Velocity = [] 

# Printing of value angular displacement and angular velocity from the result array using the for loop
for i in range(0, len(theta_results)):

    # Angular Displacement
    angular_displacement = theta_results[i, 0]
    Angular_Displacement.append(angular_displacement)

    # Angular Velocity
    angular_velocity = theta_results[i, 1]
    Angular_Velocity.append(angular_velocity)

#creating animation

ct = 0              #(counter variable)      

for j in range(0, len(t)):

    THETA = Angular_Displacement[j]

    OMEGA = Angular_Velocity[j]

    # Defining the coordinates of fixed end of string of the pendulum
    X0 = 0
    Y0 = 0

    # Defining the coordinates of free end of string of the pendulum
    X1 = L*(np.sin(THETA))
    Y1 = -L*(np.cos(THETA))

    #counter variable to store plotted images in different names.
        
    ct = ct + 1
    filename = str(ct) + '.png'

    #animation  plots
    plt.plot([-1, 1], [0, 0], color = 'black', linewidth = 10)
    plt.plot([X0, X1], [Y0, Y1], color = 'red', linewidth = 4)
    plt.plot(X1, Y1, marker = 'o', markerfacecolor = 'blue', markersize = 30)
    plt.xlim(-2, 2)
    plt.ylim(-1.5, 0.2)
    plt.close()


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

Youtube video link:

https://youtu.be/0ia3FXPvS60