Solving second order ODEs using python

AIM

To solve the differential equation representing the oscillatory motion of a damped pendulum

GOVERNING EQUATION

`(d^2theta)/dt^2 +b/m*(d theta)/dt+g/l*sintheta=0`

`theta`= angular displacement of the pendulum

b = damping coefficient

m = mass of pendulum ball (kg)

g = acceleration due to gravity (m/`s^2`)

l = length of string

OBJECTIVE

The objective of this program is to solve the given second order differential equation and use the solution to simulate the motion of the pendulum governed by this differential equation. Apart from this, we would also plot the variation of displacement and velocity with respect to time.

Using odeint function from scipy module in python, a `n^(th)` order differential equation is solved by breaking it into n linear differential equation. For example, a `2^(nd)` order ifferential equation is solved by forming 2 linear differential equations and solving them to get the final result.

Considering the given equation

`(d^2theta)/dt^2 +b/m*(d theta)/dt+g/l*sintheta=0`

let `(d theta)/dt = (d theta_1)/dt = theta_2`

`implies``(d^2theta)/dt = (d^2theta_1)/dt = d/dt((d theta_1)/dt) = (d theta_2)/dt`

`implies (d theta_2)/dt + b/m*theta_2 + g/l*sintheta_1 = 0`

`or, (d theta_2)/dt=-b/m*theta_2-g/l*sintheta_1`

PROGRAM

i) The necessary modules are imported. It includes math, matplotlib, numpy and from scipy module, odeint function is imported.

ii) A function is defined to break the second order ODE into 2 first order ODE and store the values in an array which is returned.

iii) The required inputs and the initial condition of the pendulum is defined along with the time span.

iv) odeint function is used to solve the differential equation.

v) Different positions of the pendulum are saved saved as different figures which are then stitched together to create an animation.

odeint syntax

odeint(model, theta0, t, args=(b,g,l,m)) Here,

model is the function that computes the derivative of theta at time t.

theta0 is the initial condition

t is the sequency of time points for which theta is solved. This sequence is either monotonically increasing or decreasing.

args refer to the extra arguements that are needed to be passed to the function in order to solve the equation.

# program to simulate damped oscillatng pendulum

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


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

# Defining inputs
b = 0.05
g = 9.81
l = 1
m = 1

# Initial conditions
theta0 = [0,3]
t = num.linspace(0,10,200)

theta = odeint(model,theta0,t,args = (b,g,l,m))

# plotting displacement and velocity of pendulum as function of time
plt.plot(t,theta[:,0],'b-',label=r'Displacement')
plt.plot(t,theta[:,1],'r-',label=r'Velocity')
plt.xlabel('time')
plt.ylabel('plot')
plt.legend(loc='best')
plt.grid('on')
plt.show()

# Animation
ct = 1
x0 = 0
y0 = 0
for i in range(len(t)):
	x1 = l*math.sin(theta[i,0])
	y1 = -l*math.cos(theta[i,0])
	filename = 'pendulum%04d.png' %ct

	plt.figure()
	plt.plot([-1,1],[0,0],'k-')
	plt.plot([x0,x1],[y0,y1])
	plt.plot(x1,y1,marker='o',markersize='10')
	plt.axis([-1,1,-1.2,0.1])
	plt.savefig(filename)
	ct = ct+1

OUTPUT

Animation

https://drive.google.com/file/d/19EkbXN4SK5LIPuc7BFlFJgaAo8eZGzkJ/view?usp=sharing

CONCLUSION

As per the initial conditions provided, the graph for the displacement starts from 0 whereas that for the velocity starts from 3m/s.

The velocity of pendulum is maximum at the mean position and 0 at the amplitude.

With every oscillation, the amplitude of oscillation decreases. For higher damping coefficient, the amplitude would decrease at a faster rate and the pendulum would come to rest early.