Solving second order ODEs to Simulate the Motion a Simple Pendulum using Python

OBJECTIVE: Write a program to solve the ODE provided below of the equation of motion for a simple pendulum and to simulate its transient behavior. Also, create an animation of its motion.

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

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.

ODE's are used to describe the transient behavior of a system such as motion of a simple pendulum.

The motion of a simple pendulum is governed by Newton's second law of motion.

We have a second-order ODE that gives us the position of the ball with respect to time.

1. PROGRAM FORMAT:

Let's say,

`theta = theta_1`

`(d theta)/(dt) = theta_2` ---(1)

and,

`(d^2 theta)/(dt^2) = (d^2 theta_1)/(dt^2) = d /dt((d theta_1)/(dt)) = (d (theta_2))/(dt)`

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

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

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

From eqn. (1) & (2),

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

2. INPUT PARAMETERS:

Gravity(g) = `9.8ms^-2`

Mass of the pendulum's ball(m) = 1 Kg

Length of the pendulum's string(L) = 1 m

Damping Coefficient(b) = 0.05

Time Interval = 0-20 secs

Initial Condition at time(t) = 0:

Angular Displacement = 0

Angular Velocity = 3m/s

3. PROGRAM'S APPROACH FLOW:

4. PYTHON CODE:

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


# Function Definition
def derivative_func(theta, t, b, g, l, m):
    theta1 = theta[0]
    theta2 = theta[1]

    dtheta1_dt = theta2  # Equation_1
    dtheta2_dt = -(b / m) * theta2 - (g / l) * math.sin(theta1)  # Equation_2

    dtheta_dt = [dtheta1_dt, dtheta2_dt]
    return dtheta_dt


# Input Parameters
b = 0.05
g = 9.81
l = 1
m = 1

# Initial condition
theta_0 = [0, 3]  # Displacement & Angular Velocity = 0-3 radians

# Timing points
t = np.linspace(0, 20, 200)  # 200 Time-Intervals between 0 to 20 seconds

# Solving the Ordinary Differential Equations(ODE)
theta = odeint(derivative_func, theta_0, t, args=(b, g, l, m))

# Plotting Commands
plt.plot(t, theta[:, 0], 'm-', label=r'Angular Displecement (rad)')
plt.plot(t, theta[:, 1], 'c--', label=r'Angular Velocity (rad/s)')
plt.title('Angular Displacement & Angular Velocity vs Time')
plt.xlabel('Time--->')
plt.ylabel('Displacement--->')
plt.legend()
plt.savefig('Output.png')
plt.show()

# Code for creating the animation
ct=1
x0 = 0
y0 = 0

for i in theta[:,0]:
	x1 = (1*math.sin(i))
	y1 = (-1*math.cos(i))

	filename = 'test%05d.png' % ct
	ct = ct +1

	plt.figure()
	plt.plot([-1,1],[0,0],"black",linewidth = 10)
	plt.plot([x0,x1],[y0,y1],"red",linewidth =5)
	plt.plot(x1,y1,'o',markersize = 20)
	plt.xlim([-1.5,1.5])
	plt.ylim([-1.5,1.5])
	plt.savefig(filename)

5. EXPLANATION FOR THE PYTHON CODE:

1. First, all the necessary modules including Math module, Matplotlib module, Numpy module and SciPy required for this program is imported.

2. A function definition is required to state the discretized equations and values of theta together which returns the value for 'dtheta/dt' to be used later in the program.

3. Input Parameter are stated.

3. Initial Conditions for Angular Velocity is provided for the range of 0-3 radians.

4. Overall time intended for the simulation is subdivided into 200 time segments using 'linspace' command.

5. 'odeint' function is used to create the intended solver. Here, 'args' gives us all the input arguments, 'derivative_func' stands for derivative function defined initially in the program, 'theta_0' stands for initial conditions which are displacement(`theta`) and Velocity, 't' stands for time intervals between 0 and 20secs,

6. We're then inserting the plotting commands to plot the graph between Angular Displacement & Angular Velocity VS Time using matplotlib module as 'plt'.

7. Next, we are building code to create an animation for the movement of the pendulum motion. Here, we're setting the limits for X and Y-axis from -1.5 to +1.5. We're using a for loop to go through each values stored in 'theta'. Using x1 and y1, we are setting the axis for the movement of the pendulum. Also, we're using 'ct' as a loop counter variable which will give the command for the loop to move to the next value of theta for plotting and using 'savefig' comamnd we're saving all those respective plots at our designated location.

OUTPUTS DERIVED FROM THE PYTHON CODE:

A) Plot for Angular Displacement & Angular Velocity VS Time:

B) Animation to simulate the motion of the Pendulum:

Link for the Animation on YouTube: https://www.youtube.com/watch?v=HejPHlfvXoo