Week 3 - Solving second order ODEs

Aim:

To write a program in Python

1. To solve the following ODE which represents the equation of motion of a simple pendulum with damping: `(d^2 θ)/dt^2 +b/m (dθ)/dt+g/L sinθ=0`
2. To simulate the motion between 0 - 20 seconds, for the given conditions:
   Length of pedulum, L = 1m
   Mass of the ball, m =1 kg
   Damping co-efficient, b = 0.05
   Acceleration due to gravity = 9.81 m/`sec^2`
   At time, t=0, angular displacement = 0, angular velocity = 3 rad/sec.

Procedure:

Step 1: Solving the given ODE

• The given second order ODE is broken down into two first order ODEs.
  ODE 1: Represents angular displacement as a function of time
  `(dθ_1)/dt=θ_2`
  ODE 2: Represents angular velocity as a function of time
  `(dθ_2)/dt=-b/m θ_2-g/L sinθ_1`
• Then these two first order ODEs are solved as described in the Python code shown below.

Step 2: Simulating the motion of the pendulum

• In order to simulate the motion of the pendulum we need to generate an array of images and stitch them to form an animation. We use "for" loop to generate these images which depict the position of the pendulum at a given time.
• The pivot of the pendulum is set as the origin (0,0) and the co-ordinates of the ball at the other end is calculated with repect to time based on the "sin" and "cos" relations as shown in the Python code below.
• Plots are plotted and the animation depicting the motion of the pendulum is created by using the images generated using "for" loop.

Python code:

# Title - Solving Second Order ODEs

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

#function that returns dθ/dt
def model(theta, t, b, g, l, m):

	theta1 = theta[0]
	theta2 = theta[1]
	dtheta1_dt = theta2 	#ODE1 - represents angular displacement
	dtheta2_dt = -(b/m) * theta2 - (g/l) * math.sin(theta1) 	#ODE2 -  represents angular velocity
	dtheta_dt = [dtheta1_dt, dtheta2_dt] 	#combining ODE1 and ODE2
	
	return dtheta_dt

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

#Initial condition for angular displacement and angular velocity at t = 0
theta_0 = [0,3]

#Creating time points using linspace command
n_t = 300	#number of values for time
t = np.linspace(0, 20, n_t)

#Solve the ODE using odeint
theta = odeint(model, theta_0, t, args = (b, g, l, m))

#Plot the solution of the two first order ODEs vs Time
plt.plot(t, theta[:,0], 'b-', label = r'$\frac{d\theta_1}{dt} = \theta_2$')
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.plot([-5, 25], [0, 0], 'grey', linestyle = ':')
plt.xlabel('Time')
plt.ylabel('Plot')
plt.xlim(-1, 21)
plt.ylim(-3.5,3.5)
plt.legend(loc = 'best')
plt.savefig('W3_ODE_Time')


#Co-ordinates of the pivot of the pendulum
x0 = 0
y0 = 0

ct = 1	#counter variable

for i in range(0, n_t):
	
	#Co-ordinates of the ball of the pendulum
	x1 = l * math.sin(theta[i,0])		#theta[i,0] returns the value of theta at every instance of time created using the linspace command
	y1 = -l * math.cos(theta[i,0])

	filename = 'pendulum%05d.png' % ct 	#unique file name for each png image with 5 digits. 'ct' is the counter variable in loop
	ct = ct + 1			#incrementing 'ct' after each execution of the loop
	
	#Plotting the horizontal surface, the pivot point, ball of the pendulum, and the string
	plt.figure()
	plt.plot([-0.5,0.5], [0,0], 'black')
	plt.plot([0,0], [0.25, -1.25], 'grey', linestyle = ':')
	plt.plot(x0, y0, 'o', markerfacecolor = 'black')
	plt.plot([x0,x1], [y0,y1], 'green')
	plt.plot(x1, y1, 'o', markersize = 20, markerfacecolor = 'tab:blue', markeredgecolor = 'lime')
	plt.xlim(-1.5, 1.5)
	plt.ylim(-1.5, 0.5)
	plt.title('Motion of a Simple Pendulum')
	plt.savefig(filename)		#inputting 'filename' as the argument in order to generate images with unique filenames

Results:

• Solution of the given ODE (Plot):

• Motion of the simple pendulum:

Youtube link: Motion of the simple pendulum