李群数值积分器 (二) · LGVI 与离散变分力学流派：Hamilton 变分原理、辛几何保积性与长时域守恒

Lie Group Numerical Integrators II: Lie Group Variational Integrators (LGVI), Discrete Euler-Poincaré & Symplectic Momentum Invariance
基于 Taeyoung Lee, Melvin Leok, N. Harris McClamroch (CMAME 2007, IEEE CDC, AIAA JGCD) 经典奠基之作

---

💡 导言：走向辛力学与变分积分的皇冠
在前篇《李群数值积分器 (一) · RKMK 篇》中，我们系统见证了代数回拉（Lie-algebraic Pull-back）的数学魔力：通过指数逆微分 $\mathrm{d}{\exp }^{-1}$ 算子，RKMK 完美实现了在李群流形上的任意高阶推进。然而，我们同时也指出了 RKMK 的致命死穴——保李群流形几何并不等于保相空间辛几何！在深空引力场、卫星自由翻滚以及天体轨道模拟中，非辛的 RKMK 无论阶数多高，都会因为非线性坐标拉回而在相空间不可逆地产生人造数值阻尼，导致能量持续单调泄漏。
为了彻底攻克这一百年力学难题，密歇根大学与加利福尼亚大学的 Taeyoung Lee、Melvin Leok 与 N. Harris McClamroch 于 2007 年在计算力学顶级期刊 CMAME 发表了划时代的里程碑长文，创立了 LGVI (Lie Group Variational Integrators，李群变分积分器) 理论体系。
LGVI 彻底摒弃了“离散化连续微分方程”的传统思维，而是开辟了**“直接在李群上离散化哈密顿最小作用量原理”的崭新道路。本篇专论将带您深入离散变分力学的微积分底层，推导流形上的离散欧拉-庞加莱方程 (DEP)、离散诺特定理与阴影哈密顿量**，曝光高阶经典 RK4 在全刚体引力轨道模拟中引发的“坠毁惨剧”，并揭示 LGVI 对现代空间几何最优控制 (DMOC) 的终极反哺！

---

---

目录

• 1. 经典微分力学的破局：从离散微分方程到离散变分原理
  • 1.1 传统离散微分方程的结构缺陷与能量耗散
  • 1.2 变分力学思想源头：哈密顿最小作用量原理
• 2. 连续与离散运动方程的系统推导闭环 (The Commutative Cube)
  • 2.1 构型空间 SE(3) 的几何半直积与动能构造
  • 2.2 流形构型变分 δR 与两步相对增量变分 δF_k 恒等式
  • 2.3 连续与离散力学的八面体可交换立方体
• 3. 李群变分积分器 (LGVI) 的核心方程组
  • 3.1 自伴随离散拉格朗日量与离散勒让德变换
  • 3.2 隐式姿态增量方程的罗德里格公式向量化与牛顿解算
  • 3.3 离散哈密顿形式全状态步进映射
• 4. 辛保积性、离散诺特定理与长时域零能量耗散机理
  • 4.1 相空间辛两形式的机器精度守恒
  • 4.2 离散诺特定理与空间总动量矩严格守恒
  • 4.3 后向误差分析与阴影哈密顿量 (Shadow Hamiltonian)
  • 4.4 Yoshida 对称时间步长组合与高阶拓展
• 5. 全刚体两哑铃系统 (Two-Dumbbell Problem) 基准实验全景剖析
  • 5.1 全刚体两哑铃系统的引力势能与轨道-姿态强耦合
  • 5.2 长时域仿真中的轨道坠毁惨剧：RK4 持续能量泄漏
  • 5.3 LGVI 千年稳态闭合与动量机器精度守恒验证
• 6. LGVI 与第一篇 RKMK 的全方位横向决选
• 7. LGVI 对 Chu 2026 几何最优控制与 DMOC 的技术反哺
• 8. 完整 Python / NumPy 可运行工业级实现
  • 8.1 paper_lgvi.py 姿态隐式增量向量化牛顿解算器
  • 8.2 integrators.py 相对哈密顿 LGVI 单步推进内核
• 9. 权威参考文献与延伸阅读

---

1. 经典微分力学的破局：从离散微分方程到离散变分原理

1.1 传统离散微分方程的结构缺陷与能量耗散

在经典计算力学中，几乎所有的标准常微分方程数值求解器（包括欧拉法、Adams 多步法、经典 Runge-Kutta 乃至前篇介绍的 RKMK）都遵循如下既定工作路线：

graph LR
    A["连续拉格朗日作用量 $\mathfrak{G} = \int L \mathrm{d}t$"] -->|"连续 Hamilton 原理 $\delta \mathfrak{G} = 0$"| B["连续微分方程 $\dot{x} = f(x)$"]
    B -->|"数值差分截断 (RK4, 多步法)"| C["❌ 离散轨迹 (破坏辛结构与动量守恒)"]

这一传统范式存在一个根本性的力学缺陷：力学系统深层的对称性、诺特定理守恒量以及相空间辛两形式，是变分原理自身的内在几何属性，而不是微分方程孤立数值形式的属性。当我们在连续微分方程上施加差分近似时，哪怕局部截断误差高达 $\mathcal{O}(h^4)$ 或 $\mathcal{O}(h^6)$，系统的相空间微元体积与动量曲面依然遭到破坏。在长时间推演中，这些微小的破损会无情地演化为人造能量耗散或发散。

1.2 变分力学思想源头：哈密顿最小作用量原理

变分积分器（Variational Integrators，起源于 Moser、Veselov 与 Marsden、West 2001）采取了颠覆性的**“离散先行”哲学**：

graph LR
    A["连续拉格朗日量 $L(q, \dot{q})$"] -->|"离散近似 $L_d(q_k, q_{k+1}) \approx \int_0^h L \mathrm{d}t$"| D["离散作用和 $\mathfrak{G}_d = \sum L_d$"]
    D -->|"离散 Hamilton 原理 $\delta \mathfrak{G}_d = 0$"| E["✅ 离散变分积分器 (机器精度保辛与守恒)"]

变分积分器根本不去直接离散连续微分方程，而是首先将时间轴划分为离散节点 $t_k=kh$，将连续作用量积分近似为离散拉格朗日作用和 (Discrete Action Sum)：

$$\begin{array}{r}{\mathfrak{G}}_d\triangleq \sum _{k=0}^{N-1}{L}_d(q_k,q_{k+1})\approx {\int }_{t_0}^{t_f}L(q,\dot{q})\mathrm{d}t\end{array}\text{\hspace{1em}(1-1)}$$

然后直接在离散序列上施加离散变分原理：

$$\begin{array}{r}\delta {\mathfrak{G}}_d=\delta \sum _{k=0}^{N-1}{L}_d(q_k,q_{k+1})=0\end{array}\text{\hspace{1em}(1-2)}$$

由式 (1-2) 导出的差分方程被称为离散欧拉-拉格朗日方程 (Discrete Euler-Lagrange, DEL)。由于数值方程是从纯粹的变分作用量驻值条件中“原汁原味”诞生出来的，它在离散层面完整继承了哈密顿力学的全部辛几何结构！

---

2. 连续与离散运动方程的系统推导闭环 (The Commutative Cube)

2.1 构型空间 SE(3) 的几何半直积与动能构造

对于空间中的每一个三维刚体 $\mathcal{B}$，其构型空间为特殊欧氏群 $\mathrm{SE}(3)\triangleq {\mathbb{R}}^3⋊\mathrm{SO}(3)$。质心在惯性系中的位置记为 $\mathbf{x}\in {\mathbb{R}}^3$，姿态由旋转矩阵 $\mathbf{R}\in \mathrm{SO}(3)$ 描述（体坐标系到惯性系的变换）。

刚体 $\mathcal{B}$ 的动能可展开为平动动能与转动动能之和：

$$\begin{array}{r}T=\frac{1}{2}m\Vert \dot{\mathbf{x}}{\Vert }_2^2+\frac{1}{2}\mathrm{tr}\left[\hat{\mathbf{\Omega }}{\mathbf{J}}_d{\hat{\mathbf{\Omega }}}^{\top }\right]\end{array}\text{\hspace{1em}(2-1)}$$

其中 $\mathbf{\Omega }\in {\mathbb{R}}^3$ 为体坐标系下的角速度（$\dot{\mathbf{R}}=\mathbf{R}\hat{\mathbf{\Omega }}$），$\hat{\mathbf{\Omega }}\triangleq S(\mathbf{\Omega })\in \mathfrak{so}(3)$ 为反对称叉积矩阵。

💡 标准惯量矩阵与非标准惯量矩阵的代数转换：
式中 ${\mathbf{J}}_d\triangleq {\int }_{\mathcal{B}}\mathbf{\rho }{\mathbf{\rho }}^{\top }\mathrm{d}m\in {\mathbb{R}}^{3\times 3}$ 为非标准转动惯量矩阵。它与物理教科书中的标准转动惯量矩阵 $\mathbf{J}\triangleq {\int }_{\mathcal{B}}S(\mathbf{\rho }{)}^{\top }S(\mathbf{\rho })\mathrm{d}m$ 满足极其优美的迹对偶关系：

