Week 3 - Solving second order ODEs

Aim:
Solving differential equation using ODEINT and simulating motion of a simple pendulum .
Requirements:
In the given equation,

g = gravity in m/s2,
L = length of the pendulum in m,
m = mass of the ball in kg,
b=damping coefficient.
L=1 metre,
m=1kg,
b=0.05.
g=9.81 m/s2.
Simulate the motion between 0-20 sec, for angular displacement=0,angular velocity=3 rad/sec at time t=0.

Python Code:

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

# function that returns d_theta/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


# inputs	

b=0.05
g=9.81	
l=1
m=0.1

#initial conditions
theta_0=[0,5]

# time intervels
t=np.linspace(0,20,150)

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

# Animation
THETA3=[]
for theta3 in theta[:,0]: #If you have a 2-dimensional list/matrix/array, this notation will give you all values in column 0 (from all rows).

	THETA3.append(theta3)

ct=1
for THETA4 in THETA3:

	x0=0
	y0=0
	x3=-1.5
	y3=0
	x4=1.5
	y4=0

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

	filename=str(ct)+'.png'
	ct=ct+1
	
	plt.figure()
	plt.plot([x3,x4],[y3,y4])  # Horizontal support of a  pendulum
	plt.plot([x0,x1],[y0,y1],x1,y1,'o',markersize=10,color='red')
	plt.xlim([-2,2])
	plt.ylim([-2,2])
	plt.title('oscillating Pendulum')
	plt.savefig(filename)
	


Procedure :

>In order to solve the second order diffential equation we have to convert it in to two single order differntial equation.
>To solve the ODE function,import odeint module from scipy module.
>Define a function model  with inputs,theta,t,b,g,l and m which returns dtheta/dt as a matrix.
>The first column stores angular displacment and  2nd column stores velocity.
>The time interval is defined by numpy module using linspace command for 150 values.
>To solve the ODE equation we use odeint command.The output of the diffential equation is stored in the theta.
>Then we can observe the plotted results where at t=0 velocity is 5m/s which can be depicted by orange lines.And the displacement=0 can be observed by a blue dotted lines.

>To create an animation, Using for loop we draw the values of i from theta which shows the movement of the pendulum for every value in theta which is displacment. 

Results:

> We can observe the time plot for dispalcement and velocity.
>We can also observe the movement of the simple pendulum in a given time interval.