Solving second order ODEs using Python

Simple Pendulum

A pendulum is a body suspended from a fixed support so that it swings freely back and forth under the influence of gravity. 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 it back and forth. The mathematics of pendulums are in general quite complicated. Simplifying assumptions can be made, which in the case of a simple pendulum allow the equations of motion to be solved analytically for small-angle oscillations.

Representation of the motion of a simple pendulum is given by the below equation

`(d^2theta)/dt^2+(gsintheta)/l=0`

where

g = gravity in m/`s^2`

l= length of the pendulum in m

`theta=` angle displaced

Considering a damping system is incorporated, the equation of motion of the pendulum becomes

`(d^2theta)/dt^2+b/m (d theta)/ dt+g/lsintheta=0`

where

b=Damping Coefficient

m= mass of the bob in kg

g = gravity in m/`s^2`

l= length of the pendulum in m

`theta=` angle displaced in radian

Formulation of the Program

To find the path we need to solve the second order diferential equation using the odeint function in Python we need to convert the given equation to a first order diferential equation.

`(d^2theta)/dt^2+b/m (d theta)/ dt+g/lsintheta=0`

Let `theta=theta_1` and `(d theta_1)/dt=theta_2` and substituting in the above equation, it becomes

`(d theta_2)/dt+b/m theta_2+g/lsintheta_1=0`

Transforming the above equation

`(d theta_2)/dt=-b/m theta_2-g/lsintheta_1`

To form the above equation we create a function named model where two first-order differential equation is stored in an array.

#fucntion that returns dz/dt

def model(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
	

Declaring Inputs

The initial angular displacement=0 and angular velocity=3 rad/sec at t=0 are stored in an array. The first initial condition of angular dispalcment=0 is used to solve the equation to find velocity and the angular velocity=3 rad/sec to solve the equation to find acceleration respectively.

#Inputs
b=0.05
g=9.81
l=1
m=0.1

#initial condition
theta_0=[0,3]

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

Explanation to get the plot

The size of the variable result is used to formulate the for loop according to the number of rows of theta.

Here we use only the result corresponding to the equation used to find the displacement as a function of time.

We fix the coordinates of the fixed end of the pendulum at (0,0) and the coordinates of the free end are found using the equation

`x=x_0+lsintheta`

`y=y_0-lcostheta`

`-` sign is used with the equation of y because `cos` of `-theta` is positive in the 3rd quadrant.

The for loop provides the corresponding values of `theta` and as per the index i. The coordinates are plotted based on the instantaneous values of x1 and y1.

The origin 1.5,1.5 is given a square shape with the required size.

The end position (x1,y1) of the pendulum a circle is created with the required size.

Complete Program

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

#fucntion that returns dz/dt

def model(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
g=9.81
l=1
m=0.1

#initial condition
theta_0=[0,3]

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

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

plt.figure(1)
plt.plot(t,theta[:,0],'b-',label=r'$frac{dtheta_1}{dt}=theta_2$')
plt.plot(t,theta[:,1],'r--',label=r'$frac{dtheta_2}{dt}=-frac{b}{m}theta_2-frac{g}{l}sintheta_1$')
plt.ylabel('Plot')
plt.xlabel('Time')
plt.legend(loc='best')
plt.show()

# coordinates for the pendulum
x0=1.5
y0=1.5

ct=1

for i in range (0,len(theta)):
    x1=x0+l*math.sin(theta[i,0])
    y1=y0+-l*math.cos(theta[i,0])
    filename='plot'+'%03d.png'%ct
    ct=ct+1
    plt.figure(i)
    plt.plot(x0,y0,marker="s",markersize=15)
    plt.plot([x0,x1],[y0,y1],linewidth=5)
    plt.plot(x1,y1,marker=".",markersize=90)
    plt.xlim([0,3])
    plt.ylim([0,2])
    plt.savefig(filename)

Results

The given second-order partial differential equation for a simple pendulum was successfully solved using the scipy 'odeint' function. The results can be used to analyse the behaviour of the simple pendulum by varying the parameters.

It is plot clearly represent the change in displacement and velocity with respect to time at a given damping coefficient. The damping coefficient can be varied in order to get the desired time in which the oscillations must be stopped.

Animation

The animation gives a better visualization of the pendulum. It depicts how the angular displacement and angular velocities vary with time.

Error

Colon was missing after for loop

  File "C:UsersbooksbekaDownloadsSimple_Pendulum.py", line 46
    for i in range (0,len(theta))
                                 ^
SyntaxError: invalid syntax
[Finished in 0.2s with exit code 1]
[shell_cmd: python -u "C:UsersbooksbekaDownloadsSimple_Pendulum.py"]
[dir: C:UsersbooksbekaDownloads]
[path: C:WINDOWSsystem32;C:WINDOWS;C:WINDOWSSystem32Wbem;C:WINDOWSSystem32WindowsPowerShellv1.0;C:Program FilesMATLABR2020bbin;C:WINDOWSSystem32OpenSSH;C:UsersbooksbekaAppDataLocalProgramsPythonPython39Scripts;C:UsersbooksbekaAppDataLocalProgramsPythonPython39;C:UsersbooksbekaAppDataLocalMicrosoftWindowsApps;C:UsersbooksbekaAppDataLocalProgramsPythonPython39;]