Solving second order ODEs - Motion of pendulum
Solving second order ODEs
TRANSIENT BEHAVIOUR OF A SIMPLE PENDULUM
To solve the second order ODEs that show the transient behaviour of a simple pendulum using python and also to create the animation of its motion.
ODE is used to describe the transient behaviour of a system. Simple pendulum deals with the ODE concepts. The ODE represents the equation of motion of a simple pendulum with damping.
Consider the 2nd order differential equation i.e.
This 2nd order ODE can be used for a simple pendulum and also for vehicle dynamics.
Main reason for solving ODEs is when we have higher order ODEs we break them into a set of 1st order ODEs. (because we are considering 2nd order differential equation)
Program :-
import math
import matplotlib.pyplot as plt
import numpy as np
from scipy.integrate import odeint
# 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 # two ODE\'S to single ODE
# inputs
b = 0.05 # damping coefficient
g = 9.8 # gravity in \'m/s2\'
l = 1 # length of pendulum in \'m\'
m = 1 # mass of ball in \'kg\'
# initial condition
theta_0 = [0,3] # angular displacemet = 0 & angular velocity = 3 rad/sec at t=0
# time points
t = np.linspace (0,20,150) # motion between (0-20)sec
# solving ODE
theta = odeint(model, theta_0, t, args = (b,g,l,m))
# plotting results
plt.plot(t,theta[:,0], \'b-\', label = \'Displacement\')
plt.plot(t,theta[:,1], \'r--\', label = \'Velocity\')
plt.xlabel(\'Time in sec\')
plt.ylabel(\'plot\')
plt.legend(loc = \'best\')
plt.grid()
plt.show()
Explanation :
- Initially before starting main code , \"numpy module \" is imported.
- Also import \" scipy module\".
- To plot the result, we require \" matplot module \", so import \" matplotlib.pyplot as plt \".
- Define \" def model(theta,t,b,g,l,m): \"
- Then define constants that is the given inputs as per the question. i.e. b= 0.05, g=9.8, l=1, m=1.
- Based on the condition given in question, define initial condition i.e. angular displacement and angular velocity.
- As we have 2-ODEs, 1st value is called displacement value and 2nd value is calle d velocity value. i.e according to question [0,3] here 0- angular displacement & 3- angular velocity.
- For linspace command we require numpy module.
- To solve the ODE equation we should use \' odeint command \'.
- odeint command works only after importing scipy module.
- The odeint is bascially a function that integrate ODE.
- To integrate we bascially need to define the model/system, initial condition, time and arugements are entered.
- Finally we get Theta = solution vector.
Errors faced while programming :
1.
Errors made -
- Incomplete program without defining the xlabel for the graph obtained.
- I found both curves on the graph have same lines with different curves.
Solutions -
- plt.xlabel(\' Time in sec \') is added just before 36th line.
- plt.plot(t,theta[:,0], \'b-\', label = \'Displacement\') ; plt.plot(t,theta[:,1], \'r--\', label = \'Velocity\').
2.
Errors made -
- In line 31, key word argument \'arg\' is found as error.
Solution -
- AS we have multiple arguments, we need to give the command as \'args\'.
3.
Errors made -
- In line 39, theta[:,0] is defined as angular _ displacement and in line 41 i.e starting of for loop statment it is denoted by ang_displacement.
Solution -
- Change \"angular_displacement\" to \"ang_displacement\".
Output :-
The (displacement and velocity) vs (time) plot is generated and it shows that displacement starts at point-0 and velocity starts at point-3. Hence the plot obtained is correct. It is also observed that magnitude of both velocity and displacement are droping with increasing time.
2. Write a program in Python that will simulate the pendulum motion.
Given input :
L=1 m
m=1 kg
b=0.05
g=9.81 m/s2.
Simulating the motion between 0-20 sec, for angular displacement=0,angular velocity=3 rad/sec at time t=0.
Program :-
import math
import matplotlib.pyplot as plt
import numpy as np
from scipy.integrate import odeint
# function
def oscillation(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 # two ODE\'S to single ODE
# inputs
b = 0.05 # damping coefficient
g = 9.8 # gravity in \'m/s2\'
l = 1 # length of pendulum in \'m\'
m = 1 # mass of ball in \'kg\'
# initial condition
theta_0 = [0,3] # angular displacemet = 0 & angular velocity = 3 rad/sec at t=0
# time points
t = np.linspace (0,20,150) # motion between (0-20)sec
# solving ODE
theta = odeint(oscillation, theta_0, t, args = (b,g,l,m))
# co-ordinates
x0 = 0
y0 = 0
x = []
y = []
ang_displacement = theta[:,0]
for i in ang_displacement:
x1 = l*math.sin(i)
y1 = -l*math.cos(i)
x.append(x1)
y.append(y1)
image_num = 1
for j in range(0,150):
# results
image_name = \'image%05d.png\' % image_num
plt.figure()
plt.plot([x0,x[j]], [y0,y[j]], label = \'string\')
plt.plot([-1,1],[0,0], linewidth = 2, label = \'wall\')
plt.plot(x[j],y[j],\'.\', markersize = 30, label = \'bob\')
plt.xlim(-2,2)
plt.ylim(-2,2)
plt.xlabel(\'x co-ordinate\')
plt.ylabel(\'y co-ordinate\')
plt.title(\'simple pendulum\')
plt.savefig(image_name)
image_num = image_num+1
Importanat Considerations :-
- Defining a function that is necessary for coding.
- List or define the inputs before using them .
- Mention time points i.e. time interval that is mentioned in question.
- Declair initial co-ordiantes as[0,0] and define an empty array for x & y respectively.
- Mentioning the condition for angular displacement is mandatory.
- When plotting the pendulum it is necessary to declair the limits of both the axis.
Output -
Result : The motion of a Simple Pendulum is observed on excuting the program.
Youtube link -
Projects link -
https://projects.skill-lync.com/projects/Solving-second-order-ODEs-Simple-Pendulum-10032