Week 3 - Solving second order ODEs using PYTHON
OBJECTIVE
- To write a program in python to simulate the transient behavior of a Simple pendulum.
- To plot graph Time vs Displacement.
- To create an animation of its motion.
THEORY ABOUT SIMPLE PENDULUM
The pendulum is a simple mechanical system that follows a differential equation. The pendulum is initially at rest in a vertical position. When the pendulum is displaced by an angle θ and released, the force of gravity pulls it back towards its resting position.
A pendulum is a weight suspended from a pivot so that it can swing freely. When a pendulum is displaced sideways from its resting, equilibrium position, it is subject to a restoring force due to gravity that will accelerate it back toward the equilibrium position. When released, the restoring force acting on the pendulum's mass causes it to oscillate about the equilibrium position, swinging back and forth. The time for one complete cycle, a left swing and a right swing, is called the period. The period depends on the length of the pendulum and also to a slight degree on the amplitude, the width of the pendulum's swing.
GOVERNING EQUATION
Second order differential equation for simple pendulum is
`(d^2 θ)/(dt^2 )+ b/m dθ/dt+ g/lsinθ=0`
Where,
g=gravity (`m/s^2`)
l=Length of the pendulum (m)
m=mass of the ball (kg)
b=damping coefficient
INPUT VALUES
g=9.81 `m/s^2`
l=1 m
b=0.05 m
m=1 kg
PROCEDURE
- Initially we are going to import the “matplotlib.pyplot” , “math” , “scipy” and “numpy” library for the calculations.
- Scipy library provides us a equation of the second order ODE to solve the theta value.
- Then just input the given parameters are g, b, l and m to calculate the other parameters.
- Using the calling function command we are going to call theta values.
- And then using “numpy” command we provide the steps in the time.
- Then plot the Time against Displacement graph.
- For the motion of simple pendulum, first mention the co-ordinates of the pendulum at the start(x0, y0) and end(x1, y1).
- Using the string and counter command we’ll give different png names and save all the png files.
- For the Animation we use some other software to create the animation.
PROGRAM
#second order ODE
#by Sallavoddin
import matplotlib.pyplot as plt
import math
from scipy.integrate import odeint
import numpy as np
# calling 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
g=9.81
l=1
m=1
#intial condition
theta_0=[0,3] #angular displacement abd angular velcity
#time point
t=np.linspace(0,20,250) #0 to 20 sec motion
#Solve ODE
theta=odeint(model,theta_0,t,args=(b,g,l,m))
#point results
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}{l}sin\theta_1$')
plt.xlabel('Displacement')
plt.ylabel('Time')
plt.legend(loc='best')
plt.show()
#Animation of Simple pendulum
x0=0 #this is fixed point at the wall
y0=0
ct=1
for ang in theta[:,0]: #angle
x1=l*math.sin(ang)
y1=-l*math.cos(ang)
filename=str(ct)+'.png'
ct=ct+1
plt.figure()
plt.plot([-1,1],[0,0],'red',linewidth=10)
plt.plot([x0,x1],[y0,y1],'black')
plt.plot(x1,y1,'o',markersize=20)
plt.xlim(-1.5,1.5)
plt.ylim(-1.5,1.5)
plt.savefig(filename)
RESULTS
Plot
Animation Video
Simulate the motion between 0-20 sec, for angular displacement=0,angular velocity=3 rad/sec at time t=0.
CONCLUSION
- Simulation of transient behavior of a Simple pendulum motion by solving second order differential equation.
- Animation of the Simple pendulum is simulated.
REFFERENCES