两自由度Duffing振子非线性系统映射图（Poincaré截面） 的MATLAB程序

两自由度Duffing振子非线性系统映射图（Poincaré截面） 的MATLAB程序

完整MATLAB程序

%% 两自由度Duffing振子Poincaré映射图程序
% 系统方程： 
% m1*x1'' + c1*x1' + k1*x1 + k3*(x1-x2)^3 = F1*cos(ω*t)
% m2*x2'' + c2*x2' + k2*x2 + k3*(x2-x1)^3 = F2*cos(ω*t)

clear all; close all; clc;

%% 1. 参数设置
% 质量参数
m1 = 1.0;
m2 = 1.0;

% 线性阻尼系数
c1 = 0.1;
c2 = 0.1;

% 线性刚度系数
k1 = 1.0;
k2 = 1.0;

% 非线性刚度系数（立方项）
k3 = 0.5;

% 激励参数
F1 = 0.8;      % 第一个振子的激励幅值
F2 = 0.3;      % 第二个振子的激励幅值
omega = 1.2;   % 激励频率

% 仿真参数
T_total = 5000;        % 总仿真时间（足够长以达到稳态）
dt = 0.01;            % 时间步长
tspan = 0:dt:T_total; % 时间向量
n_steps = length(tspan);

% Poincaré截面参数（选择x1=0平面）
section_tol = 1e-3;   % 截面检测容差
skip_transient = 100; % 跳过初始瞬态过程的时间

%% 2. 定义系统微分方程（状态空间形式）
% 状态变量: y = [x1; v1; x2; v2]
% 其中: v1 = dx1/dt, v2 = dx2/dt

duffing_sys = @(t, y) [
    y(2);  % dx1/dt = v1
    (F1*cos(omega*t) - c1*y(2) - k1*y(1) - k3*(y(1)-y(3))^3) / m1;  % dv1/dt
    y(4);  % dx2/dt = v2
    (F2*cos(omega*t) - c2*y(4) - k2*y(3) - k3*(y(3)-y(1))^3) / m2   % dv2/dt
];

%% 3. 数值积分求解系统响应
fprintf('正在求解系统响应...\n');
% 初始条件
y0 = [0.1; 0; 0.2; 0];  % [x1_0; v1_0; x2_0; v2_0]

% 使用Runge-Kutta方法（ode45）进行数值积分
options = odeset('RelTol', 1e-8, 'AbsTol', 1e-10);
[t, Y] = ode45(duffing_sys, tspan, y0, options);

% 提取状态变量
x1 = Y(:, 1); v1 = Y(:, 2);
x2 = Y(:, 3); v2 = Y(:, 4);

fprintf('系统响应求解完成！\n');

%% 4. Poincaré截面采样（在x1=0平面）
fprintf('正在构建Poincaré映射...\n');

% 寻找x1穿过零平面的点（从负到正）
poincare_points = [];
phase_angles = [];  % 激励相位（mod 2π）

for i = 2:n_steps
    % 检测是否穿越x1=0平面（从负到正）
    if (x1(i-1) < 0 && x1(i) >= 0) || (x1(i-1) > 0 && x1(i) <= 0)
        % 使用线性插值获得更精确的穿越点
        if abs(x1(i) - x1(i-1)) > 1e-10
            alpha = -x1(i-1) / (x1(i) - x1(i-1));
            t_cross = t(i-1) + alpha * dt;
            x2_cross = x2(i-1) + alpha * (x2(i) - x2(i-1));
            v1_cross = v1(i-1) + alpha * (v1(i) - v1(i-1));
            v2_cross = v2(i-1) + alpha * (v2(i) - v2(i-1));
            
            % 只保留瞬态过程之后的数据
            if t_cross > skip_transient
                poincare_points = [poincare_points; x2_cross, v2_cross];
                
                % 计算激励相位（对2π取模）
                phase = mod(omega * t_cross, 2*pi);
                phase_angles = [phase_angles; phase];
            end
        end
    end
end

fprintf('找到 %d 个Poincaré截面点\n', size(poincare_points, 1));

%% 5. 可视化结果
fprintf('正在绘制结果...\n');

% 5.1 时域响应
figure('Position', [100, 100, 1200, 800]);

