两自由度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/dtduffing_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/dty(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π）fori=2:n_steps% 检测是否穿越x1=0平面（从负到正）if(x1(i-1)<0&&x1(i)>=0)||(x1(i-1)>0&&x1(i)<=0)% 使用线性插值获得更精确的穿越点ifabs(x1(i)-x1(i-1))>1e-10alpha=-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));% 只保留瞬态过程之后的数据ift_cross>skip_transient poincare_points=[poincare_points;x2_cross,v2_cross];% 计算激励相位（对2π取模）phase=mod(omega*t_cross,2*pi);phase_angles=[phase_angles;phase];endendendendfprintf('找到 %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]);elsetext(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);ifsize(poincare_points,1)>10x2_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;endsgtitle('两自由度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)));% 判断周期性ifsize(poincare_points,1)<20fprintf('系统可能处于周期运动状态\n');else% 计算相邻点的距离变化distances=sqrt(diff(poincare_points(:,1)).^2+diff(poincare_points(:,2)).^2);ifstd(distances)<0.05*mean(distances)fprintf('系统可能处于准周期运动状态\n');elsefprintf('系统可能处于混沌运动状态\n');endendend%% 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');forw_idx=1:length(omega_values)omega_current=omega_values(w_idx);% 使用上一个频率的终值作为初始条件（连续性）ifw_idx==1y0_bif=y0;elsey0_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=[];fori=2:length(t_bif)if(x1_bif(i-1)<0&&x1_bif(i)>=0)&&t_bif(i)>400alpha=-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];endendbifurcation_data{w_idx}=poincare_simple;end% 绘制分岔图hold on;forw_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);endendxlabel('激励频率 \omega');ylabel('Poincaré截面上的 x_2');title('两自由度Duffing振子的分岔图');grid on;xlim([min(omega_values),max(omega_values)]);fprintf('分岔图绘制完成！\n');%% 辅助函数：带可变频率的系统方程functiondydt=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];endfprintf('\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

参数调整建议

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

参数	典型范围	对系统行为的影响
k₃	0.1~10	控制非线性强度，值越大混沌区域越广
F₁, F₂	0~2	激励幅值，影响系统响应幅值和混沌阈值
ω	0.5~2.5	激励频率，改变共振特性和分岔结构
c₁, c₂	0.01~0.5	阻尼系数，抑制混沌和振荡