Damped Harmonic Oscillator

da PowerShell -- jupyter nbconvert --to html notebook.ipynb

1. Modello fisico

Lo script Python simula numericamente il moto di un oscillatore armonico smorzato descritto dall'equazione differenziale lineare del secondo ordine:

$$\large m \ddot{x} + b \dot{x} + k x = 0 $$

dove:

• $b$ è il coefficiente di smorzamento viscoso,
• $k$ è la costante elastica,
• $x(t)$ è lo spostamento.

2. Riduzione a sistema del primo ordine

Introducendo la velocità:

$$\large v = \dot{x} $$

l'equazione diventa il sistema:

$$\large \begin{cases} \dot{x} = v \\ \dot{v} = -\dfrac{b}{m} v - \dfrac{k}{m} x \end{cases} $$

3. Discretizzazione numerica (Metodo di Eulero)

Lo script utilizza il metodo di Eulero esplicito con passo temporale $\Delta t$:

$$\large a_n = -\frac{b v_n + k x_n}{m} $$$$\large v_{n+1} = v_n + a_n \Delta t $$$$\large x_{n+1} = x_n + v_n \Delta t $$

Questa procedura approssima la soluzione continua tramite un avanzamento iterativo nel tempo.

4. Condizioni iniziali

La simulazione parte da:

$$\large x(0) = x_0, \qquad v(0) = v_0 $$

Nel codice:

$$\large x_0 = 1, \quad v_0 = 0 $$

5. Regime dinamico

Il comportamento dipende dal discriminante:

$$\large b^2 - 4mk $$

Nel caso dello script:

$$\large b^2 < 4mk $$

quindi il sistema è $\textbf{sottosmorzato}$ e la soluzione analitica è:

$$\large x(t) = A e^{-\frac{b}{2m}t} \cos(\omega_d t + \phi) $$

con

$$\large \omega_d = \sqrt{\frac{k}{m} - \frac{b^2}{4m^2}} $$

6. Significato della simulazione

Il grafico prodotto mostra:

• oscillazioni armoniche,
• ampiezza che decade esponenzialmente,
• dissipazione di energia dovuta allo smorzamento.

La simulazione numerica riproduce quindi qualitativamente la dinamica prevista dalla teoria.

import numpy as np
import matplotlib.pyplot as plt 

def damped_oscillator_equation(m, b, k, x0, v0, t_max, dt):
    """
    Numerical simulation of a damped harmonic oscillator.
        :param m: Mass of the oscillator.
        :param b: Damping coefficient.
        :param k: Spring constant.
        :param xO: Initial position.
        :param vO: Initial velocity.
        :param t_max: End time for simulation.
        :param dt: Time step for numerical integration.
    :return: 
        Arrays of time, position, and velocity.
    """
    # Initialize time and state variable arrays
    t_values = np.arange(0, t_max, dt)
    x_values = np.zeros_like(t_values)
    v_values = np.zeros_like(t_values)

    # Set initial conditions
    x_values[0] = x0
    v_values[0] = v0

    # Time-stepping solution using Euler's method
    i = 0
    for j in range(len(t_values)-1) :
        i += 1
        a = -(b * v_values[i-1] + k * x_values[i-1]) / m
        v_values[i] = v_values[i-1] + a * dt
        x_values[i] = x_values[i-1] + v_values[i-1] * dt

    return t_values, x_values, v_values 
    
def plot_damped_oscillator(t, x, title='Damped Harmonic Oscillator'):
    """
    Plot the position of the damped harmonic oscillator over time.
        :param t: Time array.
        :param x: Position array.
        :param title: Title for the plot.
    """
    plt.figure(figsize=(8, 3))
    plt.plot(t, x, label='Position (x)')
    plt.title(title)
    plt.xlabel('Time (s)')
    plt.ylabel('Displacement (m)')
    plt.grid()
    plt.legend()
    plt.show()
    
# Define parameters
mass = 1.0 # Mass (kg)
damping = 0.5 # Damping coefficient (Ns/m)
stiffness = 1.0 # Spring constant (N/m)
initial_position =1.0 # Initial displacement (m)
initial_velocity = 0.0 # Initial velocity (m/s)
time_end = 10.0 # Duration of simulation (s)
time_step = 0.01 # Time step (s)

# Run the simulation
time, position, velocity = damped_oscillator_equation(mass, damping, stiffness, initial_position, initial_velocity, time_end, time_step) 

# Plot the results
plot_damped_oscillator(time, position)