Week 3 - Solving second order ODEs

Aim: To solves the ODE which represents the equation of motion of a simple pendulum with damping.

Objective:

 To write a program that solves the following ODE which represents the equation of motion of a simple pendulum with damping and create an animated video of output obtains by solving this ODE.

Theory:

The motion of a simple pendulum is a basic classical example of simple harmonic motion, consisting of a small bob and a massless string. Movement of the pendulum depends on Newton's second law. When this law is written down, a second-order Ordinary Differential Equation is formed that describes the position of the end of the pendulum w.r.t time.

Equation of motion is:

In the above equation,

• g = gravity in m/s^2.
• L = length of the pendulum in the meter.
• m = mass of the ball in kg.
• b = damping coefficient.

Since the given function is in the second-order differential equation we will convert it into first order to solve the ODE.

     Let
                θ=θ1

                dθ1/dt=θ2

                dθ^2/dt^2=dθ2/dt

Input :

• L=1 metre
• m=1 kg
• b=0.05
• g=9.81 m/s^2
• Time span= 0-20 sec

Code :

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

#Funtion that returns dz/dt
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=0.1

#initial condition
theta_0 = [0,5]

#time points
t = np.linspace(0,10,150)

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

#plot results
plt.plot(t,theta[:,0],'b-',label=r'$\frac{d\theta_1}{dt}=\theta2$')
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.xlabel('plot')
plt.ylabel('time')
plt.legend(loc='best')
plt.show()

Detail explanation each step of the program:

1.First step to inputs the python liberary.

2.Write the funtion dz/dt.

3. Inputs the constant values (b,g,l,m).

4. Putting initial condition for theta_0.

5. Using the odient funtion to solve this equation.

6. Gives the ploting command.

Output:

1. The output obtained is graphical representation velocity curve and displacement curve w.r.t Time.

•   Angular Displacement and Angular Velocity plots were generated for a time range of 20 secs.