Week 3 - Solving second order ODEs

AIM:

To write a Matlab program to simulate the motion of a simple pendulum with damping.

INPUTS:

g = 9.81 m/s2, gravity

L = 1 metre, length of the pendulum

m = 1 kg, mass of the ball

B = 0.05, damping coefficient

OBJECTIVE:

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

ODE represents the equation of motion of a simple pendulum with damping

PROGRAM:

import numpy as np

from scipy.integrate import odeint

import math

import matplotlib.pyplot as plt

#function 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

#Given data

b=0.5

g=9.81

l=1

m=1

# initial conditions

theta_0 = [0,3]

# input for time

t = np.linspace(0,20,200)

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}{1}sin\theta_1$')

plt.xlabel('Plot')

plt.ylabel('Time')

plt.legend()

plt.savefig('output.png')

plt.show()

# To create Animation

ct = 1

x0 = 0

y0 = 0

for i in theta[:,0]:

x1=x0+(1*math.sin(i))

y1=y0+(-1*math.cos(i))

filename = str(ct)+'.png'

ct = ct+1

plt.figure()

plt.plot([-1,1],[0,0],"black")

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)

plt.plot()

Results:

Plot

Animation: