Week 3 - Solving second order ODEs

WEEK 3

ORDINARY DIFFERENTIAL EQUATIONS – THE SIMPLE PENDULUM

Objectives

• To solve the governing equation of a simple pendulum in the form of a second-order ordinary differential equation.
• To obtain a plot showing the variation of angular displacement and velocity with respect to time.
• To develop an animation of an oscillating simple pendulum.

Problem statement

According to Newton’s second law of motion, the governing equation of simple pendulum is as follows.

Here, b = damping co-efficient

m = mass of the ball(kg)

L = length of the string(m)

g = acceleration due to gravity(m/s^2)

To solve a higher order differential equation in python, it is initially broken down into a series of first-order differential equations by considering few assumptions. They are as follows.

By substituting these assumptions in the governing equations, two first-order ODEs are obtained as follows.

In this case, the first equation is used to obtain the values of angular displacement and the second equation is used to obtain the values of angular velocity. The input in the case of python is provided in a matrix format, in this case, consisting of 2 columns (angular displacement and velocity respectively) and multiple rows. Odeint is the command used to solve such problems which require a function with solving steps, initial conditions, and time as inputs. Matplotlib animation commands are used to obtain the animation later.

Solution Procedure

• Initially, the required modules are imported that would be used for calculations, plotting graphs, and generating animations.
• A function is defined that contains information regarding the ordinary differential equation to be solved. It includes defining arrays for storing results and statements that calculate necessary equations.
• Later, certain parameters – the damping co-efficient, length of the string, mass of the bob, and acceleration due to gravity are defined and initiated.
• Initial conditions for solving the ordinary differential equation – values of angular displacement and velocity at 0sec are defined.
• An array is created using the linspace command that creates and stores the values of the time.
• Now, the ode solution statement is written which uses the function defined above, initial conditions, and time as inputs.
• From the results, the plotting statements are written including plot indicators and legends.
• Now, the matplotlib animation module is used to obtain the animation of oscillating pendulum.
• In this case, initially, the overall plot size and subplot sizes are fixed by using axes class commands.
• Now, a for loop is used to determine different co-ordinates of the pendulum bob that varies with respect to time. Trigonometric relations are used in this case. Since the results of ode are stored in an array, len () command is used to access elements inside it and later carry out scalar operations.
• The color and line weight of the string are also defined along with marker command that defines the shape and size of the pendulum bob.
• Axes class is used again to define the aspect ratio and grid.
• Finally, an animation is created by using ArtistAnimation class which also includes defining several intervals between each frame.
• Later, FFmpeg is used to extract and save the video in .mp4 format.

The code is as follows.

Code

# The following code is used to solve the governing equation of a simple pendulum
import math
import matplotlib.pyplot as plt 
from scipy.integrate import odeint
import numpy as np 
import matplotlib.animation as animation


# A function is defined to solve the ode
def pendulum(theta,time,l,m,b,g):
	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
	
# intial parameters
l = 1
m = 1
b = 0.05
g = 9.81 

# initial conditions 
theta_0 = [0,3]

# time condition
time = np.linspace(0,20,150)

# solving the ode 
theta= odeint(pendulum,theta_0,time,args = (l,m,b,g))

# plotting results
plt.plot(time,theta[:,0],'b-',label = r'$\frac{d\theta_1}{dt} = \theta2$')
plt.plot(time,theta[:,1],'r--',label = r'$\frac{d\theta_2}{dt} = -\frac{b}{m}\theta_2 - \frac{g}{L}sin\theta_1$')
plt.xlabel('time')
plt.ylabel('plot')
plt.legend()
plt.show()
plt.savefig('Graph for corresponding ODE solution')

# Simple pendulum animation
# Figure sizes and subplot sizes are defined 
fig = plt.figure(figsize=(5, 5))
ax = fig.add_subplot(1, 1, 1)

# An empty array is defined to store the co-ordinates of all points during oscillation 
pdlm = []
for i in range(len(theta)):
    x, = ax.plot([0, np.sin(theta[i, 0])], [0, -np.cos(theta[i, 0])],color='k', lw=3,marker='o',markersize=30)
    
    pdlm.append([x])
ax.set_aspect('equal', 'datalim')
ax.grid()
ani = animation.ArtistAnimation(fig, pdlm, interval=20)


# Saving the animation
fn = 'simple_pendulum_animation.mp4'
writer = animation.FFMpegWriter(fps = 30)
ani.save(fn, writer = writer)

Errors

• TclError: invalid command name ".!canvas"

The error occurred is as follows

This error occurs when the code involves plt.show() after plotting and creating an animation. The error is rectified by eliminating this statement which leads to directly saving the animation before viewing it first.

• ValueError: x and y must have same first dimension, but have shapes (150,) and (0, 2)

The error occurred is as follows

This error because of a missing syntax. In this case, consider the code theta [:,0]. This statement indicates the use of individual elements in the array of theta and its corresponding use in plotting the co-ordinates. The error is rectified by correcting the syntax and using appropriate symbols.

• NameError: name 'animation' is not defined

The error obtained is as follows

Matplotlib is a collection of many modules including some regarding plotting and animations. The error occurs by not importing matplotlib.animation and is rectified by adding the import statement in the beginning along with other import statements.

Results

The plot showing the variation of angular displacement and velocity with respect to time is obtained as follows

The YouTube link for the animation is https://www.youtube.com/watch?v=c09snxC99T0

Discussions

Although the damping co-efficient is involved in this case, the effect is minimal since the value is very small in the problem statement (b = 0.05). From the graph above, the curve in blue represents angular displacement and the curve in red represents the angular velocity. Upon observation, the amplitude of the waves decreases slightly towards the end as damping plays a role in this case. Predominant change can be observed in the amplitudes of these curves by increasing the value of b.