Symmetry Reduction¶
da PowerShell -- jupyter nbconvert --to html notebook.ipynb
In [8]:
import numpy as np
def noether_conserved_quantity(L, q, q_dot, eta):
"""
Calculate the conserved quantity from Noether's theorem.
:param L: Lagrangian function.
:param q: Generalized coordinates.
:param q_dot: Generalized velocities.
:param eta: Symmetry transformation function.
:return:
Conserved quantity Q.
"""
Q = np.sum([np.dot(np.gradient(L, q_dot_i), eta_i(q)) for q_dot_i, eta_i in zip(q_dot, eta)])
return Q
def hamiltonian_symmetry_reduction(H, q, p, G):
"""
Perform symmetry reduction on a Hamiltonian system.
:param H: Hamiltonian function.
:param q: Generalized coordinates.
:param p: Conjugate momenta.
:param G: Lie group representing the symmetries.
:return:
Reduced Hamiltonian.
"""
# Placeholder for transformation (syrrplectomorphism)
Q, P = G.act_on_space(q, p)
# Calculate reduced Hamiltonian
H_reduced = H(Q, P)
return H_reduced
def poisson_bracket(f, g, J):
'''
Compute the Poisson bracket of two functions f and g.
:param f: Function f in phase space.
:param g: Function g in phase space.
:param J: Symplectic matrix.
:return:
Poisson bracket {f, g}.
'''
grad_f = np.gradient(f)
grad_g = np.gradient(g)
return np.dot(grad_f, np.dot(J, grad_g))
def central_force_reduction(H, r, L, mu):
'''
Reduce Hamiltonian for central force problem using conserved
angular momentum.
:param H: Hamiltonian function.
:param r: Radial coordinate.
:param L: Angular momentum.
:param mu: Reduced mass.
:return:
Effective one-dimensional Hamiltonian.
'''
V_eff = lambda r: H(r) + L**2 / (2 * mu * r**2)
return V_eff
def symmetry_based_reduction_algorithm(H, initial_conditions, G):
'''
Perform symmetry reduction on a dynamical system.
:param H: Hamiltonian function.
:param initial_conditions: Initial state of the system.
:param G: Symmetry group acting on the system.
:return:
Reduced dynamics.
'''
# Identify invariant manifolds and conserved quantities
invariant_manifolds = G.invariant_manifolds(initial_conditions)
conserved_quantities = [Noether_conserved_quantity(H, q, p, eta) for q, p, eta in invariant_manifolds]
# Perform reduction
reduced_trajectory = []
for quantities in conserved_quantities:
reduced_state = hamiltonian_symmetry_reduction(H, quantities['q'] , quantities ['p'] , G)
reduced_trajectory.append(reduced_state)
return reduced_trajectory
# Demonstration
L = lambda q, q_dot: np.sum(q_dot**2) / 2 - np.sum(q**2) / 2
H = lambda q, p: np.sum(p**2) / 2 + np.sum(q**2) / 2
J = np.array([[0, -1], [1, 0]])
eta = [lambda q: q]
G = type('LieGroup', (object,), {'act_on_space': lambda self, q, p: (q, p)})()
q = np.array([1.0, 0.0])
q_dot = np.array([0.0, 1.0])
p = np.array([0.0, 1.0])
r = 1.0
L_val =1.0
mu = 1.0 # Compute conserved quantity, Poisson bracket, and reduced dynamics
'''
Q_conserved = noether_conserved_quantity(L, q, q_dot, eta)
poisson_val = poisson_bracket(H(q, p), L(q, q_dot), J)
reduced_H = central_force_reduction(H, r, L_val, mu)
reduced_trajectory = symmetry_based_reduction_algorithm(H, q, G)
print("Conserved Quantity:", Q_conserved)
print("Poisson Bracket:", poisson_val)
print("Reduced Hamiltonian:", reduced_H(r))
print("Reduced Trajectory:", reduced_trajectory)
'''
a=1
In [12]:
import numpy as np
# -----------------------------
# Lagrangiana e Hamiltoniana
# -----------------------------
def L(q, q_dot):
return 0.5 * np.dot(q_dot, q_dot) - 0.5 * np.dot(q, q)
def H(q, p):
return 0.5 * np.dot(p, p) + 0.5 * np.dot(q, q)
# -----------------------------
# Noether (simmetria lineare)
# -----------------------------
def noether_charge(q, q_dot, eta):
"""
Q = p · eta(q), con p = ∂L/∂q_dot
"""
p = q_dot
return np.dot(p, eta(q))
# -----------------------------
# Poisson bracket
# -----------------------------
def poisson_bracket(f_grad, g_grad, J):
"""
{f, g} = ∇f^T J ∇g
"""
return f_grad @ J @ g_grad
# -----------------------------
# Riduzione per forza centrale
# -----------------------------
def effective_radial_hamiltonian(r, p_r, L, mu, V):
"""
H_eff = p_r^2/(2μ) + L^2/(2μ r^2) + V(r)
"""
return p_r**2 / (2 * mu) + L**2 / (2 * mu * r**2) + V(r)
# -----------------------------
# ESEMPIO CONCRETO
# -----------------------------
q = np.array([1.0, 0.0])
q_dot = np.array([0.0, 1.0])
p = q_dot
# Simmetria di rotazione infinitesima
eta = lambda q: np.array([-q[1], q[0]])
# Matrice simplettica standard
J = np.block([
[np.zeros((2,2)), np.eye(2)],
[-np.eye(2), np.zeros((2,2))]
])
print(J)
# Noether
Q = noether_charge(q, q_dot, eta)
# Poisson bracket {H, L}
grad_H = np.concatenate((q, p))
grad_L = np.concatenate((eta(q), eta(p)))
PB = poisson_bracket(grad_H, grad_L, J)
# Riduzione centrale
V = lambda r: 0.5 * r**2
H_eff = effective_radial_hamiltonian(r=1.0, p_r=0.2, L=Q, mu=1.0, V=V)
print("Noether charge (angular momentum):", Q)
print("Poisson bracket {H, L}:", PB)
print("Effective radial Hamiltonian:", H_eff)
[[ 0. 0. 1. 0.]
[ 0. 0. 0. 1.]
[-1. -0. 0. 0.]
[-0. -1. 0. 0.]]
Noether charge (angular momentum): 1.0
Poisson bracket {H, L}: -2.0
Effective radial Hamiltonian: 1.02
In [7]:
import numpy as np
def lagrangian(T, V):
"""
Calculate the Lagrangian of a system.
:param T: Kinetic energy.
:param V: Potential energy.
:return:
Lagrangian.
"""
return T - V
def rayleigh_dissipation(velocities, coefficients):
"""
Calculate the Rayleigh dissipation function.
:param velocities: Generalized velocities.
:param coefficients: Damping coefficients.
:retum:
Rayleigh dissipation.
"""
return 0.5 * np.sum(coefficients * velocities**2)
def harmonic_oscillator(m, k, gamma, x, v, dt):
"""
Simulates a damped harmonic oscillator using the velocity Verlet method.
:param m: Mass of the oscillator.
:param k: Spring constant.
:param gamma: Damping coefficient.
:param x: Initial position.
:param v: Initial velocity.
:param dt: Time step.
:return:
Simulated position and velocity over time.
"""
# Initialize arrays to store position and velocity
num_steps = 1000
positions = np.zeros(num_steps)
velocities = np.zeros(num_steps)
positions[0] = x
velocities [0] = v
for n in range(1, num_steps):
# Calculate acceleration from force terms
a = -(k / m) * positions[n-1] - (gamma / m) * velocities [n-1]
# Update position
positions[n] = positions[n-1] + velocities[n-1] * dt + 0.5 * a * dt**2
# Calculate new acceleration
a_new = -(k / m) * positions[n] - (gamma / m) * velocities[n-1]
# Update velocity
velocities[n] = velocities[n-1] + 0.5 * (a + a_new) * dt
return positions, velocities
# Simulation parameters
mass = 1.0
stiffness = 10.0
damping = 0.5
initial_position = 1.0
initial_velocity = 0.0
time_step = 0.1 # Run simulation
pos, vel = harmonic_oscillator(mass, stiffness, damping, initial_position, initial_velocity, time_step)
print("Final position:", pos[-1])
print("Final velocity:", vel[-1])
Final position: 1.4012065282144866e-11 Final velocity: 7.488868522582272e-12
In [ ]: