Transient behaviour of a simple pendulum using 2nd order ODE in python
Objective: Write a program to simulate the transient behaviour of a simple pendulum and to create an animation of its motion.
Brief: In Engineering, ODE is used to describe the transient behaviour 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.
Similarly, there are hundreds of ODEs that describe key phenomena and if you have the ability to solve these equations then you will be able to simulate any of these systems.
Your objective is to write a program that solves the following ODE. This ODE represents the equation of motion of a simple pendulum with damping.
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.
Write a program in Python that will simulate the pendulum motion, just like the one shown in the start of this challenge.
use,
L=1 metre,
m=1 kg,
b=0.05.
g=9.81 m/s2.
Simulate the motion between 0-20 sec, for angular displacement=0, angular velocity=3 rad/sec at time t=0.
When you do this, you will get the position of the pendulum as a function of time. Use this information to animate the motion of the pendulum, as shown in the below figure.
Your animation does not have to look exactly like this and here is the logic.
- For each value in time
- Create an image that represents the position of the pendulum
- You can assume the length of the string and the size of the ball
- Use the marker command to create a sphere
- One end of the pendulum is fixed
- The other end moves according to the solution of the ODE you had solved.
Theoretical explanation:
Here we are going to understand the basic understanding of 2nd order differential equation of simple pendulum analytically and how it is derived in an equation. Here\'s how the position and velocity of pendulum ball can be determined with respect to time with the below equation.
where,
b is damping coefficient
m is a mass in kg
g is acceleration due to gravity
l is length of the pendulum
this is how we proceed here as
by substituting the above eq in eq.1
By getting the value of theta1 and theta2 from here,
we have here matrix formed,
Step by step explanation of coding: Here we understand how this coding works with python etiquettes.
# 2nd order ODE of a simple pendulum
# using differential equation
import math
import matplotlib.pyplot as plt
import numpy as np
from scipy.integrate import odeint
before we begin the code there are some math fucntions as odeint and numpy from scipy.integrate as they help us to do mathematical operation on an element of an array and it also stores value in a range.
#function returning the dz/dt value
def model(theta,t,b,g,l,m):
here we have defined the function called as def model(theta,t,b,g,l,m) which stores the dtheta_dt values which shows the transient behaviour of pendulum.
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
here again the code in python returns the values of dtheta_dt, hnece mentioned return dtheta_dt.
#input
# damping coefficient
b = 0.02
g = 9.81
l = 1
m = 0.1
here we have mentioned the given values.
#initial condition
theta_0 = [math.pi/2,0]
theta_0 has values
#time points
t = np.linspace(0,10,150)
print (t)
the function we used at the beginning of the code it deploys this math opeartion where t is ranging in an array.
# ODE function definition
theta = odeint(model,theta_0,t,args = (b,g,l,m))
here we are calling the fucntion we made earlier as function definition.
#plot results
plt.plot(t,theta[:,0], \'b-\',label = r\'$\\frac{d\\theta_1}{dt}=\\theta$\')
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()
here we have plotted the the whole values to get desirable outcome.
#for loop for animation
the below code is to perfom animation of the above coding and exhibit the transient behaviour of a pendulum.
ct = 1
for i in (theta[:,0]):
x0 = 0
y0 = 0
x1 = l*math.sin(i)
y1 = -l*math.cos(i)
filename = str(ct) + \'.png\'
ct = ct +1
print(filename)
plt.figure()
plt.plot([-1,1],[0,0],\'k\',linewidth=4)
plt.plot([x0,x1],[y0,y1])
plt.plot(x1,y1,\'o\',markersize = 20)
plt.xlim([-1.5,1.5])
plt.ylim([-1.5,1.5])
plt.savefig(filename)
again we have made 150 .png files with this from where we can generate our animation by combining the whole pngs to make gif file where it shows th animation.
Original coding:
# 2nd order ODE of a simple pendulum
# using differential equation
import math
import matplotlib.pyplot as plt
import numpy as np
from scipy.integrate import odeint
#function returning the dz/dt value
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
# damping coefficient
b = 0.02
g = 9.81
l = 1
m = 0.1
#initial condition
theta_0 = [math.pi/2,0]
#time points
t = np.linspace(0,10,150)
print (t)
# ODE function definition
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}=\\theta$\')
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()
#for loop for animation
ct = 1
for i in (theta[:,0]):
x0 = 0
y0 = 0
x1 = l*math.sin(i)
y1 = -l*math.cos(i)
filename = str(ct) + \'.png\'
ct = ct +1
print(filename)
plt.figure()
plt.plot([-1,1],[0,0],\'k\',linewidth=4)
plt.plot([x0,x1],[y0,y1])
plt.plot(x1,y1,\'o\',markersize = 20)
plt.xlim([-1.5,1.5])
plt.ylim([-1.5,1.5])
plt.savefig(filename)And here again we have function program of ODE as def model(theta,t,b,g,l,m)
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_dtAnimation of output:
Output:
The angular displacement & velocity vs time has been plotted and we can see that in the below graph.
The red dotted line wave is angular velocity as it start from 0 and the green wave is angular displacement as it starts from pi/2.