$$\begin{array}{r}\mathbf{J}=\mathrm{tr}[{\mathbf{J}}_d]{\mathbf{I}}_3-{\mathbf{J}}_d,S(\mathbf{J}\mathbf{\Omega })=S(\mathbf{\Omega }){\mathbf{J}}_d+{\mathbf{J}}_{d}S(\mathbf{\Omega })\end{array}\text{\hspace{1em}(2-2)}$$

结合引力势能 $U(\mathbf{x},\mathbf{R})$，全系统的拉格朗日量为 $L=T-U$。

2.2 流形构型变分 δR 与两步相对增量变分 δF_k 恒等式

在李群流形上推导变分时，必须保证微扰始终严格切于流形：

• 连续时间姿态变分：$\delta \mathbf{R}=\mathbf{R}\mathbf{\eta }$，其中 $\mathbf{\eta }\in \mathfrak{so}(3)$ 在时间边界端点处为零；
• 角速度变分：对运动学 $\dot{\mathbf{R}}=\mathbf{R}\hat{\mathbf{\Omega }}$ 取变分，可得： $$\begin{array}{r}\widehat{\delta \mathbf{\Omega }}=\dot{\mathbf{\eta }}+[\hat{\mathbf{\Omega }},\mathbf{\eta }]=\dot{\mathbf{\eta }}+\hat{\mathbf{\Omega }}\mathbf{\eta }-\mathbf{\eta }\hat{\mathbf{\Omega }}\end{array}\text{\hspace{1em}(2-3)}$$

在离散时间步进中，定义两连续时间步之间的相对姿态变换矩阵：

$$\begin{array}{r}{\mathbf{F}}_k\triangleq {\mathbf{R}}_k^{\top }{\mathbf{R}}_{k+1}\in \mathrm{SO}(3)⟹{\mathbf{R}}_{k+1}={\mathbf{R}}_k{\mathbf{F}}_k\end{array}\text{\hspace{1em}(2-4)}$$

🏆 离散相对步进变分基本引理 (Taeyoung Lee 2007)：
对式 (2-4) 两端同时施加流形变分 $\delta {\mathbf{R}}_k={\mathbf{R}}_k{\mathbf{\eta }}_k$ 与 $\delta {\mathbf{R}}_{k+1}={\mathbf{R}}_{k+1}{\mathbf{\eta }}_{k+1}$：

$$\begin{array}{rl}\delta {\mathbf{F}}_k& =\delta ({\mathbf{R}}_k^{\top }{\mathbf{R}}_{k+1})=\delta {\mathbf{R}}_k^{\top }{\mathbf{R}}_{k+1}+{\mathbf{R}}_k^{\top }\delta {\mathbf{R}}_{k+1}\\ & =-{\mathbf{\eta }}_k{\mathbf{R}}_k^{\top }{\mathbf{R}}_{k+1}+{\mathbf{R}}_k^{\top }{\mathbf{R}}_{k+1}{\mathbf{\eta }}_{k+1}\\ & =-{\mathbf{\eta }}_k{\mathbf{F}}_k+{\mathbf{F}}_k{\mathbf{\eta }}_{k+1}\end{array}\text{\hspace{1em}(2-5)}$$

式 (2-5) 结构高度对称，它正是离散流形变分能够像连续微积分一样实施“分部求和（Summation by parts）”的核心代数钥匙！

2.3 连续与离散力学的八面体可交换立方体

Taeyoung Lee 等人在论文中系统归纳了连续/离散、拉格朗日/哈密顿、惯性系/相对系三维正交属性所构成的八套运动方程体系，其逻辑拓扑由如图所示的可交换立方体统一概括：

图 2.3.1：Lee 2007 原论文图 2，从连续 Hamilton 原理到离散 Hamilton 原理的系统化推导流线图。

图 2.3.2：原论文图 7，全身问题八套运动方程在“连续/离散、惯性/相对、拉格朗日/哈密顿”坐标轴下的可交换立方体拓扑。

---

3. 李群变分积分器 (LGVI) 的核心方程组

3.1 自伴随离散拉格朗日量与离散勒让德变换

利用中点辛规则，将线速度近似为 ${\dot{\mathbf{x}}}_k\approx ({\mathbf{x}}_{k+1}-{\mathbf{x}}_k)/h$，角速度近似为 ${\hat{\mathbf{\Omega }}}_k\approx \frac{1}{h}({\mathbf{F}}_k-{\mathbf{I}}_3)$。代入动能后，构造自伴随离散拉格朗日量 (Self-Adjoint Discrete Lagrangian)：

$$\begin{array}{r}{L}_d({\mathbf{x}}_k,{\mathbf{x}}_{k+1},{\mathbf{R}}_k,{\mathbf{F}}_k)\triangleq \frac{m}{2h}\Vert {\mathbf{x}}_{k+1}-{\mathbf{x}}_k{\Vert }_2^2+\frac{1}{h}\mathrm{tr}\left[({\mathbf{I}}_3-{\mathbf{F}}_k){\mathbf{J}}_d\right]-\frac{h}{2}U({\mathbf{x}}_k,{\mathbf{R}}_k)-\frac{h}{2}U({\mathbf{x}}_{k+1},{\mathbf{R}}_{k+1})\end{array}\text{\hspace{1em}(3-1)}$$

定义离散纤维导数（离散勒让德变换）：

$$\begin{array}{r}{\mathbf{\gamma }}_k=-D_{{\mathbf{x}}_k}{L}_d,S({\mathbf{\Pi }}_k)=-D_{{\mathbf{R}}_k}{L}_d+D_{{\mathbf{R}}_k}{L}_d^{\top }\end{array}\text{\hspace{1em}(3-2)}$$

对式 (3-1) 取变分并整理，获得线动量 ${\mathbf{\gamma }}_k=m{\mathbf{v}}_k$ 与体坐标系角动量 ${\mathbf{\Pi }}_k=\mathbf{J}{\mathbf{\Omega }}_k$ 的离散关联。

3.2 隐式姿态增量方程的罗德里格公式向量化与牛顿解算

通过对离散作用和实施变分并应用引理 (2-5)，离散变分原理 $\delta \sum L_d=0$ 最终导出了决定下一时步姿态相对旋转 ${\mathbf{F}}_k\in \mathrm{SO}(3)$ 的核心隐式矩阵方程：

$$\begin{array}{r}hS\left({\mathbf{\Pi }}_k+\frac{h}{2}{\mathbf{M}}_k\right)={\mathbf{F}}_k{\mathbf{J}}_d-{\mathbf{J}}_d{\mathbf{F}}_k^{\top }\end{array}\text{\hspace{1em}(3-3)}$$

其中 ${\mathbf{M}}_k\in {\mathbb{R}}^3$ 为外力势能诱导的力矩（$S({\mathbf{M}}_k)=\frac{\partial U_k^{\top }}{\partial {\mathbf{R}}_k}{\mathbf{R}}_k-{\mathbf{R}}_k^{\top }\frac{\partial U_k}{\partial {\mathbf{R}}_k}$）。

⚙️ 极速罗德里格向量化解法 (Rodrigues Newton Solve)：
令已知常数向量为 $\mathbf{g}\triangleq h({\mathbf{\Pi }}_k+\frac{h}{2}{\mathbf{M}}_k)\in {\mathbb{R}}^3$。将待求解的未知旋转矩阵参数化为李代数向量 $\mathbf{f}\in {\mathbb{R}}^3$：

$$\begin{array}{r}{\mathbf{F}}_k=\exp (S(\mathbf{f}))={\mathbf{I}}_3+\frac{\sin \Vert \mathbf{f}\Vert }{\Vert \mathbf{f}\Vert }S(\mathbf{f})+\frac{1-\cos \Vert \mathbf{f}\Vert }{\Vert \mathbf{f}{\Vert }^2}S(\mathbf{f}{)}^2\end{array}\text{\hspace{1em}(3-4)}$$

将式 (3-4) 代入矩阵方程 (3-3)，利用等式 $S(\mathbf{f}){\mathbf{J}}_d+{\mathbf{J}}_{d}S(\mathbf{f})=S(\mathbf{J}\mathbf{f})$，矩阵方程惊人地化简为如下平坦三维向量非线性方程：

$$\begin{array}{r}\mathbf{g}=\mathbf{G}(\mathbf{f})\triangleq \frac{\sin \Vert \mathbf{f}\Vert }{\Vert \mathbf{f}\Vert }\mathbf{J}\mathbf{f}+\frac{1-\cos \Vert \mathbf{f}\Vert }{\Vert \mathbf{f}{\Vert }^2}\mathbf{f}\times (\mathbf{J}\mathbf{f})\end{array}\text{\hspace{1em}(3-5)}$$

方程 (3-5) 的解析雅可比矩阵 $\nabla \mathbf{G}(\mathbf{f})\in {\mathbb{R}}^{3\times 3}$ 具有完全解析表达式：

