Routhian Reduction¶
da PowerShell -- jupyter nbconvert --to html notebook.ipynb
La $Routhian Reduction$ è una tecnica della meccanica analitica che consente di ridurre il numero di gradi di libertà di un sistema lagrangiano quando sono presenti coordinate cicliche. Essa può essere vista come una riduzione parziale che combina elementi del formalismo lagrangiano e hamiltoniano.
Sistema lagrangiano¶
Si consideri un sistema meccanico con coordinate generalizzate $$ (q^1, \dots, q^n), $$ descritto da una Lagrangiana $$ L(q^1, \dots, q^n, \dot{q}^1, \dots, \dot{q}^n, t). $$
Una coordinata $q^k$ è detta $ciclica$ (o ignorable) se la Lagrangiana non dipende esplicitamente da essa: $$ \frac{\partial L}{\partial q^k} = 0. $$
Momento coniugato e costante del moto¶
Se $q^k$ è ciclica, il momento coniugato associato $$ p_k = \frac{\partial L}{\partial \dot{q}^k} $$ è una costante del moto: $$ \dot{p}_k = 0. $$
Questo fatto consente di eliminare la coordinata ciclica dal problema dinamico.
Definizione del Routhiano¶
Si separino le coordinate in due insiemi: $$ (q^1, \dots, q^r, q^{r+1}, \dots, q^n), $$ dove le ultime $n-r$ coordinate sono cicliche.
Il $Routhiano$ è definito come: $$ R(q^1, \dots, q^r, \dot{q}^1, \dots, \dot{q}^r, p_{r+1}, \dots, p_n, t) = \sum_{k=r+1}^{n} p_k \dot{q}^k - L. $$
In questa espressione, le velocità $\dot{q}^k$ cicliche vengono eliminate usando le relazioni $$ p_k = \frac{\partial L}{\partial \dot{q}^k}. $$
Equazioni del moto ridotte¶
Le equazioni di Eulero–Lagrange per le coordinate non cicliche derivano dal Routhiano: $$ \frac{d}{dt} \left( \frac{\partial R}{\partial \dot{q}^i} \right)
\frac{\partial R}{\partial q^i} = 0, \quad i = 1, \dots, r. $$
Per le coordinate cicliche, invece, valgono le equazioni hamiltoniane: $$ \dot{q}^k = \frac{\partial R}{\partial p_k}, \quad \dot{p}_k = -\frac{\partial R}{\partial q^k} = 0. $$
Interpretazione¶
Il Routhiano è una trasformazione di Legendre parziale della Lagrangiana, effettuata solo rispetto alle coordinate cicliche. In questo senso:
- per le coordinate cicliche il formalismo è hamiltoniano;
- per le coordinate non cicliche il formalismo resta lagrangiano.
Esempio: particella in campo centrale¶
Consideriamo una particella in un campo centrale con Lagrangiana: $$ L = \frac{1}{2}m(\dot{r}^2 + r^2 \dot{\theta}^2) - V(r). $$
La coordinata $\theta$ è ciclica, con momento angolare: $$ p_\theta = m r^2 \dot{\theta} = \ell. $$
Il Routhiano è: $$ R(r, \dot{r}, \ell) = \ell \dot{\theta} - L = \frac{1}{2}m \dot{r}^2 - \left( V(r) + \frac{\ell^2}{2 m r^2} \right). $$
Il problema si riduce a un moto unidimensionale in un potenziale efficace: $$ V_{\text{eff}}(r) = V(r) + \frac{\ell^2}{2 m r^2}. $$
Relazione con altre riduzioni¶
- La riduzione di Routh è un caso particolare della riduzione di Marsden--Weinstein.
- È strettamente legata ai teoremi di Noether e alle simmetrie cicliche.
- Fornisce un ponte concettuale tra meccanica lagrangiana e hamiltoniana.
Osservazioni finali¶
La Routhian Reduction è particolarmente utile in problemi con simmetrie continue, come sistemi rotazionalmente invarianti, meccanica celeste e dinamica dei sistemi rigidi.
Simple Pendulum with Rotating Support¶
$$ \large interazione \; \propto\, c\, (\dot{\theta_j}-\dot{\theta_i})$$
In [ ]:
import numpy as np
import matplotlib.pyplot as plt
from matplotlib.animation import FuncAnimation
from IPython.display import HTML
# Parametri fisici
g = 9.81
L = 1.0
Omega = 1.0
R = 0.5
dt = 0.02
t_max = 20
theta0 = 0.2
omega0 = 0.0
def simulate():
t_vals = np.arange(0, t_max, dt)
theta = theta0
omega = omega0
xs, ys = [], []
for t in t_vals:
omega += -(g/L) * np.sin(theta) * dt
theta += omega * dt
xp = R * np.cos(Omega * t)
yp = R * np.sin(Omega * t)
x = xp + L * np.sin(theta)
y = yp - L * np.cos(theta)
xs.append((xp, x))
ys.append((yp, y))
return t_vals, xs, ys
t_vals, xs, ys = simulate()
fig, ax = plt.subplots()
ax.set_aspect("equal")
ax.set_xlim(-R - L - 0.5, R + L + 0.5)
ax.set_ylim(-L - 1, L + 1)
line, = ax.plot([], [], "o-", lw=2)
def update(i):
xp, x = xs[i]
yp, y = ys[i]
line.set_data([xp, x], [yp, y])
return line,
anim = FuncAnimation(fig, update, frames=len(t_vals), interval=20)
# ---- QUI la parte importante per Jupyter ----
HTML(anim.to_html5_video())
Spherical Pendulum with Runge-Kutta¶
In [1]:
# Animazione 3D del pendolo sferico (Jupyter-friendly: salva GIF e la mostra)
import numpy as np
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D # noqa: F401 (serve per la proiezione 3D)
from matplotlib.animation import FuncAnimation, PillowWriter
from IPython.display import Image, display
# ---------------- Parametri fisici ----------------
g = 9.81 # gravità [m/s^2]
l = 1.0 # lunghezza del filo [m]
# ---------------- Parametri numerici ----------------
dt = 0.04 # passo temporale (s). Aumenta per meno frame / velocità.
t_max = 6.0 # durata totale (s). Aumenta per animazione più lunga.
t_vals = np.arange(0, t_max, dt)
# ---------------- Condizioni iniziali ----------------
theta0 = 0.8 # angolo polare iniziale (rad), 0 = verso il basso
theta_dot0 = 0.0 # velocità iniziale in theta
phi0 = 0.0 # angolo azimutale iniziale (rad)
phi_dot0 = 2.2 # velocità angolare azimutale iniziale (rad/s)
# ---------------- Equazioni del pendolo sferico ----------------
# y = [theta, theta_dot, phi, phi_dot]
# theta'' = sin(theta)cos(theta)*phi_dot^2 - (g/l)*sin(theta)
# phi'' = -2*(theta_dot*phi_dot)*cot(theta) = -2*theta_dot*phi_dot*(cos(theta)/sin(theta))
def deriv(y):
theta, th_dot, phi, ph_dot = y
sin_theta = np.sin(theta)
cos_theta = np.cos(theta)
# protezione per theta ~ 0 evitando divisione per zero
if abs(sin_theta) < 1e-8:
ph_ddot = 0.0
else:
ph_ddot = -2.0 * th_dot * ph_dot * (cos_theta / sin_theta)
th_ddot = sin_theta * cos_theta * ph_dot**2 - (g / l) * sin_theta
return np.array([th_dot, th_ddot, ph_dot, ph_ddot])
# ---------------- Integrazione RK4 ----------------
y = np.array([theta0, theta_dot0, phi0, phi_dot0], dtype=float)
traj = np.zeros((len(t_vals), 4))
for i, t in enumerate(t_vals):
traj[i] = y
k1 = deriv(y)
k2 = deriv(y + 0.5 * dt * k1)
k3 = deriv(y + 0.5 * dt * k2)
k4 = deriv(y + dt * k3)
y = y + (dt / 6.0) * (k1 + 2*k2 + 2*k3 + k4)
theta_vals = traj[:,0]
phi_vals = traj[:,2]
# ---------------- Coordinate cartesiane ----------------
x = l * np.sin(theta_vals) * np.cos(phi_vals)
y = l * np.sin(theta_vals) * np.sin(phi_vals)
z = -l * np.cos(theta_vals) # z = -l quando theta=0 (massa sotto il pivot)
# ---------------- Preparazione grafica 3D ----------------
fig = plt.figure(figsize=(6,6))
ax = fig.add_subplot(111, projection='3d')
ax.set_box_aspect([1,1,1])
lim = l + 0.2
ax.set_xlim3d(-lim, lim)
ax.set_ylim3d(-lim, lim)
ax.set_zlim3d(-lim, 0.5)
# elementi grafici: filo+massa e traiettoria
line, = ax.plot([], [], [], 'o-', lw=2)
trail, = ax.plot([], [], [], lw=1)
pivot_point, = ax.plot([0], [0], [0], 'o')
ax.set_xlabel('x')
ax.set_ylabel('y')
ax.set_zlabel('z')
ax.view_init(elev=25, azim=40)
trail_len = int(1.0 / dt) # quanto della traiettoria mostrare (in secondi)
def update(i):
xi = x[i]; yi = y[i]; zi = z[i]
# filo dalla sorgente (0,0,0) alla massa
line.set_data([0, xi], [0, yi])
line.set_3d_properties([0, zi])
start = max(0, i - trail_len)
trail.set_data(x[start:i+1], y[start:i+1])
trail.set_3d_properties(z[start:i+1])
return line, trail, pivot_point
anim = FuncAnimation(fig, update, frames=len(t_vals), interval=40)
# ---------------- Salvataggio GIF (PillowWriter, niente JS/ffmpeg necessario) ----------------
gif_name = "pendolo_sferico_3d.gif"
writer = PillowWriter(fps=int(round(1/dt)))
anim.save(gif_name, writer=writer)
plt.close(fig) # evita doppia figura
# Mostra la GIF nel notebook
display(Image(filename=gif_name))
print(f"GIF creata: {gif_name}")
<IPython.core.display.Image object>
GIF creata: pendolo_sferico_3d.gif
Spherical Pendulum with Routhian-Reduction¶
In [30]:
import sympy as sp
# Define symbolic variables
q_k, q_j, p_j , T, V, m, l, g, theta, phi, Omega, t = sp.symbols('q_k q_j p_j T V m l g theta phi Omega t')
dot_q_k, dot_q_j = sp.symbols('dot_q_k dot_q_j') # Lagrangian function
L = T - V
# Example Lagrangian for a simple pendulum with rotating support
T_example = (1/2) * m * l**2 * dot_q_k**2 + (1/2) * m * l**2 * Omega**2 * sp.sin(theta)**2
V_example = -m * g * l * sp.cos(theta)
L_example = T_example - V_example # Calculate momentum for cyclic coordinate
p_phi = sp.diff(L_example, dot_q_j)
# Define Routhian
R = L_example - p_phi * dot_q_j
# Derive equations of motion from Routhian
dR_ddotqk = sp.diff(R, dot_q_k)
dR_dqk = sp.diff(R, q_k)
# Time derivative of partial derivative with respect to dot_q_k
d_dt_dR_ddotqk = sp.diff(dR_ddotqk, t)
# Equations of motion
equation_of_motion = sp.Eq(d_dt_dR_ddotqk, dR_dqk)
# Simplify equation of motion for demonstration
simplified_eq = equation_of_motion.subs({q_k: theta, dot_q_k:
sp.Function('dot.theta')(t)}).simplify()
# Outputs
print("Routhian:", R)
print("Equation of Motion:", equation_of_motion)
print("Simplified Equation:", simplified_eq)
# Example function for Routhian reduction algorithm
def routhian_reduction_algorithm(L, cyclic_coordinates):
"""
Simplify the system by identifying cyclic coordinates and using
the Routhian.
:param L: Symbolic Lagrangian.
:param cyclic_coordinates: List of cyclic coordinates.
:return:
Reduced equations symbolic.
"""
R = L
equations = []
for q_cyclic in cyclic_coordinates:
p_cyclic = sp.diff(L, sp.Function(f'dot_{q_cyclic}')(t))
R -= p_cyclic * sp.Function(f'dot_{q_cyclic}')(t)
for q in L.free_symbols - set(cyclic_coordinates):
dR_dq = sp.diff(R, q)
d_dt_dR_ddotq = sp.diff(sp.diff(R, sp.Function(f'dot_{q}')(t)), t)
equation = sp.Eq(d_dt_dR_ddotq, dR_dq)
equations.append(equation)
return equations
# Example usage
cyclic_coords = [phi]
print(cyclic_coords)
reduced_eqs = routhian_reduction_algorithm(L_example, cyclic_coords)
for eq in reduced_eqs:
print("Reduced Equation:", eq)
a=1
Routhian: 0.5*Omega**2*l**2*m*sin(theta)**2 + 0.5*dot_q_k**2*l**2*m + g*l*m*cos(theta) Equation of Motion: True Simplified Equation: True [phi] Reduced Equation: Eq(0, 1.0*Omega**2*l**2*m*sin(theta)*cos(theta) - g*l*m*sin(theta)) Reduced Equation: Eq(0, 1.0*Omega**2*l*m*sin(theta)**2 + 1.0*dot_q_k**2*l*m + g*m*cos(theta)) Reduced Equation: Eq(0, 0.5*Omega**2*l**2*sin(theta)**2 + 0.5*dot_q_k**2*l**2 + g*l*cos(theta)) Reduced Equation: Eq(0, 1.0*Omega*l**2*m*sin(theta)**2) Reduced Equation: Eq(0, l*m*cos(theta)) Reduced Equation: Eq(0, 1.0*dot_q_k*l**2*m)
Pendolo sferico con Routhian-Reduction - esempio numerico¶
In [40]:
# Pendolo sferico basato su riduzione di Routh
import numpy as np
import matplotlib.pyplot as plt
from matplotlib.animation import FuncAnimation, PillowWriter
from IPython.display import Image, display
# -----------------------
# Parametri fisici
# -----------------------
m = 1.0
l = 1.2
g = 9.81
# Momento angolare p_phi (costante perché φ è ciclica)
p_phi = 0.9 # puoi cambiarlo liberamente
# -----------------------
# Equazioni ridotte
# -----------------------
def theta_ddot(theta, theta_dot):
sin_t = np.sin(theta)
cos_t = np.cos(theta)
return (p_phi**2 / (m**2 * l**4)) * (cos_t / sin_t**3) - (g/l)*sin_t
def phi_dot(theta):
return p_phi / (m * l**2 * np.sin(theta)**2)
# -----------------------
# Integrazione RK4
# -----------------------
dt = 0.03
t_max = 7.0
t_vals = np.arange(0, t_max, dt)
theta0 = 0.7
theta_dot0 = 0.0
phi0 = 0.0
theta = theta0
theta_dot = theta_dot0
phi = phi0
θ_list = []
φ_list = []
for t in t_vals:
# salva
θ_list.append(theta)
φ_list.append(phi)
# RK4 su θ
k1 = theta_ddot(theta, theta_dot)
k2 = theta_ddot(theta + 0.5*dt*theta_dot, theta_dot + 0.5*dt*k1)
k3 = theta_ddot(theta + 0.5*dt*theta_dot, theta_dot + 0.5*dt*k2)
k4 = theta_ddot(theta + dt*theta_dot, theta_dot + dt*k3)
theta_dot += dt*(k1 + 2*k2 + 2*k3 + k4)/6
theta += dt * theta_dot
# aggiornamento φ usando p_phi costante
phi += dt * phi_dot(theta)
θ = np.array(θ_list)
φ = np.array(φ_list)
# -----------------------
# Coordinate cartesiane
# -----------------------
x = l * np.sin(θ) * np.cos(φ)
y = l * np.sin(θ) * np.sin(φ)
z = -l * np.cos(θ)
# -----------------------
# Animazione 3D
# -----------------------
fig = plt.figure(figsize=(6,6))
ax = fig.add_subplot(111, projection="3d")
ax.set_box_aspect([1,1,1])
lim = l + 0.3
ax.set_xlim(-lim, lim)
ax.set_ylim(-lim, lim)
ax.set_zlim(-lim, 0.4)
line, = ax.plot([], [], [], "o-", lw=2)
trail, = ax.plot([], [], [], lw=1)
pivot, = ax.plot([0], [0], [0], "o", color="black")
ax.set_xlabel("x")
ax.set_ylabel("y")
ax.set_zlabel("z")
trail_len = int(1.0/dt)
def update(i):
line.set_data([0, x[i]], [0, y[i]])
line.set_3d_properties([0, z[i]])
start = max(0, i - trail_len)
trail.set_data(x[start:i+1], y[start:i+1])
trail.set_3d_properties(z[start:i+1])
return line, trail, pivot
anim = FuncAnimation(fig, update, frames=len(t_vals), interval=40)
gif_name = "pendolo_routh.gif"
writer = PillowWriter(fps=int(1/dt))
anim.save(gif_name, writer=writer)
plt.close(fig)
display(Image(filename=gif_name))
print("GIF creata:", gif_name)
<IPython.core.display.Image object>
GIF creata: pendolo_routh.gif
In [ ]: