求两自由度碰撞碰撞真相
本帖最后由 lihaitao123 于 2011-12-4 22:45 编辑我做这图和陆启韶原文图差距很大，不知道怎么回事啊，希望搞这方面指教！
clc
clear all;
global m r wf
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
=ode45('xqythkfun',,y0,options);
%\fs{将每次得到的数据依次存在同一矩阵}
tout=;
yout=;
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'
=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 =events(t,y)
Q=y(1)-y(2); %\fs{当Q为0时,解微分方程终止}
value=Q;
isterminal=1;%\fs{开启判断终止功能}
direction=-1;%\fs{由Q减小的方向终止}
陆启韶的原图太大原图，不过能肯定是时程是周期1的
页:
[1]