Week 3 - Solving second order ODEs
Solving 2nd order differential equation using python programming
Aim:-
- Understand the Ordinary differential equation(ODE).
- Study pendulum motion and damping of mass.
- Break down higher order ODE into set of First order ODE.
- Write the function code for ODE
- Write code for main program
- Create report on overall program and ODE also on pendulum motion and damping.
Objective:-
- Study how the Ordinary differential equation are made for pendulum motion and spring mass damping problem.
- Study how dynamic analysis work.
- Brake down Pendulum motion 2nd order ODE into first order ODE.
- Use that ODE to write function codes to run the main program.
- Write main program with changing theta angle in range.
- Plot Displacement and angular velocity diagram.
- Create image array of pendulum motion.
- Convert pendulum motion frames into video file using Open cv library.
- Write the report explaining all codes also about ode and Pendulum mass movement and spring mass system.
An Ordinary Differential Equation (ODE) is a differential equation containing one or more function of one independent variable and the derivatives of those function. The term ordinary is used in contrast with the term partial differential equation which may be with respect to more than one independent variable.
Following is the 2nd order differential equation for pendulum mass damping system.
`(d^2theta)/(dt^2) + b/m (d theta)/(dt) + g/L sin theta = 0`
Where,
`b`= damping
`m` = mass
`g`= gravity
`L` = length
At the very beginning of the coding, import required library for this program
import numpy as np
from scipy.integrate import odeint
import matplotlib.pyplot as plt
import cv2
import mathA simple pendulum consist of a ball(point-mass) `m` hanging from a string of length L and fixed at a pivot point P. When displaces to an initial angle and released, the pendulum will swing bank and forth with periodic motion, By applying Newton’s second law for rotational system, the equation of motion for the pendulum may be obtained.
Newtons’s Second Law for rotation:-
If more than one torque act on a rigid body about a fix axis, then the sum of the torques equals the moment of inertia times the angular displacement.
`tau = I alpha` `rArr` `-mgsin theta L = mL^2 (d^2 theta)/dt^2`
`:. (d^2 theta)/dt^2 + g/L sin theta = 0`
If the amplitude of the angular displacement is small enough, so the small angle approximation holds true, then the equation of motion reduces to the equation of simple harmonic motion
`:. (d^2 theta)/(d theta^2) + g/L theta = 0`
Next, we add damping to the model. We make the simplest possible assumption about the damping force, that it is proportional to velocity. Since arc length and central angle are themselves proportional (with proportionality constant `L`), it makes no difference whether we use linear or angular velocity. Having selected `theta` as our dependent variable, we will represent the damping as proportional to angular velocity, say, `-b (d theta) / dt`. The negative sign is because the damping force has to be opposite the direction of motion. When we include this term in the model, our equation becomes.
`:. (d^2theta)/(dt^2) + b/m (d theta)/(dt) + g/L sin theta = 0`
Break down second order ODE into first order ODE
So first consider
`theta = theta_1`
And `(d theta)/(dt) = (d theta_1)/(dt) = theta_2` `…(1)`
`:. (d^2theta)/(dt^2) = (d^2theta_1)/(dt^2) = d/(dt) (theta_2)` `…(2)`
So `(d^2theta)/(dt^2) + b/m (d theta)/(dt) + g/L sin theta = 0`
Will be `(d^2theta_1)/(dt^2) + b/m (d theta_1)/(dt) + g/L sin theta_1 = 0`
`:. (d theta_2)/(dt) + b/m theta_2 + g/L sin theta_1 = 0`
`:. (d theta_2)/(dt) = - b/m theta_2 - g/L sin theta_1`
`:. d/(dt) [[theta_1] , [theta_2]] = [[theta_2],[- b/m theta_2 - g/L sin theta]]`
Now, write function program using above ordinary differential equation to solve the pendulum equation of motion problem.
# function 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
Next, import the input parameters for the problem which is…
`b`= 0.05
`m` = 1 kg
`g`= 9.81
`L` = 1m
Also, mention the angular velocity with time 0 and create the time range using linspace command.
b = 0.05
g = 9.81
l = 1
m = 1
# initial condition
theta_0 = [0,3]
# time points
t = np.linspace(0,40,350)call the function using odent from scipy,integrate library.
scipy.integrate.odeint(func, y0, t, args=(), Dfun=None, col_deriv=0, full_output=0, ml=None, mu=None, rtol=None, atol=None, tcrit=None, h0=0.0, hmax=0.0, hmin=0.0, ixpr=0, mxstep=0, mxhnil=0, mxordn=12, mxords=5, printmessg=0, tfirst=False)
Integrate a system of ordinary differential equations.
# Solve ODE
theta = odeint(model,theta_0,t,args = (b,g,l,m))After solving ODE we can create animation of pendulum using theta array and plot angular displacement and velocity plot
Firstly, we are going to create sub plots to plot all three figures in one single figure.
fig = plt.figure(figsize=(10,10))
ax1 = plt.subplot2grid((9, 9), (0, 0), colspan=9,rowspan=3)
ax2 = plt.subplot2grid((9, 9), (4, 1), colspan=3,rowspan=3)
ax3 = plt.subplot2grid((9, 9), (4, 5), colspan=3,rowspan=3)
fig.suptitle(' The Pendulum Motion',size='20',weight='bold')Secondly, we are going to plot Angular displacement diagram following way.
ax1.plot(t,theta[:,0],'b-',label=r'$\frac{d\d\theta_1}{dt}=\theta_2$')
ax1.plot(t,theta[:,1],'r--',label=r'$\frac{d\theta_2}{dt}=-\frac{b}{m}\theta_2-\frac{g}{L}sin \theta_1$')
ax1.set_title('Angular Displacement',size='16',weight='bold')
ax1.set_ylabel('Plot',size='12',weight='bold')
ax1.set_xlabel('Angular Displacement',size='12',weight='bold')
ax1.legend(loc='best')
ax1.grid(color='b',linestyle=':',linewidth='0.5')Next, angular velocity diagram using following codes.
ax3.plot(theta[:,0],theta[:,1],color='r',linewidth='1')
ax3.set_title('Angular Velocity',size='16',weight='bold')
ax3.set_ylabel('Position radious',size='12',weight='bold')
ax3.set_xlabel('Angular velocity',size='12',weight='bold')
ax3.grid(color='b',linestyle=':',linewidth='0.5')
Lastly, for animation, we fix the one end of the string by providing (x,y) value.
Then, to create figure array at each calculation we are going to first value for ct = 1, witch are going to improve one by one using for loop.
x0 = 0
y0 = 0
ct=1
img_array = []
for i in theta[:,0]:
# Link 1, Link 2 connector
x1 = l*math.sin(i)
y1 = -l*math.cos(i)
x2 = (l/2)*math.sin(i)
y2 = -(l/2)*math.cos(i)
filename = str(ct) +'.jpg'
ct=ct+1
ax2.plot([x0,x1],[y0,y1],linewidth='3',color='b')
ax2.plot(x1,y1,marker="o", markersize=20, color='r')
ax2.plot([-0.6,0.6],[0.05,0.05],linewidth='10',color='k')
ax2.plot(x0,y0,marker="o", markersize=16, color='k')
ax2.text((x1-0.08), (y1-0.05), 'm', fontsize=12)
ax2.text((x2+0.05), y2, 'L=1m', fontsize=12)
ax2.set_title('The Pendulum motion',size='16',weight='bold')
ax2.set_xlim([-1.2,1.2])
ax2.set_ylim([-1.3,0.5])
ax2.set_xlabel('X-Axis',size='12',weight='bold')
ax2.set_ylabel('Y-Axis',size='12',weight='bold')
ax2.grid(color='b',linestyle=':',linewidth='0.5')
plt.savefig(filename)Output of above codes
Random sample Figure:-
Use cv2 to read all the saved images. Then mention image shape by providing height, weight and layers. Afterwards, mention size of the video file as a image file. And store all the frame in empty img array string.
img = cv2.imread(filename)
height, width, layers = img.shape
size = (width,height)
img_array.append(img)This time we create a VideoWriter object. We should specify the output file name (eg: output.avi). Then we should specify the FourCC code (details in next paragraph). Then number of frames per second (fps) and frame size should be passed. And the last one is the is Color flag. If it is True, the encoder expect color frame, otherwise it works with grayscale frame.
FourCCis a 4-byte code used to specify the video codec. The list of available codes can be found in fourcc.org. It is platform dependent. The following codecs work fine for me.
- In Fedora: DIVX, XVID, MJPG, X264, WMV1, WMV2. (XVID is more preferable. MJPG results in high size video. X264 gives very small size video)
- In Windows: DIVX (More to be tested and added)
- In OSX: MJPG (.mp4), DIVX (.avi), X264 (.mkv).
out = cv2.VideoWriter('project.avi',cv2.VideoWriter_fourcc(*'DIVX'), 15, size)Use for loop and write all the images in video file that we created above.
for i in range(len(img_array)):
out.write(img_array[i])
out.release()Output video Link (YOUTUBE LINK)
Conclusion:-
- Understood the basic analysis of the pendulum motion and spring mass damping system.
- Generated 350 number of images by assigning them separate increment numbers of name.
- Learn basics of cv2 library in python.
- Finally, generated the video file using cv2 library in python.
References:-
[1] ODE - https://mathinsight.org/ordinary_differential_equation_introduction
[2] Oscillation of a Simple Pendulum - https://www.acs.psu.edu/drussell/Demos/Pendulum/Pendulum.html#:~:text=By%20applying%20Newton's%20secont%20law,enough%2C%20so%20the%20small%20angle
[3] Newton’s Second Law for Rotation https://openstax.org/books/university-physics-volume-1/pages/10-7-newtons-second-law-for-rotation#:~:text=If%20more%20than%20one%20torque,i%CF%84i%3DI%CE%B1.
Complete Program:-
import numpy as np
from scipy.integrate import odeint
import matplotlib.pyplot as plt
import cv2
import math
# function 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 = 1
# initial condition
theta_0 = [0,3]
# time points
t = np.linspace(0,40,350)
# Solve ODE
theta = odeint(model,theta_0,t,args = (b,g,l,m))
x0 = 0
y0 = 0
ct=1
img_array = []
for i in theta[:,0]:
fig = plt.figure(figsize=(10,10))
ax1 = plt.subplot2grid((9, 9), (0, 0), colspan=9,rowspan=3)
ax2 = plt.subplot2grid((9, 9), (4, 1), colspan=3,rowspan=3)
ax3 = plt.subplot2grid((9, 9), (4, 5), colspan=3,rowspan=3)
fig.suptitle(' The Pendulum Motion',size='20',weight='bold')
ax1.plot(t,theta[:,0],'b-',label=r'$\frac{d\d\theta_1}{dt}=\theta_2$')
ax1.plot(t,theta[:,1],'r--',label=r'$\frac{d\theta_2}{dt}=-\frac{b}{m}\theta_2-\frac{g}{L}sin \theta_1$')
ax1.set_title('Angular Displacement',size='16',weight='bold')
ax1.set_ylabel('Plot',size='12',weight='bold')
ax1.set_xlabel('Angular Displacement',size='12',weight='bold')
ax1.legend(loc='best')
ax1.grid(color='b',linestyle=':',linewidth='0.5')
ax3.plot(theta[:,0],theta[:,1],color='r',linewidth='1')
ax3.set_title('Angular Velocity',size='16',weight='bold')
ax3.set_ylabel('Position radious',size='12',weight='bold')
ax3.set_xlabel('Angular velocity',size='12',weight='bold')
ax3.grid(color='b',linestyle=':',linewidth='0.5')
# Link 1, Link 2 connector
x1 = l*math.sin(i)
y1 = -l*math.cos(i)
x2 = (l/2)*math.sin(i)
y2 = -(l/2)*math.cos(i)
filename = str(ct) +'.jpg'
ct=ct+1
ax2.plot([x0,x1],[y0,y1],linewidth='3',color='b')
ax2.plot(x1,y1,marker="o", markersize=20, color='r')
ax2.plot([-0.6,0.6],[0.05,0.05],linewidth='10',color='k')
ax2.plot(x0,y0,marker="o", markersize=16, color='k')
ax2.text((x1-0.08), (y1-0.05), 'm', fontsize=12)
ax2.text((x2+0.05), y2, 'L=1m', fontsize=12)
ax2.set_title('The Pendulum motion',size='16',weight='bold')
ax2.set_xlim([-1.2,1.2])
ax2.set_ylim([-1.3,0.5])
ax2.set_xlabel('X-Axis',size='12',weight='bold')
ax2.set_ylabel('Y-Axis',size='12',weight='bold')
ax2.grid(color='b',linestyle=':',linewidth='0.5')
#plt.show()
#plt.savefig('test.png')
plt.savefig(filename)
img = cv2.imread(filename)
height, width, layers = img.shape
size = (width,height)
img_array.append(img)
out = cv2.VideoWriter('project.avi',cv2.VideoWriter_fourcc(*'DIVX'), 15, size)
for i in range(len(img_array)):
out.write(img_array[i])
out.release()Output video:-
https://www.youtube.com/watch?v=0_EZEqcDpSc
Screenshot of Video