小球在弹簧顶端的木块上的弹性跳动问题之问题描述

小球在弹簧顶端的木块上的弹性跳动问题

                                      ——问题描述

 

最近做了个小项目，是关于用Matlab解决小球碰撞问题的，觉得这个例子对用Matlab解微分方程有启发性的作用。我将以系列课程的形式写出来与大家分享。通过该案例，你将了解用Matlab解微分方程的基本过程，同时对Matlab函数（function）功能，odefile模板，绘图技巧将有更深层次的认识。闲话不多说，今天先把这个案例的理论基础与全部代码分享给大家。

                                               

图1.实体模型

 

如图1所示，实体模型由两部分组成，一是小球，二是底面连接弹簧的木板。小球从高为h的位置自然下落后撞击在木板上，此后，木板因为受冲击往下运动，小球与木板撞击后向上运动，如果忽略一切摩擦力，则小球与木板会一直运动下去。我们首先建立坐标系，以弹簧自然伸长位置的上端作为坐标原点，y轴竖直向上，m1，m2分别表示小球与木块的质量，h为小球自由落体的高度，k为弹簧的弹性模量，如下图所示。

图2.参考坐标系建立

 

为了分析整个运动过程，把运动分为两个部分：

部分一：小球与木板没有接触时

该情况下，根据各自的受力情况，列出方程组（1）：

                     （1）

其中，y1，y2分别为小球和木块的位移，g为重力加速度。对方程进行降阶，令  y3=dy1/dt y4=dy2/dt，y3，y4的物理意义为小球与木块的速度，则方程组（1）变为方程组（2），

                                                       （2）

部分二：小球与木板碰撞时

该碰撞过程十分短暂，忽略碰撞过程中二者的高度变化，根据动量守恒与能量守恒定律有方程组（3），其中y3'与y4'分别为碰撞后小球与木板的速度，碰撞后用此速度代替原来二者的速度。

      （3）

化简后有方程组（4）

                 

                                （4）

所有问题变成了如何解方程组2与4，下节课将对此进行分析。

 

PS.今天给出这个案例的全部代码，最近将对这些代码进行剖析，敬请期待：

%%以下是主函数，当主函数1与子函数1,2,3在同一目录下，运行主

%%函数即可

h0=50;

k=60;

m1=20;m2=50;

tstart=0;

tfinal=1000;

y0=[h0;0;0;0];

tout=tstart;

yout=y0';

options=odeset('Events','on');

for i=1:25

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

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

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

    y0(1)=y(end,1);y0(2)=y(end,2);

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

    y0(3)=(-m2*v10+2*m2*v20+m1*v10)/(m2+m1);

    y0(4)=(2*m1*v10+m2*v20-v20*m1)/(m2+m1);

    tstart=t(end);

end

figure

ylabel('高度');

xlabel('时间');

hold on

plot(tout,yout(:,1),tout,yout(:,2));

legend('小球','木板');

figure

axis([-1 1 -50 h0+10]);

axis off

hold on

yt1=-45:0.3:0;

xt1=0.06*sin(yt1);

tanhuang=line(xt1,yt1,'color','k','erasemode','xor','linewidth',2);

qiu=line(0,yout(1,1)+4,'color','k','erasemode','xor','marker','.','markersize',50);

tank=line([-0.1,0.1],[yout(1,2),yout(1,2)],'color',[0.3 0.1 0.5],'erasemode','xor','linewidth',8);

ground=line([-0.5 0.5],[-50 -50],'color',[0.6 0.1 0.2],'linewidth',20);

%%

for i=1:length(tout)

    yt=-45:0.3:yout(i,2);%    xt=0.06*sin((yt-yout(i,2))*(-45)./(-45-yout(i,2)));

    set(tanhuang,'xdata',xt,'ydata',yt);

    set(qiu,'ydata',yout(i,1)+4);

    set(tank,'ydata',[yout(i,2),yout(i,2)]);

    drawnow;

end

%%以下是子函数 1

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

%%以下是子函数 2

function ydot=f(t,y)

ydot=[y(3);y(4);-9.8;-9.8-(k/m2)*y(2);];

%%以下是子函数3

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

Q=y(1)-y(2);

value=Q;

isterminal=1;

direction=-1;

Ps.一些基本的注释因为出现乱码，全部删除，故给出代码下载位置。

代码下载位置

http://www.jishulink.com/content/doc/feebd153-d319-454b-a3dd-18e6624aa9e9