Week 3 - Solving second order ODEs
Objective:
To write a program that solves a second order ODE and represents the equation of motion of a simple pendulum with damping.
Introduction:
In Engineering, ODE is used to describe the transient behavior of a system. A simple example is a pendulum.
The way the pendulum moves depends on the Newtons second law. When this law is written down, we get a second order Ordinary Differential Equation that describes the position of the "ball" w.r.t time.
Formula :
`(d^2 theta)/dt^2 + b/m (d theta)/dt +g/l sin theta = 0`
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.
Code:
import numpy as np
from scipy.integrate import odeint
import matplotlib.pyplot as plt
import math
#defining the 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
#inputs
b=0.05 #damp coefficient
g=9.81 #gravity in m/s^2
l=1 #length of the pendulum in meter
m=1 #mass of the ball in kg
#initial condition
theta_0 = [0,3]
#time points
t = np.linspace(0,20,150)
#solving ode
theta = odeint(model,theta_0,t,args = (b,g,l,m))
#animation plots
a=[]
f=1
for i in theta[:,0]:
a.append(i)
x0=0
y0=0
x1=(l*math.sin(i))
y1=-(l*math.cos(i))
filename = str(f)+'.png'
f=f+1
plt.figure()
#horizontal support for the pendulum
plt.plot([-0.2,0.2],[0,0])
#length of the string
plt.plot([x0,x1],[y0,y1])
#the bob
plt.plot(x1,y1,'-o')
plt.xlim([-1.5,1.5])
plt.ylim([-1.5,1.5])
plt.savefig(filename)
#plot results
plt.plot(t,theta[:,0],'b-')
plt.plot(t,theta[:,1],'r--')
plt.ylabel('plot')
plt.xlabel('time')
plt.show()
Code explanation:
1. The scipy module is installed to solve any ODE.
2. The second order ode is divided into 2 first order ode .The first one is with the displacement and the second one for the velocity.
3. Initial conditions are declared and defined as theta0=(0,3) and for the pendulum motion the time spacings are given as 20 sec using linspace command.
4. ODEINT is used to solve ode and the values are saved as theta.
Errors faced:
Correction done:
The syntax error was corrected from brackets to braces.
Results:
Plot:
Animation: