在分子动力学的数值方法模拟中，经常使用Leapfrog算法来“追踪”每个粒子的速度以及位置。Leapfrog的最简单总结就是，其在时间上把对速度求解和对位置求解（即对速度积分）分离开，在此刻算速度，在下一刻算位置，循环往复。
那么，为什么要这样做呢？为什么Euler方法的「速度和位置同时算」就会带来能量增长以及精度下降？一个可行的分析角度是从哈密顿力学的相空间出发。
一、Leapfrog算法简述
对于一个只考虑静电场的数值程序，标准的staggered Leapfrog模式为：
- 位置 $\mathbf{x}$ 定义在整数时间步 $t^n$；
- 速度 $\mathbf{v}$ 定义在半整数时间步 $t^{n+1/2}$；
- 电场由 $t^n$ 时刻的位置/电荷密度计算。
时间推进循环为：
$$
\mathbf{v}^{n+1/2} = \mathbf{v}^{n-1/2} + \frac{q}{m} \mathbf{E}(\mathbf{x}^n) \Delta t
$$
$$
\mathbf{x}^{n+1} = \mathbf{x}^n + \mathbf{v}^{n+1/2} \Delta t
$$
这个模式满足时间反演对称性，当把$\Delta t$换以$-\Delta t$，公式的形式完全不变。这意味着：将同一物理过程向前推进100步再向后推进100步，理论上将精确回到起点。
而在实际数值方法实现上（即 离散化后），上述时间反演对称性仍然将被以一种振荡的方式满足。为了在数学上更清晰的理解这件事情，我们需要借助相空间这一工具来看。
二、相空间面积及刘维尔定理
在由位置-动量组成的坐标系，即相空间中，每个相点$(q, p)$的轨迹称为相曲线。对于保守哈密顿系统（无耗散，$\mathcal{H}$二阶连续可微），刘维尔定理断言这个相曲线为一封闭曲线，并且其围成的面积不变。或者其规范表述：
保守哈密顿系统的相流在演化中保持相空间体积（测度）不变，等价于分布函数沿相轨迹的随体导数为零。
2.1 辛结构与体积守恒
设相空间坐标 $\mathbf{z} = (q_1, \ldots, q_n, p_1, \ldots, p_n)$，哈密顿方程为：
$$
\dot{\mathbf{z}} = \mathbf{J} \nabla \mathcal{H}
$$
其中 $\mathbf{J}$ 是辛矩阵（symplectic matrix）：
$$
\mathbf{J} = \begin{pmatrix} \mathbf{0} & \mathbf{I} \ -\mathbf{I} & \mathbf{0} \end{pmatrix}
$$
相流 $\phi_t$ 将初始点 $\mathbf{z}_0$ 映射到 $t$ 时刻的 $\mathbf{z}(t)$。其雅可比矩阵 $\mathbf{M}(t) = \partial\mathbf{z}(t)/\partial\mathbf{z}_0$ 满足辛条件：
$$
\mathbf{M}^\top \mathbf{J} \mathbf{M} = \mathbf{J}
$$
取行列式：
$$
\det(\mathbf{M}^\top) \det(\mathbf{J}) \det(\mathbf{M}) = \det(\mathbf{J})
$$
由于 $\det(\mathbf{J}) \neq 0$ 且 $\det(\mathbf{M}^\top) = \det(\mathbf{M})$：
$$
\det(\mathbf{M})^2 = 1 \quad \Rightarrow \quad \det(\mathbf{M}) = \pm 1
$$
对于从恒等映射连续演化的相流，$\det(\mathbf{M}(0)) = +1$，故对所有 $t$ 有 $\det(\mathbf{M}(t)) = +1$。这意味着相空间体积元 $dV = dq_1 \wedge dp_1 \wedge \cdots \wedge dq_n \wedge dp_n$ 在演化中严格守恒。
2.2 面积与能量的关系
对于一维简谐振子，哈密顿量为：
$$
\mathcal{H} = \frac{p^2}{2m} + \frac{1}{2}kq^2 = E
$$
等能面是相空间中的椭圆，半轴分别为：
$$
a = \sqrt{\frac{2E}{k}}, \quad b = \sqrt{2mE}
$$
椭圆面积为：
$$
A = \pi ab = \pi \sqrt{\frac{2E}{k}} \sqrt{2mE} = 2\pi E \sqrt{\frac{m}{k}} = \frac{2\pi E}{\omega}
$$
其中 $\omega = \sqrt{k/m}$ 是振子角频率。面积与能量成正比：$A \propto E$。
因此，若数值方法导致相空间面积膨胀（$\det \mathbf{M} > 1$），则等效于系统能量被人工放大；若面积收缩（$\det \mathbf{M} < 1$），则能量被人工耗散。只有 $\det \mathbf{M} = 1$ 时，能量才不会出现系统性漂移。
2.3 离散映射的体积变化
将连续相流替换为离散映射 $\mathbf{z}_{n+1} = \mathbf{F}(\mathbf{z}_n)$ 时，每一步的雅可比矩阵 $\mathbf{M}_n = \partial\mathbf{F}/\partial\mathbf{z}_n$ 决定了体积变化：
$$
\Delta V_{n+1} = \det(\mathbf{M}_n) \cdot \Delta V_n
$$
- 若 $\det(\mathbf{M}_n) > 1$（如 Euler 方法），体积指数膨胀，能量单调增长；
- 若 $\det(\mathbf{M}_n) < 1$，体积指数收缩，能量单调衰减；
- 若 $\det(\mathbf{M}_n) = 1$（如 Leapfrog 方法），体积守恒，能量围绕真实值振荡。
三、Euler方法的相空间面积不守恒
以简谐振子为例：
$$
H = \frac{p^2}{2m} + \frac{1}{2}kq^2
$$
运动方程：
$$
\dot{q} = \frac{p}{m}, \quad \dot{p} = -kq
$$
前向 Euler 的离散映射为（整步推进）：
$$
q_{n+1} = q_n + \frac{p_n}{m}\Delta t
$$
$$
p_{n+1} = p_n – kq_n\Delta t
$$
写成矩阵：
计算这个矩阵的行列式：
$$
\det(M_{\text{Euler}}) = 1 \cdot 1 – \left(\frac{\Delta t}{m}\right)(-k\Delta t) = 1 + \frac{k(\Delta t)^2}{m}
$$
因为 $k, m, \Delta t$ 都是正的，所以
$$
\det(M_{\text{Euler}}) = 1 + \underbrace{\frac{k(\Delta t)^2}{m}}_{> 0} > 1
$$
行列式 > 1 意味着，相空间中任何微小面积 $dA = dq \times dp$，经过 Euler 一步映射后：
$$
dA_{n+1} = \det(M) \cdot dA_n = \left(1 + \frac{k(\Delta t)^2}{m}\right) dA_n
$$
可见，每步面积都在膨胀。那么例如1000步后：
$$
dA_{1000} = \left(1 + \frac{k(\Delta t)^2}{m}\right)^{1000} dA_0 \approx e^{1000 \cdot \frac{k(\Delta t)^2}{m}} dA_0
$$
这就是指数增长的来源。Euler 把相空间“拉伸”了，并且只拉伸不压缩。因此在长步数下计算得到的结果难以满足能量守恒等的物理先验。
四、Leapfrog方法的相空间面积不变
Leapfrog对简谐振子模型可以写成：
$$
p_{n+1/2} = p_{n-1/2} – kq_n\Delta t
$$
$$
q_{n+1} = q_n + \frac{p_{n+1/2}}{m}\Delta t
$$
或者更清楚地看成两个正则变换的复合：
Step 1: Kick（只动量变，位置不变）：
矩阵：
$$
M_{\text{kick}} = \begin{pmatrix} 1 & 0 \ -\frac{k\Delta t}{2} & 1 \end{pmatrix}, \quad \det = 1
$$
Step 2: Drift（只位置变，动量不变）
矩阵：
$$
M_{\text{drift}} = \begin{pmatrix} 1 & \frac{\Delta t}{m} \ 0 & 1 \end{pmatrix}, \quad \det = 1
$$
Step 3: 再 Kick
完整的 Leapfrog 映射：
$$
M_{\text{Leapfrog}} = M_{\text{kick}} \cdot M_{\text{drift}} \cdot M_{\text{kick}}
$$
注意到，每个子步的映射都是辛映射，矩阵行列式都为1，3次或者有限次矩阵相乘，结果行列式仍为1：
$$
\det(M_{\text{Leapfrog}}) = 1 \times \ldots \times 1 = 1
$$
所以，Leapfrog推进的相空间面积严格守恒。
在代码实现上，需要对上述计算公式进行离散化，此时Leapfrog推进的相空间面积将不是定值，但仍然将在初始面积的附近振荡。这在一定精度下确保了能量守恒的物理先验。
如图是一个以Python程序实现的两种方法对比，分别绘制了两种方法运行10步、100步后的相空间面积。可见，Leapfrog相比于Euler算法更好地贴切了初始状态的相空间面积。
五、总结
Euler 的离散映射不是正则变换（$\det M \neq 1$），它在每一步“放大”了相空间面积，而面积 $\propto$ 能量，所以能量指数增长。
Leapfrog 的离散映射是正则变换（$\det M = 1$），它只“拧转”相空间（剪切变形），面积绝对值不变，所以能量不会系统性漂移，只会围绕真实值振荡。
附：Python程序
import numpy as np
import matplotlib.pyplot as plt
# ========== 参数设置 ==========
m = 1.0
k = 1.0
dt = 0.3
n_points = 200 # 圆盘上的点数
r = 1.0
# 生成初始圆盘
theta = np.linspace(0, 2*np.pi, n_points, endpoint=False)
q0 = r * np.cos(theta)
p0 = r * np.sin(theta)
def polygon_area(q, p):
"""鞋带公式计算多边形面积"""
return 0.5 * abs(np.sum(q[:-1]*p[1:] - q[1:]*p[:-1]) + q[-1]*p[0] - q[0]*p[-1])
def run_euler(q, p, n_steps, dt, m, k):
"""前向Euler推进"""
q_e, p_e = q.copy(), p.copy()
for _ in range(n_steps):
q_new = q_e + (p_e / m) * dt
p_new = p_e - k * q_e * dt
q_e, p_e = q_new, p_new
return q_e, p_e
def run_leapfrog(q, p, n_steps, dt, m, k):
"""Leapfrog推进，返回整数时间层的(q, p)"""
q_l = q.copy()
# 初始化半步速度: v^{-1/2} = v^0 - 0.5*a^0*dt
p_half = p.copy() - 0.5 * (-k * q_l) * dt
for _ in range(n_steps):
q_l = q_l + (p_half / m) * dt # Drift
p_half = p_half + (-k * q_l) * dt # Kick
# 最后半步kick回到整数时间层: p^N = p^{N+1/2} - 0.5*a^N*dt
p_l = p_half - 0.5 * (-k * q_l) * dt
return q_l, p_l
# ========== 运行模拟 ==========
steps_list = [10, 100]
results = {}
for n in steps_list:
q_e, p_e = run_euler(q0, p0, n, dt, m, k)
q_l, p_l = run_leapfrog(q0, p0, n, dt, m, k)
results[n] = {
'euler': (q_e, p_e, polygon_area(q_e, p_e)),
'leapfrog': (q_l, p_l, polygon_area(q_l, p_l)),
'det_theory': (1 + k*dt**2/m)**n
}
area0 = polygon_area(q0, p0)
# ========== 绘图 ==========
fig, axes = plt.subplots(2, 3, figsize=(14, 9))
fig.suptitle('Phase Space Area Conservation: Euler vs Leapfrog',
fontsize=14, fontweight='bold', y=0.98)
c_euler, c_leapfrog, c_init = '#E74C3C', '#3498DB', '#2C3E50'
# 第一行：10步
for col, (title, color, data_key) in enumerate([
('Initial State', c_init, None),
('Forward Euler (10 steps)', c_euler, 'euler'),
('Leapfrog (10 steps)', c_leapfrog, 'leapfrog')
]):
ax = axes[0, col]
if data_key is None:
q, p, area = q0, p0, area0
else:
q, p, area = results[10][data_key]
ax.fill(q, p, color=color, alpha=0.15)
ax.plot(q, p, 'o', color=color, markersize=2, alpha=0.6)
ax.set_title(f'{title}\nArea = {area:.3f}', fontsize=11,
color=color if data_key else c_init)
ax.set_ylabel('p')
ax.set_aspect('equal')
ax.grid(True, alpha=0.3)
ax.set_xlim(-4, 4)
ax.set_ylim(-4, 4)
# 第二行：100步
for col, (title, color, data_key) in enumerate([
('Initial State', c_init, None),
('Forward Euler (100 steps)', c_euler, 'euler'),
('Leapfrog (100 steps)', c_leapfrog, 'leapfrog')
]):
ax = axes[1, col]
if data_key is None:
q, p, area = q0, p0, area0
else:
q, p, area = results[100][data_key]
ax.fill(q, p, color=color, alpha=0.15)
ax.plot(q, p, 'o', color=color, markersize=1.5, alpha=0.5)
if data_key == 'euler':
det_th = results[100]['det_theory']
ax.set_title(f'{title}\nArea = {area:.0f}\n(expansion factor: {area/area0:.0f}×)',
fontsize=11, color=color)
max_c = max(np.abs(q).max(), np.abs(p).max()) * 1.1
ax.set_xlim(-max_c, max_c)
ax.set_ylim(-max_c, max_c)
else:
ax.set_title(f'{title}\nArea = {area:.3f}', fontsize=11,
color=color if data_key else c_init)
ax.set_xlim(-4, 4)
ax.set_ylim(-4, 4)
ax.set_xlabel('q')
ax.set_ylabel('p')
ax.set_aspect('equal')
ax.grid(True, alpha=0.3)
plt.tight_layout(rect=[0, 0, 1, 0.96])
plt.savefig('phase_space_demo.png', dpi=200, bbox_inches='tight',
facecolor='white', edgecolor='none')
plt.show()
# ========== 打印总结 ==========
print("=" * 60)
print("相空间面积守恒演示结果")
print("=" * 60)
print(f"参数: m={m}, k={k}, dt={dt}, 初始圆盘半径 r={r}")
print(f"Euler行列式: det = 1 + k·dt²/m = {1 + k*dt**2/m:.4f}")
for n in steps_list:
_, _, area_e = results[n]['euler']
_, _, area_l = results[n]['leapfrog']
print(f"\n{n} 步后:")
print(f" Euler: 面积 = {area_e:10.2f} (膨胀 {area_e/area0:.0f}×)")
print(f" Leapfrog: 面积 = {area_l:10.4f} (守恒)")备注：①代码由Kimi生成。②实务中，考虑计算速度，通常使用C++程序实现，此处只做效果演示。