2. 非线性振子数学模型建立

非线性振子的数值模拟，第一步不是写代码，而是把真实物理系统“翻译”成可计算的数学对象。这一步决定了后续所有仿真结果是否有意义。本节以单摆和Duffing振子为主线，完整走一遍从物理建模、无量纲化、状态空间表示，到平衡点求解、稳定性判断与雅可比矩阵计算的流程。

一、从物理系统到运动方程

考虑长度为 L、质量为 m 的单摆，摆角 θ 从竖直向下位置量起，逆时针为正。对摆球做受力分析：重力切向分量为 mgsinθ，方向指向平衡位置；阻尼力矩取与角速度成正比，系数为 c。由牛顿第二定律的转动形式得：

m L² θ̈ = -m g L sinθ - c θ̇

两边除以 mL²，得到标准形式：

θ̈ + (c/(m L²)) θ̇ + (g/L) sinθ = 0

这里的关键是：非线性来自 sinθ 这一项，而不是来自阻尼或惯性。当 θ 很小时 sinθ ≈ θ，方程退化为线性谐振子；一旦摆幅增大，sinθ 的高阶项不可忽略，系统就表现出振幅依赖的周期、倍频成分等非线性特征。

另一类常见模型是Duffing振子，用于描述硬化或软化弹簧、梁的大挠度振动等：

ẍ + δ ẋ + α x + β x³ = F cos(ω t)

其中 δ 为阻尼系数，α 为线性刚度，β 为非线性刚度系数（β>0 为硬化，β<0 为软化），F、ω 为外激励幅值与频率。与单摆不同，Duffing振子的非线性直接写在恢复力中，且常带周期外激励，属于非自治系统。

二、无量纲化：减少参数、暴露本质

直接对含 m、L、g、c 的方程做数值积分并非不可，但参数多、量纲混杂，不利于分析。无量纲化的目的是把方程化为只含少数无量纲参数的形式。对单摆，取时间尺度 τ = t√(g/L)，则 d²θ/dt² = (g/L) d²θ/dτ²，代入后得：

θ'' + 2ζ θ' + sinθ = 0,   2ζ = c/(m√(gL))

这里撇号表示对无量纲时间 τ 求导，ζ 是无量纲阻尼比。原方程四个参数被压缩为一个 ζ，物理意义更清晰：ζ 越小，系统越接近无阻尼振荡；ζ≥1 时线性化系统不再振荡。

对Duffing振子，取 τ = ω₀t（ω₀=√α），并令 γ = δ/ω₀、ε = β/α、f = F/α、Ω = ω/ω₀，可得：

x'' + γ x' + x + ε x³ = f cos(Ω τ)

无量纲化的操作要点是：先列出所有参数的量纲，选取由系统本身决定的时间尺度（如 √(L/g) 或 1/ω₀），再把每个变量除以它的特征尺度。这样得到的参数组合才是真正控制系统行为的“旋钮”。

三、状态空间表示

数值求解二阶方程的标准做法是降阶为一阶方程组。对单摆，令 x₁ = θ，x₂ = θ'，则：

x₁' = x₂
x₂' = -2ζ x₂ - sin x₁

写成向量形式 x' = f(x)，其中 x = [x₁, x₂]ᵀ。状态空间的好处是：任意时刻系统的完整状态由二维平面上的一个点表示，轨迹就是相平面中的曲线，平衡点、闭轨、分界线等结构一目了然。对带外激励的Duffing振子，还需引入相位变量 φ = Ωτ，状态变为三维 x = [x, x', φ]ᵀ，系统成为自治的三维流。

四、平衡点与稳定性分析基础

平衡点满足 f(x*) = 0。对单摆，由 x₂=0 和 sin x₁=0 得 x₁ = kπ。物理上，k 为偶数对应摆球在最低点，k 为奇数对应倒立位置。稳定性判断的直观依据是：最低点附近受扰动会来回摆动，是稳定的；倒立位置受微小扰动就会远离，是不稳定的。但“稳定”的严格判断需要线性化。

设 x = x* + ξ，对 f 在平衡点处做一阶泰勒展开：

ξ' ≈ J ξ,   J = ∂f/∂x |_(x=x*)

J 即雅可比矩阵。对单摆系统：

J = [ 0        1      ]
    [ -cos x₁  -2ζ    ]

在 (0,0) 处，cos0=1，J = [[0,1],[-1,-2ζ]]，特征方程 λ² + 2ζλ + 1 = 0，特征值 λ = -ζ ± √(ζ²-1)。当 0<ζ<1 时，特征值为一对共轭复数且实部为负，平衡点是稳定焦点，轨迹螺旋收敛；ζ=0 时为纯虚根，是中心，轨迹为闭合椭圆；ζ>1 时为两个负实根，是稳定结点。

在 (±π, 0) 处，cos(±π) = -1，J = [[0,1],[1,-2ζ]]，特征方程 λ² + 2ζλ - 1 = 0，特征值一正一负（因常数项为 -1），故为鞍点，必然不稳定。这与物理直觉一致。

五、雅可比矩阵的实际计算

对Duffing振子（无激励时），状态为 x₁=x、x₂=x'，方程为 x₁'=x₂、x₂'=-x₁-εx₁³-γx₂。平衡点由 x₂=0、x₁(1+εx₁²)=0 给出。当 ε>0 时只有原点一个实平衡点；当 ε<0 时还有 x₁=±√(-1/ε) 两个平衡点。雅可比矩阵为：

J = [ 0            1      ]
    [ -(1+3εx₁²)  -γ     ]

以 ε=-0.5、γ=0.2 为例，非零平衡点为 x₁=±√2 ≈ ±1.414。在原点处，J 的特征方程为 λ²+0.2λ+1=0，特征值实部为 -0.1，原点为稳定焦点。在 x₁=√2 处，1+3εx₁² = 1-3 = -2，特征方程 λ²+0.2λ-2=0，判别式为正且常数项为负，特征值一正一负，该平衡点为鞍点。这说明软化Duffing系统存在两个稳定焦点和一个鞍点，相图中会出现典型的双井结构。

六、建模流程的检查要点

完成上述步骤后，建议按以下清单自查，避免把错误带入数值积分阶段：

• 运动方程中每一项的量纲是否一致，符号约定（角度正方向、阻尼方向）是否全文统一；
• 无量纲化后参数个数是否减少，每个无量纲参数的物理含义是否明确；
• 状态变量是否覆盖了系统的全部自由度，非自治系统是否补充了相位变量；
• 平衡点是否求全（包括周期解对应的平衡点），是否遗漏复数或高次方程的根；
• 雅可比矩阵是否在平衡点处求值，特征值实部符号是否与稳定性结论对应；
• 线性化结论是否仅用于局部判断，大范围行为是否仍需依赖相图或数值积分验证。

需要强调的是，线性化只能回答“平衡点附近的小扰动如何演化”，不能直接给出大振幅下的周期、分岔或混沌行为。雅可比矩阵给出的是局部稳定性判据，而全局结构（如分界线、闭轨、吸引子）必须通过相平面分析或数值积分来揭示。因此，建模阶段的输出应当是一组清晰的方程、无量纲参数、状态向量、平衡点坐标和对应的雅可比矩阵，这些正是后续编写数值模拟程序时最直接的输入。