哈密顿相空间与辛积分 notebook

这是手机端静态预览：保留说明和代码；图表请在电脑端 Jupyter 中运行 notebook 查看。

??????????

Audience:

• ????????????????????????????????

Prerequisites:

• ?????? p = partial L / partial qdot?
• ??????????? H = p^2/(2m) + kq^2/2?
• ???? Python?NumPy ? Matplotlib ???

Learning goals:

• ?????????????????
• ?????? Euler?RK4?symplectic Euler?Velocity Verlet ??????
• ???????????????????????

Outline

1. ?????????????
2. ?????????????
3. ??????????????
4. ?????????????
5. ???????????????????

from __future__ import annotations

from dataclasses import dataclass
from typing import Callable

import matplotlib.pyplot as plt
import numpy as np

plt.rcParams["figure.figsize"] = (9, 4.8)
plt.rcParams["axes.grid"] = True
plt.rcParams["grid.alpha"] = 0.25

1. ?????

??????????

??????

???????????????????? (qdot, pdot)?

@dataclass(frozen=True)
class Oscillator:
    mass: float = 1.0
    k: float = 1.0


def hamiltonian(q: np.ndarray, p: np.ndarray, system: Oscillator) -> np.ndarray:
    return p**2 / (2 * system.mass) + 0.5 * system.k * q**2


def vector_field(y: np.ndarray, system: Oscillator) -> np.ndarray:
    q, p = y
    q_dot = p / system.mass
    p_dot = -system.k * q
    return np.array([q_dot, p_dot], dtype=float)


system = Oscillator(mass=1.0, k=1.0)
q0, p0 = 1.0, 0.0
hamiltonian(np.array(q0), np.array(p0), system).item()

2. ??????????

????????????????????????????????????

q_grid, p_grid = np.meshgrid(np.linspace(-2.0, 2.0, 25), np.linspace(-2.0, 2.0, 25))
q_dot = p_grid / system.mass
p_dot = -system.k * q_grid
energy_grid = hamiltonian(q_grid, p_grid, system)

fig, ax = plt.subplots(figsize=(6.5, 6.0))
ax.contour(q_grid, p_grid, energy_grid, levels=10, colors="#66706d", alpha=0.55)
ax.streamplot(q_grid, p_grid, q_dot, p_dot, color="#087c75", density=1.1, linewidth=1.1, arrowsize=1.0)
ax.scatter([q0], [p0], color="#c35a2b", s=60, zorder=3, label="initial state")
ax.set_title("Hamiltonian phase flow")
ax.set_xlabel("q")
ax.set_ylabel("p")
ax.set_aspect("equal", adjustable="box")
ax.legend()
plt.show()

3. ??????

?????????

• explicit Euler?????????????????
• RK4???????????????????
• symplectic Euler?????????????
• Velocity Verlet?????????

?????????????????????????????????

StepFn = Callable[[np.ndarray, float, Oscillator], np.ndarray]


def explicit_euler_step(y: np.ndarray, dt: float, system: Oscillator) -> np.ndarray:
    return y + dt * vector_field(y, system)


def rk4_step(y: np.ndarray, dt: float, system: Oscillator) -> np.ndarray:
    k1 = vector_field(y, system)
    k2 = vector_field(y + 0.5 * dt * k1, system)
    k3 = vector_field(y + 0.5 * dt * k2, system)
    k4 = vector_field(y + dt * k3, system)
    return y + (dt / 6.0) * (k1 + 2 * k2 + 2 * k3 + k4)


def symplectic_euler_step(y: np.ndarray, dt: float, system: Oscillator) -> np.ndarray:
    q, p = y
    p_next = p - dt * system.k * q
    q_next = q + dt * p_next / system.mass
    return np.array([q_next, p_next], dtype=float)


def velocity_verlet_step(y: np.ndarray, dt: float, system: Oscillator) -> np.ndarray:
    q, p = y
    p_half = p - 0.5 * dt * system.k * q
    q_next = q + dt * p_half / system.mass
    p_next = p_half - 0.5 * dt * system.k * q_next
    return np.array([q_next, p_next], dtype=float)


def simulate(step_fn: StepFn, y0: np.ndarray, dt: float, steps: int, system: Oscillator) -> dict[str, np.ndarray]:
    y = np.zeros((steps + 1, 2), dtype=float)
    y[0] = y0
    for i in range(steps):
        y[i + 1] = step_fn(y[i], dt, system)
    t = np.arange(steps + 1) * dt
    q = y[:, 0]
    p = y[:, 1]
    return {"t": t, "q": q, "p": p, "H": hamiltonian(q, p, system)}


