双吸引子与双稳态：基于双势阱杜芬振子的非线性动力学仿真与可视化

“他有两个吸引子。”第一次听到这句话很多人会以为是在聊人格分裂。但在非线性系统语境里这句话描述的是一种真正具有“选择能力”的动力学结构同一个微分方程同一组参数仅仅因为初始状态不同系统最终会收敛到两个完全不同的稳定终态。生活中双稳态电路、相变存储器的“写入/擦除”、大脑联想记忆的某个状态、神经网络训练时从不同随机初始化落到不同局部最优背后都可以抽象成“有两个吸引子”的问题。我的判断是如果你想真正理解动力系统不要一开始就去啃 Lyapunov 指数、分岔图、奇异吸引子这些高阶概念而是先从一个能直接看出来、能亲手仿出来的“两个稳定不动点”开始。双吸引子是最低成本的“记忆单元”也是相空间结构里第一个真正的分水岭。谁掌握了双吸引子谁就理解了“状态收敛”和“初值敏感”到底是怎么回事。这篇文章会用一个非常经典、物理图像清晰的系统——双势阱杜芬振子带你把“两个吸引子”从数学定义、稳定性分析、数值模拟到盆地可视化完整走一遍。整篇文章不需要高性能服务器不需要 GPU只需要一台装了 Python 的电脑。读完以后你不仅能解释“他为什么有两个吸引子”还能用代码验证它的吸引子位置、稳定性以及不同初始点到底“属于”哪一个吸引子。1. 这篇文章真正要解决的问题很多初学者看动力系统教材最容易卡住的地方是这样的书上定义吸引子用“不变集”“拓扑传递”“极限集”这些概念看起来非常严谨但你回到代码里不知道怎么写一个最小系统去验证。尤其是“两个吸引子”或“多个吸引子”的场景如果不亲手画出相图很难建立直觉。这篇文章要解决的问题并不是单纯解释“吸引子是什么”而是回答四个工程里更常遇到的疑问。第一一个系统凭什么说它有“两个吸引子”判断标准是什么第二如果我想复现一个有两个吸引子的系统最少的代码量是多少第三两个吸引子的边界在哪里为什么初始条件稍微差一点最终落点就完全不同第四当我要给同事或团队讲透这个概念时应该用什么实验来佐证而不是只贴公式把这些问号逐个拆开你会发现“两个吸引子”并不是一个抽象哲学概念而是可以被数值模拟直接观察到的几何事实。用双势阱杜芬振子做实验平台是因为它在保持极低复杂度的同时具备了稳定不动点、鞍点、盆域边界这三个要素。对入门者来说这个系统是“刚刚好”的复杂度太简单则看不到两个选择太复杂则会淹没在混沌和分岔里。从工程价值看理解两个吸引子还有一层现实原因造一个“有记忆”的系统通常就是制造两个吸引子。一个只会收敛到唯一状态的系统不可能存住任何信息只有当存在两个稳定状态时这个系统才能像 1 比特存储器、双稳态触发器、联想记忆节点一样在输入结束后继续保持状态。你在处理强化学习策略多峰、优化损失函数多个极小值、神经形态计算器件状态切换时本质上都是在跟“多吸引子”打交道。所以这篇内容适合所有想进入非线性动力学、物理信息机器学习、控制理论或类脑计算方向的开发者。2. 基础概念与核心原理2.1 什么是吸引子吸引子是一个动力学系统在长时间演化后最终停留或无限逼近的状态集合。可以把它理解为系统状态在相空间里的“归宿”。一个系统通常由微分方程或差分方程描述系统的状态会随时间变化。初始状态确定之后状态按照规则演化。如果演化到足够长的时间状态不再离开某个集合并且倾向于稳定在该集合附近那么这个集合就是一个吸引子。最直观的例子是小球滚进碗底无论你从碗口的哪个位置轻轻放下小球只要碗里有摩擦它最后都会停在碗底最低点。这个最低点就是一个稳定不动点吸引子。技术一点的说法是吸引子是一个不变集合并且在集合的某邻域内几乎所有的轨道最终都会进入并保持在它的任意小邻域内。这个定义强调了两件事其一是不变性一旦状态进入吸引子后续状态不会跑出去其二是吸引性附近轨道都有向它汇拢的趋势。吸引子的种类并非只有一种。常见的有稳定不动点、极限环和奇异吸引子。稳定不动点是系统最终静止在某一点比如本节的小球停在碗底极限环是系统最终进入周期性振荡比如范德波尔振荡器奇异吸引子则是混沌系统在相空间中的有界非周期集合比如洛伦兹系统著名的蝶形吸引子。对本文来说双势阱杜芬振子产生的是两个稳定不动点吸引子。2.2 两个吸引子意味着双稳态一个系统只有一个吸引子时它的行为比较单调不管从哪出发最后都去往同一个终点。这虽然稳定但没有“选择”。当系统存在两个吸引子时情况发生了本质变化最终状态依赖于初始条件落在哪个吸引子的“势力范围”也就是落在哪个盆域里。双稳态在工程里极其常见。一个金属球放在倒扣的碗面上可能向左滚也可能向右滚但不会停在顶部D 触发器的输出要么是高电平要么是低电平不会稳定在中间值忆阻器的一个单元可能处于高阻态或低阻态用来表示逻辑 0 和逻辑 1。这些系统的共同点都是存在两个稳定状态中间隔着不稳定的鞍点或势垒。用两个吸引子的语言来描述这些系统的状态空间有两个“碗底”。系统没有外界输入时状态会停留在其中一个碗底外界以足够强的扰动让状态越过中间势垒后系统便可能切换到另一个碗底。这个过程就是写入信息、切换状态或者完成计算。理解了这一点你就明白了为什么研究“两个吸引子”不只是数学爱好者的趣味而是一切状态型器件的动力学本质。2.3 双势阱杜芬振子最经典的“两个吸引子”实验平台在众多存在两个吸引子的系统里双势阱杜芬振子是最适合教学的系统之一。它的物理图像是一个质量为 1 的质点在一个 W 形的势能曲线中运动同时受到线性弹簧、非线性回复力和阻尼的共同作用。这里使用的无量纲方程是x y y x - x^3 - δ * y其中 x 表示位移y 表示速度δ 表示阻尼系数。为了产生明显且稳定的两个吸引子本文统一取 δ 0.2。对应的势能函数是V(x) -0.5 * x^2 0.25 * x^4画出这个函数会得到两个局部极小值点分别在 x -1 和 x 1它们的势能都是 -0.25。在 x 0 处势能取到局部极大值 0相当于两个谷之间的“山脊”。质点在这个 W 形轨道上运动时如果没有阻尼它会在谷内来回振荡甚至可能越过中间山脊加入阻尼后动能会被不断消耗最终质点会停在某个谷底附近。这个系统的状态不止一个维度。x 和 y 一起组成的二维状态空间称为相平面。一个稳定状态在相平面里对应一个点因此两个吸引子对应相平面对称位置的两个点左侧稳定点 (-1, 0) 和右侧稳定点 (1, 0)。质点停在哪一边取决于初始位移和初始速度的组合。2.4 平衡点与稳定性谁能成为吸引子并不是所有平衡点都能成为吸引子。平衡点只是状态不再变化的点即 y 0 且 x - x^3 0解出来有三个A- (-1, 0) S (0, 0) A (1, 0)为了判断它们是否稳定可以计算雅可比矩阵J [[0, 1], [1 - 3*x^2, -δ]]代入各个平衡点之后求矩阵的特征值在 (-1, 0)矩阵的局部线性化特征值实部为负说明邻域内的轨道会螺旋收敛到这个点。在 (1, 0)情况完全对称同样是一个稳定的螺旋吸引子。在 (0, 0)两个特征值一个为正、一个为负因此这是一个鞍点。鞍点不是稳定吸引子但它就像分水岭上的“山口”是决定状态最终流向哪一侧的关键。所以这个系统里真正有资格被称为吸引子的只有两个也就是 A- 和 A。中间那个平衡点 S 不是吸引子而是一条分界线的源头。一个重要细节是阻尼必须大于 0。如果 δ 0系统变成保守系统没有能量耗散A- 和 A 附近的轨道会一直绕圈不会收敛也就不满足吸引子的严格定义。实际工程中要让系统具有“记忆力”往往需要引入合适的耗散或误差衰减机制这正是优化算法里阻尼、动量项、正则项在动力学层面扮演的角色。2.5 稳定流形与盆域边界两个吸引子各自有一个“控制范围”人们把这个范围称为吸引盆或盆域。落进左侧盆域的初始状态最终会去往 A-落进右侧盆域的初始状态最终会去往 A。盆域之间的边界并不是随意画的一条竖线而是鞍点 S 的稳定流形。在二维相平面里鞍点有一条稳定流形和一条不稳定流形。稳定流形上的轨道会沿时间正向趋近鞍点但在实际数值计算中初始状态很难恰好落在稳定流形上一点微小的扰动就会让状态偏向其中一侧最终被两个吸引子之一捕获。因此靠近鞍点的初始状态对“选左还是选右”极其敏感。初始状态越靠近鞍点系统判断归属需要的时间就越长任何数值误差或物理噪声都会影响最终选择。这种现象用下面的初值扫描实验可以非常直观地看到。下表总结了本节的核心内容名称坐标局部线性化角色A-(-1, 0)复特征值实部为负左侧吸引子S(0, 0)一个正特征值、一个负特征值鞍点、盆域边界A(1, 0)复特征值实部为负右侧吸引子3. 环境准备与前置条件本文的代码基于 Python主要使用 NumPy、SciPy 和 Matplotlib 三个库。SciPy 的 solve_ivp 用来做常微分方程数值积分NumPy 负责数组运算Matplotlib 用于绘制相图和盆域图。建议使用 Python 3.8 或更高版本。版本号并不过分敏感只要保证依赖库能正常安装即可。安装命令如下pip install numpy scipy matplotlib如果你的电脑同时存在多个 Python 环境建议先在虚拟环境中执行上述命令。如果不希望通过命令行安装也可以在 PyCharm 或 VS Code 的 Python 解释器设置里直接安装这三个依赖。实际操作中你会用到 4 个脚本duffing_phase.py、potential_well.py、stability_check.py、basin_map.py。前两个脚本用来建立直观印象第三个脚本做严格的局部稳定性验证第四个脚本绘制两个吸引子的盆域图。建议按顺序运行因为每个脚本都从不同角度回答同一个核心问题这个系统为什么有两个吸引子。4. 核心流程拆解整个实验的流程并不复杂但每一步都有明确目的。完整流程可以拆成六个阶段。第一阶段是建立模型。我们把二阶运动方程改写成一阶常微分方程组。改写方法很简单令位移为 x速度为 y那么位移的时间导数是 y速度的时间导数是外力项也就是 x - x^3 减去阻尼项 δ*y。这里的核心思想是引入状态变量把一个二阶方程变成两个一阶方程方便用数值积分器处理。第二阶段是找到系统的全部平衡点。令 x 0 且 y 0可以得到 y 0以及 x - x^3 0。解出平衡点后不要急于下结论因为平衡点不一定是吸引子还要看局部稳定性。第三阶段是计算雅可比矩阵和特征值。雅可比矩阵本质上是把非线性系统在每个平衡点附近线性化特征值的实部决定了状态会不会收敛到该点。如果所有特征值实部都为负那么该平衡点在非线性系统里是一个局部吸引子如果一个特征值实部为正、一个为负那么它是鞍点。第四阶段是从两个不同的初始条件出发对完整非线性系统做数值积分。模拟起点选在小球分别位于左侧势阱和右侧势阱的位置例如 (-0.8, 0.0) 和 (0.8, 0.0)时间跨度设为 0 到 40。观察两个轨迹最终到达的位置。第五阶段是绘制相图。相图的横轴是 x纵轴是 y轨迹在相平面中画出一圈圈向吸引子靠近的螺旋线。你可以同时把 A-、A 和鞍点 S 标记在图上这样就能直观看到两条轨道分别在不同区域画圈最终停在不同星标位置。第六阶段是扫描初值并绘制盆域图。在相平面内取一个矩形网格每个网格点作为一个初始条件独立做数值积分然后用最终落点是左侧还是右侧给网格点着色。这样得到的就是两个吸引子的盆域可视化盆地分界也会自动显现出来。这套流程的核心价值是它可复用。你以后遇到任意二维自治系统都可以按照同样的方法完成“找平衡点、做局部稳定性分析、全局模拟、画盆域”四步而不只是看懂本文里的特例。5. 完整示例与代码实现5.1 示例 1从两个初始条件看相轨迹收敛先创建一个文件duffing_phase.py代码如下# duffing_phase.py import numpy as np import matplotlib.pyplot as plt from scipy.integrate import solve_ivp DELTA 0.2 T_MAX 40.0 def duffing(t, state): x, y state return [y, x - x * x * x - DELTA * y] def integrate(x0, y0): sol solve_ivp( duffing, [0, T_MAX], [x0, y0], dense_outputTrue, rtol1e-8, atol1e-10, ) x_final, y_final sol.y[:, -1] print(f初值 ({x0:.2f}, {y0:.2f}) - 终值 ({x_final:.8f}, {y_final:.8f})) return sol def plot_phase(sol_a, sol_b): t np.linspace(0, T_MAX, 3000) fig, ax plt.subplots(figsize(7, 6)) for sol, label, color in [ (sol_a, from left well, tab:blue), (sol_b, from right well, tab:red), ]: x_traj sol.sol(t)[0] y_traj sol.sol(t)[1] ax.plot(x_traj, y_traj, colorcolor, lw1.2, labellabel) ax.scatter([-1, 1], [0, 0], s200, marker*, colorblack, labelattractors) ax.scatter([0], [0], markers, facecolorwhite, edgecolorblack, labelsaddle, zorder3) ax.set_xlabel(x) ax.set_ylabel(y) ax.set_title(Double-well Duffing: two stable attractors) ax.legend() ax.grid(alpha0.4) plt.savefig(duffing_phase.png, dpi150) plt.show() if __name__ __main__: sol_left integrate(-0.8, 0.0) sol_right integrate(0.8, 0.0) plot_phase(sol_left, sol_right)这段代码的关键点有三个。第一solve_ivp的dense_outputTrue会让积分器额外生成一个连续可访问的解函数方便后面绘制光滑轨迹曲线。第二积分容差设置为rtol1e-8, atol1e-10可以避免因数值误差造成吸引子判断失误。第三最终状态用sol.y[:, -1]取出打印在终端供验证。运行命令python duffing_phase.py程序会打印两个结果。从左侧势阱出发的初始状态终态应该接近 (-1, 0)从右侧势阱出发的初始状态终态应该接近 (1, 0)。这个输出就是“系统有两个吸引子”最直接的数值证据。5.2 示例 2画出 W 形势能曲线为了解释为什么会有两个稳定点可以单独绘制势能函数曲线。创建potential_well.py# potential_well.py import numpy as np import matplotlib.pyplot as plt def potential(x): return -0.5 * x * x 0.25 * x**4 x_vals np.linspace(-2.0, 2.0, 1000) v_vals potential(x_vals) plt.figure(figsize(8, 5)) plt.plot(x_vals, v_vals, lw2, colorsteelblue) for point_x in [-1, 0, 1]: point_v potential(point_x) plt.plot(point_x, point_v, o, colorblack) plt.annotate( f({point_x:.0f}, {point_v:.2f}), xy(point_x, point_v), xytext(point_x 0.15, point_v 0.15), fontsize10, ) plt.axhline(0, colorgray, lw0.8, ls--) plt.xlabel(x) plt.ylabel(V(x)) plt.title(V(x) -0.5*x^2 0.25*x^4) plt.grid(alpha0.3) plt.savefig(potential_well.png, dpi150) plt.show()这张图的意义在于解释“选择”为什么会发生。两个势阱底分别位于 -1 和 1中间 x 0 是势能局部极大值。小球如果有足够能量越过中间凸起它就能从一个势阱跑到另一个势阱如果没有足够能量就只能待在初始所在的势阱里。势能曲线给出了系统“两个吸引子”的物理来源但要注意真正吸引子的位置必须结合速度维度看否则容易把势能最低点直接等同于吸引子这个理解在相空间中并不完整。运行命令python potential_well.py你会看到一条清晰的 W 形曲线左侧最低点、右侧最低点和中间顶点分别对应 A-、A 与平衡点 S。5.3 示例 3用雅可比矩阵验证稳定性前面的相图只是数值现象为了从理论上确认三个平衡点的类型需要计算雅可比矩阵的特征值。创建stability_check.py# stability_check.py import numpy as np DELTA 0.2 def jacobian(x, y): return np.array( [ [0.0, 1.0], [1.0 - 3.0 * x * x, -DELTA], ] ) for name, x, y in [(A-, -1.0, 0.0), (S, 0.0, 0.0), (A, 1.0, 0.0)]: eigenvalues np.linalg.eigvals(jacobian(x, y)) print(f{name}: coordinate({x:.1f}, {y:.1f}), eigenvalues{eigenvalues})运行命令python stability_check.py预期输出类似A-: coordinate(-1.0, 0.0), eigenvalues[-0.11.4106736j -0.1-1.4106736j] S: coordinate(0.0, 0.0), eigenvalues[