回复 15 # 雨人 的帖子
%% 输入系统的结(3个状态变量)
R=c1*Y(:,1)+c2*Y(:,2)+c3*Y(:,3)-d;%庞加莱平面(R=0)
RES=zeros(floor(length(R)/2),4);
k=0;
% Use intersection from - to +
for i=2:length(R)-2
if sign(R(i))<sign(R(i+1)) %%-1->>0->>1,随时间变大,正穿越
t1=abs(R(i)/R(i+1))*(T(i+1)-T(i))/(1+abs(R(i)/R(i+1)));
nT=T(i)+t1;
%% 插值计算
nD=[spline(T(i-1:i+2),Y(i-1:i+2,1),nT),spline(T(i-1:i+2),Y(i-1:i+2,2),nT), spline(T(i-1:i+2),Y(i-1:i+2,3),nT)];
k=k+1;
RES(k,1)=nT; RES(k,2:end)=nD; %% 庞加莱截面上的点
end
end
%% plot(RES(:,1),RES(:,2),RES(:,3))
|