两自由度碰振程序不知道哪出问题了？

lihaitao123

clc
clear all
global delta a b c d  m r w f
delta=0;
m=0.5;
r=0.8;
w=4.1;
f=10;
w2=1;
delta=0;
a=(m-r)/(m+1);
b=(1-r)/(m+1);
c=m*(1+r)/(m+1);
d=(1-m*r)/m+1;
tstart=0;
tfinal=10;
  u0=[1,-1,-1,1];
tout=tstart;
  uout=u0.';
options=odeset('Event','on');
for i=1:15
     [t,u,event]=ode45(@pengzhen,[tstart,0.03,tfinal],u0,options);
     tout=[tout;t(2:end)];
     uout=[uout;u(2:end,:)];
     u0(1)=u(end,1);
     u0(2)=u(end,2);
       v10=u(end,3);
       v20=u(end,4);
     u0(3)=a*v10+b*v20;
     u0(4)=c*v10+d*v20;
     tstart=t(end);
end
figure
ylalel('位移');xlabel('时间');
  subplot(2,1,1)
  plot(tout,uout(:,1))
  subplot(2,1,2)
  plot(tout,uout(:,2))

lihaitao123

回复 1 # lihaitao123 的帖子

function vararout=pengzhen(t,u,flag) %%判断何时运用时间
switch flag
    case''
        varargout1=f(t,u);
    case'events'
        [varargout1:3]=events(t,u);
    otherwise
        error{'Unkown flag ' ''flag''  '.'};
end

lihaitao123

回复 2 # lihaitao123 的帖子

function dy=f(t,u)    %%子函数
global a b c d  m r w f
m=0.5;
r=0.8;
w=4.1;
f=10;
w2=1;

du=[u(3);
    u(4);
   -u(1)+cos(w.*t);
   -w2^2*u(2)+f*cos(w.*t)];

function[value,isterminaldirection]=events(t,u)
Q=u(1)-u(2);value=Q;
isterminal=1;
direction=-1