PYTHON Project: Solving Second Order ODEs (Simple Pendulum Animation)

Aim
To write a PYTHON program to solve the second order differential equation and to simulate the pendulum motion.

Theory
A differential equation is an equation with a function and one or more of its derivatives [1].

`y+(dy)/(dx)=5x`: It is a differential equation with the function `y` and its derivative `(dy)/(dx)`.

The Order of the differential equation is the highest derivative.

A simple representation of Second order differential equation is,
`(d^2y)/(dx^2)+P(x)*(dy)/(dx)+Q(x)*y=f(x)`
where, `P(x)`, `Q(x)` and `f(x)` are functions of `x`.

For the function to be homogeneous, `f(x)=0`.
i.e., `(d^2y)/(dx^2)+P(x)*(dy)/(dx)+Q(x)*y=0`

Also, if the functions `P(x)` and `Q(x)` are constants `p` and `q`, then,
`(d^2y)/(dx^2)+p*(dy)/(dx)+q*y=0`

In engineering, ODE is used to describe the transient behaviour of a system. An example is a simple pendulum.

The Equation of Motion of a Simple Pendulum
A pendulum is body suspended from a fixed support so that it swings freely back and forth under the influence of gravity [2]. So, when a pendulum is displaced sideways from its resting or the equilibrium position, it is subject to a restoring force due to gravity that will accelerate it back towards the equilibrium position. So, the pendulum will swing back and forth with periodic motion, if it displaced to an initial angle and released [3].

Figure 1: Representation of Simple Pendulum

The simple pendulum consists of a mass `m` hanged from a straight string of length `L` and fixed at a pivot point. The way the pendulum moves depends on the Newton's second law. The differential equation which represents the motion of a simple pendulum is,
`(d^2theta)/(dt^2)=b/m*(d theta)/(dt)+g/l*sintheta=0` ... ... ... ...(i)

where, `g` is the acceleration due to gravity, `b` is the damping coefficient, `l` is the length of the pendulum and `theta` is the angular displacement.

This system is equivalent to the mass attached to spring and damper system. Is has various applications such as, in vehicle dynamics, suspension system can be modelled as a set of mass, spring and damper using the equation of motion of this system. This reveals the significance of the solution to this equation.

Generally, the higher order ODEs are solved by breaking it into first order ODEs. For solving the second order ODE, same approach can be used in PYTHON by making use of the odeint command.

For solving Equation (i),
Let `theta_2=(d theta_1)/(dt)`... ... ... ... ... ... ... ... ... ... ... ... ... ...(ii)

when equation (ii) is equated with equation (i), we get,
`(d theta_2)/(dt)=-b/m*theta_2-g/l*sintheta_1`

i.e., `d/(dt)[[theta_1],[theta2]]=[[theta_2],[-b/mtheta_2-g/lsintheta_1]]`
where the matrix on the LHS is the solution vector because it is this values we are finding out using the odeint function.

Objectives

1. Understanding the theoretical concept of Simple pendulum and the solving of the Ordinary differential equation.
2. Understanding the odeint command in PYTHON that to be used for solving the ODE in the PYTHON program.
3. Writing a program in PYTHON to solve the equation of motion of simple pendulum (second order ODE) and to plot the angular displacement and angular velocity graph with respect to time.
4. Modifying the program to simulate the pendulum motion and to create an animation of it.

Program to solve ODE

The programming is done in PYTHON.

# Solving 2nd order Ordinary differential equation

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

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

# Initial Condition
i_c = [0,3]

# Time Points
t_span = np.linspace(0,20,1000)		#motion from 0 to 20 seconds

