Solving second order ODEs (Python)

Solving Second Order ODEs in Python

1. AIM : To simulate the transient behaviour (motion) of a simple pendulum in Python.

2.OBJECTIVE:

• To write a program in MATLAB to solve ODE (ordinary differential equation) for equation of motion of a simple pendulum with damping.
• To simulate the motion of the pendulum between time period of 0 – 20 sec.
• To create the animation of the motion of pendulum by getting its position as a function of time.

3.THEORY :

3.1 Simple Pendulum: - A simple pendulum is one which can be considered to be a point mass suspended from a string or rod of negligible mass. It is a resonant system with a single resonant frequency.

3.2 Pendulum Motion: - The motion of a simple pendulum is like ‘simple harmonic motion’ .

fig.1: Animation of the motion of a Simple Pendulum

3.3 Period of Simple Pendulum: - A point mass hanging on a massless string is an idealized example of a simple pendulum. When displaced from its equilibrium point, the restoring force which brings it back to the center.

3.4 Pendulum Geometry: - The period of a simple pendulum for small amplitudes 𝚹 is dependent only on the pendulum length and gravity.

fig.2: Simple Pendulum (free body diagram)

3.5 Pendulum Equation: -

• 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 2^nd order Ordinary Differential Equation that describes the position of the "ball" w.r.t time.

4. CALCULATION :

The ODE which represents the equation of motion of simple pendulum with damping is –

In the above equation,

g = gravity in m/s²,

L = length of the pendulum in m,

m = mass of the ball in kg,

b = damping coefficient.

Now,

Solving 2^nd Order differential equation –:

Let’s say,

𝚹 = 𝚹₁

Then,

`(dθ) /dt = (dθ_1)/dt = θ_2` (assumed) ------------------- ①

`(d^2 θ)/(dt^2 ) = (d^2 θ_1)/(dt^2 ) = d/dt ((dθ_1)/dt) = (dθ_2 )/dt` -------------------- ②

Now, going back to our original 2^nd order ODE equation,

`(d^2 θ)/(dt^2 ) + "b" /"m" (dθ)/dt + "g" /"L" sin⁡θ = 0`

`=>``(d^2 θ_1)/(dt^2 ) + "b" /"m" (dθ_1)/dt + "g" /"L" sin⁡θ_1 =0`

`=>``(d^2 θ_2)/(dt^2 ) + "b" /"m" (dθ_2)/dt + "g" /"L" sin⁡θ_1 =0`

So, in original words to solve this 2^nd order ODE; what we have done is,

we have created two ODEs –

1st ODE	2nd ODE
`(dθ_1)/dt = θ_2`	`(dθ_2)/dt = -"b" /"m" θ_2 - "g" /"L" sin⁡θ_1`

So, this is what I meant as breaking down a 2^nd order ODE into a series of 1^st order ODEs, and that is why we did this transformation.

• So, the idea is if we have a 3^rd order ODE, then we can break it down into a system of three 1^st order ODEs.

Then, we are going to say that –

`d/dt [ [θ_1], [θ_2] ] = [ [θ_2], [-"b" /"m" θ_2 - "g" /"L" sin⁡θ_1] ]`

`θ`

Above one is a solution vector. [We are calling it as a solution vector because it contains two quantities].

So, now when we are trying to solve it with programming; it can be any of the ODE solver through Python, we need to provide the equation on this particular format. And that is what we are going to do at next –

5. BODY OF CONTENT :

Python code for solving ODEs and animation -

"""
	A program for solving second order ODEs and simulate its motion
"""
import math
import matplotlib.pyplot as plt
import numpy as np
from scipy.integrate import odeint

def ode_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

# Inputs
# damping coefficient
b = 0.05
# acceleration due to gravity(m/s^2)
g = 9.81
# length of string (m)
l = 1
# mass of bob (kg)
m = 1

# Initial Conditions
theta_0 = [0,3]
# time points
t_span = np.linspace(0,20,100)

# Solve ODE
theta = odeint(ode_function,theta_0,t_span,args = (b,g,l,m))

# Plotting the graph
plt.plot(t_span,theta[:,0],'b-',label=r'$frac{dtheta1}{dt}=theta2$')
plt.plot(t_span,theta[:,1],'r--',label=r'$frac{dtheta2}{dt}=-frac{b}{m}theta2-frac{g}{l}sintheta1$')
plt.xlabel('Time (second)')
plt.ylabel('Plot')
plt.title('Damped Harmonic motion of a Simple Pendulum')
plt.legend(loc='best')
plt.grid()
plt.show()

# Iteration using for loop to create animation of motion of pendulum
loop_counter =1
tmp = theta[:,0]
for i in tmp:
	x_start = 0
	y_start = 0
	x_end = l*math.sin(i)
	y_end = -l*math.cos(i)

	# Plotting for animation
	plt.figure()
	# Plotting String of pendulum
	plt.plot([x_start,x_end],[y_start,y_end],linewidth=3)
	# Plotting Bob of pendulum
	plt.plot(x_end,y_end,'o',markersize=20,color='r')
	# Plotting the rigid support through which string of pendulum is attached
	plt.plot([-0.5,0.5],[0,0],linewidth=8,color='k')
	# Plotting hatched vertical line for refernce
	plt.plot([0,0],[0,-1.5],'--',color='y')
	plt.xlim([-1.25,1.25])
	plt.ylim([-2,1])
	plt.title('Animation of transient behaviour of a Simple Pendulum')
	filename = 'test%05d.png' % loop_counter
	plt.savefig(filename)
	loop_counter = loop_counter + 1

Code overview : -

• Initial inputs are taken according to the provided data as follows –
  1. L = 1 metre
  2. m = 1 kg
  3. b = 0.05
  4. g = 9.81 m/s²
  5. time period = 0-20 sec
  6. angular displacement at time (t=0) = 0
  7. angular velocity at time (t=0) = 3 rad/sec

• At first math, matplotlib, numpy and scipy modules were imported for their various functions.
• Initial inputs such as b,g,l,m are provided according to data in the question.
• An array of initial theta was created to define initial angular displacement and initial angular velocity.
• An array of a time span was created to store 100 equally spaced intervals between time 0 to 20 seconds.
• To solve the second-order ODE, a function named ode_function which stores the above matrix data in it.
• The given second-order ODE was solved using 'odeint' command followed by the proper syntax of the calling function and the result obtained was used to plot the graph of variation of angular displacement and angular velocity wrt time.
• Now for creating the animation of the obtained harmonic motion of the pendulum, for loop was used and a loop counter variables were initiated, the values of angular position is stored in a separate variable called 'tmp' and is iterated in for loop. In this execution, the one end of the string of pendulum was fixed and the other end moves according to the solution of the ODE solved and was plotted for each time interval of the time span.

6. RESULTS :

The plot generated as the output of code is -

• Animation of Pendulum Motion :-

7. CONCLUSION : The above plot of angular displacement and angular velocity shows the pendulum's motion to be gradually decreasing. This is due to the damping effect which exponentially decreases the amplitude and the velocity of motion.