行星齿轮非线性程序，能出相图，庞加莱，分叉图。
行星齿轮系统这玩意儿搞动力学分析，满屏都是非线性项。刚接触那会儿看着微分方程头皮发麻，后来发现直接数值求解才是真香。咱今天就拿Python撸个能出相图、庞加莱截面和分叉图的脚本，手把手看看这非线性系统能整出什么活。
先上核心代码骨架：
import numpy as np from scipy.integrate import odeint import matplotlib.pyplot as plt mu = 0.1 # 阻尼比 k = 10.0 # 刚度 omega_n = np.sqrt(k) # 固有频率 Omega = 0.8 * omega_n # 激励频率 F = 0.5 # 激励幅值 def system(y, t): x, dx = y ddx = -mu*dx - k*x + F*np.cos(Omega*t) + 0.1*x**3 # 含三次非线性项 return [dx, ddx] t = np.linspace(0, 1000, 50000) # 长时间仿真 sol = odeint(system, [0.1, 0], t)这段代码把行星齿轮常见的非线性因素——三次刚度项给塞进去了。注意时间序列取得特别长，这是为了后面取庞加莱截面时能过滤掉瞬态响应。不过直接画相图的话得截取后面稳定段：
# 截取稳定段 x, dx = sol[-2000:,0], sol[-2000:,1] plt.plot(x, dx, ',k', alpha=0.5) plt.xlabel('Displacement') plt.ylabel('Velocity') plt.title('Phase Portrait')跑出来的相图要是出现个极限环，说明系统在做周期性振动。但要是看到轨迹开始乱窜，那可能进入混沌状态了——这时候就该庞加莱截面出场了。
取庞加莱截面的骚操作在于按激励周期采样：
# 计算激励周期对应的采样间隔 T = 2*np.pi / Omega dt = t[1] - t[0] stride = int(T/dt) # 取截面点 poincare_x = x[::stride] poincare_dx = dx[::stride] plt.figure() plt.plot(poincare_x, poincare_dx, '.r', markersize=2) plt.title('Poincare Section')当庞加莱截面呈现孤立点时是周期运动，出现成片的点云可能就混沌了。不过最直观的还是分叉图——看参数变化时系统如何从老实变癫狂：
bifurcation = [] omega_range = np.linspace(0.5, 2.0, 300) for Om in omega_range: # 每次微调参数重新仿真 def local_system(y, t): x, dx = y ddx = -mu*dx -k*x + F*np.cos(Om*t) + 0.1*x**3 return [dx, ddx] sol = odeint(local_system, [0.1,0], t) x = sol[-2000:,0] # 记录局部极值 peaks = (x[1:-1] > x[:-2]) & (x[1:-1] > x[2:]) bifurcation.append(x[1:-1][peaks]) plt.figure() for i, Om in enumerate(omega_range): plt.plot([Om]*len(bifurcation[i]), bifurcation[i], ',k', markersize=0.1) plt.title('Bifurcation Diagram')这段分叉图代码有个坑：参数循环时每次都要重新定义微分方程。如果直接修改外层变量会导致闭包问题，所以得在循环内部重新定义local_system。当分叉图从单支分裂成多支，说明系统开始出现倍周期分叉——这是进入混沌的前戏。
跑完这三板斧，基本上能把行星齿轮的非线性特性摸个大概。不过实际工程中参数选择更讲究，比如阻尼比μ别超过0.3，否则啥非线性现象都给你压没了。另外激励频率Ω在固有频率附近最容易出活，调参时可以重点照顾这个区间。