Solving 2nd order ODE in Python using ODEINT

OBJECTIVE: Solving differential equation using ODEINT and simulating motion of pendulum using solution from ODE.

Given ODE of 2nd order.

PROGRAMMING EXPLAINED:

1.To solve ODE of 2nd order, we have to convert above equation into two, single order ODE.

`(d(theta1))/(dt)=(theta2)`---->1

`(d(theta2))/(dt)+(b/m)*(theta2)+gsin(theta1)/(l)=0`------2

2. Define all constants like b,m,l and g.

3.Initial condition: diplacement and veocity>>[0,5]

4. Time was defined using numpy module for 150 values.

5. Differential equation was solved usind ODEINT command.

6.To solve the ODE we need to define our system that is model, initial conditions, time points and argumentts like m,g,etc.

7. Model is the function which define our ODE. The two single order ode was defined under this model function.

8. The two single order ODE was stored in single in variable dtheta_dt.

9. The output of the differential equation was stored in variable theta.

10. Output has two columns and multiple rows.

11.1st column stores displacement and 2nd stores velocity.

As we can see from figure at t=0 velocity =5 , indicated by orange line and displacement =0,indicated by blue line.

The motion was simulated by saving each plot for time between[1-150]. ImageMagick software was used to get GIF file which is shown below:

#Solving ODE in python using ODEINT 
#Luv Kumar

import numpy as np
from scipy.integrate import odeint
import matplotlib.pyplot as plt
import math
#function that return dtheta/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.02
g = 9.81
l = 1
m = 0.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])
plt.plot(t,theta[:,1])
plt.ylabel(\'plot\')
plt.xlabel(\'time\')
plt.legend()
plt.show()

ct = 1
for i in theta[:,0]:
	x0 = 0
	y0 = 0
	x3 = -1.5
	y3 = 0
	x4 = 1.5
	y4 = 0
	x1 = l*math.sin(i)
	y1 = -l*math.cos(i)
	filename = str(ct)+\'.png\'
	ct = ct +1
	plt.figure()
	plt.plot([x3,x4],[y3,y4])
	plt.plot([x0,x1],[y0,y1],x1,y1,\'o\',markersize=20,color=\'red\')
	plt.xlim([-2,2])
	plt.ylim([-2,2])
	plt.title(\'oscillating pendulum\')
	plt.savefig(filename)

CONCLUSION:

1. By solving ODE, we can understand the system by visualizing its motion.

2. For solving ODE , ODEINT command was used.

3. To solve higher order ODE, set of single order ODE was defined.