$$\begin{array}{rl}\nabla \mathbf{G}(\mathbf{f})& =\frac{\cos \Vert \mathbf{f}\Vert \Vert \mathbf{f}\Vert -\sin \Vert \mathbf{f}\Vert }{\Vert \mathbf{f}{\Vert }^3}(\mathbf{J}\mathbf{f}){\mathbf{f}}^{\top }+\frac{\sin \Vert \mathbf{f}\Vert }{\Vert \mathbf{f}\Vert }\mathbf{J}\\ & +\frac{\sin \Vert \mathbf{f}\Vert \Vert \mathbf{f}\Vert -2(1-\cos \Vert \mathbf{f}\Vert )}{\Vert \mathbf{f}{\Vert }^4}(\mathbf{f}\times \mathbf{J}\mathbf{f}){\mathbf{f}}^{\top }+\frac{1-\cos \Vert \mathbf{f}\Vert }{\Vert \mathbf{f}{\Vert }^2}\left(-S(\mathbf{J}\mathbf{f})+S(\mathbf{f})\mathbf{J}\right)\end{array}\text{\hspace{1em}(3-6)}$$

以 ${\mathbf{f}}_0={\mathbf{J}}^{-1}\mathbf{g}$ 作为初始猜测，标准 Newton-Raphson 迭代 ${\mathbf{f}}_{i+1}={\mathbf{f}}_i+[\nabla \mathbf{G}({\mathbf{f}}_i){]}^{-1}(\mathbf{g}-\mathbf{G}({\mathbf{f}}_i))$ 仅需 3 次迭代即可收敛至 $10^{-15}$ 机器浮点精度！

3.3 离散哈密顿形式全状态步进映射

一旦求得局部相对增量 ${\mathbf{F}}_k=\exp (S(\mathbf{f}))\in \mathrm{SO}(3)$，全系统状态按如下闭式顺序显式更新：

$$\begin{array}{r}\{\begin{array}{l}{\mathbf{x}}_{k+1}={\mathbf{x}}_k+\frac{h}{m}{\mathbf{\gamma }}_k-\frac{h^2}{2m}\frac{\partial U_k}{\partial {\mathbf{x}}_k}\\ {\mathbf{R}}_{k+1}={\mathbf{R}}_k{\mathbf{F}}_k\\ {\mathbf{\gamma }}_{k+1}={\mathbf{\gamma }}_k-\frac{h}{2}\frac{\partial U_k}{\partial {\mathbf{x}}_k}-\frac{h}{2}\frac{\partial U_{k+1}}{\partial {\mathbf{x}}_{k+1}}\\ {\mathbf{\Pi }}_{k+1}={\mathbf{F}}_k^{\top }{\mathbf{\Pi }}_k+\frac{h}{2}{\mathbf{F}}_k^{\top }{\mathbf{M}}_k+\frac{h}{2}{\mathbf{M}}_{k+1}\end{array}\end{array}\text{\hspace{1em}(3-7)}$$

这构成了严格保流形且保辛的离散哈密顿流映射 $({\mathbf{x}}_k,{\mathbf{\gamma }}_k,{\mathbf{R}}_k,{\mathbf{\Pi }}_k)\mapsto ({\mathbf{x}}_{k+1},{\mathbf{\gamma }}_{k+1},{\mathbf{R}}_{k+1},{\mathbf{\Pi }}_{k+1})$。

---

4. 辛保积性、离散诺特定理与长时域零能量耗散机理

4.1 相空间辛两形式的机器精度守恒

对于离散系统，定义相空间 $T^{\ast }Q$ 上的离散辛形式 ${\omega }_{L_d}\triangleq -\mathrm{d}{\theta }_{L_d}$。
根据变分力学的基本定理（Marsden & West 2001）：

$$\begin{array}{r}{F}_{L_d}^{\ast }{\omega }_{L_d}={\omega }_{L_d}⟺\mathrm{d}{\mathbf{x}}_{k+1}\wedge \mathrm{d}{\mathbf{\gamma }}_{k+1}+\mathrm{d}{\mathbf{R}}_{k+1}\wedge \mathrm{d}{\mathbf{\Pi }}_{k+1}=\mathrm{d}{\mathbf{x}}_k\wedge \mathrm{d}{\mathbf{\gamma }}_k+\mathrm{d}{\mathbf{R}}_k\wedge \mathrm{d}{\mathbf{\Pi }}_k\end{array}\text{\hspace{1em}(4-1)}$$

这证明 LGVI 在离散层面上严格保持相空间微分体积不变，绝对不存在欧氏算法中由加法引起的相空间人工膨胀或收缩。

4.2 离散诺特定理与空间总动量矩严格守恒

引力势能 $U(\mathbf{x},\mathbf{R})$ 仅取决于天体之间的相对几何构型，具有全空间刚体变换群 $\mathrm{SE}(3)$ 的整体平移与旋转不变性。

根据离散诺特定理 (Discrete Noether’s Theorem)，这一连续对称性在离散拉格朗日量中完好保留，直接诱导出：

1. 系统惯性系总线动量严格守恒： $$\begin{array}{r}{\mathbf{\gamma }}_T\triangleq \sum _{i=1}^n{\mathbf{\gamma }}_{i,k}=\text{const}(\text{误差}\le 10^{-15})\end{array}\text{\hspace{1em}(4-2)}$$
2. 系统空间总角动量矩严格守恒： $$\begin{array}{r}{\mathbf{\pi }}_T\triangleq \sum _{i=1}^n\left({\mathbf{x}}_{i,k}\times {\mathbf{\gamma }}_{i,k}+{\mathbf{R}}_{i,k}{\mathbf{\Pi }}_{i,k}\right)=\text{const}(\text{误差}\le 10^{-15})\end{array}\text{\hspace{1em}(4-3)}$$

4.3 后向误差分析与阴影哈密顿量 (Shadow Hamiltonian)

许多读者会困惑：既然 LGVI 是辛积分器，那它是不是每一步都严格保证系统总能量 $E$ 恒等于初值？

严格的数学答案是：并不严格守恒，但其能量误差具有惊人的界限！ 根据后向误差分析（Backward Error Analysis，Hairer et al. 2006）：

• LGVI 的离散序列并不是原哈密顿系统 $H$ 的带漂移近似解；
• 它严格等价于某个极其贴近原系统的‘修正哈密顿系统（阴影哈密顿量 Shadow Hamiltonian）’ $\tilde{H}=H+h^2{H}_2+h^4{H}_4+\dots$ 的真解！
• 只要步长 $h$ 位于分析收敛半径内，系统的真实能量 $E(t)$ 将被严格约束在初值附近以微小的振幅（$\mathcal{O}(h^2)$）做零均值高频振荡，历经数万个轨道周期绝对不发生任何单调能量漂移！

4.4 Yoshida 对称时间步长组合与高阶拓展

由于离散拉格朗日量式 (3-1) 是自伴随的（时间反转不变），其固有精度阶数恒为偶数（至少 2 阶）。
利用 Yoshida (1990) 对称步长组合技术，可以将基础二阶辛映射 ${\Phi }_h$ 串联构造出严格四阶的辛积分器：

$$\begin{array}{r}{\Psi }_h={\Phi }_{{\lambda }_{1}h}\circ {\Phi }_{{\lambda }_{2}h}\circ {\Phi }_{{\lambda }_{1}h}\end{array}\text{\hspace{1em}(4-4)}$$

其中组合权重满足代数方程 $2{\lambda }_1+{\lambda }_2=1,2{\lambda }_1^3+{\lambda }_2^3=0$，其解析解为：

$$\begin{array}{r}{\lambda }_1=\frac{1}{2-2^{1/3}}\approx 1.3512,{\lambda }_2=-\frac{2^{1/3}}{2-2^{1/3}}\approx -1.7024\end{array}\text{\hspace{1em}(4-5)}$$

通过中间引入的一个精妙负时间子步 ${\lambda }_2<0$，系统在保持辛保积性与动量守恒的同时，局部精度阶数跃升至严格四阶！

---

5. 全刚体两哑铃系统 (Two-Dumbbell Problem) 基准实验全景剖析

5.1 全刚体两哑铃系统的引力势能与轨道-姿态强耦合

为了最严苛地检验数值积分器的长时域保结构性能，论文设计了双哑铃全刚体引力耦合动力学（Full Body Two-Dumbbell Problem）。每个哑铃由两个相距为 $l_i$ 的刚性球体与无质量刚性连接杆构成：

图 5.1.1：Lee 2007 原论文图 3，全刚体引力耦合双哑铃天体力学物理模型。

全系统的归一化万有引力势能为四个球体之间的引力势叠加：

$$\begin{array}{r}U(\mathbf{X},\mathbf{R})=-\frac{1}{4}\sum _{p=1}^2\sum _{q=1}^2\frac{1}{\Vert \mathbf{X}+{\mathbf{\rho }}_{2p}+\mathbf{R}{\mathbf{\rho }}_{1q}{\Vert }_2}\end{array}\text{\hspace{1em}(5-1)}$$

势能对相对姿态 $\mathbf{R}$ 与相对位置 $\mathbf{X}$ 的非线性偏导数直接产生相互耦合的引力梯度力矩与轨道摄动力。

5.2 长时域仿真中的轨道坠毁惨剧：RK4 持续能量泄漏

在相同的物理参数与初始大离心率倾斜轨道下，对比四阶经典 Runge-Kutta (RK4) 与二阶 LGVI 的长时域表现：

图 5.2.1：原论文图 5，LGVI 算法下两个哑铃体在惯性系中长期闭合的稳定轨道与平稳姿态演化。

图 5.2.2：原论文图 6，经典高阶 RK4 算法下由于能量持续单调泄漏导致的轨道快速内卷螺旋与最终灾难性相撞坠毁！

图 5.2.3：本地工程精确复现生成的 RK4 能量演化曲线，系统总能量发生不可逆的线性雪崩式暴跌。

图 5.2.4：本地工程复现生成的 LGVI 能量演化曲线，总能量在无偏水平基准线上做严格有界的微幅高频弹性振荡。

图 5.2.5：本地复现总动量矩守恒性对比：LGVI 的总动量矩在机器浮点精度级（$10^{-14}$）恒定不变，而 RK4 的角动量彻底发散漂移。

图 5.2.6：LGVI 算法理论 2 阶收敛斜率验证。

---

6. LGVI 与第一篇 RKMK 的全方位横向决选

至此，我们可以对李群几何积分领域的两大开山流派给出决定性的权威裁决：

评测维度	第一篇 · RKMK (代数回拉流派)	第二篇 · LGVI (离散变分力学流派)
理论出发点	微分方程切向量回拉 $\mathrm{d}{\exp }^{-1}$	作用量积分离散变分 $\delta \sum L_d=0$
算法求解属性	纯显式单步推进（免非线性求解）	隐式步进（每步需解三维牛顿方程）
流形正交性保真度	机器精度严格保持（$\mathbf{R}\in \mathrm{SO}(3)$）	机器精度严格保持（$\mathbf{R}\in \mathrm{SO}(3)$）
相空间辛形式守恒	❌ 不保辛（代数回拉破坏辛两形式）	✅ 严格辛保积（$F_{L_d}^{\ast }\omega =\omega$）
第一积分动量守恒	一般不守恒（产生数值漂移）	✅ 离散诺特定理机器精度守恒
长时域能量演化	单调人造耗散/发散（导致轨道塌缩）	零长期漂移（阴影哈密顿量高频振荡）
单步计算开销	极低（纯矩阵乘法与指数映射）	中等（需 3 次局部牛顿迭代解 $g=G(f)$）
最优适用工程场景	机器人轨迹优化与短时域 NMPC 闭环	深空天体探测、卫星长期自由伴飞与轨道推演

---

7. LGVI 对 Chu 2026 几何最优控制与 DMOC 的技术反哺

变分力学在最优控制领域的最高结晶是 离散力学与最优控制 (Discrete Mechanics and Optimal Control, DMOC)（Kobilarov 2010, Leok 2007）：

1. 伴随系统的辛对偶一致性：将 LGVI 导出的离散欧拉-庞加莱方程作为约束，构建拉格朗日乘子系统时，对偶伴随方程在离散层面上严格满足离散辛差分格式，杜绝了“连续求导再离散”导致的伴随漂移；
2. 为 Chu 2026 提供空间伴飞真值与安全边界：在空间非合作翻滚航天器交会的远距保持（Hold）与安全停靠走廊中，航天器处于长时间无控或微推力漂移状态。LGVI 提供了绝不失真的保辛物理基准，彻底消除了由于数值假阻尼导致的虚假冲量消耗与避障误判！

---

8. 完整 Python / NumPy 可运行工业级实现

以下为依循 Taeyoung Lee (CMAME 2007) 理论方程与本地工程库提炼的高效、轻量化工业级 LGVI 实现：

8.1 paper_lgvi.py 姿态隐式增量向量化牛顿解算器

import numpy as np

def hat(w: np.ndarray) -> np.ndarray:
    """三维向量到反对称矩阵映射"""
    return np.array([
        [0.0, -w[2], w[1]],
        [w[2], 0.0, -w[0]],
        [-w[1], w[0], 0.0]
    ])

def exp_so3(w: np.ndarray) -> np.ndarray:
    """Rodrigues 旋转矩阵指数映射"""
    th = np.linalg.norm(w)
    if th < 1e-12:
        return np.eye(3) + hat(w)
    K = hat(w) / th
    return np.eye(3) + np.sin(th)*K + (1.0 - np.cos(th))*(K @ K)

def G_map(f: np.ndarray, J: np.ndarray) -> np.ndarray:
    """向量形式隐式姿态方程 G(f, J) (Lee 2007, Eq. 130)"""
    th = np.linalg.norm(f)
    Jf = J @ f
    if th < 1e-10:
        a = 1.0 - th**2 / 6.0
        b = 0.5 - th**2 / 24.0
    else:
        a = np.sin(th) / th
        b = (1.0 - np.cos(th)) / (th**2)
    return a * Jf + b * np.cross(f, Jf)

def jac_G(f: np.ndarray, J: np.ndarray) -> np.ndarray:
    """G_map 的精确解析雅可比矩阵 (Lee 2007, Eq. 131)"""
    th = np.linalg.norm(f)
    Jf = J @ f
    if th < 1e-7:
        eps = 1e-7
        A = np.zeros((3, 3))
        for i in range(3):
            e = np.zeros(3); e[i] = eps
            A[:, i] = (G_map(f + e, J) - G_map(f - e, J)) / (2 * eps)
        return A
    
    term1 = ((np.cos(th)*th - np.sin(th))/(th**3)) * np.outer(Jf, f)
    term2 = (np.sin(th)/th) * J
    term3 = ((np.sin(th)*th - 2.0*(1.0 - np.cos(th)))/(th**4)) * np.outer(np.cross(f, Jf), f)
    term4 = ((1.0 - np.cos(th))/(th**2)) * (-hat(Jf) + hat(f) @ J)
    return term1 + term2 + term3 + term4

def solve_increment(g: np.ndarray, J: np.ndarray, tol: float = 1e-14, max_iter: int = 15):
    """三维牛顿-拉夫逊极速求解 g = G(f)"""
    f = np.linalg.solve(J, g)
    for it in range(max_iter):
        res = g - G_map(f, J)
        if np.linalg.norm(res) < tol:
            return exp_so3(f), f, it
        A = jac_G(f, J)
        f += np.linalg.solve(A, res)
    return exp_so3(f), f, max_iter

8.2 integrators.py 相对哈密顿 LGVI 单步推进内核

def lgvi_rigid_body_step(R_k: np.ndarray, Pi_k: np.ndarray, h: float, J: np.ndarray, torque_func, t_k: float):
    """
    单刚体在力矩作用下的经典哈密顿 LGVI 单步推进 (Lee 2007, Eqs. 93-97)
    
    参数:
        R_k: 当前姿态矩阵 in SO(3)
        Pi_k: 当前体坐标系角动量 J * Omega in R^3
        h: 积分步长
        J: 标准转动惯量矩阵
        torque_func: 力矩函数 M(t, R)
    """
    M_k = torque_func(t_k, R_k)
    
    # 1. 组装已知向量 g = h * (Pi_k + 0.5 * h * M_k)
    g = h * (Pi_k + 0.5 * h * M_k)
    
    # 2. 求解隐式相对姿态步进矩阵 F_k in SO(3)
    F_k, f_vec, _ = solve_increment(g, J)
    
    # 3. 严格流形前推更新姿态
    R_next = R_k @ F_k
    
    # 4. 计算下一时刻外力矩并更新角动量 (保辛离散动量映射)
    M_next = torque_func(t_k + h, R_next)
    Pi_next = F_k.T @ Pi_k + 0.5 * h * (F_k.T @ M_k) + 0.5 * h * M_next
    
    return R_next, Pi_next

---

9. 权威参考文献与延伸阅读

1. Lee, Taeyoung, Leok, Melvin, and McClamroch, N. Harris. (2007). “Lie group variational integrators for the full body problem.” Computer Methods in Applied Mechanics and Engineering (CMAME), 196(29–32):2907–2924. (本文核心奠基文献)
2. Lee, Taeyoung, Leok, Melvin, and McClamroch, N. Harris. (2006). “Attitude dynamics on $\mathrm{SO}(3)$ using Lie group variational integrators.” In IEEE Conference on Decision and Control (CDC), pp. 6006–6011.
3. Leok, Melvin. (2004). “Foundations of Numerical Geometric Mechanics.” PhD thesis, California Institute of Technology. (离散变分几何基础博士论文)
4. Marsden, Jerrold E., and West, Matthew. (2001). “Discrete mechanics and optimal control.” Acta Numerica, 10:357–514. (离散力学与变分积分器世界奠基作)
5. Hairer, Ernst, Lubich, Christian, and Wanner, Gerhard. (2006). “Geometric Numerical Integration: Structure-Preserving Algorithms for Ordinary Differential Equations.” Springer.
6. Kobilarov, Marin, et al. (2010). “Discrete geometric optimal control on Lie groups.” IEEE Transactions on Robotics (T-RO), 26(2):306–319.
7. Yoshida, Haruo. (1990). “Construction of higher order symplectic integrators.” Physics Letters A, 150(5–7):262–268.