Week 3 - Solving second order ODEs

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{d\theta_1}{dt}=\theta2$')
plt.plot(t,theta[:,1],'b--',label=r'$\frac{d\theta2}{dt}=-\frac{b}{m}\theta_2-\frac{g}{l}sin\theta_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

Results: -

Output: - Angular displacement, angular velocity with respect to time.

Video link: - https://drive.google.com/file/d/1yobgSFVf096akfKe-LZ-rNuBCyVvrjYg/view?usp=sharing

Conclusion: - The above graph shows the angular displacement and angular velocity and it changes over the time due to the effect of the damping, the values of reduce angular displacement and angular velocity with time as is clearly visible in above plot.