Objective: To solve the second-order ODE equation and create an animated gif video of a simple pendulum.
Solution:
The second order ode represents the motion of the simple pendulum and the equation is as follows:
`(d^2 theta)/dt + b/m (dtheta)/dt + g/l sin theta = 0` ....... (1)
Python uses scipy integration to solve ODE but it solves only simple single-order ODEs
Equation (1) can be modified as follows:
let `theta = theta_1` and
`(d theta_1)/dt = theta_2` ------- (2)
`(d theta_2)/dt = b/m theta_2 - g/l sin(theta_1)` --- (3)
Procedure:
import math to solve mathematical functions
import numpy to create and use arrays elements
import odeint from scipy.integration to solve ODE
import os to get the folder details and save figures
import matplotlib.pyplot to create the graph figures
Create a function model with equations (2) and (3) to solve eqn(1) and get the values for theta_1 and theta_2 which are angular displacements and angular velocities.
Pass arguments
b bulk coffeficient = 0.01
g acc. due to gravity = 9.8 m/s^2
l length = 1 m
m mass of bob = 0.1 kg
After plotting and saving the plot figures in os,
the figures are created as gif video using imageio package, which is must be installed to access.
The code is as follows in detail:
Code to solve the motion of simple pendulum equation and create an animated gif file
import numpy as np
import math
from scipy.integrate import odeint
import matplotlib.pyplot as plt
import imageio
import os
function that returns dtheta/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
#parameters
b= 0.02
g= 9.81
l= 1
m= 0.1
#initial condition
theta_0 = (0,3)
time points
t = np.linspace(0, 20, 150)
#solve ode
theta = odeint(model, theta_0, t, args =(b, g, l, m))
#plot results
plt.plot(t, theta[:, 0], 'b-', label=r'$\frac{d\theta_1}{dt} = \theta2$')
plt.plot(t, theta[:,1], 'r--', label=r'$\frac{d\theta_2}{dt}= -\frac{b}{m}\theta_2-\frac{g}{t}sin\theta_1$')
plt.ylabel('plot')
plt.xlabel('time')
plt.legend(loc='best')
plt.show()
ct=1
filenames=[]
for Theta in theta[:,0]:
#Co-ordinates of Fixed Point of String
x0 = 0
y0 = 0
#Co-ordinates of End Point of String
x1 = l*math.sin(Theta)
y1 = -l*math.cos(Theta)
filename = str(ct) + '.png'
filenames.append(filename)
ct = ct + 1
plt.figure()
Plotting Origin, string and ball of pendulum
plt.plot([-0.5,0.5],[0,0],'g',linewidth=5)
plt.plot([x0,x1],[y0,y1],'r')
plt.plot(x1,y1,'o',markersize=25)
#Limiting X-axis and Y-axis
plt.xlim([-1.5,1.5])
plt.ylim([-1.5,1.5])
#Title
plt.title('Animation of Simple Pendulum')
#To save the file name
plt.savefig(filename)
#imageio module to generate animated gif video using saved files
with imageio.get_writer('Pendulum.mp4', mode='I', fps=10) as writer:
for filename in filenames:
image = imageio.imread(filename)
writer.append_data(image)
for filename in set(filenames):
os.remove(filename) #removes the png files saved
Conclusion:
The simple pendulum motion equation, a second-order ODE has been solved and a animated video describing motion of simple pendulum with angular velocity 3 rad/sec has been created using python.