振动联盟浪子 发表于 2012-5-1 10:26

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

本人能力有限,对庞加莱截面没有很好的理解,从论坛找了一个类似的程序修改了一下,但感觉相平面图和相应的庞加莱截面对应不起来,或者说时间历程图和相平面和庞加莱都对应不起来,不知道是哪里的错误???程序如下
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 =;
options=odeset('reltol',1e-8);
tic
=ode45(odefun,,q0,options);
toc
figure('position',);
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',);%绘制相平面图
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;
=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(,,1);
X1(k)=polyval(Q,(k-1)*T-T0);
Q=polyfit(,,1);
X2(k)=polyval(Q,(k-1)*T-T0);
end
end
figure('position',);
plot(X1,X2,'.');
grid on
xlabel('η');
ylabel('dη/dτ');

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

振动联盟浪子 发表于 2012-5-3 15:44

继续等待...

振动联盟浪子 发表于 2012-5-7 16:49


继续等待...

微风breeze 发表于 2012-5-13 10:27

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

振动联盟浪子 发表于 2012-5-13 10:53

回复 4 # 微风breeze 的帖子

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

微风breeze 发表于 2012-5-13 14:31

你先把程序看懂。。。

sl0412 发表于 2012-5-17 12:26

这个要好好学习一下

振动联盟浪子 发表于 2012-5-17 15:52

回复 7 # sl0412 的帖子

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

1634259766 发表于 2012-5-17 22:03

你这太专业了吧!

张某某zyl 发表于 2015-3-4 20:31

就是不知道截面怎么取?
页: [1]
查看完整版本: 求高人指导庞加莱截面程序是否正确!