Ordinary Differential Equations : Practical work on the harmonic oscillator¶

In this example, you will simulate an harmonic oscillator and compare the numerical solution to the closed form one.

Theory¶

Read about the theory of harmonic oscillators on Wikipedia

Mechanical oscillator¶

The case of the one dimensional mechanical oscillator leads to the following equation:

\[m \ddot x + \mu \dot x + k x = m \ddot x_d\]

Where:

\(x\) is the position,
\(\dot x\) and \(\ddot x\) are respectively the speed and acceleration,
\(m\) is the mass,
\(\mu\) the
\(k\) the stiffness,
and \(\ddot x_d\) the driving acceleration which is null if the oscillator is free.

Canonical equation¶

Most 1D oscilators follow the same canonical equation:

\[\ddot x + 2 \zeta \omega_0 \dot x + \omega_0^2 x = \ddot x_d\]

Where:

\(\omega_0\) is the undamped pulsation,
\(\zeta\) is damping ratio,
\(\ddot x_d\) is the imposed acceleration.

In the case of the mechanical oscillator:

\[\omega_0 = \sqrt{\dfrac{k}{m}}\]

\[\zeta = \dfrac{\mu}{2\sqrt{mk}}\]

Undampened oscillator¶

First, you will focus on the case of an undamped free oscillator (\(\zeta = 0\), \(\ddot x_d = 0\)) with the following initial conditions:

\[\begin{split}\left \lbrace \begin{split} x(t = 0) = 1 \\ \dot x(t = 0) = 0 \end{split}\right.\end{split}\]

The closed form solution is:

\[x(t) = a\cos \omega_0 t\]

# Setup
%matplotlib inline
import numpy as np
import matplotlib.pyplot as plt
from scipy.integrate import odeint

# Setup
f0     = 1.
omega0 = 2. * np.pi * f0
a      = 1.

Part 1: theoretical solution¶

Plot the closed form solution of the undamped free oscillator for 5 periods.

Steps:

Create an array \(t\) reprenting time,
Create a function \(x_{th}\) representing the amplitude of the closed form solution,
Plot \(x_{th}\) vs \(t\).

# Complete here
#t =
#xth =

Part 2: Numerical solution with Euler integrator¶

Solve the problem introduced in question 1 with the Euler integrator.

Steps:

Rewrite the canonical equation as a system of first order ODEs depending of the variable \(X = [x, \dot x]\),
Code the derivative function \(f(X,t) = \dot X\),
Define initial conditions \(X_0\),
Solve the problem.
Plot the position \(x\) along and compare it with the theoretical solution.

Part 3: Energies an errors¶

Calculate and plot the kinetic energy \(E_c\), the potential energy \(E_p\) and the total energy \(E_t = E_c + E_p\), comment the result.

Steps:

Calculate \(E_c\),
Calculate \(E_p\),
Calculate \(E_t\),
Plot the evolution of the 3 energies. You can use plt.fill_between instead of plt.plot,
Use the results to define a relative error estimator base on energies.

Part 4: Numerical solution convergence¶

Plot the effect of the number time steps \(n_t\) on the error \(e\).

Steps:

Create an array containing the different number of time steps from 100 to 100000,
Loop over this array and calculate the the error for each configuration,
Plot the error as a function of \(n_t\).

Part 5: integrator benchmark¶

Rewrite the code of part 4 in order to compare the RK4 and ODEint solvers with the Euler solver. Comment the efficiency of each solver.

Part 6: Error vs. time¶

Modify the code of part 5 in order to measure the computing time of each method in each case. Plot the error vs. the computing time.