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:

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

The closed form solution is: $$ x(t) = a\cos \omega_0 t $$

Part 1: theoretical solution

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

Steps:

1. Create an array $t$ reprenting time,
2. Create a function $x_{th}$ representing the amplitude of the closed form solution,
3. Plot $x_{th}$ vs $t$.

Part 2: Numerical solution with Euler integrator

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

Steps:

1. Rewrite the canonical equation as a system of first order ODEs depending of the variable $X = [x, \dot x]$,
2. Code the derivative function $f(X,t) = \dot X$,
3. Define initial conditions $X_0$,
4. Solve the problem.
5. 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:

1. Calculate $E_c$,
2. Calculate $E_p$,
3. Calculate $E_t$,
4. Plot the evolution of the 3 energies. You can use plt.fill_between instead of plt.plot,
5. 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:

1. Create an array containing the different number of time steps from 100 to 100000,
2. Loop over this array and calculate the the error for each configuration,
3. 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.