subplot(3, 3, 1);
plot(t, x1, 'b', 'LineWidth', 1);
xlabel('时间 t'); ylabel('位移 x_1');
title('第一个振子的时域响应');
grid on; xlim([T_total-100, T_total]); % 显示最后100秒

subplot(3, 3, 2);
plot(t, x2, 'r', 'LineWidth', 1);
xlabel('时间 t'); ylabel('位移 x_2');
title('第二个振子的时域响应');
grid on; xlim([T_total-100, T_total]);

subplot(3, 3, 3);
plot(x1, v1, 'b', 'LineWidth', 0.5);
xlabel('位移 x_1'); ylabel('速度 v_1');
title('第一个振子的相轨迹');
grid on; axis equal;

% 5.2 Poincaré映射图
subplot(3, 3, [4, 5, 6]);
if ~isempty(poincare_points)
    scatter(poincare_points(:,1), poincare_points(:,2), 10, 'filled', ...
        'MarkerFaceColor', 'b', 'MarkerEdgeColor', 'b', 'MarkerFaceAlpha', 0.6);
    xlabel('位移 x_2'); ylabel('速度 v_2');
    title('Poincaré映射图 (截面: x_1=0)');
    grid on;
    
    % 根据点的聚集程度自动调整坐标轴
    x_center = mean(poincare_points(:,1));
    y_center = mean(poincare_points(:,2));
    x_range = max(abs(poincare_points(:,1) - x_center));
    y_range = max(abs(poincare_points(:,2) - y_center));
    axis([x_center-1.2*x_range, x_center+1.2*x_range, ...
          y_center-1.2*y_range, y_center+1.2*y_range]);
else
    text(0.5, 0.5, '未找到Poincaré截面点', 'HorizontalAlignment', 'center');
end

% 5.3 相位与Poincaré点的关系
subplot(3, 3, 7);
if ~isempty(phase_angles) && ~isempty(poincare_points)
    scatter(phase_angles, poincare_points(:,1), 10, 'filled', ...
        'MarkerFaceColor', 'r', 'MarkerEdgeColor', 'r');
    xlabel('激励相位 (rad)'); ylabel('x_2');
    title('Poincaré点 vs 激励相位');
    grid on; xlim([0, 2*pi]);
end

% 5.4 位移关系图
subplot(3, 3, 8);
plot(x1, x2, 'g', 'LineWidth', 0.5);
xlabel('位移 x_1'); ylabel('位移 x_2');
title('两个振子的位移关系');
grid on; axis equal;

% 5.5 庞加莱点的返回映射（一维映射）
subplot(3, 3, 9);
if size(poincare_points, 1) > 10
    x2_points = poincare_points(:,1);
    plot(x2_points(1:end-1), x2_points(2:end), 'k.', 'MarkerSize', 10);
    hold on;
    plot([min(x2_points), max(x2_points)], [min(x2_points), max(x2_points)], 'r--', 'LineWidth', 1);
    xlabel('x_2^{(n)}'); ylabel('x_2^{(n+1)}');
    title('Poincaré映射的一维返回映射');
    grid on; axis equal;
end

sgtitle('两自由度Duffing振子非线性动力学分析', 'FontSize', 14, 'FontWeight', 'bold');

%% 6. 计算统计特征（可选）
if ~isempty(poincare_points)
    fprintf('\n=== Poincaré映射统计特征 ===\n');
    fprintf('截面点数量: %d\n', size(poincare_points, 1));
    fprintf('x2的均值: %.4f\n', mean(poincare_points(:,1)));
    fprintf('x2的标准差: %.4f\n', std(poincare_points(:,1)));
    fprintf('v2的均值: %.4f\n', mean(poincare_points(:,2)));
    fprintf('v2的标准差: %.4f\n', std(poincare_points(:,2)));
    
    % 判断周期性
    if size(poincare_points, 1) < 20
        fprintf('系统可能处于周期运动状态\n');
    else
        % 计算相邻点的距离变化
        distances = sqrt(diff(poincare_points(:,1)).^2 + diff(poincare_points(:,2)).^2);
        if std(distances) < 0.05 * mean(distances)
            fprintf('系统可能处于准周期运动状态\n');
        else
            fprintf('系统可能处于混沌运动状态\n');
        end
    end
end

