Driven Harmonic Oscillator

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

1. Modello fisico

Lo script simula numericamente un oscillatore armonico smorzato soggetto a una forza periodica esterna:

$$\large m \ddot{x} + b \dot{x} + k x = F_0 \cos(\omega t) $$

dove:

• $m$ = massa
• $b$ = coefficiente di smorzamento viscoso
• $k$ = costante elastica
• $F_0$ = ampiezza della forza esterna
• $\omega$ = frequenza della forzante

---

2. Sistema equivalente del primo ordine

Definendo:

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

si ottiene:

$$\large \begin{cases} \dot{x} = v \\ \dot{v} = \dfrac{F_0 \cos(\omega t) - b v - k x}{m} \end{cases} $$

Questo sistema è implementato nella funzione:

$$\large f(\mathbf{y}, t) = \begin{pmatrix} v \\ \dfrac{F_0 \cos(\omega t) - b v - k x}{m} \end{pmatrix} $$

con $\mathbf{y} = (x,v)$.

---

3. Metodo Runge--Kutta del 4° ordine

Lo script integra numericamente usando:

$$\large y_{n+1} = y_n + \frac{\Delta t}{6}(k_1 + 2k_2 + 2k_3 + k_4) $$

dove:

$$\large k_1 = f(y_n, t_n) $$$$\large k_2 = f\left(y_n + \frac{\Delta t}{2}k_1, t_n + \frac{\Delta t}{2}\right) $$$$\large k_3 = f\left(y_n + \frac{\Delta t}{2}k_2, t_n + \frac{\Delta t}{2}\right) $$$$\large k_4 = f\left(y_n + \Delta t\, k_3, t_n + \Delta t\right) $$

Metodo molto più accurato dell'Eulero (errore $\mathcal{O}(\Delta t^4)$).

---

4. Parametri numerici usati

$$\large m = 1, \quad b = 0.2, \quad k = 1 $$$$\large F_0 = 1, \quad \omega = 1 $$

Condizioni iniziali:

$$\large x(0) = 0, \quad v(0) = 0 $$

Intervallo temporale:

$$\large t \in [0,50] $$

---

5. Comportamento fisico atteso

Il sistema mostra:

• fase transiente iniziale (dipende dalle condizioni iniziali)
• successivo regime stazionario forzato
• oscillazioni alla frequenza della forzante $\omega$
• ampiezza limitata dallo smorzamento
• possibile risonanza se $\omega \approx \sqrt{k/m}$

Nel caso numerico scelto:

$$\large \omega = \omega_0 = \sqrt{k/m} = 1 $$

⇒ il sistema lavora in condizioni di quasi risonanza smorzata.

---

6. Output dello script

Il codice produce:

• integrazione numerica della traiettoria $(x(t), v(t))$
• grafico dello spostamento nel tempo
• osservazione del passaggio dal transiente al regime forzato stazionario

---

Conclusione

Lo script implementa correttamente un integratore Runge--Kutta 4 per studiare la risposta dinamica di un oscillatore armonico smorzato forzato, evidenziando i fenomeni di risonanza e stabilizzazione dovuti allo smorzamento.

import numpy as np
import matplotlib.pyplot as plt 

def runge_kutta_4(f, y0, t):
    """
    Solve ODE using ^th order Runge-Kutta method.
        :param f: Derivative function describing the system of DDEs.
        :param yO: Initial state vector.
        :param t: Time array.
    :retum: 
        Array with the solution for each time step.
    """
    n = len(t)
    y = np.zeros((n, len(y0)))
    y[0] = y0
    for i in range(n - 1):
        dt = t[i+1] - t[i]
        k1 = f(y[i], t[i])
        k2 = f(y[i] + dt/2.0 * k1, t[i] + dt/2.0)
        k3 = f(y[i] + dt/2.0 * k2, t[i] + dt/2.0)
        k4 = f(y[i] + dt * k3, t[i] + dt)
        y[i+1] = y[i] + (dt/6.0) * (k1 + 2*k2 + 2*k3 + k4)
    return y 

def driven_oscillator_system(state, t, m, b, k, F0, omega):
    """
    Define the system of equations for the driven harmonic
    oscillator.
        :param state: Current state vector [x, v].
        :param t: Current time.
        :param m: Mass.
        :param b: Damping coefficient.
        :param k: Spring stiffness.
        :param FO: Driving amplitude.
        :param omega: Driving frequency.
    :return: 
        Derivatives [dx/dt, dv/dt],
    """
    x, v = state
    dxdt = v
    dvdt = (F0 * np.cos(omega *t) -b*v-k*x) / m
    return np.array([dxdt, dvdt]) 
    
# Parameters for the driven oscillator
m = 1.0 # mass
b = 0.2 # damping coefficient
k = 1.0 # spring constant
F0 = 1.0 # driving force amplitude
omega =1.0 # driving frequency
initial_state = [0.0, 0.0] # initial position and velocity 

# Time array
t = np.linspace(0, 50, 1000) 

# Solve ODE
solution = runge_kutta_4(lambda state, t: driven_oscillator_system(state, t, m, b, k, F0, omega), initial_state, t) 

# Extract position and velocity
x, v = solution.T

# Plot results
plt.figure(figsize=(8, 3))
plt.subplot(1, 2, 1)
plt.plot(t, x, label='Position')
plt.xlabel('Time')

Text(0.5, 0, 'Time')