Solving second order ODEs Using 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

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()

Explanation of code:

• Initially we are importing required modules using 'import' command.
• Function model is created with required argumeents, inside which the term which are broken to solve second order ODE are calculated.
• Inputs are provided for the function.
• theta_0 ranges between 0 to 3 m/s^2.
• odeint command is used to solve ode function.
• Then results are plotted using plotting command.
• For creating animation we open a for loop, ranging from 0 to length of t.
• Images are saved using savefig command.

Output:

Animation link:

https://youtu.be/TjBXo0hP0Ck

since image gif format was not supported by youtube,so it is recorded and uploaded.

Errors :

Basic syntax errors are made during first attempt, those are sooner rectified and corrected