Liouville Equation

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

Un esempio per studiare l’evoluzione temporale generata da una Hamiltoniana

tramite equazioni di Hamilton (flusso di Liouville).

Scelgo come esempio l’oscillatore armonico 1D, perché:

• è semplice ma non banale
• conserva il volume di fase (Liouville)

è ideale per visualizzare traiettorie in spazio delle fasi

$$\large H(q,p)=\frac{p^2}{2m}+\frac{1}{2}m \omega^2q^2$$

Oscillatore armonico 1D

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

def Hamiltonian(q, p, m=1.0, omega=1.0):
    """
    Hamiltoniana dell'oscillatore armonico 1D
    """
    return p**2 / (2*m) + 0.5 * m * omega**2 * q**2

def hamiltons_equations(state, t, m=1.0, omega=1.0):
    """
    Equazioni di Hamilton:
    state = [q, p]
    """
    q, p = state

    dqdt = p / m
    dpdt = -m * omega**2 * q

    return [dqdt, dpdt]

def liouville_flow(initial_conditions, t_span, m=1.0, omega=1.0):
    """
    Calcola le traiettorie nello spazio delle fasi
    """
    trajectories = odeint(
        hamiltons_equations,
        initial_conditions,
        t_span,
        args=(m, omega)
    )
    return trajectories

# Condizioni iniziali
q0 = 1.0     # posizione iniziale
p0 = 0.0     # impulso iniziale
initial_conditions = [q0, p0]

# Intervallo temporale
t_span = np.linspace(0, 20, 2000)

traj = liouville_flow(initial_conditions, t_span)

q = traj[:, 0]
p = traj[:, 1]

H = Hamiltonian(q, p)
plt.figure()
plt.plot(t_span, H)
plt.xlabel("t")
plt.ylabel("H")
plt.ylim(0.4, 0.6)
plt.title("Conservazione dell'energia")
plt.grid()
plt.show()

print(H.min(), H.max())

0.4999999700647026 0.5000003305453562

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