Understanding Oscillators (Python)

From the simple harmonic oscillation intuition to the numerical solution of differential equations.

The simple harmonic oscillator is one of the most fundamental phenomena in Physics. It essentially describes the motion of a mass attached to a spring. The equation describing its behaviour is born as a consequence of combining Newton’s second law and Hooke’s law, both definitions of the concept of force.

Newton’s 2nd Law. Eq(1)

Hooke’s Law. Eq(2)

Plugging Newton’s law into Hooke’s law results in the simple harmonic oscillator differential equation:

Where omega squared = √k/m. Eq(3)

Voilà! Such a simple formula, yet so useful in Newtonian mechanics.

Differential equations can be solved either analytically (if the solution is exact) or numerically (if the solution has no exact form). The solution to this equation can be found analytically, and the general form is:

Linear superposition of trigonometric functions, where A and B are amplitude constants. Eq(4)

I recommend playing around with these functions in websites that provide tools to graph them, as Desmos. This is really helpful to gain visual insight so that the learning process is enhanced. Alternatively, you can always plot them using Python or your preferred programming language.

With this brief explanation in mind, we are going to jump to the coding side, where we are going to use the SciPy library to find a numerical solution. The great deal of this methodology is its versatility since it can be generalized for multiple types of differential equations.

In Python

First, we import the relevant libraries.

import numpy as np 
from scipy.integrate import solve_ivp
import matplotlib.pyplot as plt
import seaborn as sns
sns.set()

Then we define a range of integration for time, a value for omega squared and the initial values for the position (x) and velocity (v), which in this case are 0 and 2 respectively as shown in the variable y

t = np.linspace(0,15,1000)
omega_sq = 1
y = [0,2] #y[0]=x and y[1]=v

Now is where the fun begins. First, we decompose the second-order ODE (Ordinary Differential Equation) into two first-order differential equations.

Velocity is the rate of change of position with respect to time. Eq(5)

Then the process is as follows, we define a vector with the variables being differentiated, as you can see above those are x and v, both varying as a function of time (Eq(5)).

General definition. Eq(6)

If the understanding of Eq(5) is correct, then it is only natural to insert our variables in the general definition (Eq(6)).

The general definition used for our example (from Eq(3))

Once this point has been reached, we only have to create a function that returns the vector dy/dt, as shown in the code below

def harmonic(t,y):
    solution = [y[1],-omega_sq*y[0]]
    return solution
sho = solve_ivp(harmonic, [0,1000], y0 = y, t_eval = t)

Solve_ivp is a function of the SciPy library that solves initial value problems for systems of ODEs, which happens to be this case. For extra information, visit the online documentation.

The solution to any initial value problem is a particular solution, and for this case given the initial conditions stated before, we reach a trigonometric function, in agreement with the general solution, Eq(4). Now we just use the matplotlib library to visualize this particular solution.

plt.plot(t,sho.y[0])
plt.ylabel("Position")
plt.xlabel("Time")
plt.title('SHO', fontsize = 20)

The oscillatory pattern x(t) given certain initial conditions.

This methodology works similarly for damped vibrations. Since in real life there is always energy dissipation (eternal motion machines don’t work in real life), we have to introduce a dragging force term in Eq(3) to obtain a more realistic approach to oscillations:

Damped Oscillator where the term dx/dt represents the dragging force. Eq(7)

The way to solve this is the same as for the harmonic oscillator:

t = np.linspace(0,15,1000)
y = [0,1]
gamma = 1
omega_sqr = 100
def sho(t,y):
    solution = (y[1],(-gamma*y[1]-omega_sqr*y[0]))
    return solution
solution = solve_ivp(sho, [0,1000], y0 = y, t_eval = t)
plt.plot(t,solution.y[0])
plt.ylabel("Position")
plt.xlabel("Time")
plt.title('Damped Oscillator', fontsize = 20)

Damped vibration function x(t)

We can see from this particular numerical solution how the oscillation exponentially decays into equilibrium, where x(t) = 0. In the previous analysis for the SHO, such a decay does not appear and the system remains in motion forever.

There are entire books dedicated to studying oscillations and waves, there is a lot more to explore on this, and it sets a basis for several physical phenomena in various fields, even quantum mechanics. This post intended to give a basic overview of vibrations with the addition of handy numerical solutions. Feel free to copy the code and play around with various values and explore the solutions, possibly one of the best ways to gain insight from equations.

Have fun and ride the wave!