1. 非线性振子数值模拟导论

线性振子——弹簧质量块、小角度单摆、RLC 电路——之所以好解，是因为它们的恢复力与位移成正比，叠加原理成立，解析解可以写成正弦与余弦的组合。一旦恢复力偏离线性，或者系统中出现与速度相关的非线性阻尼，叠加原理失效，解析解往往不存在闭式表达，振子的行为也随之变得丰富：振幅与频率耦合、周期随能量改变、出现自持振荡、跳跃与迟滞、对初值极端敏感。数值模拟正是研究这类现象的主要手段。

为什么必须数值模拟

非线性振子的运动方程一般写作 ẍ = f(x, ẋ, t)，其中 f 对 x 或 ẋ 非线性。以带三次刚度的 Duffing 振子为例：

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

当 β ≠ 0 时，方程不再有通用的解析解。我们真正关心的问题——给定初值后振幅如何演化、稳态响应在哪个频段发生跳跃、极限环的周期是多少、两个几乎相同的初值何时分道扬镳——都只能通过逐步积分得到数值答案。数值模拟的价值不只是"算出曲线"，而是让我们能够系统地扫描参数、观察相图与庞加莱截面、估计李雅普诺夫指数，从而把定性理论与定量结果对上。

需要强调的是，数值模拟不是"随便跑一下"。非线性系统对积分步长、方法阶数和总积分时长都很敏感：步长过大会引入数值阻尼，把极限环"抹平"；积分时间过短会把暂态误判为稳态；混沌系统还会放大舍入误差。因此从第一节课起，就要建立"参数—方法—验证"三位一体的意识。

三个典型非线性振子

Duffing 振子描述硬弹簧或软弹簧系统，恢复力含三次项。取 δ=0.2, α=-1, β=1, γ=0.3, ω=1.2 时，它表现出双势阱、混沌与周期窗口交替的复杂行为，是研究分岔与混沌的标准模型。它的相图在无阻尼无激励时呈现同宿轨道与两个稳定焦点，加入周期激励后出现马蹄映射式的拉伸折叠。

Van der Pol 振子的特征是负阻尼项 μ(1-x²)ẋ：小振幅时系统吸能、振幅增大，大振幅时耗能、振幅减小，最终收敛到唯一的极限环。取 μ=1 时极限环接近圆形，周期约 6.66；取 μ=5 时波形出现明显的张弛振荡特征——缓慢爬升后快速翻转，此时固定步长的显式方法容易失稳，需要更小步长或刚性求解器。

单摆是最容易被低估的非线性振子。方程 θ̈ + (g/L)·sinθ = 0 中，sinθ 的泰勒展开首项是 θ，所以小角度下近似线性，周期为 2π√(L/g)；但振幅增大后周期随之变长，θ₀=90° 时实际周期比线性近似长约 18%。当能量足够高时，摆会越过顶点进入转动状态，相图上分界线（separatrix）把摆动与转动两类轨道分开。这个模型非常适合用来验证数值方法的能量守恒性。

课程项目总览

本课程围绕一个贯穿始终的项目展开：为上述三类振子建立统一的数值实验框架，并逐步加入更复杂的分析能力。项目按能力递进分为四个阶段，每个阶段都产出可运行的代码与可解释的图。

工具链配置

Python 与 MATLAB 都能胜任，课程示例以 Python 为主、MATLAB 作为对照。Python 侧建议使用 numpy 做向量化运算、scipy.integrate.solve_ivp 调用 RK45 与 Radau 求解器、matplotlib 绘图。环境准备命令如下：

python -m venv venv
source venv/bin/activate        # Windows: venv\Scripts\activate
pip install numpy scipy matplotlib

一个最小可运行骨架如下，它把二阶方程降为一阶方程组，并用 solve_ivp 积分。注意 rtol 与 atol 的设定直接决定结果可信度，建议先取 1e-8 量级做基准，再逐步放宽观察差异。

import numpy as np
from scipy.integrate import solve_ivp

def duffing(t, s, d=0.2, a=-1.0, b=1.0, g=0.3, w=1.2):
    x, v = s
    return [v, -d*v - a*x - b*x**3 + g*np.cos(w*t)]

sol = solve_ivp(duffing, [0, 200], [0.1, 0.0],
                method='RK45', rtol=1e-8, atol=1e-10,
                dense_output=True, max_step=0.01)
print(sol.y.shape, sol.t[-1])

MATLAB 侧对应写法是 ode45(@duffing,[0 200],[0.1 0])，并用 odeset('RelTol',1e-8,'AbsTol',1e-10) 控制精度。两者结果应当一致到绘图分辨率以内，这本身就是一次有价值的交叉验证。