Week 3 - Solving second order ODEs

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

#Creating a function to return the required variables to be solved
def solve(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

#Inputs(Constants present in the equation to be solved) 
b=0.05
g=9.81
l=1
m=1

#Initial Condition(Given initial conditions for the equation) 
theta_0=[0,3]

#TimePoints(Creates a set of points from 0 to 20 with a total of 150 values )
t=np.linspace(0,20,150)


#SolvingODE
theta=odeint(solve,theta_0,t,args=(b,g,l,m))

#Plot Results(Plotting Displacement and Velocty obtained by solving the equation against time)
#First column of Theta contains all the displacement values, second column contains all the velocity values
plt.plot(t,theta[:,0])
plt.plot(t,theta[:,1])	
plt.xlabel("Time")
plt.ylabel("Plot")
plt.savefig("DispVelPlot")
plt.show()

#Animation-Code to plot the displacement points of the pendulum and combine the images to create a gif file
a=theta[:,0]
th=np.array(a)
count=1
for t in th:
	xpos=l*math.sin(t)
	ypos=l*math.cos(t)
	filename=str(count)+'.png'
	count=count+1
	plt.figure()
	plt.plot([0,xpos],[0,-ypos])
	plt.scatter(xpos,-ypos)
	plt.xlim([-2,2])
	plt.ylim([-2,2])
	plt.savefig(filename)

This code is used to solve the governing equation of- the damping motion of a pendulum and obtain the displacement and velocity at different instances of time. We first create a function to store the unknown quantities. We code in these quantities and the function returns them. The constant parameters in the equation such as-b,g,l,m are hard coded.The initial conditions given are stored in an array. The 'odeint' function in scipy module is used to solve the given equation. The equation,initial conditions,and the constant values are passed to the function as per the syntax. The equation is solved and the computed values of displacement and velocity are stored in the theta[] array. The first column theta[:,0] contains all the displacement values. The second column theta[:,1] contains all the velocity values.The values are extracted and a graph of Displacement and Velocity against Time is plotted as shown.

The next part of the code is to plot the displacements of the pendulum,and thus to create an animation of the motion.For each value of displacement , we take the position of the bob of the pendulum in terms of the sin and cos component. A counter is initialized to run through the loop and is used to create filenames to store the plots created ,in sequential order.The scatter() function is used at the xpos and ypos cooridnates ,to create a dot to represent the bob of the pendulum.An empty figure is create using figure(),the value is plotted and then saved with the appropriate sequential filename. Here the ypos representing y coordinate of the bob is taken as negative, because the pendulum has to be plotted downwards,ie, in the fourth quadrant. xlim() and ylim() functions are used to limit the x and y value range.

Errors Faced-

1.Initially took the y coordinate of the pendulum bob as positive. Got all the plots of the displacements upside down, as shown in the figure. Corrected this by taking negative of ypos and then plotting it.

2.Used the marker argument in the plot() function to create a dot (bob) of the pendulum. A dot was created at both ends of the line , instead of only one end, representing the bob. Corrected this by removing the marker command and using the scatter() function.