Solution of ODE for pendulum using Python

Aim: To evaluate the motion of simple pendulum by solving ODE using Python

Theory:

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.

Consider a pendulum with ball of mass `m` attached to the end string of length `L` with other end being fixed.

where,

`theta` = Angular Displacement

`g` = Gravitational Acceleration

`t` = time

`b` = Damping Coefficient

For solving ODE using python, higher order differential equations needs to converted into a set of single order differential equations.

A second order differential equation should be convereted into two single order differential equations. This is done by considering,

`theta = theta_1`

`(d theta_1)/dt = theta_2`

and this leads to following differential equations,

`(d theta_1)/dt = theta_2`

`(d theta_2)/dt = -(b/m)*theta_2 - (g/l)*sintheta_1`

Now these equations can be solved using ode45 function of MATLAB.

The parameters choosen for study are as follows:

`L=1 m`

`m=1 kg`

`b=0.05`

`g=9.81 m/s^2`

with initial conditions,

`theta = 0`

`(d theta)/dt = 3`

Python Program:

"""
    Solving 2nd order ODE
"""

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

def model(theta,t,b,l,g,m):
    theta1 = theta[0]
    theta2 = theta[1]
    dtheta_1 = theta2
    dtheta_2 = -(b/m)*theta2 - (g/l)*math.sin(theta1)
    dtheta = [dtheta_1, dtheta_2]

    return dtheta

#input parameters
b = 0.05
g = 9.81
m = 1
l = 1

# initial values
theta_0 =  [0,3]

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

# Solve ODE
theta = odeint(model, theta_0, t, args = (b,l,g,m))

# plot result curve
plt.figure()
plt.plot(t,theta[:,0],'b-',label=r'$\frac{d\theta_1}{dt}=\theta_2$')
plt.plot(t,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(loc='best')
plt.savefig('ODE Solution')

# capture pendulum position
x0 = 0
y0 = 0
ct = 0
for i in range(0,len(t)):
    x1 = l*math.sin(theta[i,0])
    y1 = -l*math.cos(theta[i,0])
    filename = str(ct) + '.png'
    plt.figure()
    plt.plot([-1, 1],[0,0],'b',linewidth=5)
    plt.plot([x0, x1],[y0,y1],'g',linewidth=5)
    plt.plot(x1,y1,'o',markersize=30,markerfacecolor = 'r')
    plt.xlim([-2,2])
    plt.ylim([-2,2])
    plt.savefig(filename)
    ct = ct + 1

odeint is imported scipy,integrate to solve differential equation. A function model is defined to convert 2nd order differential equation into two single order differential equation. It return an array of single of differential equation.

Input parameters and initial conditions are defined by declaring appropriate variables. np.linspace is used to creates 200 equally spaced time points between 0 and 20 seconds. odeint function is used to solve differential equation, it takes model function, initial values, time points and required equation constants as arguments and returns angular dispacement and angular velocity at each time point in form of a 2D array. Results curve is plotted for dispacement and velocity with respect to time using plt.plot with appropriate legend and axis labels.

A for loop is defined for a range 0 to length of t to capture pendulum location at each time point. The coordinates are calculated using geometric relations and plotted. A counter variable ct is defined to keep track of the configuration at given time plot and it is used as filename to save the configuration at each time point.

Results:

Conclusion:

It can be seen from the plot that with time the displacement and velocity decreases, it is because damping is provided. If damping value is increased than displacement and velocity would come to zero in lesser time.

Animation: