# 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 npimport 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_ = self.omega_0        self.theta_ = 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_ = self.omega_0        self.theta_ = 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_ = self.theta_0*np.pi/180.0        self.time_= self.eta        self.theta_ = self.theta_+self.omega_0*self.eta +0.5* (self.eta**2)*alpha(self.theta_)                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.

--

--