如何用Python基于哈密顿量正确编写二次微分方程求解器绘制相空间轨迹

现有代码核心问题 #

• 坐标更新逻辑错误：你当前使用x = x + p*t_step更新坐标，其中t_step是累计总时长，会导致坐标增量计算完全偏离真实物理过程，正确的单步增量应为p * dt。
• 动量符号丢失：你通过直接对(E - U(x))（你代码中的H(x,p0) - p0**2本质就是势能项U(x)）开平方计算动量，默认仅取正根，完全丢失了动量的方向信息。

简易调整方案（保留能量约束计算逻辑） #

如果你希望继续基于能量守恒约束计算p，可以通过保留历史符号+转折点自动判断的方式处理符号，不需要手动干预：

核心逻辑：正常运动过程中动量符号和上一步保持一致，当粒子运动到势场转折点（即E - U(x)接近0，粒子速度接近0即将反向）时自动翻转符号。
调整后的代码示例：

import numpy as np

x = x0
p = p0
E = p**2 + U(x)
dt = 0.01
t_step = 0
eps = 1e-6 # 数值判断阈值，避免浮点误差导致判断失效

while t_step < tmax:
    # 先更新坐标
    x = x + p * dt
    # 计算当前位置的势能
    U_x = U(x)
    # 判断是否到达转折点，是则翻转动量符号
    if E - U_x < eps:
        p = -p
    # 计算动量大小，乘当前符号得到正确动量
    p = np.sqrt(np.abs(E - U_x)) * np.sign(p)
    t_step += dt

更优数值求解方案（推荐） #

对于哈密顿系统的相空间轨迹模拟，更推荐使用辛积分算法（比如辛欧拉法、蛙跳法），这类算法天生适配哈密顿系统的守恒特性，不仅会自动处理动量符号，还能保证长时间模拟下能量不漂移，非常适合相空间绘图场景。
你不需要再手动通过能量约束计算p，直接基于哈密顿正则方程迭代即可：
对于你给出的哈密顿量E = p**2 + U(x)（已做归一化，p=dx/dt），正则方程为：
dp/dt = -dU(x)/dx
对应的辛欧拉迭代代码：

x = x0
p = p0
dt = 0.01
t_step = 0

# 先定义势能的导数函数，例如谐振子U(x)=0.5*k*x²的话，导数为k*x
def dU_dx(x):
    # 替换为你的势能导数计算逻辑
    return k * x

while t_step < tmax:
    # 辛欧拉迭代顺序：先更新动量，再更新坐标
    p = p - dU_dx(x) * dt
    x = x + p * dt
    t_step += dt

内容的提问来源于stack exchange，提问作者Zachary