积分300周期去掉前50周期的瞬态解

bohua1208

clear;close all;clc;
% Set integration time
a0=0.05;f1=5;f2=10;mu=0.1;w0=1;beta=0.1;ks=2;alafa=1.1;
b0=(f1-(f1^2-4*mu^2*w0^2*a0^2)^0.5)/(2*mu*w0);
omega=1.1;
t=2*pi/omega;
tspan=[0 300*t];
% Set initial conditions for the equations
x0=[a0;b0];
[t,x]=ode45('poincarem',tspan,x0);
N1=length(x(:,1))-50;
ZZ=x(N1:end,1);
YY=x(N1:end,2);
XX=t(N1:end);
figure(1)
plot(x(N1:end,1),x(N1:end,2),'k','Markersize',14)
figure(2)
plot(t(N1:end),x(N1:end,1),'k','Markersize',14)
i=1;
for h=2:length(ZZ);
if ZZ(h)>0 & ZZ(h-1)<0
ZZ1(i)=1;
a=(0-ZZ(h-1))/ZZ(h);
XX1(i)=(XX(h-1)+a*XX(h))/(1+a);
YY1(i)=(YY(h-1)+a*YY(h))/(1+a);
i=i+1;
end
end
T=XX1;
XX1(1)=[];T(end)=[];
T1=XX1-T;
YY1(end)=[];
figure(3)
plot(T1,YY1,'k.','Markersize',14)
xlabel('x');ylabel('dx/dt');
请问这是画poincare截面的程序吗？我想做的是积分300周期去掉前50周期的瞬态解，怎么改下程序啊？N1=length(x(:,1))-50;
这是去掉多少周期瞬态？啊

meiyongyuandeze

回复 1 # bohua1208 的帖子

N1=length(x(:,1))-50;
是去掉了前50个积分步长的瞬态解，不一定是前50个周期！

bohua1208

回复 2 # meiyongyuandeze 的帖子

怎么改下，才能去掉前50个周期呢

meiyongyuandeze

回复 3 # bohua1208 的帖子

你可以在调用ode45的时候采用定步长积分，把步长定为周期的多少分之一就可以了！
[t,x]=ode45('poincarem',[0:t/100:300*t],x0);
N1=length(x(:,1))-50*100