Solving second order ODEs Using PYTHON

Objective:

Write a program that solves the following ODE. This ODE represents the equation of motion of a simple pendulum with damping. And also write a program to simulate the transient behavior of a simple pendulum and to create an animation of its motion.

Theory:

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.

Similarly, there are hundreds of ODEs that describe key phenomena and if you have the ability to solve these equations then you will be able to simulate any of these systems.

Solution: -

In the above equation,

• g = gravity in m/s2,
• L = length of the pendulum in m,
• m = mass of the ball in kg,
• b = damping coefficient.

Given,

• L=1 meter,
• m=1 kg,
• b=0.05.
• g=9.81 m/s2.

Simulate the motion between 0-20 sec

• for angular displacement=0
• angular velocity=3 rad/sec
• at time t=0.

Coding: -

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


#funcion that returns dtheta/dt

def function(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.05
g = 9.81
l = 1
m = 1

#initial conditions
theta_0 = [0,3]


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

#solve ode

theta = odeint(function,theta_0,t,args = (b,g,l,m))

#plot result
plt.plot(t,theta[:,0],'g-',label=r'$frac{dtheta_1}{dt}=theta2$')
plt.plot(t,theta[:,1],'b--',label=r'$frac{dtheta2}{dt}=-frac{b}{m}theta_2-frac{g}{l}sintheta_1$')
plt.xlabel('Plot')
plt.ylabel('Time')
plt.legend()
plt.savefig('Output.png')
plt.show()

#animation
ct = 1

x0 = 0
y0 = 0

for i in theta[:,0]:
	
	x1 = (1*math.sin(i))
	y1 = (-1*math.cos(i))

	filename = 'test%05d.png' % ct
	ct = ct +1
	
	plt.figure()
	plt.plot([-1,1],[0,0],"black",linewidth=10)
	plt.plot([x0,x1],[y0,y1],"red")
	plt.plot(x1,y1,'o',markersize = 20)
	plt.xlim([-1.5,1.5])
	plt.ylim([-1.5,1.0])
	plt.savefig(filename)

Steps used for coding:

• Import module numpy as np, scipy.intedrate import odeint for second ODE, math and for plots matplotlib.
• Define function that returns the value of dtheta/dt.
• Initialize the input values of b, g, m, l.
• Define initial condition form theta i.e. 0 and 3 rad/sec.
• Use np.linspace for creating numeric sequence; np.linspace(start point, last point, number of values take in between start and end point).
• Solve ODE by use of odeint command.
• Plot result by using plt.plot command to plot graph.
• For animation: - use for loop
• Name the filename
• Use plt.plot command to create every single point and after that stiches all the pic and create a gif

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.

Results / Output:

Angular displacement, angular velocity with respect to time.

Animated Vedio Youtube Link