import sympy as sp

def poisson_bracket(A, B, q, p):
    """
    Calculate the Poisson bracket of two functions A and B.
        :param A: Function A in terms of q and p.
        :param B: Function B in terms of q and p.
        :param q: List of position variables.
        :param p: List of momentum variables.
    :return: 
        Poisson bracket {A, B}.
    """
    pb = 0
    for qi, pi in zip(q, p):
        pb += (A.diff(qi) * B.diff(pi)) - (A.diff(pi) * B.diff(qi))
    return pb

def verify_bilinearity(A, B, C, q, p):
    """
    Verify the bilinearity property of Poisson brackets.
        :param A, B, C: Functions A, B, C.
        :param q, p: Lists of canonical variables.
    """
    a, b = sp.symbols('a b')
    lhs = poisson_bracket(a*A + b*B, C, q, p)
    rhs = a*poisson_bracket(A, C, q, p) + b * poisson_bracket(B, C, q, p)
    assert sp.simplify(lhs - rhs) == 0, "Bilinearity does not hold"

def verify_antisymmetry(A, B, q, p):
    """
    Verify the antisymmetry property of Poisson brackets.
        :param A, B: Functions A, B.
        :param q, p: Lists of canonical variables.
    """
    assert poisson_bracket(A, B, q, p) == -poisson_bracket(B, A, q, p) , "Antisymmetry does not hold" 
    
def verify_jacobi_identity(A, B, C, q, p):
    """
    Verify the Jacobi identity for Poisson brackets.
        :param A, B, C: Functions A, B, C.
        :param q, p: Lists of canonical variables.
    """
    lhs = poisson_bracket(A, poisson_bracket(B, C, q, p), q, p) + \
          poisson_bracket(B, poisson_bracket(C, A, q, p), q, p) + \
          poisson_bracket(C, poisson_bracket(A, B, q, p), q, p)
    assert sp.simplify(lhs) == 0, "Jacobi identity does not hold" 

def hamiltons_equations(H, q, p):
    """
    Derive Hamilton's equations of motion using Poisson brackets for
    a given Hamiltonian.
        :param H: Hamiltonian of the system.
        :param q, p: Lists of canonical variables.
        :retum: Equations of motion.
    """
    dq_dt = [poisson_bracket(qi, H, q, p) for qi in q]
    dp_dt = [poisson_bracket(pi, H, q, p) for pi in p]
    return dq_dt, dp_dt

# Define symbols
q1, p1 = sp.symbols('q1 p1')
q = [q1]
p = [p1]
# Define functions
A = q1**2 * p1
B = sp.sin(q1 + p1)
C = sp.exp(-q1*p1)
H = (p1**2)/(2*1) + (1/2)*1*1**2*q1**2 # Harmonic oscillator Hamiltonian

# Verify properties
verify_bilinearity(A, B, C, q, p)
verify_antisymmetry(A, B, q, p)
verify_jacobi_identity(A, B, C, q, p)

# Derive Hamilton's equations
dq_dt, dp_dt = hamiltons_equations(H, q, p)

# Example outputs
print("Hamilton s equations:")
print(f"dq/dt: {dq_dt}")
print(f"dp/dt: {dp_dt}")

Hamilton s equations:
dq/dt: [p1]
dp/dt: [-1.0*q1]

Poisson Brackets¶

Date tre funzioni A(q, p), B(q, p), C(q, p), calcolare le loro Poisson Bracket - fare verifica di simpletticità¶

Calcolare Hamiltoniana di oscillatore armonico¶