%% 7. 绘制分岔图（作为扩展功能）
figure('Position', [100, 100, 800, 600]);

% 改变激励频率，观察分岔现象
omega_values = linspace(0.5, 2.0, 100);  % 频率扫描范围
bifurcation_data = cell(length(omega_values), 1);

fprintf('\n正在计算分岔图（这可能需要一些时间）...\n');
for w_idx = 1:length(omega_values)
    omega_current = omega_values(w_idx);
    
    % 使用上一个频率的终值作为初始条件（连续性）
    if w_idx == 1
        y0_bif = y0;
    else
        y0_bif = Y_end;
    end
    
    % 短时间仿真获取稳态
    t_bif = 0:0.1:500;
    [~, Y_bif] = ode45(@(t,y) duffing_sys_omega(t, y, omega_current), ...
                       t_bif, y0_bif, options);
    
    % 保存最后一个周期的Poincaré点
    Y_end = Y_bif(end, :)';
    
    % 寻找Poincaré点（简化版本）
    x1_bif = Y_bif(:, 1);
    x2_bif = Y_bif(:, 2);
    poincare_simple = [];
    
    for i = 2:length(t_bif)
        if (x1_bif(i-1) < 0 && x1_bif(i) >= 0) && t_bif(i) > 400
            alpha = -x1_bif(i-1) / (x1_bif(i) - x1_bif(i-1));
            x2_cross = x2_bif(i-1) + alpha * (x2_bif(i) - x2_bif(i-1));
            poincare_simple = [poincare_simple; x2_cross];
        end
    end
    
    bifurcation_data{w_idx} = poincare_simple;
end

% 绘制分岔图
hold on;
for w_idx = 1:length(omega_values)
    points = bifurcation_data{w_idx};
    if ~isempty(points)
        plot(omega_values(w_idx)*ones(size(points)), points, 'b.', ...
            'MarkerSize', 4);
    end
end

xlabel('激励频率 \omega'); ylabel('Poincaré截面上的 x_2');
title('两自由度Duffing振子的分岔图');
grid on;
xlim([min(omega_values), max(omega_values)]);

fprintf('分岔图绘制完成！\n');

%% 辅助函数：带可变频率的系统方程
function dydt = duffing_sys_omega(t, y, omega_var)
    % 参数（使用主程序中的参数，这里需要重新定义或传入）
    m1 = 1.0; m2 = 1.0;
    c1 = 0.1; c2 = 0.1;
    k1 = 1.0; k2 = 1.0;
    k3 = 0.5;
    F1 = 0.8; F2 = 0.3;
    
    dydt = [
        y(2);
        (F1*cos(omega_var*t) - c1*y(2) - k1*y(1) - k3*(y(1)-y(3))^3) / m1;
        y(4);
        (F2*cos(omega_var*t) - c2*y(4) - k2*y(3) - k3*(y(3)-y(1))^3) / m2
    ];
end

fprintf('\n程序运行完毕！\n');

程序核心原理说明

1. Poincaré映射的核心思想

Poincaré映射通过记录系统轨迹与某个选定的超平面（这里选择x₁=0平面）的交点，将连续动力系统转换为离散映射。这些交点构成的图形能清晰地展示系统的周期、准周期或混沌行为。

2. 两自由度Duffing系统的特殊性

• 耦合非线性：两个振子通过立方非线性项k₃*(x₁-x₂)³相互耦合
• 双频激励：每个振子可以有不同幅值的激励
• 四维相空间：状态变量为[x₁, v₁, x₂, v₂]，Poincaré截面将其降维展示

3. 关键实现技术

• 穿越检测算法：精确检测轨迹穿越x₁=0平面的时刻
• 线性插值：提高截面点位置的精度
• 瞬态跳过：忽略初始瞬态过程，只分析稳态行为

参考代码 两自由度duffing振子非线性系统的映射图程序 www.3dddown.com/csa/96561.html

参数调整建议

你可以通过修改以下参数探索不同的动力学行为：

结果解读指南

程序输出的图形中：

1. Poincaré映射图：

   • 有限个离散点 → 周期运动
   • 闭合曲线 → 准周期运动
   • 复杂分形结构 → 混沌运动

2. 分岔图：展示系统行为随参数（如频率ω）变化的全局图像，可以清晰看到周期倍化通向混沌的路径。