Python实现FPU模拟：非线性振子链能量递归的数值复现-北京尧图网络

简介费米-帕斯塔-乌拉姆FPU问题模拟器是一份基于Python的完整仿真实现面向物理非线性系统爱好者、计算物理初学者及教学场景。程序在真实弹性基础上引入二次方与三次方非线性项通过调节alpha与beta参数即可直观观察非谐波力对能量分配的影响配合numpy与matplotlib运行并支持动画输出为gif或mp4视频。资源包共6个文件以Python仿真脚本为核心辅以参数示例的gif/mp4演示动画、Markdown说明文档、开源协议及版本控制配置整体约15.34MB。已有262人学习下载适合需要快速理解FPU问题机理、动手复现经典物理实验的Python开发者与物理专业学生。该资料提供参数可调的完整代码及可视化示例可直接运行验证并可作为非线性系统教学演示的参考模板。1. FPU 模拟第一步先跑出能量「赖着不走」的现象FPUFermi-Pasta-Ulam问题的模拟至今仍然是非线性动力学入门最值得动手复现的题目之一一条由非线性弹簧串起来的振子链把能量全部放在最低频模式上按经典统计物理的直觉能量应该逐渐在所有模式间均匀分布。然而模拟跑出来的结果完全相反——能量先在少数几个模式间来回振荡形成周期性递归热化迟迟不发生。这说明「非线性导致热化」这个直觉本身就有隐藏前提而让结论反转的关键就在参数和数值格式的选择上。这篇笔记从 FPU 链的建模出发给你可直接运行的 Python 代码、参数含义和边界条件并记录我在复现过程中踩过的坑。2. 物理模型与数值方法FPU-α 方程、初始模态激发与跳跃积分器2.1 一维链的力和势能怎么写FPU 问题的标准实验箱是一个一维振子链设位移为u_i(t)其中i1..N两端固定即u_0 u_{N1} 0。相邻粒子之间的势能不是纯二次型的而是带有一个非线性修正项。常见有三种写法纯三次项、纯四次项以及两者都保留。实做中最常用的是 FPU-α 模型和 FPU-β 模型。FPU-α 模型的势能写为V Σ [ (1/2)(u_{i1}-u_i)^2 (α/3)(u_{i1}-u_i)^3 ]其中 $\alpha$ 控制三次非线性强度。对位移求负梯度可得到作用在第i个粒子上的力f_i (u_{i1} - 2u_i u_{i-1}) α * [ (u_{i1}-u_i)^2 - (u_i-u_{i-1})^2 ]FPU-β 模型则对应四次势能V Σ [ (1/2)(u_{i1}-u_i)^2 (β/4)(u_{i1}-u_i)^4 ]对应力为f_i (u_{i1} - 2u_i u_{i-1}) β * [ (u_{i1}-u_i)^3 - (u_i-u_{i-1})^3 ]这里的物理含义是当相邻粒子间距较大时恢复力不再是线性的位移越大非线性项贡献越显著。初始条件通常取最低频本征模式——即让所有粒子按照u_i(0) A * sin(π i / (N1))排列初始速度为零。这样能量一开始完全集中在该模式上我们随后观察它如何向其他模式泄漏。在写代码时我建议把力函数单独抽出方便以后切换 α 模型和 β 模型。固定边界条件在数组下标上体现为左右两端各留一个虚拟粒子它们的位移始终为零。使用 NumPy 做向量化计算时u[1:] - u[:-1]就是相邻位移差注意边界下标处理即可。2.2 为什么不用 Runge-Kutta而用跳跃积分器做 FPU 模拟时我见过不少人直接用scipy.integrate.solve_ivp配合 RK45。这个问题虽然也能跑出近似结果但长期模拟会出现一个很棘手的现象RK 这类非辛格式会让总能量缓慢漂移。FPU 问题的核心观察对象正是能量在不同模式间的转移如果总能量本身就在数值上不断增长或衰减模态能量的归一化会完全失真递归图也会变得无法解释。更合适的选择是二阶辛跳跃法Störmer–Verlet也叫蛙跳法。它的核心优势是保证相空间体积守恒且长时间模拟下能量误差有界而不会漂移。时间离散格式为v_{i1/2} v_i (dt/2) * a(x_i) x_{i1} x_i dt * v_{i1/2} v_{i1} v_{i1/2} (dt/2) * a(x_{i1})其中a(x)是加速度。这种 kick-drift-kick 结构实现起来非常简单而且每步只需要两次力计算。另一个好处是蛙跳法在非线性扰动较弱时几乎等价于线性体系的精确解这让我们在验证程序时可以直接对照线性链的解析模态频率。当初我用 RK45 跑 FPU-α 模型能量在 20 万步之后漂移了约 3%换成蛙跳法同样步数下能量误差维持在 1e-10 量级——这个对比强烈推荐你自己做一次感受非常直观。2.3 模态能量怎么提取DST 而不是 FFT从位移数组提取物理模式能量是 FPU 模拟中被低估的一步。很多人直接用numpy.fft.fft对u_i做变换但这在物理上是错的FPU 链的两端是固定边界系统的本征模式是站波sin(k π i / (N1))不是周期性的平面波。用 FFT 变换相当于隐含了周期边界条件结果会出现频谱泄漏而且模式编号的对应关系也不准确。正确的做法是离散正弦变换DST-I。位移和速度分别做 DST然后每个模式的总能量按下式算出E_k 0.5 * ( v_k_hat^2 ω_k^2 * u_k_hat^2 )其中ω_k 2 sin(k π / (2(N1)))是第 k 个模式的线性频率u_k_hat和v_k_hat是位移和速度的 DST 系数。在 Python 中DST 可以使用 SciPy 的scipy.fft.dst函数并设置normortho保证正交归一化。这里有一个容易犯的错误normortho会让变换矩阵乘以一个归一化因子如果后续再乘上ω_k^2总能量的绝对数值会随着 N 变化不同尺寸链之间就不能直接做绝对能量对比。好在大多数场景我们只看「相对模态能量」——即每个模式能量占总能量的比例归一化因子在除法中自动抵消。所以只要保持一致normortho和normbackward并不会影响最终结论但绝对值的含义会不同这点在输出数据时要写清楚。我习惯把 DST 封装成一个独立函数这样在画模态能量时序图时每次只调用一次变换即可不用反复手写。3. Python 复现FPU 模拟的核心代码与运行流程3.1 你可以直接用的完整代码下面这份代码是我自己整理后常用的版本核心逻辑完整去掉无关的可视化样式。运行环境为 Python 3.10依赖 NumPy 和 SciPy。import numpy as np from scipy.fft import dst def fpu_alpha_force(u, alpha): 计算 FPU-α 模型每个粒子的力两端固定 # 为边界补零u[0] u[N1] 0 u_padded np.zeros(len(u) 2) u_padded[1:-1] u diff u_padded[1:] - u_padded[:-1] # 长度为 N1 的相邻差 # 线性项: 标准的三对角线形式 force u_padded[1:-1] * 2.0 force[1:] - u_padded[1:-2] force[:-1] - u_padded[2:-1] # 或者直接由相邻差二次方组成等价形式 linear diff[1:] - diff[:-1] nonlinear alpha * (diff[1:]**2 - diff[:-1]**2) return linear nonlinear def verlet_step(u, v, dt, alpha): 一个完整的时间步kick-drift-kick 形式 a fpu_alpha_force(u, alpha) v_half v 0.5 * dt * a u_new u dt * v_half a_new fpu_alpha_force(u_new, alpha) v_new v_half 0.5 * dt * a_new return u_new, v_new def mode_energies(u, v): 计算每个模式的能量使用 DST-I n len(u) k np.arange(1, n 1) omega 2.0 * np.sin(k * np.pi / (2 * (n 1))) u_hat dst(u, type1, normortho) v_hat dst(v, type1, normortho) E_k 0.5 * (v_hat**2 omega**2 * u_hat**2) return E_k这段代码里有个值得说明的细节力函数我先写了直观的差分形式又保留了从diff构建的等价形式。如果你仔细推导这两种写法在数学上完全一致但diff版本对边界更友好因为diff数组长度为N1diff[1:] - diff[:-1]直接得到内部 N 个点的力不需要额外处理左右边界。3.2 代码的关键参数说明fpu_alpha_force函数中的alpha是非线性系数。在经典文献中α 的取值从 0.01 到 0.5 都有人用过。一个关键的物理判断是alpha * (diff[1:]**2 - diff[:-1]**2)这一项的量级相对于线性项diff[1:] - diff[:-1]是否足够大。如果 α 太小非线性效应需要极长时间才能显现低温极限下甚至接近线性链的周期运动如果太大体系会进入强非线性区能量均分的现象又会变得不明显甚至出现混沌。dt的选择主要看最高模式频率。线性链的最高频率为ω_max ≈ 2当N较大时对应的振荡周期约π。因此dt原则上要小于最高频率周期的十分之一对于 N32 的链取dt0.05是安全值如果你把 N 提高到 128最高频率依然接近 2所以dt不需要等比例缩小太多但为保险起见建议dt 0.02。mode_energies函数里的omega数组直接由解析公式给出不依赖数值解后再拟合频率。这是和 FFT 方法最大的区别利用固定边界条件频率表达式是精确的即使非线性出现后瞬时频率会偏移但作为能量分解的基函数依旧有效。3.3 运行后你该看到什么用 N32α0.25初始激发模式 1即最正弦模式初始振幅 A0.1跑 5000 个时间步后你应该看到模式 1 的能量占比从接近 1 下降到某个低谷然后回升——形成周期性递归。具体的最低点取决于 α 和非线性耦合强度通常在 0.6 到 0.8 之间波动。模式 2、3 的能量会交替出现峰值形成动人的级联景象。这些图直接宣告了「非线性系统迅速热化」这一直觉破产。如果你看到的是能量在所有模式中很快均匀展开、且没有回来的趋势那么检查两个地方一是初始振幅是否太大比如 A 接近 1二是alpha是否过大导致直接进入强非线性区。把 A 降到 0.05 再跑一次递归通常会出现。4. 把实验跑起来参数设定、输出对比与调试路径4.1 参数扫描非线性强度如何影响等分时间CRUCIAL 的参数是 α、N 和初始振幅 A。三者的关系可以用一个无量纲参数衡量初始模式中单粒子位移的最大值与相邻差的乘积。当这个量远小于 1 时系统处于弱非线性区能量递归现象显著当它接近或大于 1 时系统进入强非线性区能量均分速度变快。你可以写一个简单循环扫描 α 从 0.01 到 0.5记录模式 1 能量首次降到 0.2 以下的时间点称为等分时间。你会发现等分时间随 α 的增大而减小但并非线性关系。在弱非线性极限下等分时间与 α 的平方近似成反比进一步提高 α 后等分时间降到某个平台期后基本不再变化。这个「平台期」对应着强非线性区物理上说明系统的随机化已经由非线性本身主导而非弱非线性共振。我建议你至少对比三组 α0.05、0.25、0.5并把三张模态能量图画在同一坐标下。这样能直观看到递归周期从长到短、从清晰到模糊的过程。4.2 模态能量演化图的判读方法画图时横轴是时间纵轴是每个模式的能量占比用不同颜色画前几个模式。这里有一个需要注意的点能量占比的计算前提是总能量守恒。如果总和随时间明显变化说明数值积分出了问题。在前面蛙跳的保证下这个和应该平直如一条直线。判读递归时有几个要点。首先看模式 1 的第一个低谷是否小于 0.5如果连 0.5 都没跌破说明非线性太弱需要跑更长时间或增大 α。其次看模式 2 和 3 的峰值到达时间是否与模式 1 低谷时间一致——如果有明显错位检查 DST 的omega数组是否写对了。最后看模式编号大于 10 的高频模式能量占比在弱非线性区它们几乎保持为零只有在强非线性区才会被显著激发。这是用来判断「系统是否真正进入大范围热化」的快速信号。4.3 典型参数组合与预期结果NαA初始振幅dt递归周期约模式 1 最低占比160.100.050.05约 700 时间单位0.52320.250.100.05约 400 时间单位0.61320.500.100.02约 180 时间单位0.38640.250.050.02约 600 时间单位0.55这张表仅供参考因为递归周期对初始振幅极其敏感。振幅增加一倍递归周期往往会缩到原来的四分之一到三分之一。实做中建议先从小到大扫描 A观察递归周期的变化趋势再确定你感兴趣的参数区域。调试路径上我习惯先用 N16 的小链跑通全流程确认能量守恒之后再切到大 N。小链跑得快迭代次数少打开日志就能看清每个时间步的能量分配便于快速定位问题。5. FPU 模拟常见问题与排查5.1 模态能量算出来总和随时间增长现象运行 2000 步后所有模式的能量总和逐渐上升到 5000 步已经涨了几倍。原因这种情况几乎都出在 DST 的归一化与omega公式不匹配上。当你使用normortho时DST 系数已经做了正交归一化此时omega公式里如果用2 sin(kπ/2(N1))是对的但如果用了2 sin(kπ/(N1))或弄错了分母等效频率和系数产生错配模式能量就会随时间出现「假增长」。解决打印初始时刻各模式能量的和与直接用位移速度计算的物理总能量做对比。两者必须相等误差小于 1e-10。如果不等调整omega的正弦自变量分母。这个报错是最容易自查的却也是最常被忽略的。5.2 dt 过大导致能量爆炸现象程序跑起来初期正常几百步后位移数值急剧增大最终溢出为 NaN。原因是非线性力在蛙跳格式下受到稳定性条件的限制。FPU-α 模型中非线性项包含位移差的平方当局部位移差稍微变大力就会猛增反过来让位移变得更大形成正反馈爆炸。线性格式的稳定性条件此时不再适用dt0.1在线性链没问题但配上强非线性直接崩。解决把 dt 减小到 0.02 或 0.01同时观察能量守恒。如果减小 dt 后仍然爆炸就要检查初位移是否有奇异值或者 alpha 是否过大导致初始时刻就有局部强拉伸。5.3 高频模式能量出现规律的假振荡现象模式 8、9、10 的能量一开始不为零而是出现周期性的小幅抖动幅度虽然不大但完全不衰减。原因这是边界条件或初始条件的「别名」效应。固定边界链的站波模式只有在整数 k 上才有定义如果你的初始位移不是严格的正弦站波而是任意离散数据点那么高频模式会被初始条件直接注入能量。此外在时间离散过程中跳跃格式对高频模式会产生截断误差表现为规律抖动。解决确保初始位移用精确的A * sin(p π i/(N1))其中 p 是整数模式号这是消除初始高频分量的前提。若抖动仍然存在使用低通滤波看是物理现象还是数值噪声——物理现象会随 α 变化数值噪声对 α 不敏感。5.4 递归周期与理论值差了几倍现象跑出来的递归周期总是比文献上的参考值短或长 3 到 5 倍。原因你或许把「时间单位」算错了。链中质量通常取 1刚度取 1这是无量纲化处理。如果你的初始位移中包含了质量或刚度的比例因子时间单位会随之缩放。另一个隐含问题是 N 不同时线性频率也变递归周期会显著变化。解决先将代码写成完全无量纲形式即质量、刚度均为 1。确认 A 是位移量纲而不是速度量纲。然后先用线性链验证把 α 设为 0模式 1 应该永远不向模式 2 转移能量这时递归周期无穷大。再慢慢增大 α观察递归周期的连续缩短就不会出现跳变。6. 进阶如何定量判定能量均分并把模拟扩展到 FPU-β 模型6.1 等分时间的定量判定方法递归虽然有趣但科学上更关心的是「什么时候系统真正趋近于能量均分」。你需要一个定量指标来分辨递归与热化。常用的是谱熵spectral entropy定义为S(t) -Σ p_k(t) ln p_k(t)其中p_k(t) E_k(t) / Σ E_k(t)是模式能量占比。当能量完全均分时p_k 1/N考虑有 N 个模式此时谱熵达到最大值ln N。实践中通常把熵达到 0.9 × ln N 的时刻定义为等分时间。这个指标比肉眼观察递归图可靠得多也方便做参数扫描时自动记录。你可以把这段函数加到代码中每个时间步计算一次def spectral_entropy(E_k): E_total E_k.sum() p E_k / E_total if E_total 0 else np.zeros_like(E_k) # 忽略零占比项避免 log(0) return -np.sum(p[p 0] * np.log(p[p 0]))从守恒角度你将看到α 较小时谱熵在很长时间内保持低值突然快速上升α 较大时谱熵从一开始就较快增长最终饱和在某个平台不再回升。这个「平台-陡增」的转变过程正是弱非线性共振与强非线性混沌之间过渡的标志。6.2 把代码改造成 FPU-β 模型前面力函数是针对三次非线性的 FPU-α改造为四次非线性的 FPU-β 并不复杂。你只需要把fpu_alpha_force中非线性项改为def fpu_beta_force(u, beta): linear (u[1:] u[:-1] - 2*u[1:-1]) # 注意索引 # 直接用相邻差构造 diff np.diff(np.concatenate(([0], u, [0]))) nonlinear beta * (diff[1:]**3 - diff[:-1]**3) return linear nonlinear你会发现 FPU-β 模型里递归现象比 α 模型弱得多而且在某些参数下系统更早进入热化。这是因为四次势能的正负不对称性更低非线性回旋共振的条件更难以满足。常用于实验验证的结论是β 模型在弱非线性极限下能量均分时间与 β 的平方成反比这一规律比 α 模型更好拟合。6.3 从一维链走向通用工具当你的 FPU 代码已经稳定不要急着扔掉。把位移和速度数组封装成类把力函数作为参数传入积分器它就能变成一个通用的晶格动力学模拟器。你可以把它扩展到二维三角格子、引入质量无序即某个粒子的质量不同或者加上阻尼项研究能量耗散。从那以后我每跑一个新的非线性动力学课题都会先拿 FPU 这个熟知的体系做基准测试用谱熵和能量守恒两个指标验证代码无误才敢继续做更复杂的模型。这套方法帮我避免了很多次在新模型里调试半天、最后发现是积分器写错的白费功夫。希望帮到你。本文还有配套的精品资源点击获取

