Modern Physics
Published in

Modern Physics

Simple Pendulum ODESolver using python

Illustration of a simple pendulum. Source: https://phet.colorado.edu/sims/html/pendulum-lab/latest/pendulum-lab_en.html

Introduction

The simple pendulum is an example of a classical oscillating system. Classical harmonic motion and its quantum analogue represent one of the most fundamental physical model. The harmonic oscillator model can be used to describe and model phenomena such as heat, molecular and crystal bonding, lattice vibration, electromagnetic waves, vibrational spectroscopy, water waves, shock absorbers, sound waves, acoustics, earthquakes, etc.

In this article, we describe 3 basic methods that can be used for solving the second-order ODE (ordinary differential equation) for a simple harmonic oscillating system. We then implement the 3 basic methods using a python solver.

General Formalism

We remark here that the 3 methods described above could be extended to include other external forces such as damping or frictional forces.

Python ODESolver for the Simple Pendulum

Inport Necessary Libraries

import numpy as np
import matplotlib.pyplot as plt

ODE Solver

class ODESolver(object):
"""Second-order ODE Solver.
Parameters
------------
omega_0 : float
initial angular velocity
theta_0 : float
initial angular displacement
eta : float
time step size
n_iter : int
number of steps

Attributes
-----------
time_ : 1d-array
Stores time values for each time step.
omega_ : 1d-array
Stores angular velocity values for each time step.
theta_ : 1d-arra
Stores angular displacement values for each time step.

Methods
-----------
euler(alpha): Implements the Euler algorithm for the acceleration function alpha.

midpoint(alpha): Implements the Midpoint algorithm for the acceleration function alpha.

verlet(alpha): Implements the Verlet algorithm for the acceleration function alpha.
"""
def __init__(self, omega_0 = 0, theta_0 = 10, eta=0.01, n_iter=10):
self.omega_0 = omega_0
self.theta_0 = theta_0
self.eta = eta
self.n_iter = n_iter

def euler(self,alpha):
"""Implements Euler Method.

Parameters
----------
alpha : acceleration function

Returns
-------
self : object
"""
self.time_ = np.zeros(self.n_iter)
self.omega_ = np.zeros(self.n_iter)
self.theta_ = np.zeros(self.n_iter)
self.omega_[0] = self.omega_0
self.theta_[0] = self.theta_0*np.pi/180.0

for i in range(self.n_iter-1):
self.time_[i+1] = self.time_[i] + self.eta
self.omega_[i+1] = self.omega_[i] + self.eta*alpha(self.theta_[i])
self.theta_[i+1] = self.theta_[i] + self.eta*self.omega_[i]
return self

def midpoint(self,alpha):
"""Implement Midpoint Method.

Parameters
----------
alpha : acceleration function
Returns
-------
self : object
"""
self.time_ = np.zeros(self.n_iter)
self.omega_ = np.zeros(self.n_iter)
self.theta_ = np.zeros(self.n_iter)
self.omega_[0] = self.omega_0
self.theta_[0] = self.theta_0*np.pi/180.0

for i in range(self.n_iter-1):
self.time_[i+1] = self.time_[i] + self.eta
self.omega_[i+1] = self.omega_[i] + self.eta*alpha(self.theta_[i])
self.theta_[i+1] = self.theta_[i] + 0.5*self.eta*(self.omega_[i]+self.omega_[i+1])
return self

def verlet(self,alpha):
"""Implement Verlet Method.

Parameters
----------
alpha : acceleration function
Returns
-------
self : object
"""
self.time_ = np.zeros(self.n_iter)
self.theta_ = np.zeros(self.n_iter)
self.theta_[0] = self.theta_0*np.pi/180.0
self.time_[1]= self.eta
self.theta_[1] = self.theta_[0]+self.omega_0*self.eta +0.5* (self.eta**2)*alpha(self.theta_[0])

for i in range(self.n_iter-2):
self.time_[i+2] = self.time_[i+1] + self.eta
self.theta_[i+2] = 2.0*self.theta_[i+1] -self.theta_[i] + (self.eta**2)*alpha(self.theta_[i+1])
return self

Define Angular Acceleration Function

def alpha(x):
return -np.sin(x)

Example 1: Euler Method

time=ODESolver(omega_0 = 0, theta_0 = 10, eta=0.1, n_iter=300).euler(alpha).time_
theta=ODESolver(omega_0 = 0, theta_0 = 10, eta=0.1, n_iter=300).euler(alpha).theta_
plt.plot(time,theta*180/np.pi,lw=3,color='red')
plt.xlabel('time(s)',size=13)
plt.ylabel('angle (deg)',size=13)
plt.title('Euler Method',size=13)
plt.show()

We observe that with a time step of 0.1, the Euler method gives a solution that is not stable. This problem can be solved by decreasing the time step to smaller value, for example 0.001.

Example 2: Midpoint Method

time=ODESolver(omega_0 = 0, theta_0 = 10, eta=0.1, n_iter=300).midpoint(alpha).time_
theta=ODESolver(omega_0 = 0, theta_0 = 10, eta=0.1, n_iter=300).midpoint(alpha).theta_
plt.plot(time,theta*180/np.pi,lw=3,color='green')
plt.xlabel('time(s)',size=13)
plt.ylabel('angle (deg)',size=13)
plt.title('Midpoint Method',size=13)
plt.show()

We observe that with a time step of 0.1, the Midpoint method gives a solution that is not stable, but relatively better when compared to the Euler method. This problem can be solved by decreasing the time step to smaller value, for example 0.001.

Example 3: Verlet Method

time=ODESolver(omega_0 = 0, theta_0 = 10, eta=0.1, n_iter=300).verlet(alpha).time_
theta=ODESolver(omega_0 = 0, theta_0 = 10, eta=0.1, n_iter=300).verlet(alpha).theta_
plt.plot(time,theta*180/np.pi,lw=3,color='blue')
plt.xlabel('time(s)',size=13)
plt.ylabel('angle (deg)',size=13)
plt.title('Verlet Method',size=13)
plt.show()

We observe that with a time step of 0.1, the Verlet method gives a reasonable solution that is stable. Because the Verlet method is based on the centered derivative while Euler and Midpoint uses the forward derivative, the error in the Verlet method is quite minimal.

Summary

In summary, we’ve shown how a python object can be built for implementing the 3 basic methods for solving second-order ODE’s. We also provided some sample outputs from the code developed. Based on this analysis, we observe that the Verlet method is computationally the most efficient method since it uses the centered derivative which is a more symmetric definition of a derivative compared to the forward derivative.

--

--

Get the Medium app

A button that says 'Download on the App Store', and if clicked it will lead you to the iOS App store
A button that says 'Get it on, Google Play', and if clicked it will lead you to the Google Play store