# Copyright 2026 Philippe Billet assisted by LLMs in free mode: chatGPT, Qwen, Deepseek, Gemini, Claude, le chat Mistral.
#
# Licensed under the Apache License, Version 2.0 (the "License");
# you may not use this file except in compliance with the License.
# You may obtain a copy of the License at
#
# http://www.apache.org/licenses/LICENSE-2.0
#
# Unless required by applicable law or agreed to in writing, software
# distributed under the License is distributed on an "AS IS" BASIS,
# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
# See the License for the specific language governing permissions and
# limitations under the License..
"""
cotangent_bundle.py — Symplectic geometry of cotangent bundles: from canonical structure to quantization
=====================================================================================================================
Overview
--------
This module provides a symbolic‑numerical laboratory for exploring the
symplectic geometry of cotangent bundles **before** a Hamiltonian is chosen.
It demonstrates how the canonical symplectic form ``ω = dξ ∧ dx`` emerges
directly from the geometry of T*M, and how it alone determines the Poisson
bracket, Hamilton’s equations, and ultimately the link to quantum mechanics.
The package implements:
* Construction of the canonical Liouville 1‑form ``θ`` and the symplectic
form ``ω = dθ``.
* Verification that ``ω`` is closed (``dω = 0``) and non‑degenerate
(``det ω ≠ 0``).
* Derivation of Hamiltonian vector fields ``X_H`` from any symbolic
Hamiltonian ``H(x,ξ)`` via the intrinsic relation ``ι_{X_H} ω = dH``.
* Numerical integration of the Hamiltonian flow with preservation checks
(energy, symplectic area, Jacobian determinant).
* Symbolic and numerical study of Poisson brackets, Jacobi identity, and
canonical commutation relations.
* Visualisation of phase‑space trajectories, energy surfaces, and the
symplectic area preserved by the flow.
* Canonical transformations and their characterisation (preservation of ω).
* A gentle introduction to the Dirac quantisation rule:
``{f,g} ↦ (1/iℏ)[f̂,ĝ]``, illustrated with pseudo‑differential operators
(through an external `psiop` module).
All computations are performed symbolically with **SymPy** and then
lambdified for efficient numerical integration and visualisation with
**NumPy/SciPy** and **Matplotlib**.
Mathematical background
-----------------------
Let `M` be an `n`‑dimensional manifold (here taken as ℝⁿ for simplicity).
The cotangent bundle `T*M` carries a **canonical 1‑form** (Liouville form)
defined in local coordinates ``(x¹,…,xⁿ, ξ₁,…,ξₙ)`` by
θ = ξ₁ dx¹ + … + ξₙ dxⁿ.
Its exterior derivative
ω = dθ = dξ₁ ∧ dx¹ + … + dξₙ ∧ dxⁿ
is a closed (``dω = 0``) and non‑degenerate 2‑form – the **symplectic form**.
Crucially, ω is **independent of any Hamiltonian**; it is a geometric
property of the bundle itself.
Given a smooth function `H : T*M → ℝ` (the Hamiltonian), the Hamiltonian
vector field `X_H` is uniquely defined by
ι_{X_H} ω = dH.
In coordinates this yields the familiar Hamilton equations
˙xⁱ = ∂H/∂ξᵢ , ˙ξᵢ = −∂H/∂xⁱ.
The flow `Φ_t` generated by `X_H` preserves ω (Liouville’s theorem); in
particular the symplectic volume (area in 1D) is conserved. The Poisson
bracket of two observables `f,g` is defined by
{f,g} = ω(X_f, X_g) = Σᵢ (∂f/∂xⁱ ∂g/∂ξᵢ − ∂f/∂ξᵢ ∂g/∂xⁱ),
and satisfies the Jacobi identity. The canonical brackets
{xⁱ, ξⱼ} = δⁱⱼ , {xⁱ, xʲ} = {ξᵢ, ξⱼ} = 0
are direct consequences of ω.
**Quantisation.** The Dirac correspondence replaces Poisson brackets by
commutators:
[f̂,ĝ] = iℏ {f,g} .
For the fundamental pair this gives the Heisenberg relation
``[x̂ⁱ, p̂ⱼ] = iℏ δⁱⱼ``. The module illustrates this by quantising a
classical Hamiltonian into a pseudo‑differential operator (using a
companion module `psiop`) and verifying the correspondence.
**Canonical transformations.** A change of coordinates `(x,ξ) → (X,Ξ)` is
canonical iff it preserves ω, i.e. `dΞ ∧ dX = dξ ∧ dx`. Equivalently, the
Jacobian matrix `J` satisfies `Jᵀ ω J = ω`. Several examples are provided,
together with generating functions.
**Higher dimensions.** The code handles arbitrary dimension by building the
block‑diagonal symplectic matrix
ω = [[0, -I], [I, 0]] .
References
----------
.. [1] Abraham, R. & Marsden, J. E. *Foundations of Mechanics*,
Addison‑Wesley, 1978 (2nd ed.). Chapter 3: Hamiltonian Systems.
.. [2] Arnold, V. I. *Mathematical Methods of Classical Mechanics*,
Springer‑Verlag, 1989 (2nd ed.). §37: Symplectic Structure.
.. [3] Cannas da Silva, A. *Lectures on Symplectic Geometry*,
Springer‑Verlag, 2001. Chapter 1: Symplectic Forms.
.. [4] Dirac, P. A. M. *The Principles of Quantum Mechanics*,
Oxford University Press, 1930 (4th ed.). §21: Poisson Brackets.
.. [5] McDuff, D. & Salamon, D. *Introduction to Symplectic Topology*,
Oxford University Press, 2017 (3rd ed.). Chapter 1: Symplectic
Manifolds.
.. [6] Weyl, H. *The Theory of Groups and Quantum Mechanics*,
Dover Publications, 1931 (translated from German). §II.7:
Quantisation and the Symplectic Form.
"""
from sympy import symbols, Matrix, simplify
from sympy.diffgeom import Manifold, Patch, CoordSystem
import numpy as np
[docs]
def verify_symplectic_structure(dim=1):
"""Verify that ω is closed and non-degenerate"""
if dim == 1:
x, xi = symbols('x xi', real=True)
# Symplectic form ω = dξ ∧ dx
# In coordinates: ω = [[0, -1], [1, 0]]
omega = Matrix([[0, -1], [1, 0]])
print("Symplectic form ω:")
print(omega)
# Non-degeneracy: det(ω) ≠ 0
det_omega = omega.det()
print(f"\ndet(ω) = {det_omega}")
assert det_omega != 0, "ω is degenerate!"
# Antisymmetry: ω = -ωᵀ
print(f"\nω + ωᵀ = {simplify(omega + omega.T)}")
assert (omega + omega.T).is_zero_matrix
elif dim == 2:
x, y, xi, eta = symbols('x y xi eta', real=True)
# ω = dξ∧dx + dη∧dy
# Matrix: [[0, 0, -1, 0], [0, 0, 0, -1], [1, 0, 0, 0], [0, 1, 0, 0]]
omega = Matrix([
[0, 0, -1, 0],
[0, 0, 0, -1],
[1, 0, 0, 0],
[0, 1, 0, 0]
])
print("2D symplectic form ω:")
print(omega)
print(f"\ndet(ω) = {omega.det()}")
return omega
# Test
omega = verify_symplectic_structure(dim=1)
[docs]
def hamiltonian_vector_field(H_expr, vars_x):
"""Computes the Hamiltonian vector field X_H from an arbitrary H"""
from sympy import symbols, diff
dim = len(vars_x)
if dim == 1:
x, = vars_x
xi = symbols('xi', real=True)
# Hamilton's equations: ẋ = ∂H/∂ξ, ξ̇ = -∂H/∂x
dx_dt = diff(H_expr, xi)
dxi_dt = -diff(H_expr, x)
print("Hamiltonian vector field X_H:")
print(f" dx/dt = ∂H/∂ξ = {dx_dt}")
print(f" dξ/dt = -∂H/∂x = {dxi_dt}")
return {'dx/dt': dx_dt, 'dxi/dt': dxi_dt}
elif dim == 2:
x, y = vars_x
xi, eta = symbols('xi eta', real=True)
# X_H = (∂H/∂ξ, ∂H/∂η, -∂H/∂x, -∂H/∂y)
dx_dt = diff(H_expr, xi)
dy_dt = diff(H_expr, eta)
dxi_dt = -diff(H_expr, x)
deta_dt = -diff(H_expr, y)
print("2D Hamiltonian vector field X_H:")
print(f" dx/dt = {dx_dt}")
print(f" dy/dt = {dy_dt}")
print(f" dξ/dt = {dxi_dt}")
print(f" dη/dt = {deta_dt}")
return {
'dx/dt': dx_dt,
'dy/dt': dy_dt,
'dxi/dt': dxi_dt,
'deta/dt': deta_dt
}
# Example: arbitrary H
from sympy import symbols, exp, sin, cos
x, xi = symbols('x xi', real=True)
# Hamiltonians of increasing complexity
hamiltonians = {
"Free": xi**2 / 2,
"Harmonic": xi**2/2 + x**2/2,
"Pendulum": xi**2/2 - cos(x),
"Non-polynomial": xi**2/2 + x**4/4 + sin(xi*x),
"Transcendental": exp(xi**2) * sin(x)
}
for name, H in hamiltonians.items():
print(f"\n{'='*60}")
print(f"Hamiltonian: {name}")
print(f"H(x,ξ) = {H}")
print('='*60)
X_H = hamiltonian_vector_field(H, [x])
[docs]
def verify_symplectic_preservation(H_expr, vars_x, t_max=1.0, n_times=10):
"""
Simplified and more robust version
"""
from scipy.integrate import solve_ivp
from sympy import lambdify, symbols
import numpy as np
print("\n" + "="*70)
print("VERIFICATION OF SYMPLECTIC PRESERVATION")
print("="*70)
dim = len(vars_x)
if dim == 1:
x, = vars_x
xi = symbols('xi', real=True)
# Hamiltonian system
X_H = hamiltonian_vector_field(H_expr, vars_x)
f_x = lambdify((x, xi), X_H['dx/dt'], 'numpy')
f_xi = lambdify((x, xi), X_H['dxi/dt'], 'numpy')
def system(t, y):
return [f_x(y[0], y[1]), f_xi(y[0], y[1])]
# Initial condition
y0 = np.array([1.0, 0.5])
print(f"\nInitial condition: (x₀, ξ₀) = ({y0[0]}, {y0[1]})")
print(f"Hamiltonian: H = {H_expr}")
# Integration
t_eval = np.linspace(0, t_max, 1000)
sol = solve_ivp(system, [0, t_max], y0, t_eval=t_eval,
method='DOP853', rtol=1e-10, atol=1e-12)
# 1. Verify conservation of H
H_func = lambdify((x, xi), H_expr, 'numpy')
H_values = H_func(sol.y[0], sol.y[1])
print(f"\n[1] Conservation of energy H:")
print(f" H(t=0) = {H_values[0]:.12f}")
print(f" H(t={t_max:.1f}) = {H_values[-1]:.12f}")
print(f" |ΔH|/H₀ = {abs(H_values[-1]-H_values[0])/abs(H_values[0]):.2e}")
# 2. Verify det(Jacobian) = 1 via variational method
epsilon = 1e-7
times_to_check = np.linspace(0, t_max, n_times)
print(f"\n[2] Conservation of symplectic structure (det J = 1):")
print(f" Method: finite differences on the Jacobian")
max_det_error = 0.0
for t_check in times_to_check:
# Integrate 4 perturbed trajectories
sol_ref = solve_ivp(system, [0, t_check], y0, method='DOP853', rtol=1e-10)
sol_dx = solve_ivp(system, [0, t_check], y0 + [epsilon, 0],
method='DOP853', rtol=1e-10)
sol_mx = solve_ivp(system, [0, t_check], y0 + [-epsilon, 0],
method='DOP853', rtol=1e-10)
sol_dxi = solve_ivp(system, [0, t_check], y0 + [0, epsilon],
method='DOP853', rtol=1e-10)
sol_mxi = solve_ivp(system, [0, t_check], y0 + [0, -epsilon],
method='DOP853', rtol=1e-10)
# Final states
final_ref = sol_ref.y[:, -1]
final_dx = sol_dx.y[:, -1]
final_mx = sol_mx.y[:, -1]
final_dxi = sol_dxi.y[:, -1]
final_mxi = sol_mxi.y[:, -1]
# Jacobian
J = np.array([
[(final_dx[0] - final_mx[0])/(2*epsilon),
(final_dxi[0] - final_mxi[0])/(2*epsilon)],
[(final_dx[1] - final_mx[1])/(2*epsilon),
(final_dxi[1] - final_mxi[1])/(2*epsilon)]
])
det_J = np.linalg.det(J)
det_error = abs(det_J - 1.0)
max_det_error = max(max_det_error, det_error)
if t_check == 0 or t_check >= t_max - 0.01 or det_error > 1e-6:
print(f" t={t_check:6.2f}: det(J) = {det_J:.12f} (|error| = {det_error:.2e})")
print(f"\n Maximum error: {max_det_error:.2e}")
# 3. Verify via symplectic form ω
print(f"\n[3] Direct verification: Jᵀ ω J = ω")
omega = np.array([[0, -1], [1, 0]])
omega_transformed = J.T @ omega @ J
omega_error = np.linalg.norm(omega_transformed - omega)
print(f" ω =")
print(f" {omega}")
print(f"\n Jᵀ ω J =")
print(f" {omega_transformed}")
print(f"\n ||Jᵀ ω J - ω|| = {omega_error:.2e}")
# Conclusion
print("\n" + "="*70)
print("CONCLUSION:")
print("="*70)
all_good = (abs(H_values[-1]-H_values[0])/abs(H_values[0]) < 1e-8 and
max_det_error < 1e-6)
if all_good:
print("✓ The Hamiltonian flow PRESERVES the symplectic structure")
print("✓ Conservation of H to numerical precision")
print("✓ det(Jacobian) = 1 verified")
print("\n➜ Confirmation of Liouville's Theorem!")
else:
print("⚠ Numerical errors detected (probably due to integration)")
print(f" - Error on H: {abs(H_values[-1]-H_values[0])/abs(H_values[0]):.2e}")
print(f" - Error on det(J): {max_det_error:.2e}")
return sol, H_values
# Test with different Hamiltonians
print("\n" + "="*70)
print("TEST 1: Harmonic oscillator H = (ξ² + x²)/2")
print("="*70)
verify_symplectic_preservation(xi**2/2 + x**2/2, [x], t_max=10.0)
print("\n" + "="*70)
print("TEST 2: Nonlinear Hamiltonian H = ξ²/2 + x⁴/4")
print("="*70)
verify_symplectic_preservation(xi**2/2 + x**4/4, [x], t_max=5.0)
[docs]
def poisson_bracket(f_expr, g_expr, vars_x):
"""Poisson bracket {f, g} induced by ω"""
from sympy import symbols, diff, simplify
dim = len(vars_x)
if dim == 1:
x, = vars_x
xi = symbols('xi', real=True)
# {f, g} = ∂f/∂x · ∂g/∂ξ - ∂f/∂ξ · ∂g/∂x
bracket = (
diff(f_expr, x) * diff(g_expr, xi) -
diff(f_expr, xi) * diff(g_expr, x)
)
elif dim == 2:
x, y = vars_x
xi, eta = symbols('xi eta', real=True)
bracket = (
diff(f_expr, x) * diff(g_expr, xi) - diff(f_expr, xi) * diff(g_expr, x) +
diff(f_expr, y) * diff(g_expr, eta) - diff(f_expr, eta) * diff(g_expr, y)
)
return simplify(bracket)
# Verify Jacobi identity {f, {g, h}} + cyclic permutations = 0
from sympy import symbols
x, xi = symbols('x xi', real=True)
f = x**2
g = xi**2
h = x*xi
bracket_fg_h = poisson_bracket(poisson_bracket(f, g, [x]), h, [x])
bracket_gh_f = poisson_bracket(poisson_bracket(g, h, [x]), f, [x])
bracket_hf_g = poisson_bracket(poisson_bracket(h, f, [x]), g, [x])
jacobi_identity = simplify(bracket_fg_h + bracket_gh_f + bracket_hf_g)
print(f"f = {f}, g = {g}, h = {h}")
print(f"\n{{f, {{g, h}}}} + {{g, {{h, f}}}} + {{h, {{f, g}}}} = {jacobi_identity}")
print(f"\nJacobi identity verified: {jacobi_identity == 0}")
[docs]
def verify_canonical_relations(vars_x):
"""Verifies the canonical commutation relations"""
from sympy import symbols
dim = len(vars_x)
if dim == 1:
x, = vars_x
xi = symbols('xi', real=True)
# Fundamental brackets
bracket_x_xi = poisson_bracket(x, xi, [x])
bracket_x_x = poisson_bracket(x, x, [x])
bracket_xi_xi = poisson_bracket(xi, xi, [x])
print("1D canonical relations:")
print(f" {{x, ξ}} = {bracket_x_xi} (must be 1)")
print(f" {{x, x}} = {bracket_x_x} (must be 0)")
print(f" {{ξ, ξ}} = {bracket_xi_xi} (must be 0)")
assert bracket_x_xi == 1
assert bracket_x_x == 0
assert bracket_xi_xi == 0
elif dim == 2:
x, y = vars_x
xi, eta = symbols('xi eta', real=True)
print("2D canonical relations:")
print(f" {{x, ξ}} = {poisson_bracket(x, xi, [x, y])}")
print(f" {{y, η}} = {poisson_bracket(y, eta, [x, y])}")
print(f" {{x, η}} = {poisson_bracket(x, eta, [x, y])}")
print(f" {{y, ξ}} = {poisson_bracket(y, xi, [x, y])}")
verify_canonical_relations([x])
[docs]
def visualize_symplectic_geometry(H_expr, vars_x, xlim=(-2, 2), klim=(-2, 2)):
"""Visualizes the symplectic geometry for an arbitrary H"""
import matplotlib.pyplot as plt
from matplotlib.patches import FancyArrowPatch
from mpl_toolkits.mplot3d import proj3d
x, = vars_x
xi = symbols('xi', real=True)
# Hamiltonian vector field
X_H = hamiltonian_vector_field(H_expr, vars_x)
# Grid
x_vals = np.linspace(*xlim, 30)
xi_vals = np.linspace(*klim, 30)
X, XI = np.meshgrid(x_vals, xi_vals, indexing='ij')
# Evaluate the field
from sympy import lambdify
f_x = lambdify((x, xi), X_H['dx/dt'], 'numpy')
f_xi = lambdify((x, xi), X_H['dxi/dt'], 'numpy')
U = f_x(X, XI)
V = f_xi(X, XI)
# Evaluate H
H_func = lambdify((x, xi), H_expr, 'numpy')
H_vals = H_func(X, XI)
fig, axes = plt.subplots(1, 2, figsize=(16, 6))
# Plot 1: Hamiltonian field + level curves of H
ax1 = axes[0]
# Level curves (energy surfaces)
contours = ax1.contour(X, XI, H_vals, levels=15, colors='blue', alpha=0.4)
ax1.clabel(contours, inline=True, fontsize=8)
# Vector field (always tangent to level curves!)
ax1.quiver(X[::2, ::2], XI[::2, ::2], U[::2, ::2], V[::2, ::2],
color='red', alpha=0.7, scale=20, width=0.003)
ax1.set_xlabel('x (position)', fontsize=12)
ax1.set_ylabel('ξ (momentum)', fontsize=12)
ax1.set_title(f'Hamiltonian vector field X_H\nH = {H_expr}', fontsize=13)
ax1.grid(True, alpha=0.3)
ax1.set_aspect('equal')
# Plot 2: Symplectic form (abstract visualization)
ax2 = axes[1]
# Visualize ω as a 2-form via its action on vectors
# ω(v, w) = v_x·w_ξ - v_ξ·w_x
# Choose a few points
test_points = [(0, 0), (1, 0.5), (-0.5, -1)]
for (x0, xi0) in test_points:
# Tangent vector to the flow
v = np.array([f_x(x0, xi0), f_xi(x0, xi0)])
# Symplectic orthogonal vector: J·v where J = [[0, -1], [1, 0]]
J = np.array([[0, -1], [1, 0]])
w = J @ v
# Verify: ω(v, w) = det([v, w])
omega_vw = np.linalg.det(np.column_stack([v, w]))
# Draw
ax2.arrow(x0, xi0, v[0]*0.2, v[1]*0.2,
head_width=0.1, head_length=0.05, fc='red', ec='red', alpha=0.7)
ax2.arrow(x0, xi0, w[0]*0.2, w[1]*0.2,
head_width=0.1, head_length=0.05, fc='blue', ec='blue', alpha=0.7)
ax2.text(x0+0.3, xi0+0.3, f'ω(v,w)={omega_vw:.2f}', fontsize=9)
ax2.contour(X, XI, H_vals, levels=15, colors='gray', alpha=0.2, linewidths=0.5)
ax2.set_xlabel('x', fontsize=12)
ax2.set_ylabel('ξ', fontsize=12)
ax2.set_title('Symplectic structure ω = dξ∧dx', fontsize=13)
ax2.grid(True, alpha=0.3)
ax2.set_xlim(xlim)
ax2.set_ylim(klim)
ax2.legend(['X_H (red)', 'J·X_H (blue)'], loc='upper right')
plt.tight_layout()
plt.show()
# Examples
import sympy as sp
visualize_symplectic_geometry(xi**2/2 + x**2/2, [x]) # Harmonic
visualize_symplectic_geometry(xi**2/2 - sp.cos(x), [x]) # Pendulum
[docs]
def quantize_hamiltonian(H_classical, vars_x, mode='weyl'):
"""
Quantizes a classical Hamiltonian into a pseudo-differential operator
H_classical(x, ξ) → Ĥ operator
mode: 'weyl' (symmetric) or 'kn' (Kohn-Nirenberg)
"""
from psiop import PseudoDifferentialOperator
# Create the operator
H_op = PseudoDifferentialOperator(
expr=H_classical,
vars_x=vars_x,
mode='symbol'
)
print(f"Classical Hamiltonian: H = {H_classical}")
print(f"Quantization ({mode}):")
# Properties of the quantized operator
print(f" - Order: {H_op.symbol_order()}")
print(f" - Self-adjoint: {H_op.is_self_adjoint()}")
# Commutator with x and ξ (Heisenberg relations)
x_op = PseudoDifferentialOperator(vars_x[0], vars_x, mode='symbol')
xi_op = PseudoDifferentialOperator(symbols('xi', real=True), vars_x, mode='symbol')
comm_H_x = H_op.commutator_symbolic(x_op, order=1, mode=mode)
comm_H_xi = H_op.commutator_symbolic(xi_op, order=1, mode=mode)
print(f"\n [Ĥ, x̂] = {simplify(comm_H_x)}")
print(f" [Ĥ, p̂] = {simplify(comm_H_xi)}")
# Classical-quantum correspondence
print(f"\n Correspondence principle:")
print(f" [Ĥ, x̂] ≈ iℏ {{H, x}} = iℏ (∂H/∂ξ)")
print(f" [Ĥ, p̂] ≈ iℏ {{H, ξ}} = -iℏ (∂H/∂x)")
return H_op
# Test
from sympy import symbols
x = symbols('x', real=True)
xi = symbols('xi', real=True)
print("="*70)
print("QUANTIZATION: Harmonic oscillator")
print("="*70)
H_harm = xi**2/2 + x**2/2
H_op_harm = quantize_hamiltonian(H_harm, [x], mode='weyl')
print("\n" + "="*70)
print("QUANTIZATION: Anharmonic Hamiltonian")
print("="*70)
H_anharm = xi**2/2 + x**4/4
H_op_anharm = quantize_hamiltonian(H_anharm, [x], mode='weyl')
omega = canonical_symplectic_form_construction()
[docs]
def from_omega_to_flow(H_expr, vars_x):
"""
Shows how the symplectic structure ω allows passage
from H (function) to X_H (vector field) to Φₜ (flow)
"""
from sympy import symbols, diff, simplify, lambdify
from scipy.integrate import solve_ivp
import numpy as np
import matplotlib.pyplot as plt
dim = len(vars_x)
if dim == 1:
x, = vars_x
xi = symbols('xi', real=True)
print("="*70)
print("FROM ω TO HAMILTONIAN FLOW: COMPLETE MECHANISM")
print("="*70)
# Step 1: The symplectic structure ω
print("\n[Step 1] Symplectic structure ω = dξ ∧ dx")
print("-"*70)
print("In matrix coordinates:")
print(" ω = [[0, -1], [1, 0]]")
omega = np.array([[0, -1], [1, 0]])
omega_inv = np.linalg.inv(omega)
print(f"\nω⁻¹ = {omega_inv}")
# Step 2: The Hamiltonian H
print(f"\n[Step 2] Hamiltonian H(x, ξ) = {H_expr}")
print("-"*70)
# Step 3: Gradient of H
print("\n[Step 3] Gradient dH = (∂H/∂x, ∂H/∂ξ)ᵀ")
print("-"*70)
dH_dx = diff(H_expr, x)
dH_dxi = diff(H_expr, xi)
print(f" ∂H/∂x = {dH_dx}")
print(f" ∂H/∂ξ = {dH_dxi}")
print(f"\n dH = ({dH_dx}, {dH_dxi})ᵀ")
# Step 4: Hamiltonian field X_H = ω⁻¹ · dH
print("\n[Step 4] Hamiltonian field X_H via ω⁻¹")
print("-"*70)
print("Definition: ι_{X_H} ω = dH")
print("In coordinates: X_H = ω⁻¹ · dH")
print(f"\nX_H = ω⁻¹ · dH")
print(f" = [[0, 1], [-1, 0]] · ({dH_dx}, {dH_dxi})ᵀ")
print(f" = ({dH_dxi}, -{dH_dx})ᵀ")
dx_dt = dH_dxi
dxi_dt = -dH_dx
print(f"\n dx/dt = ∂H/∂ξ = {dx_dt}")
print(f" dξ/dt = -∂H/∂x = {dxi_dt}")
# Step 5: Flow integration
print("\n[Step 5] Hamiltonian flow Φₜ by integrating X_H")
print("-"*70)
f_x = lambdify((x, xi), dx_dt, 'numpy')
f_xi = lambdify((x, xi), dxi_dt, 'numpy')
def hamiltonian_system(t, y):
return [f_x(y[0], y[1]), f_xi(y[0], y[1])]
# Initial condition
y0 = [1.0, 0.5]
t_span = [0, 10]
t_eval = np.linspace(*t_span, 200)
sol = solve_ivp(hamiltonian_system, t_span, y0,
t_eval=t_eval, method='DOP853')
print(f" Initial condition: (x₀, ξ₀) = {y0}")
print(f" Integration from t=0 to t={t_span[1]}")
# Step 6: Conservation verification
print("\n[Step 6] Verifications")
print("-"*70)
# Conservation of H
H_func = lambdify((x, xi), H_expr, 'numpy')
H_values = H_func(sol.y[0], sol.y[1])
print(f"\n✓ Conservation of H:")
print(f" H(t=0) = {H_values[0]:.10f}")
print(f" H(t={t_span[1]}) = {H_values[-1]:.10f}")
print(f" |ΔH|/H₀ = {abs(H_values[-1] - H_values[0])/abs(H_values[0]):.2e}")
# Conservation of area (symplectic volume)
# Numerically compute the Jacobian of the flow
epsilon = 1e-6
t_test = 5.0
idx = np.argmin(np.abs(sol.t - t_test))
x0, xi0 = sol.y[0, idx], sol.y[1, idx]
# Perturbations
sol_x_plus = solve_ivp(hamiltonian_system, [0, t_test],
[y0[0] + epsilon, y0[1]], dense_output=True)
sol_x_minus = solve_ivp(hamiltonian_system, [0, t_test],
[y0[0] - epsilon, y0[1]], dense_output=True)
sol_xi_plus = solve_ivp(hamiltonian_system, [0, t_test],
[y0[0], y0[1] + epsilon], dense_output=True)
sol_xi_minus = solve_ivp(hamiltonian_system, [0, t_test],
[y0[0], y0[1] - epsilon], dense_output=True)
# Jacobian
dx_dx0 = (sol_x_plus.sol(t_test)[0] - sol_x_minus.sol(t_test)[0]) / (2*epsilon)
dx_dxi0 = (sol_xi_plus.sol(t_test)[0] - sol_xi_minus.sol(t_test)[0]) / (2*epsilon)
dxi_dx0 = (sol_x_plus.sol(t_test)[1] - sol_x_minus.sol(t_test)[1]) / (2*epsilon)
dxi_dxi0 = (sol_xi_plus.sol(t_test)[1] - sol_xi_minus.sol(t_test)[1]) / (2*epsilon)
J = np.array([[dx_dx0, dx_dxi0], [dxi_dx0, dxi_dxi0]])
det_J = np.linalg.det(J)
print(f"\n✓ Conservation of symplectic structure (Liouville's Theorem):")
print(f" det(DΦₜ) at t={t_test}: {det_J:.10f}")
print(f" (must be 1 to preserve ω)")
# Visualization
fig, axes = plt.subplots(2, 2, figsize=(14, 12))
# Plot 1: Trajectory in phase space
ax1 = axes[0, 0]
# Level curves of H
x_grid = np.linspace(-2, 2, 100)
xi_grid = np.linspace(-2, 2, 100)
X, XI = np.meshgrid(x_grid, xi_grid)
H_grid = H_func(X, XI)
contours = ax1.contour(X, XI, H_grid, levels=15, colors='lightblue', alpha=0.5)
ax1.clabel(contours, inline=True, fontsize=8)
# Trajectory
ax1.plot(sol.y[0], sol.y[1], 'r-', linewidth=2, label='Trajectory')
ax1.plot(y0[0], y0[1], 'go', markersize=10, label='Initial point')
ax1.plot(sol.y[0, -1], sol.y[1, -1], 'rs', markersize=10, label='Final point')
ax1.set_xlabel('x (position)', fontsize=11)
ax1.set_ylabel('ξ (momentum)', fontsize=11)
ax1.set_title('Trajectory in phase space', fontsize=12)
ax1.legend()
ax1.grid(True, alpha=0.3)
ax1.set_aspect('equal')
# Plot 2: Conservation of H
ax2 = axes[0, 1]
ax2.plot(sol.t, H_values, 'b-', linewidth=2)
ax2.axhline(H_values[0], color='r', linestyle='--',
label=f'H₀ = {H_values[0]:.6f}')
ax2.set_xlabel('Time t', fontsize=11)
ax2.set_ylabel('H(x(t), ξ(t))', fontsize=11)
ax2.set_title('Energy conservation', fontsize=12)
ax2.legend()
ax2.grid(True, alpha=0.3)
# Plot 3: Coordinate evolution
ax3 = axes[1, 0]
ax3.plot(sol.t, sol.y[0], 'b-', linewidth=2, label='x(t)')
ax3.plot(sol.t, sol.y[1], 'r-', linewidth=2, label='ξ(t)')
ax3.set_xlabel('Time t', fontsize=11)
ax3.set_ylabel('Coordinates', fontsize=11)
ax3.set_title('Time evolution', fontsize=12)
ax3.legend()
ax3.grid(True, alpha=0.3)
# Plot 4: Hamiltonian vector field
ax4 = axes[1, 1]
x_vals = np.linspace(-2, 2, 20)
xi_vals = np.linspace(-2, 2, 20)
X_field, XI_field = np.meshgrid(x_vals, xi_vals, indexing='ij')
U = f_x(X_field, XI_field)
V = f_xi(X_field, XI_field)
ax4.quiver(X_field, XI_field, U, V, alpha=0.6, scale=20, width=0.003)
ax4.contour(X, XI, H_grid, levels=15, colors='lightblue', alpha=0.3)
ax4.plot(sol.y[0], sol.y[1], 'r-', linewidth=2, alpha=0.8)
ax4.set_xlabel('x', fontsize=11)
ax4.set_ylabel('ξ', fontsize=11)
ax4.set_title('Hamiltonian field X_H', fontsize=12)
ax4.grid(True, alpha=0.3)
ax4.set_aspect('equal')
plt.tight_layout()
plt.show()
return sol
# Test with different Hamiltonians
from sympy import symbols, cos
x = symbols('x', real=True)
xi = symbols('xi', real=True)
print("\n" + "="*70)
print("TEST 1: HARMONIC OSCILLATOR")
print("="*70)
sol_harm = from_omega_to_flow(xi**2/2 + x**2/2, [x])
print("\n" + "="*70)
print("TEST 2: PENDULUM")
print("="*70)
sol_pendulum = from_omega_to_flow(xi**2/2 - cos(x), [x])
[docs]
def interactive_symplectic_dynamics(H_options):
"""
Interactive dashboard to explore different Hamiltonians
with the SAME symplectic structure
"""
from ipywidgets import interact, Dropdown, FloatSlider
import matplotlib.pyplot as plt
@interact(
H_choice=Dropdown(
options=list(H_options.keys()),
value=list(H_options.keys())[0],
description='Hamiltonian:'
),
x0=FloatSlider(min=-2, max=2, step=0.1, value=1.0, description='x₀'),
xi0=FloatSlider(min=-2, max=2, step=0.1, value=0.5, description='ξ₀'),
t_max=FloatSlider(min=1, max=20, step=1, value=10, description='Max time')
)
def plot_dynamics(H_choice, x0, xi0, t_max):
from sympy import lambdify, symbols
from scipy.integrate import solve_ivp
import numpy as np
x = symbols('x', real=True)
xi = symbols('xi', real=True)
H_expr = H_options[H_choice]
# Hamiltonian field
from sympy import diff
dx_dt = diff(H_expr, xi)
dxi_dt = -diff(H_expr, x)
f_x = lambdify((x, xi), dx_dt, 'numpy')
f_xi = lambdify((x, xi), dxi_dt, 'numpy')
def system(t, y):
return [f_x(y[0], y[1]), f_xi(y[0], y[1])]
# Integration
sol = solve_ivp(system, [0, t_max], [x0, xi0],
dense_output=True, method='DOP853')
t_eval = np.linspace(0, t_max, 500)
trajectory = sol.sol(t_eval)
# Visualization
fig, axes = plt.subplots(1, 2, figsize=(16, 6))
# Phase space
x_grid = np.linspace(-3, 3, 100)
xi_grid = np.linspace(-3, 3, 100)
X, XI = np.meshgrid(x_grid, xi_grid)
H_func = lambdify((x, xi), H_expr, 'numpy')
H_grid = H_func(X, XI)
axes[0].contour(X, XI, H_grid, levels=20, colors='lightgray', alpha=0.5)
axes[0].plot(trajectory[0], trajectory[1], 'b-', linewidth=2)
axes[0].plot(x0, xi0, 'go', markersize=12, label='Initial')
axes[0].plot(trajectory[0, -1], trajectory[1, -1], 'ro',
markersize=12, label='Final')
axes[0].set_xlabel('x', fontsize=12)
axes[0].set_ylabel('ξ', fontsize=12)
axes[0].set_title(f'Phase space\nH = {H_expr}', fontsize=13)
axes[0].legend()
axes[0].grid(True, alpha=0.3)
axes[0].set_aspect('equal')
# Time evolution
axes[1].plot(t_eval, trajectory[0], 'b-', linewidth=2, label='x(t)')
axes[1].plot(t_eval, trajectory[1], 'r-', linewidth=2, label='ξ(t)')
axes[1].set_xlabel('Time t', fontsize=12)
axes[1].set_ylabel('Coordinates', fontsize=12)
axes[1].set_title('Time evolution', fontsize=13)
axes[1].legend()
axes[1].grid(True, alpha=0.3)
plt.tight_layout()
plt.show()
# Conservation
H_values = H_func(trajectory[0], trajectory[1])
print(f"\nConservation of H:")
print(f" H(t=0) = {H_values[0]:.8f}")
print(f" H(t={t_max}) = {H_values[-1]:.8f}")
print(f" |ΔH|/H₀ = {abs(H_values[-1] - H_values[0])/abs(H_values[0]):.2e}")
# Examples
from sympy import symbols, cos, sin, exp
x, xi = symbols('x xi', real=True)
H_library = {
"Free": xi**2/2,
"Harmonic": xi**2/2 + x**2/2,
"Anharmonic": xi**2/2 + x**4/4,
"Pendulum": xi**2/2 - cos(x),
"Double well": xi**2/2 - x**2/2 + x**4/4,
"Morse": xi**2/2 + (1 - exp(-x))**2,
"Hénon-Heiles (simplified)": xi**2/2 + x**2/2 + x**3,
}
interactive_symplectic_dynamics(H_library)
[docs]
def explicit_derivation_hamilton_equations():
"""
Explicitly derives Hamilton's equations
from only the data of ω and H
"""
from sympy import symbols, Matrix, solve, simplify
print("="*70)
print("EXPLICIT DERIVATION: FROM ω + H → HAMILTON'S EQUATIONS")
print("="*70)
x, xi = symbols('x xi', real=True)
dx_dt, dxi_dt = symbols('dx/dt dxi/dt', real=True)
# Input 1: Symplectic structure
print("\n[Input 1] Symplectic structure ω")
print("-"*70)
print("ω = dξ ∧ dx")
print("\nIn matrix form (basis (x, ξ)):")
omega = Matrix([[0, -1], [1, 0]])
print(omega)
# Input 2: Hamiltonian
print("\n[Input 2] Hamiltonian H(x, ξ)")
print("-"*70)
# Generic symbol
from sympy import Function
H = Function('H')(x, xi)
print(f"H = H(x, ξ) (arbitrary function)")
# Derivation
print("\n[Derivation] Equation defining X_H")
print("-"*70)
print("The Hamiltonian vector field X_H = (dx/dt, dξ/dt) is defined by:")
print(" ι_{X_H} ω = dH")
print("\nCoordinate translation:")
# Vector field
X_H = Matrix([dx_dt, dxi_dt])
print(f"\nX_H = {X_H.T}")
# Gradient of H
from sympy import diff
dH_dx = diff(H, x)
dH_dxi = diff(H, xi)
grad_H = Matrix([dH_dx, dH_dxi])
print(f"\ndH = ({dH_dx}, {dH_dxi})ᵀ")
# Contraction ι_{X_H} ω
print("\nContraction ι_{X_H} ω:")
print(" ι_{X_H} ω(·) = ω(X_H, ·)")
# In matrix form: ω · X_H
contraction = omega * X_H
print(f"\n ω · X_H = {omega} · {X_H.T}")
print(f" = {contraction.T}")
# Equation ι_{X_H} ω = dH
print("\nEquation to solve: ω · X_H = dH")
print(f" {contraction.T} = {grad_H.T}")
# System of equations
equations = [contraction[i] - grad_H[i] for i in range(2)]
print("\nSystem:")
for i, eq in enumerate(equations):
print(f" Equation {i+1} : {eq} = 0")
# Solve
solution = solve(equations, [dx_dt, dxi_dt])
print("\n[Solution] Hamilton's equations")
print("-"*70)
print(f" dx/dt = {solution[dx_dt]}")
print(f" dξ/dt = {solution[dxi_dt]}")
print("\n✨ UNIVERSAL RESULT ✨")
print("-"*70)
print("For ANY Hamiltonian H(x, ξ):")
print(" • dx/dt = ∂H/∂ξ")
print(" • dξ/dt = -∂H/∂x")
print("\nThese equations emerge SOLELY from:")
print(" 1. The symplectic structure ω = dξ ∧ dx (independent of H)")
print(" 2. The geometric equation ι_{X_H} ω = dH")
# Verification with concrete examples
print("\n" + "="*70)
print("VERIFICATION WITH CONCRETE EXAMPLES")
print("="*70)
hamiltonians = {
"Free": xi**2/2,
"Harmonic": xi**2/2 + x**2/2,
"Pendulum": xi**2/2 - sp.cos(x),
}
for name, H_concrete in hamiltonians.items():
print(f"\n{name} : H = {H_concrete}")
print("-"*60)
dH_dx = diff(H_concrete, x)
dH_dxi = diff(H_concrete, xi)
dx_dt_result = dH_dxi
dxi_dt_result = -dH_dx
print(f" dx/dt = ∂H/∂ξ = {dx_dt_result}")
print(f" dξ/dt = -∂H/∂x = {dxi_dt_result}")
return omega, solution
omega, hamilton_eqs = explicit_derivation_hamilton_equations()
[docs]
def symplectic_structure_higher_dimensions(n=2):
"""
Shows how ω generalizes to arbitrary dimension
"""
from sympy import symbols, Matrix, simplify
print("="*70)
print(f"SYMPLECTIC STRUCTURE IN DIMENSION {n}")
print("="*70)
# Variables
x_vars = symbols(f'x1:{n+1}', real=True)
xi_vars = symbols(f'xi1:{n+1}', real=True)
print(f"\n[Variables] T*ℝ^{n} has dimension {2*n}")
print("-"*70)
print(f"Spatial variables : {x_vars}")
print(f"Momentum variables : {xi_vars}")
# Symplectic form
print("\n[Canonical symplectic form]")
print("-"*70)
print("ω = Σᵢ dξᵢ ∧ dxⁱ")
# Symplectic matrix
omega = Matrix.zeros(2*n, 2*n)
for i in range(n):
# Block [[0, -I], [I, 0]]
omega[i, n+i] = -1
omega[n+i, i] = 1
print(f"\nMatrix ω (size {2*n}×{2*n}) :")
print(omega)
# Properties
print("\n[Properties]")
print("-"*70)
# Non-degenerate
det_omega = omega.det()
print(f"✓ det(ω) = {det_omega} ≠ 0 (non-degenerate)")
# Antisymmetry
is_antisymmetric = (omega + omega.T).is_zero_matrix
print(f"✓ ω = -ωᵀ : {is_antisymmetric} (antisymmetric)")
# Inverse
omega_inv = omega.inv()
print(f"\n✓ ω⁻¹ exists :")
print(omega_inv)
# Hamilton's equations
print("\n[General Hamilton's equations]")
print("-"*70)
print("For H(x₁, ..., xₙ, ξ₁, ..., ξₙ) :")
print()
for i in range(n):
print(f" dxᵢ/dt = ∂H/∂ξᵢ")
print()
for i in range(n):
print(f" dξᵢ/dt = -∂H/∂xᵢ")
# Block structure
print("\n[Block structure]")
print("-"*70)
print("ω and ω⁻¹ have canonical block form:")
print()
print(" ⎡ 0 -I ⎤ ⎡ 0 I ⎤")
print(" ω = ⎢ ⎥, ω⁻¹ = ⎢ ⎥")
print(" ⎣ I 0 ⎦ ⎣ -I 0 ⎦")
print()
print("where I is the n×n identity matrix")
return omega
# Test in dimensions 2 and 3
omega_2d = symplectic_structure_higher_dimensions(n=2)
print("\n" + "="*70 + "\n")
omega_3d = symplectic_structure_higher_dimensions(n=3)
canonical_transformations_preserve_omega()
[docs]
def from_poisson_to_commutator():
"""
Shows the fundamental link between Poisson bracket (classical)
and commutator (quantum) via correspondence
"""
from sympy import symbols, diff, simplify, I
print("="*70)
print("FROM POISSON BRACKET TO QUANTUM COMMUTATOR")
print("="*70)
x, xi = symbols('x xi', real=True)
print("\n[Classical level: Poisson bracket]")
print("-"*70)
print("For two observables f, g : T*M → ℝ :")
print(" {f, g} = ∂f/∂x · ∂g/∂ξ - ∂f/∂ξ · ∂g/∂x")
print("\nThis bracket is INDUCED by the symplectic structure ω")
# Examples of brackets
print("\n[Fundamental brackets]")
print("-"*60)
observables = {
'x and ξ': (x, xi),
'x and x': (x, x),
'ξ and ξ': (xi, xi),
'x² and ξ': (x**2, xi),
'x and ξ²': (x, xi**2),
}
for name, (f, g) in observables.items():
bracket = diff(f, x)*diff(g, xi) - diff(f, xi)*diff(g, x)
bracket = simplify(bracket)
print(f" {{{name}}} = {bracket}")
print("\n[Quantum level: Commutator]")
print("-"*70)
print("In quantum mechanics, observables become operators")
print("The Poisson bracket becomes the commutator via:")
print()
print(" DIRAC CORRESPONDENCE PRINCIPLE")
print(" ═══════════════════════════════")
print(" {f, g}_classical ↦ (1/iℏ) [f̂, ĝ]_quantum")
print()
print("Equivalently:")
print(" [f̂, ĝ] = iℏ {f, g}")
print("\n[Canonical commutation relations]")
print("-"*60)
print("From the Poisson bracket {x, ξ} = 1, we deduce:")
print(" [x̂, p̂] = iℏ")
print("\nThis is the foundation of quantum mechanics!")
print("\n[Correspondence table]")
print("-"*60)
correspondence_table = [
("Classical", "Quantum"),
("─"*30, "─"*30),
("{x, ξ} = 1", "[x̂, p̂] = iℏ"),
("{x, x} = 0", "[x̂, x̂] = 0"),
("{ξ, ξ} = 0", "[p̂, p̂] = 0"),
("{x², ξ} = 2x", "[x̂², p̂] = 2iℏx̂"),
("{H, f} = df/dt", "[Ĥ, f̂] = iℏ df̂/dt"),
]
print()
for classical, quantum in correspondence_table:
print(f" {classical:30s} ←→ {quantum}")
# Weyl quantization
print("\n[Weyl quantization]")
print("-"*70)
print("To quantize a classical Hamiltonian H(x, ξ):")
print()
print("1. STANDARD QUANTIZATION (Kohn-Nirenberg):")
print(" H(x, ξ) ↦ H(x̂, p̂) (operator ordering: position then momentum)")
print()
print("2. WEYL QUANTIZATION (symmetric):")
print(" Operator Ĥ defined via its Weyl symbol")
print(" Better preserves symplectic structure")
print(" Self-adjoint if H is real")
print("\n[Example: Harmonic oscillator]")
print("-"*60)
H_classical = xi**2/2 + x**2/2
print(f"Classical Hamiltonian : H = {H_classical}")
print()
print("Standard quantization:")
print(" Ĥ = p̂²/2 + x̂²/2")
print()
print("Commutation relations used:")
print(" [x̂, p̂] = iℏ (derived from {x, ξ} = 1)")
print()
print("Energy spectrum:")
print(" Eₙ = ℏω(n + 1/2), n = 0, 1, 2, ...")
# Conceptual diagram
print("\n[CONCEPTUAL DIAGRAM]")
print("="*70)
print()
print(" T*M with symplectic structure ω")
print(" │")
print(" │ Induces")
print(" ↓")
print(" Poisson bracket {·, ·}")
print(" │")
print(" │ Quantization")
print(" │ (ℏ → 0 : classical limit)")
print(" ↓")
print(" Commutator [·, ·] = iℏ{·, ·}")
print(" │")
print(" │ Defines")
print(" ↓")
print(" Algebra of quantum observables")
print()
print("="*70)
return correspondence_table
correspondence = from_poisson_to_commutator()
[docs]
def conceptual_architecture():
"""
High-level overview of the conceptual architecture
"""
diagram = """
╔═══════════════════════════════════════════════════════════════════╗
║ CONCEPTUAL ARCHITECTURE ║
║ Emergence of the symplectic structure ║
╚═══════════════════════════════════════════════════════════════════╝
LEVEL 0: GEOMETRY OF THE COTANGENT BUNDLE
══════════════════════════════════════════
Base manifold M (physical space)
│
│ Geometric construction
│
↓
Cotangent bundle T*M (phase space)
│
│ Intrinsic structure
│
↓
Liouville 1-form θ = Σᵢ ξᵢ dxⁱ (CANONICAL)
│
│ Exterior derivative d
│
↓
Symplectic form ω = dθ (UNIVERSAL)
✨ ω exists BEFORE any Hamiltonian H is chosen
LEVEL 1: HAMILTONIAN DYNAMICS
═════════════════════════════
Function H : T*M → ℝ (arbitrary Hamiltonian)
│
│ Equation: ι_{X_H} ω = dH
│
↓
Hamiltonian vector field X_H (tangent to T*M)
│
│ Integration
│
↓
Flow Φₜᴴ : T*M → T*M
✨ Φₜᴴ PRESERVES ω (Liouville's theorem)
LEVEL 2: ALGEBRAIC STRUCTURE
═════════════════════════════
Algebra C∞(T*M) of smooth functions
│
│ ω induces a bracket
│
↓
Poisson bracket {f, g}
│
│ Properties:
│ • Bilinear
│ • Antisymmetric
│ • Jacobi identity
│ • Leibniz derivation
│
↓
Lie algebra structure
✨ {f, g} is COMPLETELY determined by ω
LEVEL 3: QUANTIZATION
═══════════════════════
Poisson bracket {·, ·} (classical)
│
│ Correspondence:
│ {f, g} ↦ (1/iℏ)[f̂, ĝ]
│
↓
Commutator [·, ·] (quantum)
│
│ Example: [x̂, p̂] = iℏ
│ (from {x, ξ} = 1)
│
↓
Pseudo-differential operators
✨ Quantum structure INHERITS from ω
UNIVERSAL PROPERTIES
═════════════════════
✓ ω is INDEPENDENT of the choice of H
✓ ω is CANONICAL (coordinate-independent)
✓ ω is CLOSED: dω = 0
✓ ω is NON-DEGENERATE: det(ω) ≠ 0
✓ All H → same ω → same conservation laws
✓ Quantization preserves symplectic structure
KEY THEOREMS
═════════════
• Darboux’s Theorem:
Locally, every symplectic structure is isomorphic to
the canonical form ω = Σᵢ dξᵢ ∧ dxⁱ
• Liouville’s Theorem (geometric):
Hamiltonian flow preserves ω
• Noether’s Theorem (symplectic):
Symmetries ↔ Conserved quantities via ω
"""
print(diagram)
# Concise summary
print("\n" + "="*70)
print("SUMMARY: HOW DOES ω EMERGE?")
print("="*70)
print()
print("The symplectic structure ω is NOT a consequence of the Hamiltonian H,")
print("but an INTRINSIC GEOMETRIC PROPERTY of the cotangent bundle T*M.")
print()
print("Emergence sequence:")
print(" 1. M (manifold) → T*M (cotangent bundle) [geometric construction]")
print(" 2. T*M → θ (Liouville form) [unique canonical structure]")
print(" 3. θ → ω = dθ [exterior derivative]")
print(" 4. ω + H → X_H [equation ι_{X_H}ω = dH]")
print(" 5. X_H → Φₜ [flow integration]")
print()
print("AT EACH STEP, the structure is DETERMINED by geometry,")
print("NOT by an arbitrary choice of Hamiltonian.")
print()
print("="*70)
conceptual_architecture()
[docs]
def complete_symplectic_example_unified():
"""
Complete example showing the ENTIRE pipeline:
H → ω → X_H → Φₜ → Conserved quantities
Uses classes from psiop.py for integration
"""
from sympy import symbols, cos, sin, diff, lambdify
from scipy.integrate import solve_ivp
import numpy as np
import matplotlib.pyplot as plt
from psiop import PseudoDifferentialOperator
print("="*70)
print("COMPLETE EXAMPLE: PENDULUM WITH SYMPLECTIC STRUCTURE")
print("="*70)
# Step 1: Define the Hamiltonian
print("\n[STEP 1] Hamiltonian definition")
print("-"*70)
x, xi = symbols('x xi', real=True)
g, L, m = 1.0, 1.0, 1.0 # Physical constants
H_pendulum = m*L**2*xi**2/2 - m*g*L*cos(x)
print(f"Pendulum Hamiltonian: H = {H_pendulum}")
print(f" Kinetic term: T = m·L²·ξ²/2")
print(f" Potential term: V = -m·g·L·cos(x)")
# Step 2: Symplectic structure (independent of H!)
print("\n[STEP 2] Canonical symplectic structure")
print("-"*70)
from sympy import Matrix
omega = Matrix([[0, -1], [1, 0]])
print("ω = dξ ∧ dx")
print("\nSymplectic matrix:")
print(omega)
print("\nThis structure exists BEFORE we choose H!")
# Step 3: Hamiltonian vector field via ω
print("\n[STEP 3] Hamiltonian vector field X_H = ω⁻¹·dH")
print("-"*70)
dH_dx = diff(H_pendulum, x)
dH_dxi = diff(H_pendulum, xi)
dx_dt = dH_dxi
dxi_dt = -dH_dx
print(f"dH = ({dH_dx}, {dH_dxi})ᵀ")
print(f"\nHamiltonian vector field:")
print(f" dx/dt = ∂H/∂ξ = {dx_dt}")
print(f" dξ/dt = -∂H/∂x = {dxi_dt}")
# Step 4: Flow integration
print("\n[STEP 4] Integration of the Hamiltonian flow")
print("-"*70)
f_x = lambdify((x, xi), dx_dt, 'numpy')
f_xi = lambdify((x, xi), dxi_dt, 'numpy')
def pendulum_system(t, y):
return [f_x(y[0], y[1]), f_xi(y[0], y[1])]
# Multiple initial conditions
initial_conditions = [
(np.pi/6, 0.0, "Small oscillation"),
(np.pi/2, 0.0, "Large oscillation"),
(np.pi, 0.5, "Near separatrix"),
(0.0, 2.0, "Rotation"),
]
t_span = [0, 20]
t_eval = np.linspace(*t_span, 1000)
solutions = []
for x0, xi0, label in initial_conditions:
sol = solve_ivp(
pendulum_system, t_span, [x0, xi0],
t_eval=t_eval, method='DOP853', rtol=1e-10, atol=1e-12
)
solutions.append((sol, label))
print(f" {label:20s} : (x₀, ξ₀) = ({x0:.2f}, {xi0:.2f})")
# Step 5: Conservation checks
print("\n[STEP 5] Verifications (Liouville's theorems)")
print("-"*70)
H_func = lambdify((x, xi), H_pendulum, 'numpy')
for (sol, label), (x0, xi0, _) in zip(solutions, initial_conditions):
H_values = H_func(sol.y[0], sol.y[1])
H_conservation = abs(H_values[-1] - H_values[0]) / abs(H_values[0])
# Compute Jacobian to verify det = 1
epsilon = 1e-7
t_test = 10.0
# Perturbations
sol_x_p = solve_ivp(pendulum_system, [0, t_test],
[x0 + epsilon, xi0], method='DOP853', rtol=1e-10)
sol_x_m = solve_ivp(pendulum_system, [0, t_test],
[x0 - epsilon, xi0], method='DOP853', rtol=1e-10)
sol_xi_p = solve_ivp(pendulum_system, [0, t_test],
[x0, xi0 + epsilon], method='DOP853', rtol=1e-10)
sol_xi_m = solve_ivp(pendulum_system, [0, t_test],
[x0, xi0 - epsilon], method='DOP853', rtol=1e-10)
# Jacobian at t_test
J11 = (sol_x_p.y[0, -1] - sol_x_m.y[0, -1]) / (2*epsilon)
J12 = (sol_xi_p.y[0, -1] - sol_xi_m.y[0, -1]) / (2*epsilon)
J21 = (sol_x_p.y[1, -1] - sol_x_m.y[1, -1]) / (2*epsilon)
J22 = (sol_xi_p.y[1, -1] - sol_xi_m.y[1, -1]) / (2*epsilon)
det_J = J11*J22 - J12*J21
print(f"\n {label:20s} :")
print(f" Energy conservation: |ΔH|/H₀ = {H_conservation:.2e}")
print(f" Symplectic conservation: det(DΦₜ) = {det_J:.10f} (≈ 1)")
# Step 6: Full visualization
print("\n[STEP 6] Visualization of symplectic geometry")
print("-"*70)
fig = plt.figure(figsize=(18, 12))
gs = fig.add_gridspec(3, 3, hspace=0.3, wspace=0.3)
# Plot 1: Full phase portrait
ax1 = fig.add_subplot(gs[0:2, 0:2])
# Grid for contour lines
x_grid = np.linspace(-np.pi, 2*np.pi, 300)
xi_grid = np.linspace(-3, 3, 300)
X, XI = np.meshgrid(x_grid, xi_grid, indexing='ij')
H_grid = H_func(X, XI)
# Energy level curves
contours = ax1.contour(X, XI, H_grid, levels=30, colors='lightblue',
alpha=0.4, linewidths=1)
ax1.clabel(contours, inline=True, fontsize=7, fmt='%.1f')
# Separatrix (critical energy)
E_sep = H_func(np.pi, 0)
ax1.contour(X, XI, H_grid, levels=[E_sep], colors='red',
linewidths=3, linestyles='--')
# Trajectories
colors = ['green', 'blue', 'orange', 'purple']
for (sol, label), color in zip(solutions, colors):
ax1.plot(sol.y[0], sol.y[1], color=color, linewidth=2,
label=label, alpha=0.8)
# Initial point
ax1.plot(sol.y[0, 0], sol.y[1, 0], 'o', color=color,
markersize=8, markeredgecolor='black', markeredgewidth=1)
ax1.set_xlabel('x (angle)', fontsize=12)
ax1.set_ylabel('ξ (angular momentum)', fontsize=12)
ax1.set_title('Pendulum phase portrait\n(Contour lines = constant energy surfaces)',
fontsize=13)
ax1.legend(loc='upper right', fontsize=9)
ax1.grid(True, alpha=0.3)
ax1.set_xlim(-np.pi, 2*np.pi)
ax1.set_ylim(-3, 3)
# Mark equilibrium points
ax1.plot(0, 0, 'k*', markersize=15, label='Stable equilibrium')
ax1.plot(np.pi, 0, 'rx', markersize=15, markeredgewidth=3,
label='Unstable equilibrium')
ax1.legend(loc='upper right', fontsize=9)
# Plot 2: Energy conservation
ax2 = fig.add_subplot(gs[0, 2])
for (sol, label), color in zip(solutions, colors):
H_traj = H_func(sol.y[0], sol.y[1])
H0 = H_traj[0]
relative_error = np.abs((H_traj - H0) / H0)
ax2.semilogy(sol.t, relative_error, color=color, linewidth=2,
label=label, alpha=0.8)
ax2.set_xlabel('Time t', fontsize=10)
ax2.set_ylabel('|ΔH|/H₀', fontsize=10)
ax2.set_title('Energy conservation', fontsize=11)
ax2.legend(fontsize=8)
ax2.grid(True, alpha=0.3)
# Plot 3: Time evolution (example)
ax3 = fig.add_subplot(gs[1, 2])
sol_example, label_example = solutions[1] # Large oscillation
ax3.plot(sol_example.t, sol_example.y[0], 'b-', linewidth=2, label='x(t)')
ax3.plot(sol_example.t, sol_example.y[1], 'r-', linewidth=2, label='ξ(t)')
ax3.set_xlabel('Time t', fontsize=10)
ax3.set_ylabel('Coordinates', fontsize=10)
ax3.set_title(f'Evolution: {label_example}', fontsize=11)
ax3.legend(fontsize=9)
ax3.grid(True, alpha=0.3)
# Plot 4: Hamiltonian vector field
ax4 = fig.add_subplot(gs[2, 0])
x_vec = np.linspace(-np.pi, 2*np.pi, 20)
xi_vec = np.linspace(-3, 3, 20)
X_vec, XI_vec = np.meshgrid(x_vec, xi_vec, indexing='ij')
U = f_x(X_vec, XI_vec)
V = f_xi(X_vec, XI_vec)
# Normalize for visualization
magnitude = np.sqrt(U**2 + V**2)
U_norm = U / (magnitude + 1e-10)
V_norm = V / (magnitude + 1e-10)
ax4.quiver(X_vec, XI_vec, U_norm, V_norm, magnitude,
cmap='viridis', alpha=0.6, scale=25, width=0.003)
ax4.contour(X, XI, H_grid, levels=15, colors='gray',
alpha=0.3, linewidths=0.5)
ax4.set_xlabel('x', fontsize=10)
ax4.set_ylabel('ξ', fontsize=10)
ax4.set_title('Hamiltonian vector field X_H', fontsize=11)
ax4.set_xlim(-np.pi, 2*np.pi)
ax4.set_ylim(-3, 3)
# Plot 5: Symplectic form (abstract visualization)
ax5 = fig.add_subplot(gs[2, 1])
# Visualize area preserved by the flow
# Take a small square and track its evolution
square_corners = np.array([
[np.pi/4, 0.5],
[np.pi/4 + 0.3, 0.5],
[np.pi/4 + 0.3, 0.7],
[np.pi/4, 0.7],
[np.pi/4, 0.5] # Close loop
])
# Initial area
area_init = 0.3 * 0.2 # Rectangle
# Evolution at different times
times_to_plot = [0, 2, 5, 10]
colors_time = ['red', 'orange', 'yellow', 'green']
for t_plot, c_time in zip(times_to_plot, colors_time):
evolved_corners = []
for corner in square_corners:
sol_corner = solve_ivp(
pendulum_system, [0, t_plot], corner,
method='DOP853', rtol=1e-10
)
evolved_corners.append([sol_corner.y[0, -1], sol_corner.y[1, -1]])
evolved_corners = np.array(evolved_corners)
# Compute area (shoelace formula)
area = 0.5 * np.abs(np.sum(
evolved_corners[:-1, 0] * evolved_corners[1:, 1] -
evolved_corners[1:, 0] * evolved_corners[:-1, 1]
))
ax5.plot(evolved_corners[:, 0], evolved_corners[:, 1],
color=c_time, linewidth=2, label=f't={t_plot}, A={area:.4f}')
ax5.fill(evolved_corners[:, 0], evolved_corners[:, 1],
color=c_time, alpha=0.2)
ax5.contour(X, XI, H_grid, levels=15, colors='gray',
alpha=0.2, linewidths=0.5)
ax5.set_xlabel('x', fontsize=10)
ax5.set_ylabel('ξ', fontsize=10)
ax5.set_title(f'Symplectic area conservation\n(A₀={area_init:.4f})',
fontsize=11)
ax5.legend(fontsize=7)
ax5.grid(True, alpha=0.3)
# Plot 6: Spectrum and quantization (connection)
ax6 = fig.add_subplot(gs[2, 2])
# Classical vs quantum energy levels
E_classical = np.linspace(-2, 2, 50)
# For each energy, compute the action (enclosed area)
# Approximation: period ~ enclosed area
actions_classical = []
for E in E_classical:
if E > -m*g*L: # Above minimum
# Approximation: for small oscillations, action ~ E/ω
omega_0 = np.sqrt(g/L)
action_approx = 2*np.pi*E / omega_0 if E > 0 else 0
actions_classical.append(action_approx)
else:
actions_classical.append(0)
actions_classical = np.array(actions_classical)
ax6.plot(E_classical, actions_classical, 'b-', linewidth=2,
label='Classical action')
# Bohr-Sommerfeld quantization: I = n·ℏ
hbar = 0.3 # Fictitious parameter for visualization
n_levels = np.arange(0, 10)
I_quantum = n_levels * hbar
# Mark quantum levels
for n, I_n in zip(n_levels, I_quantum):
ax6.axhline(I_n, color='red', linestyle='--', alpha=0.5, linewidth=1)
if I_n < max(actions_classical):
ax6.text(1.8, I_n, f'n={n}', fontsize=8, color='red')
ax6.set_xlabel('Energy E', fontsize=10)
ax6.set_ylabel('Action I', fontsize=10)
ax6.set_title('Quantization (Bohr-Sommerfeld)\nI = n·ℏ', fontsize=11)
ax6.legend(fontsize=9)
ax6.grid(True, alpha=0.3)
plt.suptitle('PENDULUM: Complete symplectic geometry',
fontsize=15, fontweight='bold', y=0.995)
plt.savefig('complete_symplectic_pendulum.png', dpi=300, bbox_inches='tight')
plt.show()
# Step 7: Connection with Ψ-DOs
print("\n[STEP 7] Connection with pseudo-differential operators")
print("-"*70)
# Create operator via psiop.py
H_op = PseudoDifferentialOperator(H_pendulum, [x], mode='symbol')
print(f"Hamiltonian operator created via psiop.py")
print(f" Symbol: {H_op.symbol}")
print(f" Order: {H_op.symbol_order()}")
# Hamiltonian flow via operator method
flow = H_op.symplectic_flow()
print(f"\nSymplectic flow (via operator method):")
print(f" dx/dt = {flow['dx/dt']}")
print(f" dξ/dt = {flow['dxi/dt']}")
print("\n✨ The symplectic structure is ENCODED in the operator!")
# Final summary
print("\n" + "="*70)
print("SUMMARY OF THE COMPLETE PIPELINE")
print("="*70)
print()
print("1. GEOMETRY: T*M → ω (canonical symplectic structure)")
print("2. DYNAMICS: H + ω → X_H (Hamiltonian vector field)")
print("3. EVOLUTION: X_H → Φₜ (flow)")
print("4. CONSERVATIONS: Φₜ preserves H and ω (Liouville)")
print("5. QUANTIZATION: {·,·} → [·,·] (correspondence principle)")
print("6. OPERATORS: Symbols → Ψ-DOs (microlocal analysis)")
print()
print("At EVERY STEP, the symplectic structure ω plays the central role,")
print("and it is INDEPENDENT of the particular choice of H.")
print()
print("="*70)
# Execution
complete_symplectic_example_unified()