辛积分法

辛积分法（symplectic integrator）是一类保持哈密顿系统辛结构的数值时间积分方法。它关注位置与动量之间的几何关系，而不只是某一步的数值误差。在长期振动和轨道计算中，保持这种结构常能避免某些持续的人为膨胀或收缩；但辛方法通常并不逐步精确保留原始能量，也不保证任意大步长都稳定。

振子的数值轨迹

考虑无量纲简谐振子，位置为 $q$ ，动量为 $p$ ： $q^{'} = p, p^{'} = - q, H (q, p) = \frac{q^{2} + p^{2}}{2} .$ 沿精确轨迹， $H^{'} = q q^{'} + p p^{'} = q p - p q = 0$ 。初态 $(q, p) = (1, 0)$ 的解为 $(\cos t, - \sin t)$ ，相图是一条圆。

显式 Euler 同时使用旧位置与旧动量： $q_{n + 1} = q_{n} + h p_{n}, p_{n + 1} = p_{n} - h q_{n} .$ 把两式平方相加，交叉项相消，得到 $H_{n + 1} = (1 + h^{2}) H_{n} .$ 无论步长多小，只要 $h > 0$ ，每一步都会增加能量；在固定步长下长期重复会向外螺旋。这里能明确指出问题，不能仅说图画得“不够圆”。

一种只改变更新顺序的方法是先更新动量，再用新动量更新位置： $p_{n + 1} = p_{n} - h q_{n}, q_{n + 1} = q_{n} + h p_{n + 1} .$ 它叫辛 Euler 法的一种形式。其精度仍为一阶，却具有不同的长期结构。对 $h = 0.2$ ，首步为 $p_{1} = - 0.2, q_{1} = 0.96$ ，真实能量 $H_{1} = 0.4808$ ；能量没有精确守恒，后面也不是每步都下降。

显式Euler在相平面向外螺旋，偏离单位圆 — 显式 Euler 取 h=0.2、初态 (1,0)，推进60步。虚线为精确圆轨道；图的坐标范围为正负3.4。

辛Euler在单位圆附近保持有界的相轨迹 — 相同步长、初态与步数的辛 Euler，坐标范围为正负1.15。两图尺度不同，均明确标注；比较的是更新规则的几何性质，而非提高了方法阶数。

辛结构的定义

对 $d$ 个位置与动量，按 $z = (q, p) \in ℝ^{2 d}$ 排列坐标，令 $J = (\begin{matrix} 0 & I_{d} \\ - I_{d} & 0 \end{matrix}), z^{'} = J \nabla H (z) .$ 一个可微的一步映射 $Ψ_{h}$ 称为辛映射，如果其导数矩阵 $M = D Ψ_{h} (z)$ 满足 $M^{𝖳} J M = J .$ 这一条件保留的双线性量是 $u^{𝖳} J v$ ：对两条微小方向先作映射再计算，结果与先计算相同。辛性因此描述整个微分映射，并不只是某条已画出的轨迹。

精确哈密顿流本身是辛的。设其导数矩阵为 $M (t)$ ，沿轨迹的线性化矩阵为 $A = J \nabla^{2} H$ 。Hessian 对称且 $J^{𝖳} = - J, J^{2} = - I$ ，所以 $A^{𝖳} J + J A = 0$ 。由乘积法则， $\frac{d}{d t} (M^{𝖳} J M) = M^{𝖳} (A^{𝖳} J + J A) M = 0.$ 初时刻 $M (0) = I$ ，因而等式始终等于 $J$ 。这里假设 $H$ 有足够连续二阶导数并且相应流存在。辛积分法以离散映射保留的正是这个结构。定义及其与能量的关系见 Hairer、Lubich 与 Wanner 的几何积分论文。^[1]

剪切的组合与面积

对本例，辛 Euler 的矩阵为 $(\begin{matrix} q_{n + 1} \\ p_{n + 1} \end{matrix}) = M_{h} (\begin{matrix} q_{n} \\ p_{n} \end{matrix}), M_{h} = (\begin{matrix} 1 - h^{2} & h \\ - h & 1 \end{matrix}) .$ 其行列式为 $(1 - h^{2}) + h^{2} = 1$ 。二维中直接相乘可验证 $M^{𝖳} J M = (det M) J$ ，所以此处保面积恰好等价于辛性。先做一次动量剪切，再做一次位置剪切，也分别可验证为辛映射，组合仍为辛映射。

高维中仅保总体积则不足够。例如按 $(q_{1}, q_{2}, p_{1}, p_{2})$ 排列的矩阵 $diag (2, 1, 1, 1 / 2)$ 行列式为一，却把 $(q_{1}, p_{1})$ 面积乘以二、把 $(q_{2}, p_{2})$ 面积乘以一半，不满足 $M^{𝖳} J M = J$ 。因此不能用一个行列式检查替代一般辛性检查。

对于可分离 Hamiltonian $H (q, p) = T (p) + V (q)$ ，同样的先动量、后位置更新为 $p_{n + 1} = p_{n} - h \nabla V (q_{n}), q_{n + 1} = q_{n} + h \nabla T (p_{n + 1}) .$ 这两个剪切在相应光滑条件下都是辛的。更对称的半步—整步—半步组合给出斯托默—韦尔莱法等二阶方法。一般非可分离 Hamiltonian 不能未经推导就把上述显式两式照搬过去。

一个真正精确保留的二次量

对谐振子，定义依赖步长的修正量 $\tilde{H} (q, p) = \frac{1}{2} (q^{2} + p^{2} - h q p) .$ 为把抵消过程写清，暂记旧值为 $q, p$ ，新值为 $Q, P$ 。由 $P = p - h q, Q = q + h P$ ， $Q^{2} + P^{2} - h Q P = (q + h P)^{2} + P^{2} - h (q + h P) P = q^{2} + h q P + P^{2} .$ 再代入 $P = p - h q$ ，右端成为 $q^{2} + h q (p - h q) + (p - h q)^{2} = q^{2} + p^{2} - h q p$ 。因此 $q_{n + 1}^{2} + p_{n + 1}^{2} - h q_{n + 1} p_{n + 1} = q_{n}^{2} + p_{n}^{2} - h q_{n} p_{n} .$ 因此 $\tilde{H}$ 逐步精确保留，而原始 $H$ 通常上下振荡。两者的差来自交叉项，并没有把“能量”换个名字就证明原能量相等。

由 $| q p | \leq (q^{2} + p^{2}) / 2 = H$ ，当 $| h | < 2$ 时有 $(1 - | h | / 2) H \leq \tilde{H} \leq (1 + | h | / 2) H .$ 所以修正量正定并控制真实能量，轨迹始终有界。首态 $(1, 0)$ 时 $\tilde{H} = 1 / 2$ ，例如 $h = 0.2$ 给出 $\frac{5}{11} \leq H_{n} \leq \frac{5}{9} .$ 这是可以验证的全步数界，而不是只画若干步得到的经验结论。

辛Euler的原始能量在固定带内振荡，修正二次量保持不变 — 取 h=0.2。原始 H 的界为 5/11 与 5/9，修正量始终为 1/2；曲线显示的三者并不是同一个量。

步长限制仍不可丢。 $M_{h}$ 的迹为 $2 - h^{2}$ ，特征方程为 $r^{2} - (2 - h^{2}) r + 1 = 0$ 。当 $0 < h < 2$ 时，两个根在单位圆上； $h > 2$ 时有一个根的模大于一。 $h = 2$ 也出现非平凡 Jordan 情形，通常导致增长。这个方法对所有步长都是辛的，但并非对所有步长都稳定。

可运行的振子积分

代码只实现本文的无量纲谐振子与先动量、后位置的辛 Euler 法。输入为有限初态、正步长和非负整数步数，返回包含初态的状态二元组列表，不修改输入。每步常数次运算，总时间 $O (N)$ ，保存全轨迹占 $O (N)$ ，滚动推进只需常数额外空间。代码不替用户判断步长稳定性。

Python 3

import math


def symplectic_oscillator(q0, p0, h, steps):
    if not isinstance(steps, int) or steps < 0:
        raise ValueError("steps must be a nonnegative integer")
    q, p, h = map(float, (q0, p0, h))
    if not all(math.isfinite(v) for v in (q, p, h)) or h <= 0:
        raise ValueError("finite state and positive h required")
    states = [(q, p)]
    for _ in range(steps):
        p = p - h * q
        q = q + h * p
        if not math.isfinite(q) or not math.isfinite(p):
            raise OverflowError("intermediate or result outside floating-point range")
        states.append((q, p))
    return states


h = 0.2
states = symplectic_oscillator(1, 0, h, 1000)
q, p = states[1]
print(round(q, 6), round(p, 6))  # 0.96 -0.2
invariant_error = max(abs((q*q + p*p - h*q*p)/2 - 0.5) for q, p in states)
print(invariant_error)  # near roundoff, not exactly zero in floating point

变式取 $h = 0.1$ ，首步应为 $(q_{1}, p_{1}) = (0.99, - 0.1)$ ，原能量为 $0.49505$ ，修正量仍为 $1 / 2$ 。如果先更新位置再把同一行代码中的新位置用于动量，那是另一种辛 Euler 次序，交叉项符号也会改变，不能把两种代码与同一个修正量混用。

长期结论的边界

对一般非线性 Hamiltonian，往往只能在光滑性、小步长、远离奇性等条件下，用修正方程解释长时间较小的能量误差；不能把谐振子精确修正量的计算直接推广成所有 Hamiltonian 都存在同一个精确保留函数。轨道的相位仍会积累误差，非辛方法也可能在短时间或给定容限内更准确。几何数值方法的价值在于与目标问题相匹配的结构保持，而非一种不分场景的精度保证。^[1]

参考资料

↑ ^1.0 ^1.1 Ernst Hairer, Christian Lubich, Gerhard Wanner. “Geometric numerical integration illustrated by the Störmer–Verlet method”. Acta Numerica 12 (2003), 399–450. DOI:10.1017/S0962492902000144。大学课程提供的原论文，§1.2、1.4、2.1—2.2 及 §5.1。

[hairer-1] 1.0 ^1.1 Ernst Hairer, Christian Lubich, Gerhard Wanner. “Geometric numerical integration illustrated by the Störmer–Verlet method”. Acta Numerica 12 (2003), 399–450. DOI:10.1017/S0962492902000144。大学课程提供的原论文，§1.2、1.4、2.1—2.2 及 §5.1。

[1]