Python实现FPU模拟：非线性振子链能量递归的数值复现

相关推荐

从理论到落地：构建企业级信息安全防线

【Web 渗透】10 款 Burp Suite 神器插件合集，直接拉满渗透测试效率

塔吉克语论文写到一半想求助？这份工具榜先按环节挑，别再只问谁第一 [特殊字符]

更多资讯

openHiTLS 社区贡献指南：CLA 签署、Issue 讨论与本地构建验证实战

视频播放异常排查：非APP层导致的黑屏卡顿根因定位

玲珑轮胎登上山东台：一条轮胎凭什么跑出“山东好品”的底气

PRAXIST 直接 CLI 运维指南：从校验、启动到观察、停止与恢复的完整生命周期

ATF转换PNG/JPG工具实战：从识别文件头到DXT1解压的完整指南

GuardianShell主动式AI野生动物哨兵：边缘计算架构与部署实战

今日更新

小波包能量谱原理与MATLAB实现：从频带分解到故障特征提取

基于LSTM的网易云音乐评论情感分析：数据清洗与模型调优实战

Alpha-Beta剪枝与极大极小值算法：博弈树搜索原理与实战优化

本周精选

交通拥堵预测毕设：基于Python+LSTM的车辆流量预测系统全流程解析

MATLAB/Simulink搭建10机39节点电力系统暂态稳定仿真实战

无人机航拍三维重建全流程：从SfM到网格生成的避坑指南

本月热门

我发现了一个新思路：用 Remotion + Claude Code 像写代码一样自动化生成短视频

Windows下 Codex 中 Chrome 和 Computer Use 插件不可用问题排查及解决参考方式：TaoToken 统一 Key 配置与验证

2026 大模型集体涨价：用 Python 做企业 Token 成本测算与选型避坑（附配置）

读完文章，想聊聊您的网站？