Week 3 - Solving second order ODEs

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.

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.

Solution: -

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.

Given,

• L=1 meter,
• 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.

Equation

(d^2 θ/dt^2) +(dθ/dt)*b/m+(g/l)*sin⁡θ=0

θ=θ1

(d^2 θ1/dt^2 )+(dθ1/dt)*b/m+(g/l)*sin⁡θ1=0

dθ1/dt= θ2

d^2 θ1/dt^2=dθ2/dt

θ2/dt+(b/m)*θ2+(g/l)*sin⁡θ1=0

θ2/dt=-(b/m)*θ2-(g/l)*sin⁡θ1

dθ1/dt=θ2

Coding: -

import numpy as np 
from scipy.integrate import odeint
import math
import matplotlib.pyplot as plt 


#funcion that returns dtheta/dt

def function(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 conditions
theta_0 = [0,3]


#time points
t = np.linspace(0,20,200)

#solve ode

theta = odeint(function,theta_0,t,args = (b,g,l,m))

#plot result
plt.plot(t,theta[:,0],'r-',label=r'$\frac{d\theta_1}{dt}=\theta2$')
plt.plot(t,theta[:,1],'g--',label=r'$\frac{d\theta2}{dt}=-\frac{b}{m}\theta_2-\frac{g}{l}sin\theta_1$')
plt.xlabel('Plot')
plt.ylabel('Time')
plt.legend()
plt.savefig('Output.png')
plt.show()

#animation
ct = 1

x0 = 0
y0 = 0

for i in theta[:,0]:
	
	x1 = (1*math.sin(i))
	y1 = (-1*math.cos(i))

	filename = 'test%05d.png' % ct
	ct = ct +1
	
	plt.figure()
	plt.plot([-1,1],[0,0],"black",linewidth=10)
	plt.plot([x0,x1],[y0,y1],"red")
	plt.plot(x1,y1,'o',markersize = 20)
	plt.xlim([-1.5,1.5])
	plt.ylim([-1.5,1.0])
	plt.savefig(filename)

Explanation of your code.

• Import module numpy as np, scipy.intedrate import odeint for second ODE, math and for plots matplotlib.
• Define function that returns the value of dtheta/dt.
• Initialize the input values of b, g, m, l.
• Define initial condition form theta i.e. 0 and 3 rad/sec.
• Use np.linspace for creating numeric sequence; np.linspace(start point, last point, number of values take in between start and end point).
• Solve ODE by use of odeint command.
• Plot result by using plt.plot command to plot graph.
• For animation: - use for loop
• Name the filename
• Use plt.plot command to create every single point and after that stiches all the pic and create a gif

Results: -

Upload an animation on You tube and provide the You tube link for an animation as given below

https://youtu.be/4xb7l6ZCeoc

Output: - Angular displacement, angular velocity with respect to time.

Main Window