小球在彈簧頂端的木塊上的彈性跳動(dòng)問題

                                      ——問題描述

 

最近做了個(gè)小項(xiàng)目，是關(guān)于用Matlab解決小球碰撞問題的，覺得這個(gè)例子對(duì)用Matlab解微分方程有啟發(fā)性的作用。我將以系列課程的形式寫出來(lái)與大家分享。通過(guò)該案例，你將了解用Matlab解微分方程的基本過(guò)程，同時(shí)對(duì)Matlab函數(shù)（function）功能，odefile模板，繪圖技巧將有更深層次的認(rèn)識(shí)。閑話不多說(shuō)，今天先把這個(gè)案例的理論基礎(chǔ)與全部代碼分享給大家。

                                               

圖1.實(shí)體模型

 

如圖1所示，實(shí)體模型由兩部分組成，一是小球，二是底面連接彈簧的木板。小球從高為h的位置自然下落后撞擊在木板上，此后，木板因?yàn)槭軟_擊往下運(yùn)動(dòng)，小球與木板撞擊后向上運(yùn)動(dòng)，如果忽略一切摩擦力，則小球與木板會(huì)一直運(yùn)動(dòng)下去。我們首先建立坐標(biāo)系，以彈簧自然伸長(zhǎng)位置的上端作為坐標(biāo)原點(diǎn)，y軸豎直向上，m1，m2分別表示小球與木塊的質(zhì)量，h為小球自由落體的高度，k為彈簧的彈性模量，如下圖所示。

圖2.參考坐標(biāo)系建立

 

為了分析整個(gè)運(yùn)動(dòng)過(guò)程，把運(yùn)動(dòng)分為兩個(gè)部分：

部分一：小球與木板沒有接觸時(shí)

該情況下，根據(jù)各自的受力情況，列出方程組（1）：

                     （1）

其中，y1，y2分別為小球和木塊的位移，g為重力加速度。對(duì)方程進(jìn)行降階，令  y3=dy1/dt y4=dy2/dt，y3，y4的物理意義為小球與木塊的速度，則方程組（1）變?yōu)榉匠探M（2），

                                                       （2）

部分二：小球與木板碰撞時(shí)

該碰撞過(guò)程十分短暫，忽略碰撞過(guò)程中二者的高度變化，根據(jù)動(dòng)量守恒與能量守恒定律有方程組（3），其中y3'與y4'分別為碰撞后小球與木板的速度，碰撞后用此速度代替原來(lái)二者的速度。

      （3）

化簡(jiǎn)后有方程組（4）

                 

                                （4）

所有問題變成了如何解方程組2與4，下節(jié)課將對(duì)此進(jìn)行分析。

 

PS.今天給出這個(gè)案例的全部代碼，最近將對(duì)這些代碼進(jìn)行剖析，敬請(qǐng)期待：

%%以下是主函數(shù)，當(dāng)主函數(shù)1與子函數(shù)1,2,3在同一目錄下，運(yùn)行主

%%函數(shù)即可

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('時(shí)間');

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

%%以下是子函數(shù) 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

%%以下是子函數(shù) 2

function ydot=f(t,y)

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

%%以下是子函數(shù)3

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

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

value=Q;

isterminal=1;

direction=-1;

Ps.一些基本的注釋因?yàn)槌霈F(xiàn)亂碼，全部刪除，故給出代碼下載位置。

代碼下載位置

http://www.yqgqt.org.cn/content/doc/feebd153-d319-454b-a3dd-18e6624aa9e9