Week 3 - Solving second order ODEs

AIM:

A program that solves the ODE which represents the equation of motion of a simple pendulum with damping.

THEORY:

The below 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.

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.

Simulating the motion between 0-20 sec, for angular displacement=0,angular velocity=3 rad/sec at time t=0.

CODE EXPLANATION:

1.Declaring the input values i.e., b,m,g,l

2. Initially providing t_span (integration limits) and initial values of displacement and velocity

3. Building the function of ode

4. ODE solving using python:

(i)Here there are 2 odes so, first, we provide 2 initial conditions to solve them i.e., displacement and velocity

(ii)Then including time array using numpy module

(iii)For solving ODE odeint command is used which is in scipy module

(iv)This odeint integrates a system of an ode by taking initial condition to a particular time step and after getting some solution it keeps on integrating this time step depends on type of ODE system is used

(v)Then we do plotting

5. In figure(1):plotting angular displacament and angular velocity with respect to time.

6. In figure(2):plotting pendulum with marker 'o' using robort arm technique.

7. Next animation plot using same technique that done in robort arm.

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

# Defining constants(Damping coeff,length of pendulum,gravity and mass)
b=0.05
g=9.81
l=1
m=1
# intial condtn
theta_0=[0,3]
# time points
t_span=np.linspace(0,20,100)
# Fuction defination for solving ODE
def pendulum(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
# calling function
theta=odeint(pendulum,theta_0,t_span,args=(b,g,l,m))
# plotting
plt.figure()
plt.plot(t_span,theta[:,0])
plt.plot(t_span,theta[:,1])
plt.xlabel('Time')
plt.ylabel('angular displacement and angular velocity')
plt.legend(['angular displacement','angular velocity'])
plt.show()
angular_displacements=theta[:,0]
ct=1
for a_d in angular_displacements:
	x0=0
	y0=0
	x1=l*math.sin(a_d)
	y1=-l*math.cos(a_d)
	plt.figure()
	plt.plot([-1,1],[0,0])
	plt.plot([x0,x1],[y0,y1])
	plt.plot(x1,y1,'o',markersize=25)
	plt.xlim([-2,2])
	plt.ylim([-2,2])
	filename=str(ct)+'_.png'
	plt.savefig(filename)
	ct=ct+1

OBSERVATIONS:
PLOT Angular Displacement vs Angular Velocity:

Angular Velocity and Angular Displacement drops down slowly due to the damping coefficient.