Small Oscillations

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

Pendoli che si sincronizzano - interazione viscosa

Sì, il fenomeno esiste ed è ben documentato: pendoli nominalmente identici e presenti nella stessa stanza possono sincronizzarsi spontaneamente grazie a un’interazione debole, tipicamente attraverso il supporto comune (la trave, il muro, il tavolo) che funziona come un mezzo di accoppiamento dissipativo/viscoso.

$$ \large interazione \; \propto\, c\, (\dot{\theta_j}-\dot{\theta_i})$$

import numpy as np
import matplotlib.pyplot as plt
from scipy.integrate import solve_ivp

# ---------------------------------------------------------
# Modello: N pendoli accoppiati debolmente (accoppiamento viscoso)
# ---------------------------------------------------------

def coupled_pendula(t, y, N, omegas, c):
    """
    y = [theta1, theta2, ..., thetaN, dtheta1, ..., dthetaN]
    """
    thetas = y[:N]
    dthetas = y[N:]
    dd = np.zeros(N)

    # Equazioni del moto: ddtheta_i = -omega_i^2 * theta_i  + accoppiamento viscoso
    for i in range(N):
        # damping accoppiato con tutti gli altri pendoli
        coupling = 0
        for j in range(N):
            if j != i:
                coupling += c * (dthetas[j] - dthetas[i])

        dd[i] = - omegas[i]**2 * thetas[i] + coupling

    return np.concatenate([dthetas, dd])


# parametri
N = 5
omegas = 2*np.pi*(1 + 0.01*np.random.randn(N))   # piccole differenze nelle frequenze
c = 0.05                                         # forza di accoppiamento
t_max = 80

# condizioni iniziali
thetas0 = np.random.uniform(-0.5, 0.5, N)
dthetas0 = np.zeros(N)
y0 = np.concatenate([thetas0, dthetas0])

# integrazione
sol = solve_ivp(
    coupled_pendula, [0, t_max], y0,
    args=(N, omegas, c), t_eval=np.linspace(0, t_max, 2000)
)

# plotting
plt.figure(figsize=(10, 6))
for i in range(N):
    plt.plot(sol.t, sol.y[i], lw=1)
plt.xlabel("t")
plt.ylabel("theta(t)")
plt.title("Sincronizzazione spontanea di pendoli accoppiati")
plt.show()

Pendoli accoppiati con una molla

import numpy as np
def potential_energy_matrix(n, k, mg_l):
    """
    Constructs the stiffness matrix for a system of n coupled oscillators.
        : param n: Number of oscillators.
        :param k: Spring constant between oscillators.
        :param mg_l: Coupled weight force constant.
    :return: 
        Stiffness matrix C.
    """
    C = np.zeros((n, n))
    for i in range(n):
        C[i, i] = mg_l + k
        if i > 0:
            C[i, i - 1] = -k
            C[i - 1, i] = -k
    return C 
    
def normal_modes_analysis(mass, C):
    """
    Analyzes normal modes by solving the eigenvalue problem for coupled oscillations.
        :param mass: Mass of each oscillator.
        :param C: Potential energy matrix.
    :return: 
        Normal mode frequencies, mode shapes.
    """
    M_inv_sqrt = np.linalg.inv(np.sqrt(np.diag([mass] * len(C))))
    D = np.dot(M_inv_sqrt, np.dot(C, M_inv_sqrt))
    eigenvalues, eigenvectors = np.linalg.eigh(D)
    frequencies = np.sqrt(eigenvalues)
    mode_shapes = np.dot(M_inv_sqrt, eigenvectors)
    return frequencies, mode_shapes 
    
def coupled_pendula_example():
    """
    Example analysis of two coupled pendula with equal mass and spring connection.
    return: 
        Normal mode frequencies, mode shapes of the system.
    """
    n = 2
    mass = 1.0
    k = 0.5
    mg_l = 9.81
    C = potential_energy_matrix(n, k, mg_l)
    frequencies, mode_shapes = normal_modes_analysis(mass, C)
    return frequencies, mode_shapes # Outputs demonstration for coupled pendula

frequencies, mode_shapes = coupled_pendula_example() 
print("Normal Mode Frequencies:", frequencies)
print("Mode Shapes:\n", mode_shapes)

Normal Mode Frequencies: [3.13209195 3.28785644]
Mode Shapes:
 [[-0.70710678 -0.70710678]
 [-0.70710678  0.70710678]]

Pendoli accoppiati con una molla

import numpy as np

# -------------------------
# Parameters (SI units)
# -------------------------
m = 1.0      # mass [kg]
l = 1.0      # length [m]
k = 5.0      # spring constant [N/m]
g = 9.81     # gravity [m/s^2]

# -------------------------
# Mass and stiffness matrices
# -------------------------
M = m * l**2 * np.eye(2)

K = np.array([
    [m*g*l + k*l**2,   -k*l**2],
    [-k*l**2,          m*g*l + k*l**2]
])

# -------------------------
# Generalized eigenproblem
# -------------------------
eigvals, eigvecs = np.linalg.eig(np.linalg.inv(M) @ K)

# Angular frequencies
omega = np.sqrt(np.real(eigvals))

# Sort modes
idx = np.argsort(omega)
omega = omega[idx]
eigvecs = eigvecs[:, idx]

# Normalize mode shapes
mode_shapes = eigvecs / np.max(np.abs(eigvecs), axis=0)

# -------------------------
# Output
# -------------------------
for i in range(2):
    print(f"Mode {i+1}")
    print(f"  omega = {omega[i]:.4f} rad/s")
    print(f"  f     = {omega[i]/(2*np.pi):.4f} Hz")
    print(f"  mode shape = {mode_shapes[:, i]}")

Mode 1
  omega = 3.1321 rad/s
  f     = 0.4985 Hz
  mode shape = [1. 1.]
Mode 2
  omega = 4.4508 rad/s
  f     = 0.7084 Hz
  mode shape = [ 1. -1.]