Week 3 - Solving second order ODEs

AIM:

To write a program in Python to simulate the pendulum motion and to solve second order ODE.

OBJECTIVE:

1. To solve the second order ODE in Python and plot amplitude vs time and velocity vs time graph.
2. To simulate the motion of simple pendulum and show the animation

GIVEN DATA:

1. The value of acceleration due to gravity is (g) 9.81 m/s^2
2. The length of the pendulum is 1 meter.
3. The damping coefficient is 0.05.
4. The mass of the pendulum bob is 1 kg.

THEORY:

The motion of the pendulum can be easilty represented by the following differential equation:

`(d^2 θ)/dt^2 +b/m dθ/dt + g/L sin⁡θ = 0`

where:

`theta`= angular displacement of the bob

b = damping coefficient

g = accelration due to gravity

m = mass of the pendulum bob

L = length of the pendulum string

t = time period

`(dθ)/dt`=angular velocity

`(d^2θ)/(dt^2)` = angular acceleration

python program:

"""
Second order ODE

Raaghul V
"""

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

#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

#input parameters

b=0.05#damping coefficient
g=9.81#acceleration due to gravity in m/s^2
l=1#length of the pendulum in m
m=1#mass of the ball in kg

theta_0=[0,3]

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

theta=odeint(model,theta_0,t,args=(b,g,l,m))#inbuilt ode solver function

plt.figure()#opens a new empty figure
plt.plot(t,theta[:,0],'b-',label=r'$\frac{d\theta_1}{dt}=\theta_2$')

#
plt.plot(t,theta[:,1],'r--',label=r'$\frac{d\theta_2}{dt}=\frac{-b}{m}\theta_2-\frac{g}{l}sin\theta_1$')
plt.ylabel('Plot')#labelling X-axis
plt.xlabel('Time')#labelling Y-axis
plt.legend(loc='best')
plt.savefig('ode result')#saves plot as given name
plt.show()#display the plot result


#to save pendulum's motion images
val=theta[:,0]
count=1

for i in range(0,len(val)):
    plt.figure()
    x1=0.5-l*np.sin(val[i])
    y1=0.5-l*np.cos(val[i])
    plt.plot([0.5,x1],[0.5,y1])
    plt.plot(x1,y1,'o',markersize=15)
    plt.plot([-1,2],[0.5,0.5],linewidth=5)
    plt.xlim([-1,2])
    plt.ylim([-1,2])
    filename=str(count)+'.png'
    plt.savefig(filename)
    count=count+1

Explanation:

1.Set the input parameters of the pendulum like mass, length, gravity and damping coefficient.
2.Set the initial conditions and timespan for which you want to run the simulation.
3.To use the ode solver, use odeint sovler from scipy module and use model function with inputs as (t,theta,l,b,g,m) and t and theta_0 as the initial conditions.
4.These results are stored in an array of dtheta_dt which gets retuned to theta
5.Since results contain 2 components, angle displacement and velocity, to plot both the graphs, use the [:,0] and (:,1) property of Python and plot those graphs.
6.To create animation by using the ct counter to make up for the file name using str(ct) and a '.png' file extension.
7.Use gif to stich all the images to create an animation and convert it to video.

ERRORS:

No errors were generated while executing.

Output:

Graph:

Pendulum animation video:

https://youtu.be/ThT9koIsZNE