Week 3 - Solving second order ODEs

SOLVING 2ND ORDER ODE USING PYTHON

INTRODUCTION:

ODE's represent the transient behaviour of a system, for example as of a pendulim.

The pendulum has a string attached at one end and a bob is attachedat the other end.

If an initial velocity is given , the bobs starts swinging but the amplitudeof osciattion decreases over time.

This motion of the bob is determined by the differential equation given below:

SOLUTION:

To simulate the motion of the pendulum, the above ode has to be solved and to solve this ode, below is the python code

import numpy as np
import math
from scipy.integrate import odeint # function used for solving the ode
import matplotlib.pyplot as plt 


# function for dz/dt
def model(theta,t,b,l,g,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   
'''
 this function returns the right hand side of the matrix , which needs to be solved using the odeint function
'''
# inputs for the ode

# Damping coefficient
b = 0.05 

# Length of the pendulum
l = 1

# Gravity
g = 9.81

# Mass of the ball
m = 1

# initial_condition
theta_0 = [0,3]

"""
The index zero of this matrix gives us the initial angular displacement of 0 and initial angular velocity of 3

"""

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

"""
This function linspace generates 150 numbers in between the interval 0 to 20
"""

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

"""
The funciton ODEINT outputs the value of theta as an array
"""

# plot for animation
X = theta[:,0]
ct = 1
for theta in X:
	x1 = 0;y1 = 3
	x2 = 3*math.sin(theta)
	y2 = -(2*math.cos(theta))
	file = 'test%05d.png'%ct
	plt.figure()
	plt.plot([-3,3],[3,3],linewidth = 10,color = 'blue')
	plt.plot([x1,x2],[y1,y2],linewidth = 4, color = 'orange')
	plt.plot(x2,y2,'o',markersize =12, color = 'pink')
	plt.xlim([-3,3])
	plt.ylim([-3,3])
	plt.savefig(file)
	ct = ct+1

Now the program for plotting the output after solving the ode:

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

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


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

#initial_condition
theta_0 = [0,5]

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

#solve ODE
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}=\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.ylabel('Plot')
plt.xlabel('time')
plt.legend(loc = 'best') 
plt.show() 	

I have explained the meanings of the functions used in the code for simulating the pendulum.

OUTPUT:

THE YOUTUBE LINK FOR THE PENDULUM SIMULATION: https://www.youtube.com/watch?v=P-3Z732ID38

THE PLOT OBTAINED AFTER SOLVING THE ODE WITH THE GIVEN CONDITIONS: