简介：这份资源围绕Fermi-Pasta-Ulam问题的数值模拟展开，面向对非线性动力学、统计物理与科学计算感兴趣的Python学习者，适合具备一定编程基础、希望动手复现经典物理实验的读者。资源包共6个文件，包含1个Python脚本、1个Markdown说明、1个许可证、1个gitignore，以及1个gif和1个mp4演示文件，压缩包约15.34MB，脚本负责实现带二次与三次非线性项的弹性力模拟，演示文件直观呈现振动演化过程。程序默认k为1，alpha与beta可手动调整，将二者设为大于0即可产生非谐波力，描述中给出的alpha=0、beta=0.3示例可作为入门参考。运行依赖numpy与matplotlib 1.4.2以上版本，若需将动画导出为mp4，还需配置电影编写器，Linux下推荐使用avconv而非ffmpeg。目前已有262人学习下载，读者可借此理解FPU问题的建模思路、非线性项对能量分布的影响，并掌握动画生成与常见排错方法。
1. 从一次“能量不守恒”的翻车说起：FPU 问题到底在算什么
如果你用 Python 写过分子动力学或者非线性振动模拟，大概率遇到过这种玄学现象：明明没有外加阻尼，系统总能量却在几十步之后肉眼可见地漂移，动能谱里还冒出一堆莫名其妙的频率成分。我第一次复现 Fermi-Pasta-Ulam 问题时就栽在这上面，当时以为是积分器精度不够，把步长从 1e-3 一路压到 1e-6，结果能量该漂还是漂。后来才想明白，FPU 问题本身就是在研究“能量为什么没有按预期均分到各个模式”，它压根不是一个能量守恒的演示，而是一个非线性系统里模式耦合与回归现象的经典实验场。
这份fpu-problem资源，核心就是一套用 Python 实现的 FPU 链模拟代码。它把一维非线性晶格里的 α-FPU 和 β-FPU 两种势能形式都做了封装，支持固定边界和周期边界，输出包括各模式能量随时间演化、实时位移曲线以及相空间轨迹。适合谁用？如果你正在做非线性动力学课程设计、想验证辛积分器在哈密顿系统里的长期行为，或者单纯想搞明白“热化”这件事在数值上到底长什么样，这套代码能直接跑起来给你看。它不依赖冷门库，numpy + matplotlib 就能撑起全部计算和可视化，对新手友好，但参数没调对的话，翻车姿势也很多。
2. 拆开 FPU 链：势能形式、模式分解与积分器选型
2.1 α-FPU 与 β-FPU 的物理差异
FPU 链的本质是一串用非线性弹簧连接的质点。设第 i 个质点的位移为 u_i，相邻质点间的伸长量为 Δu_i = u_{i+1} - u_i，两种常见势能写法是：
- α-FPU：V(Δu) = (1/2)Δu² + (α/3)Δu³
- β-FPU：V(Δu) = (1/2)Δu² + (β/4)Δu⁴
α 项是非对称非线性，会让系统产生二次谐波，能量更容易往高频模式泄漏；β 项是对称非线性，短时间内的模式耦合更“温和”，但长时间看依然会破坏简正模式的独立性。代码里用potential_type参数切换，alpha和beta分别控制非线性强度。我一般建议新手从 β-FPU 入手，因为它的数值行为更稳定，不容易在第一步就炸掉。
2.2 模式能量怎么算：离散正弦变换
要观察能量是否均分，得把实空间位移转到简正模式空间。对固定边界条件，模式坐标 Q_k 就是位移序列的离散正弦变换（DST）：
import numpy as np def compute_mode_energies(u, v, N, k_max=None): """ 将实空间位移 u 和速度 v 转到简正模式，计算各模式能量。 u: shape (N,), 质点位移 v: shape (N,), 质点速度 N: 质点数 k_max: 最大模式数，默认取 N-1 返回: E_modes, shape (k_max,) """ if k_max is None: k_max = N - 1 n = np.arange(1, N + 1) E_modes = np.zeros(k_max) for k in range(1, k_max + 1): # 模式形状：sin(pi*k*n/(N+1)) phi = np.sin(np.pi * k * n / (N + 1)) # 归一化因子 norm = np.sum(phi ** 2) Qk = np.dot(u, phi) / norm Pk = np.dot(v, phi) / norm # 模式频率 omega_k = 2*sin(pi*k/(2*(N+1))) omega_k = 2.0 * np.sin(np.pi * k / (2 * (N + 1))) E_modes[k - 1] = 0.5 * (Pk ** 2 + omega_k ** 2 * Qk ** 2) return E_modes这段代码的逻辑是：先把位移投影到每个简正模式的正弦基函数上，得到模式坐标 Q_k 和模式动量 P_k，再按简谐振子的能量公式算 E_k = 0.5(P_k² + ω_k²Q_k²)。参数k_max控制你关心前多少个模式，通常取 N/2 就够看趋势了。注意归一化因子norm不能省，否则模式能量总和和实空间总能量对不上，你会误以为积分器漏能。
2.3 积分器选型：为什么我最终用了 Störmer-Verlet
FPU 系统是哈密顿系统，长时间模拟最怕数值耗散。显式欧拉和四阶 Runge-Kutta 我都试过：欧拉直接让能量指数增长，RK4 短期精度高但长期有缓慢漂移。Störmer-Verlet 是辛积分器，虽然只有二阶精度，但它能保持相空间体积，能量误差有界振荡而不发散。代码里默认用 Verlet，核心更新步骤：
def verlet_step(u, v, a, dt): """ 单步 Störmer-Verlet 更新。 u: 当前位移 v: 当前速度 a: 当前加速度（由力计算得到） dt: 时间步长 返回: 更新后的 u, v, a """ v_half = v + 0.5 * dt * a # 半步速度 u_new = u + dt * v_half # 全步位移 a_new = compute_acceleration(u_new) # 用新位移算力 v_new = v_half + 0.5 * dt * a_new # 全步速度 return u_new, v_new, a_new参数dt的选取有个血泪经验：不要只看位移变化，要看最高模式频率。ω_max ≈ 2，所以 dt 至少小于 0.1 才能稳定，我一般取 0.01 到 0.05。如果你把 β 调到 10 以上，dt 还得再降，否则第一步就给你报 overflow。
3. 从零跑通一次 β-FPU 模拟：参数配置与结果解读
3.1 初始化：正弦模态激发与边界条件
FPU 原始实验的初始条件很讲究：只激发最低频模式，看能量多久泄漏到其他模式。代码里用initial_mode指定激发模式编号，amplitude控制初始振幅。固定边界下，两端质点位移恒为零，所以实际自由度是 N-1 个。周期边界则用 roll 操作处理邻居索引。
def init_fpu(N=32, mode=1, amplitude=0.5, boundary='fixed'): """ 初始化 FPU 链。 N: 质点数 mode: 初始激发的模式编号（从1开始） amplitude: 初始振幅 boundary: 'fixed' 或 'periodic' 返回: u0, v0 """ n = np.arange(1, N + 1) if boundary == 'fixed': # 固定边界：两端位移为零，模式形状用 sin(pi*k*n/(N+1)) u0 = amplitude * np.sin(np.pi * mode * n / (N + 1)) else: # 周期边界：模式形状用 cos(2*pi*k*n/N) u0 = amplitude * np.cos(2 * np.pi * mode * n / N) v0 = np.zeros(N) # 初始速度为零 return u0, v0这里mode=1就是只激发最低频模式，amplitude=0.5在 β=1 时已经能明显看到非线性效应。如果你把 amplitude 设到 2 以上，固定边界下两端位移约束会被破坏，模拟直接发散。我一般从 0.1 开始试，确认稳定后再往上加。
3.2 主循环：时间推进与能量记录
主循环负责调用 Verlet 步进，每隔一定步数记录一次模式能量。记录频率不用太高，否则输出文件巨大，我一般每 10 到 50 步记一次。
def run_fpu(N=32, mode=1, amplitude=0.5, beta=1.0, dt=0.02, steps=100000, record_every=50, boundary='fixed'): """ 运行 FPU 模拟并返回模式能量历史。 返回: times, E_hist """ u, v = init_fpu(N, mode, amplitude, boundary) a = compute_acceleration(u, beta=beta, boundary=boundary) times = [] E_hist = [] for step in range(steps): u, v, a = verlet_step(u, v, a, dt) if step % record_every == 0: E_modes = compute_mode_energies(u, v, N) times.append(step * dt) E_hist.append(E_modes) return np.array(times), np.array(E_hist)参数steps和record_every要配合着调：总模拟时间 = steps * dt，如果你要看回归现象，总时间至少得覆盖几十个最低频周期。最低频周期约 2π/ω_1 ≈ 2π/(2sin(π/(2(N+1))))，N=32 时大约 64 个时间单位，所以 steps 取 100000、dt=0.02 对应 2000 个时间单位，足够看到多次回归。
3.3 结果图怎么看：模式能量热图与回归现象
跑完之后最直观的是画模式能量热图：横轴时间，纵轴模式编号，颜色深浅表示能量占比。正常 β-FPU 在低振幅下会看到能量先往高频扩散，然后大部分能量又回到初始模式，这就是著名的 FPU 回归。如果你把振幅调大或者 β 调大，回归会变弱，能量最终趋向均分。代码里用imshow加colorbar就能出图，注意把纵轴模式编号从 1 开始标，别用默认的 0 索引，不然看图的人会懵。
4. 避坑与排查：FPU 模拟里最容易翻车的五个地方
4.1 能量不守恒，越跑越大
现象：总能量随时间单调上升，几十步后翻倍。原因：步长 dt 太大，Verlet 虽然辛但仍有稳定性上限，超过临界步长后数值解指数发散。解决：把 dt 降到 0.01 以下，或者检查力计算里有没有漏掉负号。我一般会先跑 1000 步看能量漂移是否小于 1%，超了就砍 dt。
4.2 模式能量总和与实空间能量对不上
现象：E_modes 求和比实空间总能量少一截，或者多一截。原因：DST 归一化因子写错，或者模式数 k_max 没取够，漏掉了高频模式的贡献。解决：确认norm = np.sum(phi**2)用的是完整正弦基，k_max 至少取到 N-1。如果还不对，检查边界条件是否和模式基函数匹配——固定边界用 sin，周期边界用 cos，混用必错。
4.3 初始振幅稍大就 overflow
现象：运行几步后位移变成 inf 或 nan。原因：β 项是四次方，振幅大时力增长极快，显式积分器直接炸。解决：降低 amplitude 到 0.1 量级，或者改用自适应步长。但 FPU 问题通常研究弱非线性，振幅本来就不该太大，硬调大反而失去物理意义。
4.4 回归现象看不到，能量一直均匀分布
现象：模式能量热图从头到尾都是花的，没有明显回流。原因：要么 β 太大，要么模拟时间太短，要么初始模式不是最低频。解决：把 β 降到 0.1 到 1 之间，mode 设为 1，总模拟时间至少覆盖 50 个最低频周期。如果还不行，检查记录频率是不是太高导致热图看起来乱，试试record_every加大到 100。
4.5 周期边界下能量突然跳变
现象：周期边界模拟中，总能量在某个时刻出现台阶式跳变。原因：邻居索引用 roll 时方向搞反，导致力计算不对称，相当于引入了虚假外力。解决：周期边界下，第 i 个质点的右邻居是 (i+1)%N，左邻居是 (i-1)%N，用np.roll(u, -1)和np.roll(u, 1)分别取右和左。检查力计算里两个方向的伸长量是否都用了正确的 roll 方向。
5. 进阶技巧：用 FFT 替代 DST 与长时间模拟的稳定性检查
5.1 周期边界下用 FFT 加速模式分解
固定边界必须用 DST，但周期边界下简正模式就是离散傅里叶模式，直接用np.fft.fft比循环投影快一个数量级。我后来把compute_mode_energies改了一版，周期边界走 FFT 分支：
def compute_mode_energies_fft(u, v, N): """ 周期边界下用 FFT 计算模式能量。 """ U = np.fft.fft(u) / N V = np.fft.fft(v) / N k = np.fft.fftfreq(N, d=1.0/N) # 模式编号 omega = 2.0 * np.abs(np.sin(np.pi * k / N)) E_modes = 0.5 * (np.abs(V) ** 2 + omega ** 2 * np.abs(U) ** 2) return E_modes[:N//2] # 只取正频率部分注意 FFT 返回的是复数，取模平方才是能量。omega用2*sin(pi*k/N)是因为周期边界下色散关系不同。这版代码跑 N=128 的长时间模拟时，单步耗时从 8ms 降到 0.5ms，省下来的时间够你多跑几组参数。
5.2 长时间模拟的稳定性自检清单
跑超过 10^6 步之前，我习惯做三件事：第一，用 dt 和 dt/2 各跑 10000 步，看模式能量曲线是否重合，不重合说明步长还不够小；第二，每 10^5 步检查一次总能量相对漂移，超过 1% 就停下来查原因；第三，把 β 设为 0 跑一遍，线性极限下模式能量应该严格守恒，如果这都不守恒，那肯定是代码有 bug 而不是物理效应。
5.3 参数扫描的批量跑法
如果你要研究 β 和振幅对回归时间的影响，别手动改参数一次次跑。写个外层循环，把结果存成 npz 文件：
import itertools beta_list = [0.1, 0.5, 1.0, 2.0] amp_list = [0.1, 0.3, 0.5] for beta, amp in itertools.product(beta_list, amp_list): t, E = run_fpu(beta=beta, amplitude=amp, steps=200000) np.savez(f'fpu_beta{beta}_amp{amp}.npz', times=t, E_hist=E)这样跑完一批之后，用 matplotlib 画个回归时间随 β 变化的曲线，比单看一张热图有说服力得多。注意文件命名别用小数点，有些系统会出问题，用下划线代替。
从那以后我每次跑 FPU 模拟，都强制先跑线性极限自检，确认能量守恒再开非线性。这个习惯帮我省了至少三次通宵排查。希望帮到你。
本文还有配套的精品资源，点击获取