methods: dict[str, StepFn] = {
    "explicit Euler": explicit_euler_step,
    "RK4": rk4_step,
    "symplectic Euler": symplectic_euler_step,
    "Velocity Verlet": velocity_verlet_step,
}

runs = {name: simulate(fn, np.array([1.0, 0.0]), dt=0.08, steps=2200, system=system) for name, fn in methods.items()}
{name: float(run["H"][-1] - run["H"][0]) for name, run in runs.items()}

4. ?????

?????????????????????????????????????????????????????????

fig, axes = plt.subplots(2, 2, figsize=(10, 8))
for ax, (name, run) in zip(axes.flat, runs.items()):
    ax.plot(run["q"], run["p"], lw=1.5, color="#087c75")
    ax.scatter([run["q"][0]], [run["p"][0]], color="#c35a2b", s=30)
    ax.set_title(name)
    ax.set_xlabel("q")
    ax.set_ylabel("p")
    ax.set_aspect("equal", adjustable="box")
fig.tight_layout()
plt.show()

5. ??????

????????????????????????????????????????????

fig, ax = plt.subplots(figsize=(10, 5))
for name, run in runs.items():
    relative_error = (run["H"] - run["H"][0]) / run["H"][0]
    ax.plot(run["t"], relative_error, label=name, lw=1.8)
ax.set_title("Relative energy error")
ax.set_xlabel("t")
ax.set_ylabel("(H(t)-H(0))/H(0)")
ax.legend()
plt.show()

max_errors = {
    name: float(np.max(np.abs((run["H"] - run["H"][0]) / run["H"][0])))
    for name, run in runs.items()
}
max_errors

6. ???????

???????????????????????????????????????? Euler ? symplectic Euler ???????

def polygon_area(points: np.ndarray) -> float:
    x = points[:, 0]
    y = points[:, 1]
    return 0.5 * abs(np.dot(x, np.roll(y, -1)) - np.dot(y, np.roll(x, -1)))


def evolve_cloud(step_fn: StepFn, cloud: np.ndarray, dt: float, steps: int, system: Oscillator) -> np.ndarray:
    evolved = cloud.copy()
    for _ in range(steps):
        evolved = np.array([step_fn(point, dt, system) for point in evolved])
    return evolved


angles = np.linspace(0, 2 * np.pi, 80, endpoint=False)
cloud = np.column_stack([1.0 + 0.12 * np.cos(angles), 0.0 + 0.12 * np.sin(angles)])
initial_area = polygon_area(cloud)

euler_cloud = evolve_cloud(explicit_euler_step, cloud, dt=0.08, steps=300, system=system)
symp_cloud = evolve_cloud(symplectic_euler_step, cloud, dt=0.08, steps=300, system=system)

areas = {
    "initial": initial_area,
    "explicit Euler": polygon_area(euler_cloud),
    "symplectic Euler": polygon_area(symp_cloud),
}
areas

fig, axes = plt.subplots(1, 2, figsize=(10, 4.6))
for ax, title, evolved in [
    (axes[0], "explicit Euler", euler_cloud),
    (axes[1], "symplectic Euler", symp_cloud),
]:
    ax.plot(cloud[:, 0], cloud[:, 1], color="#66706d", lw=1.5, label="initial")
    ax.plot(evolved[:, 0], evolved[:, 1], color="#c35a2b", lw=1.8, label="after evolution")
    ax.set_title(title)
    ax.set_xlabel("q")
    ax.set_ylabel("p")
    ax.set_aspect("equal", adjustable="box")
    ax.legend()
fig.tight_layout()
plt.show()

??

1. ? dt ? 0.08 ?? 0.2????????????????
2. ?????? (q0,p0)=(0,1.5)????????????????
3. ? mass=2 ? k=4?????????????
4. ???????????RK4 ?????RK4 ???????????

# Exercise scaffold: change dt and compare final relative energy error.
def compare_dt(dt: float) -> dict[str, float]:
    local_runs = {
        name: simulate(fn, np.array([1.0, 0.0]), dt=dt, steps=int(120 / dt), system=system)
        for name, fn in methods.items()
    }
    return {
        name: float((run["H"][-1] - run["H"][0]) / run["H"][0])
        for name, run in local_runs.items()
    }

compare_dt(0.2)

??

????????????H ?????????????????????????????????????

• ???qdot = partial H / partial p?pdot = -partial H / partial q?
• ??????????????????
• ?????????????????