Week 3 - Solving second order ODEs

Objective: To solve the 2nd order ode that represents the motion of a pendulum.

Problem statement: The solvers provided in scipy library of python like 'odeint' only solves 1st oder ode.

Solution Procedure: We create a function in which we convert the 2nd order ode into two 1st order ode's and return these ode's in matrix form so that the solver can solve both the ode's at the same time.

Procedure:

Step 1: Import the numpy, scipy, matplotlib, math and sys libraries.

Step 2: Set recursion limit to 1500 to avoid recursion error.

Step 3: Define the function that returns the 2 ode's in matrix form.

Step 4: Define the physical properties of the pendulum like length, mass, damping coefficient and acceleration due to gravity.

Step 5: Define the initial position and velocity imparted to the bob.

Step 6: Define the time interval for the plot and divide the intervals into equal partitions.

Step 7: Solve the ode's with 'odeint' function.

Step 8: Plot position and velocity vs time graph.

Step 9: To make the animation we need to plot the position of bob with respect to the fixed point of pendulum. Using the angular positions stored in 1st column of solution given by 'odeint' program and length of string 'l' we can calculate the x and y coordinates of the bob using following formulae:

x = length * cos(angular position)

y = length - length * sin(angular position)

Step 10: To plot the velocity vector, initially we calculate the absolute velocity from the angular velocity stored in 2nd column of solution given by 'odeint' function using formula:

absolute velocity = length * angular velocity

using absolute velocity and angular position of the bob, we can determine the length of velocity vector in x and y directions using following formulae:

length in x direction = absolute velocity * cos(angular position of bob)

length in y direction = absolute velocity * sin(angular position of bob)

We also scale the velocity vector to 1/10th of its original size.

Step 11: We save all the plots to create a gif animation.

import numpy as np
from scipy.integrate import odeint
import matplotlib.pyplot as plt
import math
import sys
sys.setrecursionlimit(1500) # to avoid recursion error

# function that returns dz/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	   #damping coefficient
g = 9.81	#acceleration due to gravity
l = 1	# length of string
m = 1	# mass of the bob

# initial condition (position, velocity)
theta_0 = [0,3]

# time
t = np.linspace(0,20,150)	# dividing time interval into 150 equal partitions for smooth animation

# 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}=\theta_2$')	# displacement vs time plot
plt.plot(t,theta[:,1],'r--',label=r'$\frac{d\theta_2}{dt}=-\frac{b}{m}\theta_2-\frac{g}{L}sin\theta_1$')	# velocity vs time plot
plt.ylabel('Plot')
plt.xlabel('Time')
plt.legend(loc='best')
plt.show()

# simple pendulum aimation
ct = 1
for z in theta[:,0]:
	x = l * math.sin(z)	# x coordinate of position of the bob
	y = l - l * math.cos(z)	# y coordinate of position of the bob

	absolute_velocity = theta[ct-1,1] * l 	# absolute velocity = angular velocity * length of string
	v_x = absolute_velocity * math.cos(z) / 10	# length of velocity vector in x direction (scaled to 1/10th size)
	v_y = absolute_velocity * math.sin(z) / 10	   # length of velocity vector in y direction (scaled to 1/10th size)
	
	filename = str(ct)+'.png'	# assigning names to plots
	ct = ct+1

	plt.figure()
	plt.plot([0, x], [l, y], linewidth = 2, color = 'blue')		# string
	plt.plot(0,l,'o', markersize = 5, markerfacecolor = 'white', markeredgecolor = 'black')		# fixed end of pendulum
	plt.plot(x,y,'.', markersize = 25, color = 'red')		# bob of pendulum
	plt.plot([0, 0],[-.1,1.1], linestyle = 'dashdot', linewidth = 0.8, color = 'black')		# axis of motion of pendulum
	plt.arrow(x,y,v_x,v_y,head_width = 0.02, color = 'black')		# velocity vector of pendulum
	plt.xlim([-1.2, 1.2])
	plt.ylim([-0.2, 1.2])
	plt.savefig(filename)		# saving the plots to create animation

Output:

Youtube link of the animation:

https://youtu.be/7KGVk-3JyqE