# Defining ODE function
def model(theta,t_span,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

# Solving ODE
theta = odeint(model,i_c,t_span,args=(b,g,l,m))

# Plotting
plt.plot(t_span,theta[:,0], 'b-', label='Angular Dispacement, rad')
plt.plot(t_span,theta[:,1],'r--', label='Angular Velocity, rad/s')
plt.ylim([-4,4])
plt.xlim([0,20])
plt.ylabel('Plot')
plt.xlabel('Time, s')
plt.legend(loc='best')
plt.show()

Each step of the program is explained with comments.

The function model takes the value b, g, l and m as its arguments. This function has two inputs, t_span and initial boundary conditions, apart from the arguments.

Output of the Program

odeint basically uses a Runge-Kutta method of order 4 (with smart variants). Since the Runge-Kutta method is a first order ODE solver and hence, so is odeint. Thus, ODE to be integrated by odeint must be provided as a set of first order ODEs. (Refer theory section to know how the equation of motion of simple pendulum, which is 2nd order ODE is converted into a set of 2 first order ODE by using the variable `theta_2`. The initial condition would be the value of the function at t=0. Here, initial condition is provide using the variable i_c where the angular displacement and angular velocity are provided (here, 0 and 3).

Program to simulate the pendulum motion
The programming is done in PYTHON

# Creating the animation of Simple pendulul (2nd Order ODE)

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

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

# Initial Condition
i_c = [0,3]

# Time Points
t_span = np.linspace(0,20,300)		#motion from 0 to 20 seconds

# Defining ODE function
def model(theta,t_span,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

# Solving ODE
theta = odeint(model,i_c,t_span,args=(b,g,l,m))

# ct is the count of no. of plots
ct = 0

# Plotting
for i in range(0,len(theta)):
	x1 = l*math.sin(theta[i,0])
	y1 = -l*math.cos(theta[i,0])
	
	t = str(round(t_span[i],3))
	ct = ct+1

	plt.figure()
	plt.plot([0,x1],[0,y1],'m')
	plt.plot(x1,y1,'b',marker='o',markeredgewidth=6)
	plt.plot([-0.5,0.5],[0,0],'k',linewidth=3)
	plt.xlabel('x-axis')
	plt.ylabel('y-axis')
	plt.text(0.7,0,'Time = ')
	plt.text(1.1,0,t)
	plt.title('Simple Pendulum Python Animation')
	plt.xlim([-1.7,1.7])
	plt.ylim([-1.4,0.7])
	filename = str(ct)+'.png'
	plt.savefig(filename)

For creating an animation, a for loop is generated.
The variable result will have two columns after solving the differential equation using odeint. The first column represents the displacement angle for each value of time t. Similarly, the second column of the variable result contains the angular velocity for each value of t. The position of the pendulum is found from the first column which represent angular displacement and is then translated to normal plot using simple trigonometry. For representing the ball of the pendulum, again the point corresponding to the angular displacement has been plotted and the 'o' marker is used. For indicating the current time in the simple pendulum animation, the text command is used to indicate the required text at a particular point on the plot. For creating a movie, the variable ct has been used to indicate the frame number.

Output of the Program
The video output of the corresponding program is shown below:

Video link: https://youtu.be/QakzIKNB9FA

Conclusion
The PYTHON program to solve the second order ODE has been successfully created.
The governing equation of the simple pendulum has been solved and graphs showing the variation of position and angular velocity with time has been plotted.

The odeint function has been used to solve the ODE and the pendulum motion has been simulated. From the graphs and animation, it shows that the pendulum motion is more and continuous due to the low damping coefficient. Increasing the damping coefficient would increases the resistance to flow.

References
[1] https://www.mathsisfun.com/calculus/differential-equations-second-order.html#:~:text=where%20P(x)%2C%20Q,a%20wider%20range%20of%20functions.
[2] https://en.wikipedia.org/wiki/Pendulum_(mathematics)
[3] https://www.acs.psu.edu/drussell/Demos/Pendulum/Pendula.html
[4] https://www.acs.psu.edu/drussell/Demos/Pendulum/Pendulum.html#:~:text=By%20applying%20Newton's%20secont%20law,enough%2C%20so%20the%20small%20angle
[5] http://hyperphysics.phy-astr.gsu.edu/hbase/pend.html
[6] https://www.khanacademy.org/science/ap-physics-1/simple-harmonic-motion-ap/simple-pendulums-ap/a/simple-pendulum-ap1
[7] https://in.mathworks.com/help/matlab/ref/subplot.html