Equazioni del moto di Hamilton

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

Due oscillatori armonici indipendenti - Approccio simbolico

import math

# ============================
#   HARMONIC OSCILLATOR
# ============================

def harmonic_hamiltonian(q, p, m, omega):
    """
    Hamiltonian function for a simple harmonic oscillator.
        :param q: Generalized coordinate.
        :param p: Generalized momentum.
        :param m: Mass of the oscillator.
        :param omega: Angular frequency.
    :retum: 
        Value of the Hamiltonian.
    """
    return p**2 / (2*m) + 0.5 * m * omega**2 * q**2

def harmonic_equations(q, p, m, omega):
    """
    Derive Hamilton's equations of motion.
        :param q: Generalized coordinate.
        :param p: Generalized momentum.
        :param m: Mass of the oscillator.
        :param omega: Angular frequency.
    :return: 
        dq/dt and dp/dt.
    """
    dq = p / m
    dp = -m * omega**2 * q
    return dq, dp


# ============================
#           PENDULUM
# ============================

def pendulum_hamiltonian(q, p, m, l, g=9.81):
    """
    Hamiltonian for a simple pendulum.
        :param q: Generalized coordinate (angle).
        :param p: Generalized momentum.
        :param m: Mass of pendulum bob.
        :param I: Length of the pendulum.
    :return: 
        Value of the Hamiltonian.
    """
    return p**2 / (2 * m * l**2) + m * g * l * (1 - math.cos(q))

def pendulum_equations(q, p, m, l, g=9.81):
    """
    Derive Hamilton's equations of motion for a pendulum.
        :param q: Generalized coordinate (angle).
        :param p: Generalized momentum.
        :param m: Mass of pendulum bob.
        :param I: Length of the pendulum.
    :return: 
        dq/dt and dp/dt.
    """
    dq = p / (m * l**2)
    dp = -m * g * l * math.sin(q)
    return dq, dp


# ============================
#       GENERIC LEAPFROG
# ============================

def leapfrog_integrator(q0, p0, delta_t, num_steps, equations, args=()):
    """
    Generic leapfrog integrator for Hamiltonian systems.
    equations = function(q, p, *args) -> (dq, dp)
    
    Symplectic integrator using the leapfrog method for simple
    harmonic oscillator.
        :param qO: Initial coordinate.
        :param pO: Initial momentum.
        :param m: Mass of the oscillator.
        :param omega: Angular frequency.
        :param delta_t: Time step.
        :param num_steps: Number of integration steps.
    :return: 
        Arrays of q and p over time.
    """
    q = q0
    p = p0
    qs = [q]
    ps = [p]

    for _ in range(num_steps):
        # Evaluate derivatives at (q,p)
        dq, dp = equations(q, p, *args)

        # Half momentum step
        p_half = p + 0.5 * delta_t * dp

        # Full position step using p_half
        dq_half, _ = equations(q, p_half, *args)
        q = q + delta_t * dq_half

        # Evaluate dp at new position
        _, dp_new = equations(q, p_half, *args)

        # Second half momentum step
        p = p_half + 0.5 * delta_t * dp_new

        qs.append(q)
        ps.append(p)

    return qs, ps


# ============================
#       EXAMPLE USAGE
# ============================

if __name__ == "__main__":
    m = 1.0
    omega = 1.0
    q0, p0 = 1.0, 0.0
    delta_t = 0.1
    num_steps = 100

    # Harmonic oscillator
    qs_ho, ps_ho = leapfrog_integrator(
        q0, p0, delta_t, num_steps,
        harmonic_equations,
        args=(m, omega)
    )

    # Pendulum
    l = 1.0
    qs_p, ps_p = leapfrog_integrator(
        q0, p0, delta_t, num_steps,
        pendulum_equations,
        args=(m, l)
    )

    print("=== Harmonic Oscillator ===")
    print("Final Position:", qs_ho[-1])
    print("Final Momentum:", ps_ho[-1])

    print("\n=== Pendulum ===")
    print("Final Position:", qs_p[-1])
    print("Final Momentum:", ps_p[-1])

=== Harmonic Oscillator ===
Final Position: -0.8367949271103875
Final Momentum: 0.5468316142446552

=== Pendulum ===
Final Position: -0.3660020882019804
Final Momentum: 2.7487333889022048

Perché harmonic_equations(q, p, m, omega) funziona dentro leapfrog_integrator?

La funzione leapfrog_integrator accetta:

equations → la funzione che calcola le equazioni di Hamilton

args → una tupla di parametri aggiuntivi (ad esempio (m, omega))

La parte importante è questa:

dq, dp = equations(q, p, *args)