Week 3 - Solving second order ODEs
AIM:
- Write a program that solves the following ODE. This ODE represents the equation of motion of a simple pendulum with damping.
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.
2. Write a program in Python that will simulate the pendulum motion, just like the one shown in the start of this challenge.
use,
L=1 metre,
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.
THEORY:
A pendulum is a body suspended from a fixed support so that it swings freely back and forth under the influence of gravity. When a pendulum is displaced sideways from its resting, equilibrium position, it is subject to a restoring force due to gravity that will accelerate it back toward the equilibrium position. When released, the restoring force acting on the pendulum's mass causes it to oscillate about the equilibrium position, swinging back and forth.
Damping suggests that the angular amplitude decreases exponentially over time. This means that the amplitude of the pendulum will go on decreasing with time and eventually the pendulum will stop.
In Engineering, ODE is used to describe the transient behaviour of a system. A simple example is a pendulum. The way in which 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.
In this project we will simulate the motion of the simple pendulum with damping by solving the second order ODE which represents the equation of motion of a simple pendulum with damping.
In this project we solve the second order ODE by considering the following:
`theta=theta_1`
`(d theta)/dx=(d theta_1)/dx=theta_2`
`(d^2 theta)/dt^2=(d^2 theta_1)/dt^2=(d theta_2)/dt`
`(d^2 theta_1)/dt^2+b/m((d theta_1)/dt)+g/l sin theta_1=0`
`(d theta_2)/dt +b/m(theta_2)+g/l sin theta_1=0`
`(d theta_2)/dt = -b/m(theta_2)-g/l sin theta_1`
`d/dx[[theta_1],[theta_2]]=[[theta_2],[-b/m*theta_2-g/l sin theta_1]]`
This is a system of ODEs with two first order ODE in it. We use this to get our solution.
PROGRAM:
Programming language used is Python.
Code to simulate the transient behaviour of a simple pendulum and to create an animation of it's motion using second order Ordinary Differential Equation in Python.
# Program for solving second order ODEs and simulate its motion.
import math
import matplotlib.pyplot as plt
import numpy as np
from scipy.integrate import odeint
# Defining ODE function
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
b=0.05
g=9.81
l=1
m=1
theta_0=[0,3]
t = np.linspace(0,20,150)
# Solve ODE
theta = odeint(model,theta_0,t, args = (b,g,l,m))
# Plotting the graph
plt.plot(t,theta[:,0],'b-',label=r'$\frac{d\theta_1}{dt}=\theta2$')
plt.plot(t,theta[:,1],'r--',label=r'$\frac{d\theta2}{dt}=-\frac{b}{m}\theta_2-\frac{g}{l}sin\theta_1$')
plt.ylabel('Plot')
plt.xlabel('Time(second)')
plt.title('Damped Harmonic motion of a Simple Pendulum')
plt.legend(loc='best')
plt.show()
# Iteration using for loop to create animation of motion of pendulum
p=theta[:,0]
ct=1
x0=0
y0=0
for i in p:
x1= x0 + (l* math.sin(i))
y1= y0 - (l* math.cos(i))
ct =ct+1
filename = '%05d.png' % ct
plt.figure()
# Plotting the rigid support through which string of pendulum is attached
plt.plot([-.5,.5],[0,0],'b-', linewidth = 10,color= 'brown')
plt.plot([x0,x1],[y0,y1],linewidth=4) # Plotting String of pendulum
plt.plot(x1,y1,'o',markersize = 25) # Plotting Bob of pendulum
plt.xlim([-2,2])
plt.ylim([-2,1])
plt.title('Simple Pendulum')
plt.savefig(filename)EXPLANATION OF THE PROGRAM:
- Program started by importing modules required. We need math module for mathematical calculations, matplotlib module for plotting results, numpy module to evaluate time intervals, scipy odeint module to solve ODEs (Ordinary Differential equations).
- In order to solve ODEs, a function is considered which will return the value of dΘ/dt.
def ode_function(theta,t,b,g,l,m):
Above function consists of theta, time, damping co-effcient, length of pendulum, mass of bob, acceleration due to gravity. As explained in calculation part, Θ1 is displacement function and Θ2 is angular velocity function.
- Input values are:
l = 1m, m = 1kg, b = 0.05, g = 9.81 m/s^2.
- Time interval are specified with help of numpy module. Time interval of 20 seconds with 150 step size.
- The function defined in called back in order to calculate theta.
theta = odeint(ode_function,theta_0,t_span,args = (b,g,l,m))
Above function made use of odeint function of scipy module.
- Result obtained from above function are plotted by using matplotlib module. Result obtained will have 2 columns and multiple rows. First columns are the values of displacement while that of in second column are the values of angular velocity.
Therefore, dispacement results are plotted as;
plt.plot(t, theta[:,0]) = it implies theta should consider first column and all rows. It will give response of displacement curve w.r.t time with red color.
plt.plot(t, theta[:,1]) = it implies thata should consider second row and all rows. It will give response of velocity curve w.r.t time with green dashed line.
Legend function used will display annotations of curve at "best" suited location.
- After plotting response of pendulum movement against time, next objective is to plot motion of simple pendulum.
-
After this position of the pendulum is plotted at particular time in the time span and captured frame by frame to make an animation.
- Then we can use command line tool “Cygwin” and “imagemagick” program tool to combine all these positions to form an animation.
- We can also go to google and convert images into gif. and video file.
Pendulum movement against time :-
Simulation of the pendulum motion.
GIF. file :-
Youtube link for video file of motion of pendulum :-