SOLVING SECOND ORDER ODE

INTRODUCTION

ODE is used to describe the behavior of the system. Here the motion of a simple pendulum is described in second order ODE using Newton\'s second law of motion as\'

`(d^2theta)/dt^2-(b/m)*(d theta)/dt-(g/l)*sintheta`

where

g = gravity in m/s^2,

L = length of the pendulum in m,

m = mass of the ball in kg,

b = damping coefficient.

PYTHON can be used to solve this ODEs and simulate how the pendulum moves using the program.

EXPLANATION OF THE PROGRAM

First we need to solve the ODE by giving the input values to the parameters, which are length, mass, acceleration due to gravity, damping coefficients and initial conditions for the differential equation(theta and time).

first to solve the ode and plot the graphs we need to import some libraries. import command is used to import libraries numpy as np, from scipy.integrate import odeint, math, mathplotlib.pyplot as plt.

ODE AS A FUNCTION

To solve the ODE we need to first define a function to describe the ODE. Here I used model as fuction name with input parameters t, theta, b, g, l, m.

describing the ODE,

I took theta1 is equal to theta[0], here 1 in theta[1] describes index of the array which is used to input the value ODE. Similarly I took theta2 is equal to theta[2].

Here is the logic to descibe the ODE in two first order ODE

`theta = theta_1`

`(d theta_1)/dt = theta_2`

`(d theta_2)/dt = -(b/m)*theta_2 - (g/l)*sin(theta_1);`

finally they can be represented as one column matrix

`(d theta)/dt = [(d theta_1)/dt ;(d theta_2)/dt]`

SOLVING THE ODE

Using odeint command we can solve this ODE as odeint(model,theta_0,t,args = (b,g,l,m))

PLOTTING THE RESULTS

plt.plot command is used to plot displacement with time and angular velocity with time. plt.savefig is used to save the figure. 

Again plt.plot is used to plot position vs angular velocity and plt.savefig to save the figure. 

ANIMATING THE PENDULUM MOTION

For loop is used to take each value of theta or position of the pendulum. If x0 and y0 are fixed points of the pendulum. x1 and y1 are calculated which are positions of the pendulum ball which changes with respect to angular velocities and damping coefficient. For each value of theta x1 and y1 are calculated which are then plotted and saved each time the for loop is run. 

When the program is run all the figures will be saved and these figures can be used to create an animation using image Magick tool. 

Use the command windows+R to open Run and type cmd to open command prompt.

Type cd and location of all the images.

Type dir to see what are there in that file location.

The image numbers may not be in sequence, If the are in sequence type magick *.png output.mp4.

If they are not in sequence the type magick \"%d.png\"[1-200] output.mp4

\"\"\" 
	ODE
\"\"\"

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

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
l = 1
m = 0.1
b = 0.02
g = 9.81

#initial conditions 
theta_0 = [0,5]

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

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\\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.savefig(\'fig1.png\')
plt.figure()
plt.plot(theta[:,0],theta[:,1])
plt.ylabel(\'Position\')
plt.xlabel(\'Angular velocity\')
plt.savefig(\'fig2.png\')


ct = 1
for theta in theta[:,0]:	
	x0 = 0
	y0 = 0
	x1 = -l*math.sin(theta)
	y1 = -l*math.cos(theta)
	plt.figure()
	plt.plot([-1,1],[0,0],\'black\',linewidth=5)
	plt.plot([x0,x1],[y0,y1],\'red\')
	plt.plot(x1,y1,\'o\',markersize=8)
	plt.xlim(-1.5,1.5)
	plt.ylim(-1.4,0.5)
	filename = str(ct)+\'.png\'
	
	plt.savefig(filename)
	ct=ct+1

RESULTS