Solving second order ODEs using Python

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

plt.rcParams.update({\'figure.max_open_warning\': 0})

#ODE function
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
g = 9.81
l = 1
m =1

#initial condition 
theta_0 = [0,3]

#time points
n = 1200
time = np.linspace(0,20,n)
t_range = np.linspace(0,n-1,n,dtype=\'int\')


#solving ODE
theta = odeint(model,theta_0,time,args = (b,g,l,m))


ct = 1
x0 = 4
y0 = 4
displacement = []
t_range_disp = []


for i in t_range:
	x1 = x0-l*(math.sin(theta[i,0]))
	y1 = y0-l*(math.cos(theta[i,0]))
	displacement.append([x1,y1])
	vel = theta[i,1]/8
	m = (x1-x0)/(y0-y1)
	x2 = (-vel/pow((1+pow(m,2)),0.5))+x1
	y2 = ((-m*vel)/pow((1+pow(m,2)),0.5))+y1
	filename = \'test%04d.png\' %ct
	ct = ct+1
	t_range_disp.append(i)
	
	
	plt.figure()
	
	#Plotting the dotted displacement
	for j in t_range_disp:
		plt.plot(displacement[j][0],displacement[j][1],\"g\",marker=\".\",markersize = 1)
	
	#Plotting the pendulum string
	plt.plot([x0,x1],[y0,y1],\"r\")

	#Plotting the pendulum velocity
	plt.plot([x1,x2],[y1,y2],\"k\")

	#Plotting the ball
	plt.plot(x1,y1,marker=\"o\",markersize=10,markeredgewidth=2,markeredgecolor=\"b\",markerfacecolor=\"None\")
	
	plt.xlim([3,5])
	plt.ylim([2.5,4.5])
	plt.savefig(filename)
	plt.cla()
	plt.close()

First convert the 2^nd degree ODE into two 1^st order ODE.

Create the ODE function (‘model’) which defines the two 1^st order ODE.

Start with pendulum programming

Declare the inputs of g, l, m, b.

Declare the state of the pendulum at time = 0

Create a time array of time 0 to 20 sec

Use the ‘odeint’ solver to solve the ODE function as defined in ‘model’

Initialize the counter used for capturing the figures.

Declare the co-ordinates of fixed top.

Next, find the displacement of pendulum and its velocity at the given instant from the ‘theta’ matrix.

In a ‘for’ loop, plot a line between fixed point and pendulum co-ordinate for each instant calculated using geometry.

The ‘for’ loop makes the calculation for every time instant. However, we have to pass the element number to the loop for it to work. Hence, a ‘t_range’ array is made of number sequence ([0,1,2,3,4,….]) as integers.

NOTE: Since the velocity was much larger number, the length of the velocity vector was huge compared to the length of string. Hence for representation purpose, the velocity is divided by a factor of 8.

Since the velocity vector is perpendicular to the direction of the string, find the slope of the velocity vector by using the equation: m1*m2 = -1

Length of the vector is known from the ‘theta’, i.e. the velocity.

Using the two above known characteristics, an equation was formed using geometry to find the end point co-ordinate of the velocity vector.

Using the pendulum as starting point and the above found end point, a line is plotted between them to represent the velocity vector.

To plot the dotted displacement, the following steps are undertaken:

Every co-ordinate the pendulum has passed since the start of the loop till the present is marked.

For this we need to store the co-ordinates of pendulum free end (i.e. x2, y2) and also, we need an array of the size of previous iteration. (NOTE: Because I just know ‘for’ feature which works by going through every element of an array, I made a number sequence array ([0,1,2,3,4,5,6,7,….]) to achieve: [‘from i = 0 to 100’])

To store the co-ordinates an empty array ‘displacement’ is created outside the ‘for’ loop. Inside the loop, using the ‘append’ feature, the co-ordinates are stored in the array.

To make an array of number sequence of the size of previous iteration, an empty array ‘t_range_disp’ is made outside the main ‘for’ loop. Using the ‘append’ feature the value of ‘i’ is added to the array.

To plot the dotted displacement, another ‘for’ loop is incorporated inside the main ‘for’ loop. This for loop plots the co-ordinates stored earlier using the ‘t_range_disp’.

To generate a mathematically correct video of the pendulum, the length of the video had to be 20 sec. And for smoother video 60 fps was chosen. So, 60*20 = 1200 frames were generated.

The generated image sequence was stitched into a video using Adobe Premiere Pro.

Results:

Youtube link: https://youtu.be/0MxZhoaRla0

Problems faced during programming:

• The ‘for’ loop demanded integer value hence ‘t_range’ was created as integer data type.

• There was a memory overload warning when the number of figures exceeded 20. Hence, ‘plt.rcParams.update({\'figure.max_open_warning\': 0})’ was used. This was found after online search of the warning. Using this we get rid of the warning. However, the program stops when it has generated about 369 images. Hence, ‘plt.cla()’ and ‘plt.close()’ was used to release the memory used and thereby 1200 frames were generated successfully.

• The velocity vector line was too big compared to other elements. Hence it was reduced by dividing it by a factor of 8.