Week 3 - Solving second order ODEs

Question:

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,

Objective: To obtain velocity and displacement with respect to time plots by solving a 2^nd order ODE for a simple pendulum using a Python program and animate the motion of the pendulum using the information.

Solution:

Using the mass-spring-damper system method, a simple pendulum motion can be described using the equation of motion as:

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

where,

θ = angular displacement in [radians]

c = damping co-efficient [Ns/m]
m = Mass of the pendulum in [Kg]
l = Length of the string in [meters]
g = acceleration due to gravity in [m/s^2]

As the equation has a highest degree of ‘2’, it is a 2^nd order ODE equation and needs to be converted into 2 simple 1^st order ODEs to sole it using programming. We do that by substituting 1^st order derivative term with a variable ‘θ₂’ such that,

θ = θ₁,

1. `d/dt`(θ₁) = θ₂, and
2. `d/dt`(θ₂) =`d/dt(d/dt`(θ₁)) = `d^2/dt^2`(θ₁)

To solve the ODE in python (integrate) we use ‘odeint’ command from the ‘scipy’ module. Here, we define the ODE as a function, and we solve it for certain initial conditions over an array of ‘t’ time steps which are determined by the program itself. We also have to define the constants/arguments which are a part of the function in specific order.

By integrating the ODEs we get an array/list with 2 columns and multiple rows as solution. The 1^st column has the solution of 1^st ODE that is ‘θ₁’ or the angular displacement and the 2^nd column has solution for 2^nd ODE – ‘θ₂’ = θ₁^’ = ‘ω’ or the angular velocity. Using the 1^st solution, θ₁, we can find the co-ordinates of the pendulum for each solution and by assuming the pendulum to be fixed at (0,0) we can plot it as a line for every solution and get the animation.

Python Program:

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

def equation_of_motion_model(theta,t,c,g,l,m):
    theta1 = theta[0]
    theta2 = theta[1]
    dtheta1_dt = theta2
    dtheta2_dt = -(c/m) * theta2 - (g/l) * math.sin(theta1)
    dtheta_dt = [dtheta1_dt, dtheta2_dt]
    return dtheta_dt


#the motion of a pendulum is governed by following equation: d^2/dt^2(theta) + (c/m)*d/dt(theta) + g/l*sin(theta)
#Input variables:
c = 0.05 #damping co-efficient Ns/m
m = 1 #Mass of the pendulum in Kg
l = 1 #Length of the string in meters
g = 9.81 #acceleration due to gravity in m/s^2


#Initial conditions to start the simulation from: thetax=[theta_initial, omega_initial]
theta0 = [0,3]

#defining time window and time steps for calculation:
t = np.linspace(0,20,301)

#now, integrating the ODEs to get solution arrays to be plotted:
theta = odeint(equation_of_motion_model,theta0,t,args=(c,g,l,m))
#print(theta)


'''
plt.plot(t,theta[:,0],'b-',label='Angular displacement')
plt.plot(t,theta[:,1],'r--',label='Angular velocity')
plt.legend(loc='best')
plt.xlabel('Time in seconds')
plt.ylabel('Model outputs')
plt.grid()
plt.show()
'''
'''
x = equation_of_motion_model(thetax,c,m,g,l)
print(x)
'''

#pendulum animation:theta[0] represents angular displacement of the pendulum hence is to be plotted as a plot for each timestep.
#as length of the pendulum string is 1m, the angle pendulum makes with the vertical is theta,
x0 = [0]
y0 = l

x1 = theta[:,0]
x3 = []
y3 = []
ct = 1
for x2 in x1:
    x = l * math.sin(x2)
    y = -1 * (l * math.cos(x2))
    sphere = mpath.Path.unit_circle()
    plt.plot([0,x],[0,y],'g',marker = sphere, markersize = 6)
    plt.plot([-0.75,0.75],[0,0],'black',linewidth = 7)
    plt.plot([0,0],[-1.2,0.2],'b--')
    plt.ylim(-1.25,0.05)
    plt.xlim(-1.25,1.25)
    plt.grid()
    ct = ct + 1
    filename = str(ct)+'.png'
    plt.savefig(filename)

    plt.show()

    #x3.append(x)
    #y3.append(y)

#print(x3)
#print(y3)
j = []
k = []

#plt.plot([x3],[y3],'r')
#plt.show()

Procedure:

1. Initially all the necessary modules are imported and stored to be called back later.
2. As mentioned before the equation of motion and the splitting of the 2^nd order ODE into two 1^st order ODEs is done using a function option to suit the ‘odeint’ format we have used in this program. To make sure the solution contains both the 1^st order ODEs, we store ‘dtheta1_dt’ and ‘dtheta2_dt’ in a single variable and use ‘return’ command to make sure the parameters are returned after the loop is run.
3. Now, the constants form the equation – c,l,g,m are defined and their respective values are input. As the ‘odeint’ requires initial conditions for integrating the ODEs, they were defined for initial displacement and initial velocity of the pendulum as ‘θ₀’.
4. As we were to plot the results against time, an array was defined using ‘numpy’ module and ‘np’ command for t = 0 to t = t_n, with ‘n’ steps in ‘t’.
5. Using the ‘odeint’ option, and following its format the integration command was entered as equal to variable ‘θ’ which consists of 2 columns and ‘n’ rows for each time step. As stated previously the two columns represent 2 solutions each of the ODEs and are plotted against time using ‘matplotlib’ module.
6. For the plots describing the motion of pendulum, first the co-ordinates of the pendulum were calculated from the simple relation x = l sin(θ) and y = -1 (l cos(θ)), where, ‘l’ is the length of the pendulum string and the negative sign for the ‘y’ co-ordinate inverts the direction to make it look like an actual pendulum. Now, these co-ordinated were plotted for each value of ‘θ’ using a ‘for’ loop to get individual plots for each solution. Using the counter variables ‘ct’ unique file names were generated and the plots were saved as images. These images then were stitched to form a ‘.gif’ flipbook animation similar to the one done for 2R robotic arm.

Plots:

1. Example plot:

2. Challenge problem plot:

3. Flipbook animation:

4. Error faced:

Due to the incorrect order in which the arguments were input in the ‘odeint’ command, incorrect output plots were being generated. After going through the ‘odeint’ tutorials on various sites and scipy.integrate.odepack section on scipy.org as well through my IDE pycharm I realised there is specific order to it and after jumbling them around I found the right order.

Conclusion:

1. The method used in this example can be used for even higher ordered derivative functions to be solved using programming. The splitting of higher order ODEs into a bunch of 1st order ODEs is simpler to solve and with the help of 'def' and function method can be automated with relative ease.

2. This program with some changes can be used to describe any spring-mass-damper system and can be altered to output various results for change in lets say damping coefficient, or length, initial conditions, time period etc. for a comparative study in design phase.