一、系统建模与参数设置
两自由度Duffing系统动力学方程：
参数定义：
c1 = 0.1; % 阻尼系数1
c2 = 0.15; % 阻尼系数2
k1 = 1; % 线性刚度1
k2 = 1.2; % 线性刚度2
alpha = 0.5;% 非线性系数
beta = 0.3;% 耦合强度
F1 = 0.8; % 外力幅值1
F2 = 0.6; % 外力幅值2
omega = 1;% 驱动频率
phi = pi/4;% 相位差
二、数值求解与相图绘制
1. 状态方程转换
将二阶方程转换为四阶状态空间形式：
function dydt = duffing2DOF(t,y)x = y(1); dx = y(2);y_pos = y(3); dy = y(4);dxdt = dx;d2xdt2 = -c1*dx - k1*x - alpha*x^3 - beta*(x-y) + F1*cos(omega*t);dydt = dy;d2ydt2 = -c2*dy - k2*y - alpha*y^3 - beta*(y-x) + F2*cos(omega*t + phi);dydt = [dxdt; d2xdt2; dydt; d2ydt2];
end
2. 四阶龙格库塔求解
h = 0.01; % 时间步长
t_end = 500; % 总仿真时间
t = 0:h:t_end;% 初始条件
y0 = [0.1; 0; 0.2; 0]; % [x, dx/dt, y, dy/dt]% 数值积分
Y = zeros(4,length(t));
Y(:,1) = y0;
for i = 1:length(t)-1k1 = duffing2DOF(t(i), Y(:,i));k2 = duffing2DOF(t(i)+h/2, Y(:,i)+k1*h/2);k3 = duffing2DOF(t(i)+h/2, Y(:,i)+k2*h/2);k4 = duffing2DOF(t(i)+h, Y(:,i)+k3*h);Y(:,i+1) = Y(:,i) + h/6*(k1+2*k2+2*k3+k4);
end
三、映射图生成
1. 相空间轨迹图
figure;
plot(Y(1,:), Y(3,:), 'b', 'LineWidth',2);
hold on;
plot(Y(1,1:1000:end), Y(3,1:1000:end), 'ro');
xlabel('x(t)'); ylabel('y(t)');
title('两自由度Duffing系统相图');
legend('轨迹','采样点');
grid on;
2. 庞加莱截面映射
figure;
Poincare_section = Y(3,:) .* (mod(t, 2*pi/omega) < 0.01);
plot(Y(1,:), Poincare_section, 'g.');
xlabel('x(t)'); ylabel('y(t)');
title('庞加莱截面图');
grid on;
3. 分岔图（Bifurcation Diagram）
f_range = 0.5:0.05:2; % 驱动频率范围
bifurcation = zeros(length(f_range),1);for i = 1:length(f_range)omega = f_range(i);Y = solve_system(omega); % 调用求解函数bifurcation(i) = max(Y(1,end-1000:end)); % 取稳态最大值
endfigure;
plot(f_range, bifurcation, 'r-o');
xlabel('驱动频率 \omega (rad/s)');
ylabel('最大位移 x_{max}');
title('分岔图');
grid on;
四、关键结果分析
- 相图特征： 周期运动：闭合轨道（如周期1、2） 混沌吸引子：非闭合的复杂轨迹（如Lorenz型结构）
- 分岔行为： 随着频率变化，系统经历倍周期分岔→混沌→周期窗口→再分岔
- 能量交换： 耦合项β显著影响两自由度间的能量传递效率
五、高级功能扩展
1. 参数敏感性分析
% 参数扰动实验
beta_values = 0.1:0.05:0.5;
chaos_indicator = zeros(size(beta_values));for i = 1:length(beta_values)beta = beta_values(i);[~, max_lyapunov] = compute_lyapunov(Y); % 计算Lyapunov指数chaos_indicator(i) = max_lyapunov;
endfigure;
stem(beta_values, chaos_indicator);
xlabel('\beta (耦合强度)');
ylabel('最大Lyapunov指数');
title('混沌敏感性分析');
2. 三维时频分析
figure;
spectrogram(Y(1,:),256,250,256,omega);
title('x(t)时频谱图');
xlabel('时间 (s)'); ylabel('频率 (Hz)');
六、参考
- 核心文献： Jing Z. et al. (2005) 两自由度Duffing系统复杂动力学研究 Yagasaki K. (1999) 双频激励下的分岔控制
- 代码：两自由度duffing振子非线性系统的映射图程序 www.youwenfan.com/contentcno/96561.html
- MATLAB工具： Dynamical Systems Toolbox：提供分岔分析函数 Control System Toolbox：用于稳定性验证