Second Order Ode s In Python

Program to solve second order differential equation

Introduction:

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.

Now the position(theta) of the bob w.r.t to time can be determined by using the second order Ode. Now the second order ode equation is `

b = damping coeffivient

m = mass of bob (kg)

l = length of string (m)

now to solve the second order differential equation and the find the value of the theta we will have to reduce this differential equation to first order which we have done in the function:

T1= T[0]
T2 = T[1]
dT1_dt = T2
dT2_dt = -(b/m)*T2 -(g/l)*math.sin(T1)
dTh_dt = [dT1_dt, dT2_dt]

here T(0) = theta

T(1) = theta1

so we have reduced the second order differentail equation in first order.

Now Program For Animation:

Function used in this Program:

np.linspace = used in python to create numeric sequences

ode int = solve the differential Equations

lim = to set the limit on x and y

marker = to create a sphere

Now for animation we have separated the displacement value and stored the theta value in a different variable and then the various positions it is save using savefig command now to plot the various file we have used the file = \'test%05d.png\'%ct command.

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


# function 
def model(T,t,b,l,g,m):
	T1= T[0]
	T2 = T[1]
	dT1_dt = T2
	dT2_dt = -(b/m)*T2 -(g/l)*math.sin(T1)
	dTh_dt = [dT1_dt, dT2_dt]
	return dTh_dt

# inputs
b = 0.05 # damping co-efficient

l = 1 # length of pendulum in m

g = 9.81 # gravity in m/s^2

m = 1 # mass of ball in kg

# initial conditions
T_x = [0,3]

# time points
t = np.linspace(0,20,150)
# solve ode
T = odeint(model,T_x,t,args=(b,l,g,m))

# plot
X = T[:,0]
ct = 1
for theta in X:
	x1 = 0;y1 = 0
	x2 = l*math.sin(theta)
	y2 = -(l*math.cos(theta))
	file = \'test%05d.png\'%ct
	plt.figure()
	plt.plot([-1.5,1.5],[0,0],linewidth = 5)
	plt.plot([x1,x2],[y1,y2],linewidth = 4)
	plt.plot(x2,y2,\'o\',markersize =12)
	plt.xlim([-1.5,1.5])
	plt.ylim([-1.5,1.5])
	plt.savefig(file)
	ct = ct+1

Program for Plot :

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


# function 
def model(T,t,b,l,g,m):
	T1= T[0]
	T2 = T[1]
	dT1_dt = T2
	dT2_dt = -(b/m)*T2 -(g/l)*math.sin(T1)
	dTh_dt = [dT1_dt, dT2_dt]
	return dTh_dt

# inputs
b = 0.05 # damping co-efficient

l = 1 # length of pendulum in m

g = 9.81 # gravity in m/s^2

m = 1 # mass of ball in kg

# initial conditions
T_x = [0,3]

# time points
t = np.linspace(0,20,150)
# solve ode
T = odeint(model,T_x,t,args=(b,l,g,m))

# plot
plt.plot(t,T[:,0],\'-\',label = \'displacement\')
plt.plot(t,T[:,1],\'--\',label= \'Velocity\')
plt.ylabel(\'plot\')
plt.xlabel(\'time\')
plt.legend()
plt.show()

Results:

Plot:

Animation:

https://youtu.be/eWLy7wMQw14