两自由度Duffing振子非线性系统的MATLAB实现方案

一、系统建模与参数设置

两自由度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)-1
    k1 = 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)); % 取稳态最大值
end

figure;
plot(f_range, bifurcation, 'r-o');
xlabel('驱动频率 \omega (rad/s)');
ylabel('最大位移 x_{max}');
title('分岔图');
grid on;

---

四、关键结果分析

1. 相图特征： 周期运动：闭合轨道（如周期1、2） 混沌吸引子：非闭合的复杂轨迹（如Lorenz型结构）
2. 分岔行为： 随着频率变化，系统经历倍周期分岔→混沌→周期窗口→再分岔
3. 能量交换： 耦合项β显著影响两自由度间的能量传递效率

---

五、高级功能扩展

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;
end

figure;
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)');

---

六、参考

1. 核心文献： Jing Z. et al. (2005) 两自由度Duffing系统复杂动力学研究 Yagasaki K. (1999) 双频激励下的分岔控制
2. 代码：两自由度duffing振子非线性系统的映射图程序 www.youwenfan.com/contentcno/96561.html
3. MATLAB工具： Dynamical Systems Toolbox：提供分岔分析函数 Control System Toolbox：用于稳定性验证