Week 3 - Solving second order ODEs

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

#defining 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

# input parameter
b=0.05    	#damping coefficient
g=9.81 		#gravity in m/s2
l=1 		#length of pendulum
m=1 		#mass of ball

#initial condition 
theta_0 = [0,3] 

#time points
t = np.linspace(0,10,150) 

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


# plotting graphs
plt.figure(1)
plt.plot(t,theta[:,0],'b-',label=r'$\frac{d\theta_1}{dt}=\theta2$') # theta [:,0] -- displacement values
plt.plot(t,theta[:,1],'r--',label=r'$\frac{d\theta2}{dt}=-\frac{b}{m}\theta_2-\frac{g}{L}sin\theta_1$') # theta[:,1]--velocity values
plt.xlabel('Time [sec]')
plt.ylabel('Angular Velocity [rad/sec]')
plt.grid()
plt.legend(loc='best')
plt.title ('Damping of simple pendulum using second order ODE')
plt.show()



# animation plot 
x0 = 0
y0 = 0

ct = 1
for T in theta[:,0]:

	x1 = l*math.sin(T)
	y1 = -l*math.cos(T)

	filename = str(ct) + '.png'
	ct=ct+1


	#plot results 
	plt.figure()
	plt.plot([-0.5,0.5],[0,0],linewidth=5)
	plt.plot([x0,x1],[y0,y1],linewidth=2)
	plt.plot(x1,y1,'o',markersize=20)	
	plt.xlim(-2,2)
	plt.ylim(-2,2)
	plt.grid()
	plt.title('Pendulum Animation')	
	plt.savefig(filename)

Problem Statement :

In Engineering, ODE is used to describe the transient behavior of a system. A simple example is a pendulum

The way the pendulum moves depends on the Newtons second law. When this law is written down, we get a second order Ordinary Differential Equation that describes the position of the "ball" w.r.t time.

Similarly, there are hundreds of ODEs that describe key phenomena and if you have the ability to solve these equations then you will be able to simulate any of theses systems.

Your objective is to write a program that solves the following ODE. This ODE represents the equation of motion of a simple pendulum with damping

Objective :

Write a program in Python that will simulate the pendulum motion, just like the one shown in the start of this challenge.

Solution procedure :

Stepwise explaination of second order ODE pendulum problem

1. Step number one is our most important and daily habit step which is importing libraries , so for these challenge we need 4 libraries which are numpy, scipy, ploting, math all these libraries we need for these code but before this step there is one important work to do that is theory base which is to convert second order derivative into single order derevetive .
2. Our second step is creating function so by using def command we define the function and the model means the function which represents the ode function in that function we have to give the things which we do on theory base which is converted second order ODE into single order ODE and we have to do these things by keeping one thing in mind which is we have to do these all things under the function statement.
3. Now in our next step leaving the function statement we have to define the parameters which we used in function statement. our parameters are b(damping coefficient) , g(gravity), L(length of pendulum), m(mass of ball)
4. Our next step is defining the theta_0 which is 90 degree but we cant give direct value cause it don’t take that value in degree so for that we have to use pi/2 and pi comes under math libraries so we have to use statement math.pi/2 to equalize the theta_0.
5. After defining theta we have to create the time array and to create time array our numpy libraries has to come in play because to create time array we need linspace command and our linspace command lies in numpy libraries so our statement be like t=np.linspace(0,10,150) as we need 150 values between 0-10.
6. In this step we calls our ode function which we creates above in very first step so equalizing theta we calls odeint(model,theta_0,t,args=(b,g,l,m)) in this statement odeint is ode function calling command and rest of parameters we all knows as I tells in above step 2 .
7. Now moving to next step which is ploting the velocity vs displacement and for that we use sets of command which are our everyday used command but some statements are new which we elaborates below,

Statement1. plt.plot(t,theta[:,0],'b-',label=r'$\frac{d\theta_1}{dt}=\theta2$') so in this line theta[:,0] represents the value of displacement while next data given is all about give the label to graph line

Statement2. plt.plot(t,theta[:,0],'r--',label=r'$\frac{d\theta2}{dt}=-\frac{b}{m}\theta_2-\frac{g}{L}sin\theta_1$') this statement is although same as above command one little difference hera is that theta[:,1] is representing the value of velocity. Rest of the commands given below Are daily use ploting commands

plt.xlabel('Time [sec]')

plt.ylabel('Angular Velocity [rad/sec]')

plt.grid()

plt.legend(loc='best')

plt.title ('Damping of simple pendulum using second order ODE')

plt.show()

8. Our next step is creating the animation of pendulum so for that we create the for loop in which we have to define the changing theta with respect to its coordinates so as our one point is fixed our x0 and y0 coordinates are equal to zero and as our other end is moving so our equation for that is “x1 = l1*math.sin(T)” in which T is range of theta which we defines while creating for loop by giving the statement “for T in theta1” and for y cordionates y1=-l1*math.cos(T)
9. Now our next step is creating the plotings and arrange that plots in frame by frame manner and for that we gives the statement which we gives below

filename = str(ct) + '.png'

ct=ct+1

10. After that we have to create plot for those frame so we use set of ploting commands in for loop just like we used in rocker arm challenge code . in this set first we defines the axis then in another command we defines the coordinates , after that we creates the ball by using ‘ O ‘ and marker size

plt.figure()

plt.plot([-0.5,0.5],[0,0],linewidth=5)

plt.plot([x0,x1],[y0,y1],linewidth=2)

plt.plot(x1,y1,'o',markersize=20)

plt.xlim(-2,2)

plt.ylim(-2,2)

plt.grid()

plt.title('Pendulum Animation')

plt.savefig(filename)

11. And now our final step is to steech all those frames or plot in one thread means to create one animation video by using those all 150 plots so by using online method we uploads all our 150 plots their and then that site crates the animation video using those plots like wise we get our complete animation file

Errors faced while solving

1. I faced error in line no 61 where I was defining the axis for pendulum I used the statement plt.axis([-0.5,0.5],[0,0],linewidth=5). After talking to support engineer I got the right point that in the place of axis I have to use plot so our renewed statement is

plt.plot([-0.5,0.5],[0,0],linewidth=5)

Result :

as a result we create a code and animation file of defined mass and length pendulum