求两自由度碰撞碰撞真相

lihaitao123

我做这图和陆启韶原文图差距很大，不知道怎么回事啊，希望搞这方面指教！
clc
clear all;

global   m r w  f
m=0.5;
r=0.8;

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=100;

y0=[.5;-1;0;0];    tout=tstart;     yout=y0.';

options=odeset('Events','on');   %\fs{开启事件判断功能}



for i=1:6

    [t,y,event]=ode45('xqythkfun',[tstart:0.3:tfinal],y0,options);

       %\fs{将每次得到的数据依次存在同一矩阵}

   tout=[tout;t(2:end)];   

   yout=[yout;y(2:end,:)];

   y0(1)=y(end,1);   y0(2)=y(end,2);%\fs{下一次弹跳的初位移}

   v10=y(end,3);    v20=y(end,4);

   y0(3)=a*v10+b*v20;
   y0(4)=c*v10+d*v20;

   tstart=t(end);

end

%\fs{时程}

figure

ylabel('位移');
xlabel('时间');

subplot(1 ,2 ,1)
plot(tout,yout(:,1),'r')
subplot(1, 2 ,2)
plot(tout,yout(:,2),'k')

function varargout=xqythkfun(t,y,flag)

switch flag

case ''

   varargout{1}=f(t,y);

case 'events'

   [varargout{1:3}]=events(t,y);

otherwise

   error(['Unknown flag ''' flag '''.']);

end

%\fs{---------------------------------------------}

%\fs{计算微分方程的子函数}

function ydot=f(t,y)

global w f

w=4.1;
f=10;
w2=1;


ydot=[y(3);

     y(4);

   -y(1)+cos(w*t);

   -w2^2*y(2)+f*cos(w*t);];

%\fs{----------------------------------------------}

%\fs{事件判断子函数}

function [value,isterminal,direction]=events(t,y)

Q=y(1)-y(2);   %\fs{当Q为0时,解微分方程终止}

value=Q;

isterminal=1;  %\fs{开启判断终止功能}

direction=-1;  %\fs{由Q减小的方向终止}

lihaitao123

陆启韶的原图太大原图，不过能肯定是时程是周期1的