Week 3 - Solving second order ODEs

SOLVING SECOND ORDER ODE USING PYTHON

AIM – Write python program to solve ODE representing the equation of motion of a simple pendulum with damping.

OBJECTIVES-

1. Create logic to break second order differential equation into a system of first order ODEs.
2. Write a program to plot displacement and angular velocity of pendulum.
3. Create an animation to simulate pendulum for first 20seconds.

INTRODUCTION-

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 it back and forth. The mathematics of pendulums is in general quite complicated. Simplifying assumptions can be made, which in the case of a simple pendulum allow the equations of motion to be solved analytically for small-angle oscillations.

In Engineering, ODE is used to describe the transient behaviour of a system. A simple example can be of a pendulum. The way the pendulum moves depends on the Newton’s second law of motion. When this law is used considering effects of damping force too, we get a second order Ordinary Differential Equation that describes the position of the "ball" with respect to time.

Similarly, there are hundreds of ODEs that describe key phenomena and if one has the ability to solve these equations then he/she will be able to simulate any of these systems.

Our objective is to 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

Given data -

L=1 metre,

m=1 kg,

b=0.05.

g=9.81 m/s2.

Initial conditions (t=0) - angular displacement=0, angular velocity=3 rad/sec

We have to simulate the motion between 0-20 seconds.

PROCEDURE-

1. Let’s substitute Θ = Θ1 and `(dΘ1)/dt=Θ2`

So, now our second order equation can be replaced by two 1^st order equations

`(dΘ1)/dt=Θ2` and `(dΘ2)/dt=(-b)/m Θ2-g/l sin⁡〖Θ1 〗`

2. Write equation for above logic in function model in our code.
3. Import numpy, scipy, math, matplotlib libraries. Import odeint from scipy.integrate.
4. Provide the input data and initial conditions. Use linspace to generate time points for calculation. Use 300 points to generate accurate solution.
5. Solve the ODEs using odeint. Plot the results to show angular displacement and angular velocity of pendulum.
6. To generate plot for animation generate string using line joining one fixed point and other which moves in accordance with our solution calculated.
7. For ball use marker. Set size to 12.

"""
Program to solve second order differential equation governing motion of simple pendulum
by Sushant Kumar Mishra
"""

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

def model(theta,t,b,l,g,m):
	theta1= theta[0]
	theta2= theta[1]
	dtheta1_dt=theta2
	dtheta2_dt=-(b/m)*theta2-(g/l)*sin(theta1)
	dtheta_dt= [dtheta1_dt,dtheta2_dt]
	return dtheta_dt

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

#intial condition 
theta_0= [0,3]

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

#solve ODE
theta= odeint(model,theta_0,t,args=(b,l,g,m))

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}\theta2-\frac{g}{l}sin\theta_1$')
plt.title('Displacement and angular velocity')
plt.xlabel('time')
plt.grid()
plt.legend(loc='lower right')
plt.show()


#animation
a=theta[:,0]
ct=1
for i in a:
	x0=0
	y0=0
	x1=-l*math.sin(i)
	y1=-l*math.cos(i)
	filename='test%05d.png'%ct
	ct=ct+1

	plt.figure()
	plt.plot([x0,x1],[y0,y1],linewidth=2)
	plt.plot([-1.25,1.25],[0,0],linewidth=5)
	plt.plot(x1,y1,'o',markersize=12)
	plt.xlim([-1.75,1.75])
	plt.ylim([-1.75,1.75])
	plt.savefig(filename)

RESULTS-

Plot of angular velocity and displacement:

ANIMATION:

CONCLUSION AND EXPLANATIONS-

1. Also, we must limit x axis and y axis to same value to create plots with same scale. Else, it might appear the length of string in changing in animation.
2. Give proper file name to each plot else their order will get mixed up while creating animation.
3. Damping is what reduces the amplitude of displacement and velocity bringing pendulum to rest.
4. Higher the damping constant lesser is time required for oscillations to die.
5. Thus, we obtain animation of pendulum by stitching images.