求高人指导庞加莱截面程序是否正确！

振动联盟浪子

本人能力有限，对庞加莱截面没有很好的理解，从论坛找了一个类似的程序修改了一下，但感觉相平面图和相应的庞加莱截面对应不起来，或者说时间历程图和相平面和庞加莱都对应不起来，不知道是哪里的错误？？？程序如下
clc
clear
s=1;aerf=0.0001;c=0.0001;epuslg=0.001;beita=0.6;u0=4;u1=0.5;w=16.535;w1=10;g=0;h=0.1;f1=1.58815;f2=3.36218;j1=0.8309;a31=-314.77925;a33=- 0.0501547;a34=20.7204;a42=-3107.80261;a44= -0.3804141;c11=12.3025;c22=46.0489;b12=-3.3421;b21=3.3421;d11=-12.3025;d22=-46.0489;
odefun=@(t,q)[q(3)
              q(4)
              a31*q(1)+a33*q(3)+a34*q(4)+d11*q(1)*(s*(c11*q(1)^2+c22*q(2)^2)+2*s*aerf*(c11*q(1)*q(3)+c22*q(2)*q(4)))+(-2*u1*beita^0.5*b12*q(4)-2*u0*u1*d11*q(1))*epuslg*sin(w1*t)-2*d11*q(1)*epuslg*beita^0.5*u1*w1*cos(w1*t)+(h*w^2*sin(w*t)+h*c*cos(w*t))*j1;
             ( a42*q(2)-a34*q(3)+a44*q(4)+d22*q(2)*(s*(c11*q(1)^2+c22*q(2)^2)+2*s*aerf*(c11*q(1)*q(3)+c22*q(2)*q(4)))+(-2*u1*beita^0.5*b21*q(3)-2*u0*u1*d22*q(2))-0.5*epuslg*sin(w1*t)*d22*q(1)*epuslg*beita^0.5*u1*w1*cos(w1*t))];
q0 =[0;0;0;0];
options=odeset('reltol',1e-8);
tic
[t,q]=ode45(odefun,[0,20],q0,options);
toc
figure('position',[30 400 400 370]);
subplot(311)%绘制时间历程图
plot(t,f1*(q(:,1)))
title('一阶截断部分')
xlabel('无量纲时间τ')
ylabel('η1')
subplot(312)
plot(t,f2*(q(:,2)))
title('二阶截断部分')
xlabel('无量纲时间τ')
ylabel('η2')
subplot(313)
plot(t,f1*q(:,1)+f2*q(:,2))
title('前两阶之和')
xlabel('无量纲时间τ')
ylabel('η')
figure('position',[447 400 400 370]);%绘制相平面图
plot(f1*q(:,1)+f2*q(:,2),f1*q(:,3)+f2*q(:,4))
xlabel('η')
ylabel('dη/dτ')
x=zeros(length(t),3);%绘制庞加莱截面图
x(:,1)=f1*q(:,1)+f2*q(:,2);x(:,2)=f1*q(:,3)+f2*q(:,3);x(:,3)=t;
T=6.28/200;T0=T*1/2; % 选择截面
kmax=round(max(x(:,3))/T);
X1=zeros(kmax);X2=zeros(kmax);
for k=1:kmax;
d=x(:,3)-(k-1)*T-T0;
[P,K]=sort(abs(d));
x1l=x(K(1),1);
x1r=x(K(2),1);
x2l=x(K(1),2);
x2r=x(K(2),2);
x3l=x(K(1),3);
x3r=x(K(2),3);
if abs(P(1))+abs(P(2))<3e-16;
X1(k)=x1l;
X2(k)=x2l;
else
Q=polyfit([x3l,x3r],[x1l,x1r],1);
X1(k)=polyval(Q,(k-1)*T-T0);
Q=polyfit([x3l,x3r],[x2l,x2r],1);
X2(k)=polyval(Q,(k-1)*T-T0);
end
end
figure('position',[865 400 400 370]);
plot(X1,X2,'.');
grid on
xlabel('η');
ylabel('dη/dτ');

我看这个时域图，应该是周期运动吧？如果是相平面图应该是最后形成已封闭曲线，庞加莱截面上应该是一个点吧，并且相平面图和庞加莱图的坐标为什么差距那么大（当时间截面取不同的值时，有时会更大）？？？看别人的相图和庞加莱截面图的坐标范围都一致的，是程序的问题吗？ 庞加莱截面图我不确定是对的，但相平面图应该是没错的啊 但上面的时域图和相平面图也显示矛盾！！！请各位大侠 指导一下啊

振动联盟浪子

继续等待...

振动联盟浪子

继续等待...

微风breeze

我觉得这个程序有点错误，X1(k)=polyval(Q,(k-1)*T-T0);不是减T0，而应该是加T0

振动联盟浪子

回复 4 # 微风breeze 的帖子

虽然不理解为什么这么改，但改了也没有什么变化:@(

微风breeze

你先把程序看懂。。。

sl0412

这个要好好学习一下

振动联盟浪子

回复 7 # sl0412 的帖子

此程序很可能是有问题的 建议就不要学习了 呵呵

1634259766

你这太专业了吧！

张某某zyl

就是不知道截面怎么取？