Solving 2nd order ODE

Objective: The objective of this report is the simulate the transient behaviour of a simple pendulum and the create the animation of its motion for a given time interval

Theory: The following equation (Ordinary Differential Equation) for the position of the bob w.r.t. to time was used in the program.

g: gravity in m/s^2 =9.81,

L: length of the pendulum in m =1,

m: mass of the bob in kg =1,

b: damping coefficient = 0.05

On solving the equation, different values of the bob (angular position) was obtained which was used to generate the images of the bob position for various values of time. The individual plots were then stitched used image magick to generate an animation of the pendulum motion

Program:

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

#function to return 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

b=0.05 # b: Damping Coefficient = 0.05
g=9.81 # g: accelerationn due to gravity= 9.81 m/s^2
l=1    # l: length of the pendulum = 1m
m=1    # m: mass of the bob = 1kg

#initial condition (Angular displacement = 0 & angular Velocity =3 rad/s)
theta0=[0,3]

#time points (simulating the motion for 0 to 20 sec, considering 150 time steps between 0 to 20 sec)
t = np.linspace(0,20,150)

#solving the ODE
theta = odeint(model,theta0,t, args= (b,g,l,m))

#plotting 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}{L}sin\theta_1$')
plt.xlabel('time')
plt.ylabel('Plot')
plt.legend(loc='best')
plt.show()

x0=0
y0=0
cnt=1

#Creating the plots for animation of the pendulum motion
for t1 in theta[:,0]:
	x1 = l*math.sin(t1)
	y1 = -l*math.cos(t1)

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

	plt.figure()
	plt.plot([-0.5,0.5],[0,0],[-0.5,0.5],[0.05,0.05],[-0.5,0.5],[0.1,0.1], [-0.4,-0.4],[0,0.1],[-0.2,-0.2],[0,0.1],[0.2,0.2],[0,0.1],[0.4,0.4],[0,0.1])
	plt.plot([x0,x1],[y0,y1])
	plt.plot(x1,y1,'o', markersize=10)
	plt.xlim([-1.5,1.5])
	plt.ylim([-1.5,1])
	plt.savefig(filename)

The program generated 150 images of the pendulum position for each values of time.

Note: 150 is the size of the time interval array considered for this problem

Post this the images were moved to a new folder. Using command prompt navigated the folder where the images were saved and used the command magick *.png bobanimation.gif. The command stiched the inidividual images Image Magick to create to animation file

Results:

Fig: Displacement and Velocity Graph

Inference: the initial conditions provided can be checked from the graph:

simulation time 0 to 20seconds: the graph was also generated for the time interval 0 to 20 sec

angular displacement = 0 and angular velocity = 3 rad/sec. In the plots the angular displacement and angular velocity starts from 0 and 3rad/sec respectively.

Also, due to the dampling coefficient, the decrease in the amplitude of displacement and velocity can be observed in the graph

Error Generated: The error generated was a image stitching error. The code used earlier generated the file names in the numbers from 1 to 150 (length of time array) due to which when the images were stitched using image magick, the break in the pendulum motion took place due to improper stitching of the file. This was generated as the image files were not stiched in proper order.

To correct this error the command for image file generation was changed to:

filename = 'test%05d.png'%cnt

The images generated had proper filenames and were stitched together to get a smooth animation as shown below:

Video: Pendulum Motion Video(After troubleshooting)