Python: Solving Second Order ODEs

AIM: To write a program for simulating the transient behaviour of a simple pendulum

OBECTIVES:

1. Solve the second order ODE for the motion of a simple pendulum.
2. Generate a plot for its angular displacement and angular velocity wrt time.
3. Create an animation for the motion of the pendulum.

GOVERNING EQUATION:

Ordinary differential equation for oscillating motion of non-linear pendulum,

`"[(d^2 θ)/(dt^2 )]+[(b/m)*(dθ/dt)]+[(g/l)* sinθ]=0"`

Where,

g = gravity in m/s2,

l = length of the pendulum in m,

m = mass of the ball in kg,

b = damping coefficient.

INPUTS:

L = 1 m,

m = 1 kg,

b = 0.05,

g = 9.81 m/s^2,

time span = 0-20 sec,

angular displacement = 0,

angular velocity = 3 rad/sec.

THEORY:

In mathematics, a differential equation is an equation that relates one or more functions and their derivatives. In applications, the functions generally represent physical quantities, the derivatives represent their rates of change, and the differential equation defines a relationship between the two. Such relations are common, therefore differential equations play a prominent role in many disciplines including engineering, physics, economics, and biology.

In mathematics, an ordinary differential equation (ODE) is a differential equation containing one or more functions of one independent variable and the derivatives of those functions. The term ordinary is used in contrast with the term partial differential equation which may be with respect to more than one independent variable.

A pendulum is a weight suspended from a pivot so that it can swing freely. When a pendulum is displaced sideways from its resting, equilibrium position, it is subject to a restoring force due to gravity that will accelerate it back toward the equilibrium position. When released, the restoring force acting on the pendulum's mass causes it to oscillate about the equilibrium position, swinging back and forth. The time for one complete cycle, a left swing and a right swing, is called the period. The period depends on the length of the pendulum and also to a slight degree on the amplitude, the width of the pendulum's swing.

STEPS:

First, we import important libraries/modules. In this case we import math, mathplotlib.pyplot (as plt), numpy (as np) and from scipy.integrate we import ‘odeint’ which is our ODE solver. This solver is crucial for solving any ODE in python.

To solve the second order ODE a function named ‘ode_function’ is defined. This function will be helpful for storing useful data. Initial inputs such as damping co-efficient, acceleration due to gravity, length of string, mass of bob, angular displacement and angular velocity are provided. A timespan of 20 seconds with 100 intervals is also provided.

The given second-order ODE is solved using 'odeint' command followed by the proper syntax of the calling function and the result obtained. These results are used for plot which shows the angular displacement and angular velocity with respect to time.

Now for creating the animation of the obtained harmonic motion of the pendulum, for loop is used and a loop counter variable is 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.

CODE:

''' A program for solving second order ODE and simulate the motion of pendulum
By Utkarsh Palwekar'''
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
b = 0.05 # Damping coefficient
g = 9.81 # Acceleration due to gravity(m/s^2)
l = 1 # Length of string (m)
m = 1 # Mass of bob (kg)

# Initial Conditions
# Here, angular displacement=0 & angular velocity = 3
theta_0 = [0,3]
# time span
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'Angular Displacement (rad)')
plt.plot(t_span,theta[:,1],'r--',label=r'Angular Velocity (m/s^2)')
plt.xlabel('Time (sec)')
plt.ylabel('Plot')
plt.title('Damped Harmonic motion of a Simple Pendulum')
plt.legend(loc='best')
plt.grid()
plt.show()

#  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,color='y')
   # 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')
   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

 

RESULTS/OUTPUT:

Plot of Damped harmonic motion of a simple pendulum.

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.

Animation of Simple Pendulum.

CONCLUSION:

In this way, we can understand that any second-order or any higher-order ODE can easily be solved using Python, and simulation of the transient behaviour of simple pendulum was made using it to explain the effect of damping on its motion.

ERRORS & SOLUTION:

Some minor typo error and indentation errors were encountered.

No major error was encountered.